热辐射输运问题的高效蒙特卡罗模拟方法*

许育培 李树

1) (北京应用物理与计算数学研究所,北京 100094)

2) (中国工程物理研究院研究生院,北京 100088)

惯性约束聚变研究中,热辐射光子在介质中的输运以及热辐射光子与介质的相互作用是重要研究课题,蒙特卡罗方法是该类问题的重要研究手段之一.隐式蒙特卡罗方法虽然能正确地模拟热辐射在介质中的输运过程,但当模拟重介质(材料的吸收系数大)问题时,该方法花费的计算时间将变得很长,导致模拟效率很低.本文以离散扩散蒙特卡罗方法为基础,开发了“离散扩散蒙特卡罗方法辐射输运模拟程序”,可以较好地解决重介质区的计算效率问题,但是离散扩散蒙卡罗方法在模拟轻介质区时精度不够高.辐射输运问题中通常既有轻介质也有重介质,为了能同时解决蒙特卡罗方法模拟的效率和精度问题,本文研究了离散扩散蒙特卡罗方法与隐式蒙特卡罗方法相结合的模拟方法,并提出了新的扩散区与输运区界面处理方法,研制了混合蒙特卡罗方法的辐射输运模拟程序.典型辐射输运问题模拟显示:在模拟重介质问题时,该程序能大幅缩短模拟时间,且能取得与隐式蒙特卡罗方法一致的结果;在模拟轻重介质均存在的问题时,与隐式蒙特卡罗方法相比,混合蒙特卡罗方法的模拟精度与其相当且计算效率同样能够得到显著提升.

1 引 言

惯性约束聚变(inertial confinement fusion,ICF)中,热辐射(X光波段)扮演着重要的角色.间接驱动ICF中,黑腔壁产生的热辐射光子烧蚀驱动飞层对热核材料进行压缩,产生高温高密度的等离子体,实现聚变点火.因此,热辐射光子在介质(材料)中的输运以及热辐射光子与介质的相互作用是ICF的重要研究课题之一.

在辐射输运过程中,辐射光子与物质的相互作用很复杂,辐射强度和物质温度强耦合,描述辐射场时空演化的辐射输运方程组是非线性的[1],只能用数值方法求解实际的热辐射输运问题.蒙特卡罗(Monte Carlo,MC)[2,3]方法是常用的求解粒子输运问题的数值方法.传统的显式蒙特卡罗方法基于与时间相关的、非线性的辐射输运方程,通过随机抽样来确定粒子的吸收、散射与发射,以迭代的方式模拟热辐射输运过程[4].然而,当物质的吸收系数较大或系统接近热力学平衡态时,传统的显式蒙特卡罗法的模拟结果存在严重的涨落和能量不均衡等问题[5].

为解决传统蒙特卡罗方法存在的问题,Fleck和Cummings[6]提出了隐式蒙特卡罗(implicit Monte Carlo,IMC)方法.该方法天然具有无条件稳定性,较传统的显式蒙特卡罗方法更稳定,精度、效率更高.然而,当系统的吸收系数变得更大时,IMC方法会出现计算时间过长的问题,导致模拟效率低下.为了提高模拟效率,Fleck和Canfield[7]基于概率方法首先提出“随机游走法”(randomwalk,RW),Giorla和Sentis[8]基于扩散方程推导了RW方程,建立了以粒子所在位置为球心的球形子区域以及判据,当判据满足时,球形子区域内粒子的多次散射过程将由一个扩散过程代替,从而有效提高了模拟效率.然而,球形子区域的半径有最小值(其半径不能小于数倍粒子平均自由程),当粒子接近边界时,球形子区域半径与粒子平均自由程的关系将不再满足,RW方法将失效,只能转而用IMC方法输运.因此,若吸收系数不够大,或系统网格划分较小,RW模拟辐射输运的效率提高则不明显.20世纪90年代,Urbatsch等[9]首次提出离散扩散蒙特卡罗(discrete diffusion Monte Carlo,DDMC)方法.DDMC方法从离散的扩散方程出发,把IMC方法中的多次散射过程替换成一个扩散过程,在保证相同精度的前提下,很大程度上提高了辐射输运模拟的效率.之后,Evans等[10]推导出平衡态下的DDMC方法;Gentile[11]对扩散方程进行空间和时间离散,得到矩阵方程,推导出一维的时间、空间离散DDMC方法;Densmore等[12,13]对扩散方程进行空间离散,进而推导出一维的时间连续、空间离散DDMC方法;Cleveland等[14]提出多群 DDMC (原作者称其为implicit Monte Carlo diffusion,IMD)方法,以解决与频率相关的热辐射输运问题;Densmore等[15]基于频率积分扩散方程,开发了与频率有关的DDMC新方法.

国外已有基于DDMC方法开发的数值模拟程序,并与IMC方法结合以求解热辐射输运问题,但由于问题的敏感性,国外DDMC方法模拟程序实现的关键环节及程序代码并未对外公开.经过多年努力,国内已开发出较成熟的IMC方法辐射输运模拟程序[16],且已经成功用于ICF方法实验的模拟与分析[17],但是至今未有基于DDMC方法的辐射输运模拟程序.因此有必要独立开展DDMC方法研究并开发相应的数值模拟程序,应用于ICF方法相关课题研究中.本文在DDMC方法基础上,解决了DDMC粒子(用DDMC输运的粒子)输运模拟的关键技术,研制了DDMC方法辐射输运模拟程序,研究了DDMC与IMC方法的耦合方法,并提出一种新的界面处理方法,进而研制了DDMC-IMC程序.该程序能准确且高效地模拟不同光性厚度介质中的热辐射输运问题.

2 DDMC方法及数值模拟流程设计

下面以一维平板问题为例,介绍DDMC方法的主要思想并推导相应的数学方程.考虑一维、均匀、无独立外源的无限平板,热辐射输运方程及能量方程[6,18]为

其中,I为辐射强度,t为时间变量,c为真空中的光速,µ为粒子运动的方向余弦,x为空间变量,σ为介质的吸收系数,a为辐射常数,T为介质温度,Cv为介质比热.

将平板网格化,即 0=x1＜x2＜x3＜··＜xI,x1,x2,··,xI为网格边界所在的位置,每个网格中的位置变量满足xi＜x＜xi+1;将时间离散化,即0=t1＜t2＜t3＜··tN,t1,t2,··,tN为各个时间步结束的时刻,每个时间步中的时间变量满足tn＜t＜tn+1.应用IMC方法后,方程(1)改写为[6]

其中带下标n的变量表示第n个时间步中对应的变量值,σes为有效散射系数,σea为有效吸收系数,

其中,f为Fleck因子,Ur为辐射能量密度,Um为物质能量密度,Δt为时间步长.

当介质的吸收系数很大时,粒子的平均自由程将变得很小,由方程(5)可知,fn也将变小,即有效散射系数σes将接近σn,而有效吸收系数σea将接近0,因此系统中粒子的散射将占据主导地位.综上,在吸收系数很大的情况下,IMC中粒子的历史变得过于漫长,导致计算效率降低.

当散射过程占主导时,考虑对方程(2)进行扩散近似.方程(2a)等式两边对µ全角度积分得

其中,φ为粒子通量密度,F为辐射流.对方程(7)中等式左边的空间导数项进行有限差分,

其中,Δx 为网格大小,Fi+1/2和Fi —1/2分别为第i个网格左、右边界的平均辐射流.因此对于第i个网格,方程(7)改写为

其中φi为网格平均辐射强度.

1)当 1＜i＜I时,应用斐克定律[19,20],网格边界平均辐射流的表达式为

其中φi+1/2为网格边界平均辐射强度.

联立方程 (10)和 (11),解出 Fi+1/2得

类似可得 Fi —1/2的表达式.

将 Fi+1/2,Fi —1/2代入方程 (9) 中,得

式中利用了 σR,i—1=σL,i的关系,另外 σL,i和 σR,i分别表示第i个网格中粒子向左、右相邻网格扩散的截面,

方程(13)是DDMC方法中粒子输运所依据的方程,其中等式左边第二项为吸收项,等式右边为源项,右边第一项表示从右边相邻网格向左扩散到当前网格的粒子,右边第二项表示从左边相邻网格向右扩散到当前网格的粒子,右边第三项表示介质辐射的粒子.

在DDMC方法中,粒子的位置变量不再有意义,因此定义粒子平均自由时间τ,以确定粒子在随时间变化过程中发生的反应(扩散、吸收或存活至当前时间步结束),

粒子平均自由时间与粒子平均自由程类似,表示粒子发生两次反应(扩散或吸收)的平均时间间隔.

2)对于系统边界处的网格,即当i=1或i=I时,方程(13)不再适用,因此有必要推导边界处粒子输运满足的离散扩散方程.令方程(9)中i=1,得

应用渐进扩散极限(asymptotic diffusionlimit)边界条件[21]:

其中,W(µ)为超越方程,可以近似为 W(µ) ≈ µ+3/2µ2;φ1/2为系统边界平均辐射强度;λ为外推距离.

对方程(18)中的空间导数项进行有限差分,得

令方程(11)中i+1=1,i+1/2=1/2,并与方程(19)联立,解出F1/2:

将方程(20)代入方程(17)中,得

方程(21)即为i=1时边界网格满足的离散扩散方程,其中

同理可推导出当i=I时,边界网格满足的输运方程.

除了右边第二项外,方程(21)与方程(13)类似.方程(21)右边第二项中PL(µ)表示粒子穿过边界进入DDMC区域后而不被反射的概率,与粒子运动方向和边界法向之间夹角的余弦µ有关.

需要特别指出得是:当模拟的对象为轻介质(吸收系数较小)时,扩散条件不成立,DDMC方法模拟的精度将变差;另外,当吸收截面较大时,上面的PL,σL,1与理论值相比偏小,也需对它们做适当的修正[22].

本文设计了基于DDMC方法的辐射输运问题模拟流程.主要步骤如下.

1) 时空网格划分 根据系统的几何结构及介质属性,将系统划分为若干空间网格,每个网格包含同一种介质;根据辐射源及辐射场演化特性,将要模拟的时间范围离散成若干时间步.

2) 计算当前时间步的各类参数 根据各个网格的温度、密度、材料种类等计算网格比热、吸收系数、物质能量密度、辐射能量密度、有效吸收系数、有效散射系数、fleck因子、粒子扩散截面、粒子转化概率等.

3) 从辐射源中抽取粒子 根据辐射总能量及模拟的粒子数,分配三种源(边界源、网格物质辐射源、上一时间步中存活下来的粒子)的粒子权重及粒子数,依次抽取源对应的粒子并确定其时间、能量、权重.

4) 粒子随机扩散 根据计算得到的扩散截面和有效吸收系数,粒子在各个网格间随机运动,或被物质吸收,或扩散至相邻网格,或存活至当前时间步结束.若粒子未被吸收,则记录粒子的相关参数,在下一时间步继续输运.

5) 求解温度方程 由步骤4)得到的辐射能和吸收能,以及网格的比热,利用能量方程解出当前时间步结束时的物质温度.

6) 预估下一时间步的平均温度 由于需要根据网格温度来计算网格其他参数,因此有必要预估下一时间步的温度,本文采用“预估迭代”的方法,此方法既能保持较好的精度,又能相对节省时间.

7) 返回步骤2),开始下一时间步计算,直至结束时刻.

3 DDMC与IMC耦合方法

ICF辐射输运问题中通常既包含有轻介质(如推进层),也包含有重介质(如黑腔壁).对于轻介质(光性薄)区域,最好应用IMC方法模拟,对于重介质(光性厚)区域,则考虑应用DDMC方法模拟.因此,需要研究DDMC方法模拟与IMC方法模拟的耦合方法.DDMC-IMC方法耦合输运的关键是处理好IMC区域(应用IMC方法模拟的区域)和DDMC区域(应用DDMC方法模拟的区域)的交界面,即如何正确合理的处理MC粒子穿过交界面时的各种物理数学参数,下面介绍其中的两个关键环节.

第一种情况,粒子从IMC区域穿过界面进入DDMC区域.IMC粒子(用IMC方法模拟的粒子)穿过界面时,有P(µ)的概率进入DDMC区域并转化为DDMC粒子,没有转化的粒子在界面处被反射,并重新抽取方向和位置,位置在界面上均匀分布.

一般情况下,DDMC区域与IMC区域的不透明度相差很大,从DDMC区域边界离开后进入IMC区域的粒子的运动方向不再是各向同性的.实际上其运动方向的角分布与面源粒子的角分布一样,为余弦分布[23],即 f(µ)=2µ.因此抽取粒子新方向的方法为

其中ξ为在[0,1]上均匀分布的随机数.

第二种情况,粒子从DDMC区域穿过边界进入IMC区域.DDMC粒子穿过界面后,转化为IMC粒子.因为在DDMC区域中,粒子的方向和位置没有实际意义,因此粒子进入IMC区域后同样需要重新抽取位置和方向变量.与第一种情况一样,粒子方向的抽取按照方程(24)的方法,位置在界面上均匀分布.

4 数值结果

根据上文所述的DDMC方法辐射输运流程以及DDMC-IMC界面处理方法,本文设计、编写了“离散扩散蒙特卡罗方法辐射输运程序”(DDMC)及“离散扩散蒙特卡罗方法与隐式蒙特卡罗方法耦合输运程序”(DDMC-IMC).为验证程序的适用性和正确性,本文计算了若干典型问题,并与已有的IMC程序计算结果比较,计算平台为个人电脑,CPU型号为Intel(R) Core(TM) i7-3770 @ 3.4 GHz,核数为8,内存为4.00 GB,操作系统为Windows 7旗舰版,编译器为Intel Visual Fortran.

4.1 单一介质热波传播问题

模型为0.5 cm厚的一维平板,左侧边界为一温度恒为1 keV的平面源,右侧边界为真空边界,平板分为25个网格,网格大小为0.02 cm,网格比热为cv=0.1 GJ/(cm3·keV),吸收系数与温度的三次方呈反比,即 σ=σ0/T3,其中温度T的单位取keV,σ0单位为keV3/cm,可取不同的值,越大表示光性越厚.后文中吸收系数及温度单位与之相同.

网格初始物质温度为1 eV,初始辐射温度为0 eV,时间步长取0.05 ns,总时长为10 ns,每步模拟粒子数为 10000.σ0分别取 200,500,1000,2000.分别用IMC和DDMC方法模拟并比较.开始进行辐射输运时,光子从边界源出射,与附近的物质相互作用,加热物质;被加热的物质重新辐射出光子,再加热更远处的物质,因此辐射得以从边界源处逐渐向远处传播.

图1和图2分别为σ0=500时,三个时间点的物质温度和辐射温度的空间分布.从两图中可以看出,DDMC与IMC方法的模拟结果总体上符合得很好,且随时间的演化,曲线所描绘出的热波逐渐向前传播,加热远处的物质,与热辐射传播的物理过程相符.

图1 不同时刻的物质温度空间分布比较(σ0=500)Fig.1.Material temperature in different moments (σ0=500).

图2 不同时刻的辐射温度空间分布比较(σ0=500)Fig.2.Radiative temperature in different moments (σ0=500).

表1为不同σ0取值下IMC和DDMC方法的模拟时间,以及DDMC与IMC方法的效率比较.从表1中可以看出,随着σ0的增大,IMC方法的模拟时间明显变长,而DDMC方法的模拟时间变化不大,相比较下DDMC方法的效率也随之提高.因为吸收系数越大,IMC粒子在一个时间步内散射的次数越多,而IMC方法中处理散射的过程占据了整个模拟过程的大部分时间,相反DDMC粒子在一个时间步中的扩散过程随吸收截面变大而略微减少,且一个扩散过程替代了粒子的多次散射过程,因此随着初始截面变大,IMC方法模拟时间明显变长,DDMC方法模拟时间反而有所缩短,最终使得DDMC的模拟效率变高.

表1 不同σ0取值下IMC与DDMC方法的模拟时间对比Table 1.Simulation time of IMC method and DDMC method in different initial cross sections.

4.2 单一介质热平衡问题

模型为0.5 cm厚一维平板,无边界源,且外边界均为全反射边界,平板分为50个网格,网格大小为 0.01 cm,网格比热 Cv=0.1 GJ/(cm3·keV),网格初始物质温度与初始辐射温度分别为1 (0—0.1 cm),0.9 (0.1—0.2 cm),0.8 (0.2—0.3 cm),0.7(0.3—0.4 cm),0.6 (0.4—0.5 cm).时间步长为0.05 ns,总时长为40 ns,每步模拟粒子数为10000.σ0取值为 200,500,1000,2000.同样分别用IMC和DDMC方法模拟并比较.

由于外边界设置为全反射边界,故系统实为全封闭系统.随着时间演化,理论上该系统最终将达到平衡,即所有网格的辐射温度和物质温度都将趋于某固定值.由于全封闭系统的能量守恒,可以建立如下关系:

图3和图4分别是σ0=500时,不同时间点的介质中物质温度和辐射温度的空间分布.从两图中可以看出,DDMC方法的模拟结果与IMC方法同样符合得很好,随着时间变化,温度较高处的网格温度逐渐降低,温度较低处的网格温度逐渐升高,且在40 ns时基本达到了理论平衡温度0.8 keV.模拟结果显示,当系统中存在温度梯度时,能流将从温度高处流向温度低处,其间,辐射光子与物质不断相互作用(光子被吸收后再发射),系统中温度分布、辐射场分布不断发生变化,对于本例,辐射温度与物质温度逐渐趋于理论平衡温度,并在时间足够长后达到平衡.

图3 不同时刻的物质温度空间分布比较(σ0=500)Fig.3.Material temperature in different moments (σ0=500).

表2为不同σ0取值下DDMC和IMC方法的模拟时间以及效率的比较.与表1相似,随着σ0的增大,IMC方法的模拟时间明显变长,而DDMC方法的模拟时间略微缩短,相比较下的DDMC方法效率也随之提高,最高接近30倍,其原因上文已做出解释.同样可得出物质吸收系数越大,DDMC方法模拟效率越高的结论.

图4 不同时刻的辐射温度空间分布比较(σ0=500)Fig.4.Radiative temperature in different moments (σ0=500).

表2 不同σ0取值下,IMC与DDMC方法的模拟时间比较Table 2.Simulation time of IMC method and DDMC method in different initial cross sections.

图4中辐射温度曲线的波动较大,尤其是平衡后各网格间的温度看似并没有达到平衡(温度相等).初步分析其原因是蒙特卡罗方法计算的粒子数偏少、统计误差偏大所致.为此,在其他参数不变的情况下,将模拟的粒子数提高至500000后计算,得到图5和图6所示的物质温度和辐射温度的分布情况.与图3和图4对比可看出,数值波动有明显的减小,说明较大的波动确实是由样本数较少引起的,不影响本文对DDMC方法的优点及模拟结果的分析.

图5 不同时刻的物质温度空间分布比较(σ0=500,Np=500000)Fig.5.Material temperature in different moments (σ0=500,Np=500000).

图6 不同时刻的辐射温度空间分布比较(σ0=500,Np=500000)Fig.6.Radiative temperature in different moments (σ0=500,Np=500000).

4.3 双介质中热波传播问题

模型为0.5 cm厚一维平板,左边界为一温度恒为1 keV的平面源,右边界为真空边界,平板分为左右两部分,左边部分为光性厚的区域,厚度为0.1 cm,划分为20个网格,网格大小为0.005 cm,右边部分为光性薄的区域,厚度为0.4 cm,划分为20个网格,网格大小为0.02 cm,所有网格比热均为Cv=0.1 GJ/(cm3·keV),网格初始物质温度为1 eV,初始辐射温度为0 eV,时间步长为0.01 ns,总时长为50 ns,每步模拟粒子数为100000.光性薄的区域σ0取值为1,光性厚的区域σ0取值为1000.在光性厚区域用DDMC方法模拟,在光性薄区域用IMC方法模拟(其结果称为DDMC-IMC),同样与全区域用IMC方法模拟的结果相比较.

图7为10,20,50 ns时,平板中物质温度的空间分布.从图7中可知,DDMC-IMC与IMC方法的模拟结果很一致.而且,IMC方法模拟时间为16.9 h,DDMC方法模拟时间仅为1.6 h,相比IMC效率提高将近10倍.随时间演化,热辐射从左边边界源处逐渐向右传播,先将光性厚的物质加热,再加热光性薄的物质,同样与热辐射输运过程相符.

图7 不同时刻介质中物质温度的空间分布比较Fig.7.Material temperature in different moments.

实际上,DDMC-IMC混合输运时效率的提高还与模型中DDMC区域和IMC区域的粒子数比例有关,模拟过程中利用DDMC方法进行输运的粒子越多,DDMC-IMC混合输运的效率越高.在本问题中,模拟过程的大部分时间内,DDMC区域的物质温度都高于IMC区域,DDMC区域分配到的粒子比IMC区域多,因此加速效果明显.

5 结 论

虽然IMC方法可以较准确地模拟热辐射在介质材料中输运的物理过程,但在介质吸收系数较大的情况下,粒子(伪)散射过程占据主导,导致IMC方法模拟时间过长,效率过低.而在这种情况下,粒子输运过程可以很好地由扩散近似代替,从而缩短跟踪粒子的时间,提高模拟效率.

本文研究了DDMC辐射模拟方法及流程,提出了DDMC-IMC界面新的处理方法,新研制了“离散扩散蒙特卡罗程序”,并验证了该程序的适用性、准确性及高效性.结果表明:1)在吸收截面较大的单一介质系统中,DDMC方法可以得出与IMC方法一致的结果,而效率却高于IMC方法,尤其随着吸收系数变大,DDMC方法的效率提高愈加明显;2)在同时存在轻重介质的系统中,DDMC-IMC耦合输运的结果与IMC方法结果符合得很好,且效率也同样高于IMC方法,而且利用DDMC方法输运的粒子所占比例越大,DDMC-IMC混合输运的效率越高.综上所述,本文DDMC方法的程序可以高效、准确地模拟热辐射在介质中的输运过程,并且可以很好地与IMC方法耦合,高效模拟强吸收物质和弱吸收物质共存情况下的热辐射输运问题.

在DDMC-IMC耦合输运中,目前DDMCIMC界面是预先设定好的,若DDMC-IMC界面能根据物质吸收系数变化情况自动调整,模拟效率可以得到进一步提高.因此下一步将开展可变界面的研究.