熱輻射輸運問題的高效蒙特卡羅模擬方法*

許育培 李樹

1) (北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100094)

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

慣性約束聚變研究中,熱輻射光子在介質(zhì)中的輸運以及熱輻射光子與介質(zhì)的相互作用是重要研究課題,蒙特卡羅方法是該類問題的重要研究手段之一.隱式蒙特卡羅方法雖然能正確地模擬熱輻射在介質(zhì)中的輸運過程,但當模擬重介質(zhì)(材料的吸收系數(shù)大)問題時,該方法花費的計算時間將變得很長,導(dǎo)致模擬效率很低.本文以離散擴散蒙特卡羅方法為基礎(chǔ),開發(fā)了“離散擴散蒙特卡羅方法輻射輸運模擬程序”,可以較好地解決重介質(zhì)區(qū)的計算效率問題,但是離散擴散蒙卡羅方法在模擬輕介質(zhì)區(qū)時精度不夠高.輻射輸運問題中通常既有輕介質(zhì)也有重介質(zhì),為了能同時解決蒙特卡羅方法模擬的效率和精度問題,本文研究了離散擴散蒙特卡羅方法與隱式蒙特卡羅方法相結(jié)合的模擬方法,并提出了新的擴散區(qū)與輸運區(qū)界面處理方法,研制了混合蒙特卡羅方法的輻射輸運模擬程序.典型輻射輸運問題模擬顯示:在模擬重介質(zhì)問題時,該程序能大幅縮短模擬時間,且能取得與隱式蒙特卡羅方法一致的結(jié)果;在模擬輕重介質(zhì)均存在的問題時,與隱式蒙特卡羅方法相比,混合蒙特卡羅方法的模擬精度與其相當且計算效率同樣能夠得到顯著提升.

1 引 言

慣性約束聚變(inertial confinement fusion,ICF)中,熱輻射(X光波段)扮演著重要的角色.間接驅(qū)動ICF中,黑腔壁產(chǎn)生的熱輻射光子燒蝕驅(qū)動飛層對熱核材料進行壓縮,產(chǎn)生高溫高密度的等離子體,實現(xiàn)聚變點火.因此,熱輻射光子在介質(zhì)(材料)中的輸運以及熱輻射光子與介質(zhì)的相互作用是ICF的重要研究課題之一.

在輻射輸運過程中,輻射光子與物質(zhì)的相互作用很復(fù)雜,輻射強度和物質(zhì)溫度強耦合,描述輻射場時空演化的輻射輸運方程組是非線性的[1],只能用數(shù)值方法求解實際的熱輻射輸運問題.蒙特卡羅(Monte Carlo,MC)[2,3]方法是常用的求解粒子輸運問題的數(shù)值方法.傳統(tǒng)的顯式蒙特卡羅方法基于與時間相關(guān)的、非線性的輻射輸運方程,通過隨機抽樣來確定粒子的吸收、散射與發(fā)射,以迭代的方式模擬熱輻射輸運過程[4].然而,當物質(zhì)的吸收系數(shù)較大或系統(tǒng)接近熱力學(xué)平衡態(tài)時,傳統(tǒng)的顯式蒙特卡羅法的模擬結(jié)果存在嚴重的漲落和能量不均衡等問題[5].

為解決傳統(tǒng)蒙特卡羅方法存在的問題,Fleck和Cummings[6]提出了隱式蒙特卡羅(implicit Monte Carlo,IMC)方法.該方法天然具有無條件穩(wěn)定性,較傳統(tǒng)的顯式蒙特卡羅方法更穩(wěn)定,精度、效率更高.然而,當系統(tǒng)的吸收系數(shù)變得更大時,IMC方法會出現(xiàn)計算時間過長的問題,導(dǎo)致模擬效率低下.為了提高模擬效率,Fleck和Canfield[7]基于概率方法首先提出“隨機游走法”(randomwalk,RW),Giorla和Sentis[8]基于擴散方程推導(dǎo)了RW方程,建立了以粒子所在位置為球心的球形子區(qū)域以及判據(jù),當判據(jù)滿足時,球形子區(qū)域內(nèi)粒子的多次散射過程將由一個擴散過程代替,從而有效提高了模擬效率.然而,球形子區(qū)域的半徑有最小值(其半徑不能小于數(shù)倍粒子平均自由程),當粒子接近邊界時,球形子區(qū)域半徑與粒子平均自由程的關(guān)系將不再滿足,RW方法將失效,只能轉(zhuǎn)而用IMC方法輸運.因此,若吸收系數(shù)不夠大,或系統(tǒng)網(wǎng)格劃分較小,RW模擬輻射輸運的效率提高則不明顯.20世紀90年代,Urbatsch等[9]首次提出離散擴散蒙特卡羅(discrete diffusion Monte Carlo,DDMC)方法.DDMC方法從離散的擴散方程出發(fā),把IMC方法中的多次散射過程替換成一個擴散過程,在保證相同精度的前提下,很大程度上提高了輻射輸運模擬的效率.之后,Evans等[10]推導(dǎo)出平衡態(tài)下的DDMC方法;Gentile[11]對擴散方程進行空間和時間離散,得到矩陣方程,推導(dǎo)出一維的時間、空間離散DDMC方法;Densmore等[12,13]對擴散方程進行空間離散,進而推導(dǎo)出一維的時間連續(xù)、空間離散DDMC方法;Cleveland等[14]提出多群 DDMC (原作者稱其為implicit Monte Carlo diffusion,IMD)方法,以解決與頻率相關(guān)的熱輻射輸運問題;Densmore等[15]基于頻率積分擴散方程,開發(fā)了與頻率有關(guān)的DDMC新方法.

國外已有基于DDMC方法開發(fā)的數(shù)值模擬程序,并與IMC方法結(jié)合以求解熱輻射輸運問題,但由于問題的敏感性,國外DDMC方法模擬程序?qū)崿F(xiàn)的關(guān)鍵環(huán)節(jié)及程序代碼并未對外公開.經(jīng)過多年努力,國內(nèi)已開發(fā)出較成熟的IMC方法輻射輸運模擬程序[16],且已經(jīng)成功用于ICF方法實驗的模擬與分析[17],但是至今未有基于DDMC方法的輻射輸運模擬程序.因此有必要獨立開展DDMC方法研究并開發(fā)相應(yīng)的數(shù)值模擬程序,應(yīng)用于ICF方法相關(guān)課題研究中.本文在DDMC方法基礎(chǔ)上,解決了DDMC粒子(用DDMC輸運的粒子)輸運模擬的關(guān)鍵技術(shù),研制了DDMC方法輻射輸運模擬程序,研究了DDMC與IMC方法的耦合方法,并提出一種新的界面處理方法,進而研制了DDMC-IMC程序.該程序能準確且高效地模擬不同光性厚度介質(zhì)中的熱輻射輸運問題.

2 DDMC方法及數(shù)值模擬流程設(shè)計

下面以一維平板問題為例,介紹DDMC方法的主要思想并推導(dǎo)相應(yīng)的數(shù)學(xué)方程.考慮一維、均勻、無獨立外源的無限平板,熱輻射輸運方程及能量方程[6,18]為

其中,I為輻射強度,t為時間變量,c為真空中的光速,μ為粒子運動的方向余弦,x為空間變量,σ為介質(zhì)的吸收系數(shù),a為輻射常數(shù),T為介質(zhì)溫度,Cv為介質(zhì)比熱.

將平板網(wǎng)格化,即 0=x1＜x2＜x3＜··＜xI,x1,x2,··,xI為網(wǎng)格邊界所在的位置,每個網(wǎng)格中的位置變量滿足xi＜x＜xi+1;將時間離散化,即0=t1＜t2＜t3＜··tN,t1,t2,··,tN為各個時間步結(jié)束的時刻,每個時間步中的時間變量滿足tn＜t＜tn+1.應(yīng)用IMC方法后,方程(1)改寫為[6]

其中帶下標n的變量表示第n個時間步中對應(yīng)的變量值,σes為有效散射系數(shù),σea為有效吸收系數(shù),

其中,f為Fleck因子,Ur為輻射能量密度,Um為物質(zhì)能量密度,Δt為時間步長.

當介質(zhì)的吸收系數(shù)很大時,粒子的平均自由程將變得很小,由方程(5)可知,fn也將變小,即有效散射系數(shù)σes將接近σn,而有效吸收系數(shù)σea將接近0,因此系統(tǒng)中粒子的散射將占據(jù)主導(dǎo)地位.綜上,在吸收系數(shù)很大的情況下,IMC中粒子的歷史變得過于漫長,導(dǎo)致計算效率降低.

當散射過程占主導(dǎo)時,考慮對方程(2)進行擴散近似.方程(2a)等式兩邊對μ全角度積分得

其中,φ為粒子通量密度,F為輻射流.對方程(7)中等式左邊的空間導(dǎo)數(shù)項進行有限差分,

其中,Δx 為網(wǎng)格大小,Fi+1/2和Fi —1/2分別為第i個網(wǎng)格左、右邊界的平均輻射流.因此對于第i個網(wǎng)格,方程(7)改寫為

其中φi為網(wǎng)格平均輻射強度.

1)當 1＜i＜I時,應(yīng)用斐克定律[19,20],網(wǎng)格邊界平均輻射流的表達式為

其中φi+1/2為網(wǎng)格邊界平均輻射強度.

聯(lián)立方程 (10)和 (11),解出 Fi+1/2得

類似可得 Fi —1/2的表達式.

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

式中利用了 σR,i—1=σL,i的關(guān)系,另外 σL,i和 σR,i分別表示第i個網(wǎng)格中粒子向左、右相鄰網(wǎng)格擴散的截面,

方程(13)是DDMC方法中粒子輸運所依據(jù)的方程,其中等式左邊第二項為吸收項,等式右邊為源項,右邊第一項表示從右邊相鄰網(wǎng)格向左擴散到當前網(wǎng)格的粒子,右邊第二項表示從左邊相鄰網(wǎng)格向右擴散到當前網(wǎng)格的粒子,右邊第三項表示介質(zhì)輻射的粒子.

在DDMC方法中,粒子的位置變量不再有意義,因此定義粒子平均自由時間τ,以確定粒子在隨時間變化過程中發(fā)生的反應(yīng)(擴散、吸收或存活至當前時間步結(jié)束),

粒子平均自由時間與粒子平均自由程類似,表示粒子發(fā)生兩次反應(yīng)(擴散或吸收)的平均時間間隔.

2)對于系統(tǒng)邊界處的網(wǎng)格,即當i=1或i=I時,方程(13)不再適用,因此有必要推導(dǎo)邊界處粒子輸運滿足的離散擴散方程.令方程(9)中i=1,得

應(yīng)用漸進擴散極限(asymptotic diffusionlimit)邊界條件[21]:

其中,W(μ)為超越方程,可以近似為 W(μ) ≈ μ+3/2μ2;φ1/2為系統(tǒng)邊界平均輻射強度;λ為外推距離.

對方程(18)中的空間導(dǎo)數(shù)項進行有限差分,得

令方程(11)中i+1=1,i+1/2=1/2,并與方程(19)聯(lián)立,解出F1/2:

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

方程(21)即為i=1時邊界網(wǎng)格滿足的離散擴散方程,其中

同理可推導(dǎo)出當i=I時,邊界網(wǎng)格滿足的輸運方程.

除了右邊第二項外,方程(21)與方程(13)類似.方程(21)右邊第二項中PL(μ)表示粒子穿過邊界進入DDMC區(qū)域后而不被反射的概率,與粒子運動方向和邊界法向之間夾角的余弦μ有關(guān).

需要特別指出得是:當模擬的對象為輕介質(zhì)(吸收系數(shù)較小)時,擴散條件不成立,DDMC方法模擬的精度將變差;另外,當吸收截面較大時,上面的PL,σL,1與理論值相比偏小,也需對它們做適當?shù)男拚齕22].

本文設(shè)計了基于DDMC方法的輻射輸運問題模擬流程.主要步驟如下.

1) 時空網(wǎng)格劃分 根據(jù)系統(tǒng)的幾何結(jié)構(gòu)及介質(zhì)屬性,將系統(tǒng)劃分為若干空間網(wǎng)格,每個網(wǎng)格包含同一種介質(zhì);根據(jù)輻射源及輻射場演化特性,將要模擬的時間范圍離散成若干時間步.

2) 計算當前時間步的各類參數(shù) 根據(jù)各個網(wǎng)格的溫度、密度、材料種類等計算網(wǎng)格比熱、吸收系數(shù)、物質(zhì)能量密度、輻射能量密度、有效吸收系數(shù)、有效散射系數(shù)、fleck因子、粒子擴散截面、粒子轉(zhuǎn)化概率等.

3) 從輻射源中抽取粒子 根據(jù)輻射總能量及模擬的粒子數(shù),分配三種源(邊界源、網(wǎng)格物質(zhì)輻射源、上一時間步中存活下來的粒子)的粒子權(quán)重及粒子數(shù),依次抽取源對應(yīng)的粒子并確定其時間、能量、權(quán)重.

4) 粒子隨機擴散 根據(jù)計算得到的擴散截面和有效吸收系數(shù),粒子在各個網(wǎng)格間隨機運動,或被物質(zhì)吸收,或擴散至相鄰網(wǎng)格,或存活至當前時間步結(jié)束.若粒子未被吸收,則記錄粒子的相關(guān)參數(shù),在下一時間步繼續(xù)輸運.

5) 求解溫度方程 由步驟4)得到的輻射能和吸收能,以及網(wǎng)格的比熱,利用能量方程解出當前時間步結(jié)束時的物質(zhì)溫度.

6) 預(yù)估下一時間步的平均溫度 由于需要根據(jù)網(wǎng)格溫度來計算網(wǎng)格其他參數(shù),因此有必要預(yù)估下一時間步的溫度,本文采用“預(yù)估迭代”的方法,此方法既能保持較好的精度,又能相對節(jié)省時間.

7) 返回步驟2),開始下一時間步計算,直至結(jié)束時刻.

3 DDMC與IMC耦合方法

ICF輻射輸運問題中通常既包含有輕介質(zhì)(如推進層),也包含有重介質(zhì)(如黑腔壁).對于輕介質(zhì)(光性薄)區(qū)域,最好應(yīng)用IMC方法模擬,對于重介質(zhì)(光性厚)區(qū)域,則考慮應(yīng)用DDMC方法模擬.因此,需要研究DDMC方法模擬與IMC方法模擬的耦合方法.DDMC-IMC方法耦合輸運的關(guān)鍵是處理好IMC區(qū)域(應(yīng)用IMC方法模擬的區(qū)域)和DDMC區(qū)域(應(yīng)用DDMC方法模擬的區(qū)域)的交界面,即如何正確合理的處理MC粒子穿過交界面時的各種物理數(shù)學(xué)參數(shù),下面介紹其中的兩個關(guān)鍵環(huán)節(jié).

第一種情況,粒子從IMC區(qū)域穿過界面進入DDMC區(qū)域.IMC粒子(用IMC方法模擬的粒子)穿過界面時,有P(μ)的概率進入DDMC區(qū)域并轉(zhuǎn)化為DDMC粒子,沒有轉(zhuǎn)化的粒子在界面處被反射,并重新抽取方向和位置,位置在界面上均勻分布.

一般情況下,DDMC區(qū)域與IMC區(qū)域的不透明度相差很大,從DDMC區(qū)域邊界離開后進入IMC區(qū)域的粒子的運動方向不再是各向同性的.實際上其運動方向的角分布與面源粒子的角分布一樣,為余弦分布[23],即 f(μ)=2μ.因此抽取粒子新方向的方法為

其中ξ為在[0,1]上均勻分布的隨機數(shù).

第二種情況,粒子從DDMC區(qū)域穿過邊界進入IMC區(qū)域.DDMC粒子穿過界面后,轉(zhuǎn)化為IMC粒子.因為在DDMC區(qū)域中,粒子的方向和位置沒有實際意義,因此粒子進入IMC區(qū)域后同樣需要重新抽取位置和方向變量.與第一種情況一樣,粒子方向的抽取按照方程(24)的方法,位置在界面上均勻分布.

4 數(shù)值結(jié)果

根據(jù)上文所述的DDMC方法輻射輸運流程以及DDMC-IMC界面處理方法,本文設(shè)計、編寫了“離散擴散蒙特卡羅方法輻射輸運程序”(DDMC)及“離散擴散蒙特卡羅方法與隱式蒙特卡羅方法耦合輸運程序”(DDMC-IMC).為驗證程序的適用性和正確性,本文計算了若干典型問題,并與已有的IMC程序計算結(jié)果比較,計算平臺為個人電腦,CPU型號為Intel(R) Core(TM) i7-3770 @ 3.4 GHz,核數(shù)為8,內(nèi)存為4.00 GB,操作系統(tǒng)為Windows 7旗艦版,編譯器為Intel Visual Fortran.

4.1 單一介質(zhì)熱波傳播問題

模型為0.5 cm厚的一維平板,左側(cè)邊界為一溫度恒為1 keV的平面源,右側(cè)邊界為真空邊界,平板分為25個網(wǎng)格,網(wǎng)格大小為0.02 cm,網(wǎng)格比熱為cv=0.1 GJ/(cm3·keV),吸收系數(shù)與溫度的三次方呈反比,即 σ=σ0/T3,其中溫度T的單位取keV,σ0單位為keV3/cm,可取不同的值,越大表示光性越厚.后文中吸收系數(shù)及溫度單位與之相同.

網(wǎng)格初始物質(zhì)溫度為1 eV,初始輻射溫度為0 eV,時間步長取0.05 ns,總時長為10 ns,每步模擬粒子數(shù)為 10000.σ0分別取 200,500,1000,2000.分別用IMC和DDMC方法模擬并比較.開始進行輻射輸運時,光子從邊界源出射,與附近的物質(zhì)相互作用,加熱物質(zhì);被加熱的物質(zhì)重新輻射出光子,再加熱更遠處的物質(zhì),因此輻射得以從邊界源處逐漸向遠處傳播.

圖1和圖2分別為σ0=500時,三個時間點的物質(zhì)溫度和輻射溫度的空間分布.從兩圖中可以看出,DDMC與IMC方法的模擬結(jié)果總體上符合得很好,且隨時間的演化,曲線所描繪出的熱波逐漸向前傳播,加熱遠處的物質(zhì),與熱輻射傳播的物理過程相符.

圖1 不同時刻的物質(zhì)溫度空間分布比較(σ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方法的效率也隨之提高.因為吸收系數(shù)越大,IMC粒子在一個時間步內(nèi)散射的次數(shù)越多,而IMC方法中處理散射的過程占據(jù)了整個模擬過程的大部分時間,相反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 單一介質(zhì)熱平衡問題

模型為0.5 cm厚一維平板,無邊界源,且外邊界均為全反射邊界,平板分為50個網(wǎng)格,網(wǎng)格大小為 0.01 cm,網(wǎng)格比熱 Cv=0.1 GJ/(cm3·keV),網(wǎng)格初始物質(zhì)溫度與初始輻射溫度分別為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,每步模擬粒子數(shù)為10000.σ0取值為 200,500,1000,2000.同樣分別用IMC和DDMC方法模擬并比較.

由于外邊界設(shè)置為全反射邊界,故系統(tǒng)實為全封閉系統(tǒng).隨著時間演化,理論上該系統(tǒng)最終將達到平衡,即所有網(wǎng)格的輻射溫度和物質(zhì)溫度都將趨于某固定值.由于全封閉系統(tǒng)的能量守恒,可以建立如下關(guān)系:

圖3和圖4分別是σ0=500時,不同時間點的介質(zhì)中物質(zhì)溫度和輻射溫度的空間分布.從兩圖中可以看出,DDMC方法的模擬結(jié)果與IMC方法同樣符合得很好,隨著時間變化,溫度較高處的網(wǎng)格溫度逐漸降低,溫度較低處的網(wǎng)格溫度逐漸升高,且在40 ns時基本達到了理論平衡溫度0.8 keV.模擬結(jié)果顯示,當系統(tǒng)中存在溫度梯度時,能流將從溫度高處流向溫度低處,其間,輻射光子與物質(zhì)不斷相互作用(光子被吸收后再發(fā)射),系統(tǒng)中溫度分布、輻射場分布不斷發(fā)生變化,對于本例,輻射溫度與物質(zhì)溫度逐漸趨于理論平衡溫度,并在時間足夠長后達到平衡.

圖3 不同時刻的物質(zhì)溫度空間分布比較(σ0=500)Fig.3.Material temperature in different moments (σ0=500).

表2為不同σ0取值下DDMC和IMC方法的模擬時間以及效率的比較.與表1相似,隨著σ0的增大,IMC方法的模擬時間明顯變長,而DDMC方法的模擬時間略微縮短,相比較下的DDMC方法效率也隨之提高,最高接近30倍,其原因上文已做出解釋.同樣可得出物質(zhì)吸收系數(shù)越大,DDMC方法模擬效率越高的結(jié)論.

圖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中輻射溫度曲線的波動較大,尤其是平衡后各網(wǎng)格間的溫度看似并沒有達到平衡(溫度相等).初步分析其原因是蒙特卡羅方法計算的粒子數(shù)偏少、統(tǒng)計誤差偏大所致.為此,在其他參數(shù)不變的情況下,將模擬的粒子數(shù)提高至500000后計算,得到圖5和圖6所示的物質(zhì)溫度和輻射溫度的分布情況.與圖3和圖4對比可看出,數(shù)值波動有明顯的減小,說明較大的波動確實是由樣本數(shù)較少引起的,不影響本文對DDMC方法的優(yōu)點及模擬結(jié)果的分析.

圖5 不同時刻的物質(zhì)溫度空間分布比較(σ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 雙介質(zhì)中熱波傳播問題

模型為0.5 cm厚一維平板,左邊界為一溫度恒為1 keV的平面源,右邊界為真空邊界,平板分為左右兩部分,左邊部分為光性厚的區(qū)域,厚度為0.1 cm,劃分為20個網(wǎng)格,網(wǎng)格大小為0.005 cm,右邊部分為光性薄的區(qū)域,厚度為0.4 cm,劃分為20個網(wǎng)格,網(wǎng)格大小為0.02 cm,所有網(wǎng)格比熱均為Cv=0.1 GJ/(cm3·keV),網(wǎng)格初始物質(zhì)溫度為1 eV,初始輻射溫度為0 eV,時間步長為0.01 ns,總時長為50 ns,每步模擬粒子數(shù)為100000.光性薄的區(qū)域σ0取值為1,光性厚的區(qū)域σ0取值為1000.在光性厚區(qū)域用DDMC方法模擬,在光性薄區(qū)域用IMC方法模擬(其結(jié)果稱為DDMC-IMC),同樣與全區(qū)域用IMC方法模擬的結(jié)果相比較.

圖7為10,20,50 ns時,平板中物質(zhì)溫度的空間分布.從圖7中可知,DDMC-IMC與IMC方法的模擬結(jié)果很一致.而且,IMC方法模擬時間為16.9 h,DDMC方法模擬時間僅為1.6 h,相比IMC效率提高將近10倍.隨時間演化,熱輻射從左邊邊界源處逐漸向右傳播,先將光性厚的物質(zhì)加熱,再加熱光性薄的物質(zhì),同樣與熱輻射輸運過程相符.

圖7 不同時刻介質(zhì)中物質(zhì)溫度的空間分布比較Fig.7.Material temperature in different moments.

實際上,DDMC-IMC混合輸運時效率的提高還與模型中DDMC區(qū)域和IMC區(qū)域的粒子數(shù)比例有關(guān),模擬過程中利用DDMC方法進行輸運的粒子越多,DDMC-IMC混合輸運的效率越高.在本問題中,模擬過程的大部分時間內(nèi),DDMC區(qū)域的物質(zhì)溫度都高于IMC區(qū)域,DDMC區(qū)域分配到的粒子比IMC區(qū)域多,因此加速效果明顯.

5 結(jié) 論

雖然IMC方法可以較準確地模擬熱輻射在介質(zhì)材料中輸運的物理過程,但在介質(zhì)吸收系數(shù)較大的情況下,粒子(偽)散射過程占據(jù)主導(dǎo),導(dǎo)致IMC方法模擬時間過長,效率過低.而在這種情況下,粒子輸運過程可以很好地由擴散近似代替,從而縮短跟蹤粒子的時間,提高模擬效率.

本文研究了DDMC輻射模擬方法及流程,提出了DDMC-IMC界面新的處理方法,新研制了“離散擴散蒙特卡羅程序”,并驗證了該程序的適用性、準確性及高效性.結(jié)果表明:1)在吸收截面較大的單一介質(zhì)系統(tǒng)中,DDMC方法可以得出與IMC方法一致的結(jié)果,而效率卻高于IMC方法,尤其隨著吸收系數(shù)變大,DDMC方法的效率提高愈加明顯;2)在同時存在輕重介質(zhì)的系統(tǒng)中,DDMC-IMC耦合輸運的結(jié)果與IMC方法結(jié)果符合得很好,且效率也同樣高于IMC方法,而且利用DDMC方法輸運的粒子所占比例越大,DDMC-IMC混合輸運的效率越高.綜上所述,本文DDMC方法的程序可以高效、準確地模擬熱輻射在介質(zhì)中的輸運過程,并且可以很好地與IMC方法耦合,高效模擬強吸收物質(zhì)和弱吸收物質(zhì)共存情況下的熱輻射輸運問題.

在DDMC-IMC耦合輸運中,目前DDMCIMC界面是預(yù)先設(shè)定好的,若DDMC-IMC界面能根據(jù)物質(zhì)吸收系數(shù)變化情況自動調(diào)整,模擬效率可以得到進一步提高.因此下一步將開展可變界面的研究.