輻射流體力學(xué)數(shù)值模擬中的隱式蒙特卡羅方法

閆 凱

(1.北京大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，北京 100084；2.西北核技術(shù)研究院，陜西 西安 710024)

輻射流體力學(xué)廣泛應(yīng)用于天體物理、慣性約束聚變和武器物理等領(lǐng)域中，這些領(lǐng)域均涉及到高能量密度物理，輻射是主要的能量和動量傳輸方式之一，由于實驗環(huán)境一般較極端，數(shù)值模擬是重要的研究手段之一。輻射流體力學(xué)的數(shù)值方法大致可分為確定論方法和隨機(jī)模擬方法。在確定論方法中，由于輻射強(qiáng)度是時間、空間、角度和頻率的函數(shù)，在數(shù)值計算中一般要對輻射輸運方程做一定程度的近似，進(jìn)而得到不同輻射輸運近似模型。按照對頻率離散方式的不同，輻射輸運近似模型可分為灰體近似、多群近似和多帶近似等；按照對角度離散方式的不同，輻射輸運近似模型包括矩近似、Sn近似和Pn近似等[1]。矩近似是對輻射輸運方程進(jìn)行角度積分，然后求解相應(yīng)的矩方程組。目前應(yīng)用較多的是0階矩近似(經(jīng)典擴(kuò)散近似和限流擴(kuò)散近似)和1階矩近似。但矩近似方程組并不封閉，因此需給出封閉條件，而這種封閉條件往往是非物理的[2]。Sn近似將粒子運動方向近似為某些特定的方向，然后在每條方向上分別進(jìn)行求解，但在高維問題中Sn近似方法受到射線效應(yīng)的困擾[3]。Pn近似是將輻射強(qiáng)度按照角分布展開為一個基本函數(shù)全集，但Pn近似在數(shù)值模擬中受到波效應(yīng)的影響，應(yīng)用相對狹窄[1]?？傮w來說，確定論方法計算效率相對較高，但由于對輻射輸運方程做了某種程度的近似，因此在復(fù)雜輻射環(huán)境下其計算結(jié)果往往需驗證。

隨機(jī)模擬方法又稱為蒙特卡羅方法，是一種基于概率統(tǒng)計的數(shù)值模擬方法。在20世紀(jì)40年代，蒙特卡羅方法就被應(yīng)用到數(shù)值計算中，但受限于計算機(jī)的計算能力，直到1971年Fleck等[4]提出隱式蒙特卡羅(IMC)方法后，蒙特卡羅方法才逐漸大規(guī)模地應(yīng)用到熱輻射輸運問題中[5-7]。近十幾年，輻射流體力學(xué)的蒙特卡羅方法發(fā)展較迅速，Nayakshin等[8]采用蒙特卡羅方法模擬輻射輸運過程，通過統(tǒng)計得到輻射壓力，用于光滑粒子流體動力學(xué)(SPH)的計算。Harries等[9-11]討論其開發(fā)的輻射輸運蒙特卡羅程序TORUS，在后續(xù)的工作中，該程序用來模擬輻射流體力學(xué)問題[12-14]。Noebauer等[15]給出了一種輻射與物質(zhì)完全耦合的蒙特卡羅方法，其在輻射輸運過程中提出對光子的能量沉積和動量沉積進(jìn)行統(tǒng)計，該方法能保證守恒性。Roth等[16]首次給出了輻射流體力學(xué)的IMC方法，但對IMC方法的實施細(xì)節(jié)討論較少。另外，其在蒙特卡羅模擬中統(tǒng)計的物理量為輻射能密度和輻射能流，然后通過一個關(guān)于時間的向前歐拉格式來更新流體動量和能量，但這種做法不能保證守恒性。本工作結(jié)合Noebauer和Roth數(shù)值方法中的優(yōu)點，采用IMC方法模擬輻射輸運問題，在模擬過程中通過統(tǒng)計能量和動量沉積來更新流體能量和動量。

1 輻射流體力學(xué)方程組

不考慮黏性和體積力，實驗室坐標(biāo)系下的流體力學(xué)方程可寫為：

(1)

(2)

(3)

省略相關(guān)物理量的空間和時間變量，實驗室坐標(biāo)系下輻射輸運方程為：

(4)

其中：Iν為單頻輻射強(qiáng)度；n為光子運動方向的單位向量；ni為n在i方向的分量；ην為單頻發(fā)射率；χν為單頻總不透明度。為方便介紹，引入單頻輻射能密度Eν、單頻輻射能流Fν和單頻輻射壓力pν，其表達(dá)式分別為：

Fν=∮nIνdn

(5)

由于輻射與物質(zhì)相互作用的復(fù)雜性，為方便計算，做如下假設(shè)：1) 物質(zhì)滿足局部熱力學(xué)平衡條件；2) 在流體靜止系中散射是彈性碰撞；3) 散射碰撞后的方向在流體靜止系下是各向同性的。為區(qū)分兩種坐標(biāo)系，用下標(biāo)0表示流體靜止系，若無下標(biāo)0則表示實驗室坐標(biāo)系下的物理量。流體靜止系下不透明度和發(fā)射率可表示為：

(6)

(7)

(8)

給出3種常用的平均不透明度系數(shù)，分別為平均輻射能密度不透明度χ0E、平均普朗克不透明度χ0P和平均輻射能流不透明度χ0F，其表達(dá)式分別為：

(9)

四元組表達(dá)式可簡化為：

(10)

(11)

其中，γ=(1-uiui/c2)-1/2為洛倫茲因子。

2 數(shù)值方法

采用算子分裂思想將計算分為3步：第1步計算不含輻射與物質(zhì)相互作用項的流體力學(xué)方程組；第2步計算輻射輸運方程，采用IMC方法進(jìn)行計算；第3步通過統(tǒng)計沉積能量和動量來更新流體能量和動量。流體力學(xué)的數(shù)值方法已相對成熟，這里主要介紹考慮相對論效應(yīng)的IMC方法。

2.1 熱輻射輸運問題的IMC方法

IMC方法最初應(yīng)用于熱輻射輸運問題，即只考慮輻射與物質(zhì)能量的耦合，而不考慮流體的運動過程。以一維灰體為例，不考慮散射情況下，即χ=χt，熱輻射輸運方程組可表示為：

(12)

(13)

其中：I為輻射強(qiáng)度；μ為光子方向與x軸正方向夾角的余弦；S為獨立外源項。

引入平衡輻射能密度ur為：

ur=aT4

(14)

定義β為：

(15)

在每個時間層[tn,tn+1]中用tn時刻的值來近似β、χ和S，即：

β(t)≈βn

χ(t)≈χn

S(t)≈Snt∈[tn,tn + 1]

(16)

則熱輻射輸運方程為：

(17)

(18)

(19)

其中，α為給定的常數(shù)，α∈[0.5,1]。對式(18)關(guān)于時間[tn,tn+1]積分，可得到：

(20)

化簡后得到：

(21)

定義Fleck因子f為：

(22)

則式(21)可寫為：

(23)

t∈[tn,tn+1]

(24)

將式(24)代入到式(17)、(18)可得到：

(25)

(26)

利用內(nèi)能和平衡輻射能的轉(zhuǎn)化關(guān)系，并對能量平衡方程關(guān)于時間積分，最終得到內(nèi)能的時間離散形式為：

(27)

(28)

2.2 輻射流體力學(xué)問題中的IMC方法

在輻射流體力學(xué)問題中，當(dāng)流體速度較大時，經(jīng)典IMC方法不再適用，本文推導(dǎo)了考慮相對論效應(yīng)的IMC方法。其中時間和空間的離散均在實驗室坐標(biāo)系下進(jìn)行，粒子的運動也是在實驗室坐標(biāo)系下的時空中，但輻射與物質(zhì)的相互作用即蒙特卡羅粒子的碰撞事件發(fā)生在局部的(每個單元)流體靜止系下。通常所說的不透明度均是指其在流體靜止系下的值，散射是各向同性的假設(shè)也是基于流體靜止系的。兩種坐標(biāo)系下光子頻率和角度滿足的關(guān)系為：

(29)

(30)

在流體靜止系下，假設(shè)滿足局部熱力學(xué)平衡，根據(jù)霍爾基夫定律，每個單元(體積為ΔV0)在單位時間(Δt0)內(nèi)發(fā)射的能量ε0為：

ε0=χ0Pcu0rΔV0Δt0≈χ0PcurΔVΔt

(31)

在IMC方法中，吸收系數(shù)被有效吸收系數(shù)代替，則：

ε0≈fχ0PcurΔVΔt

(32)

假設(shè)流體靜止系下蒙特卡羅粒子的能量權(quán)為w0，其中w0是任意給定的，則該單元上蒙特卡羅粒子的個數(shù)Nemit為：

(33)

每個蒙特卡羅粒子代表一束具有相同位置、發(fā)射時間和頻率的光子，因此粒子的能量權(quán)在兩種坐標(biāo)系下的轉(zhuǎn)化關(guān)系為：

(34)

兩種坐標(biāo)系下蒙特卡羅粒子的數(shù)目是一致的，通過式(29)、(30)可得到蒙特卡羅粒子在實驗室坐標(biāo)系下的角度、頻率及能量權(quán)。蒙特卡羅粒子的空間位置和發(fā)射時間的抽樣與經(jīng)典IMC方法一致。抽樣給出熱發(fā)射源的粒子狀態(tài)后，需統(tǒng)計實驗室坐標(biāo)系下因熱發(fā)射而產(chǎn)生的每個單元上的物質(zhì)能量損失ΔE′m和動量損失Δp′，則：

(35)

(36)

其中，下標(biāo)i表示單元上的第i個粒子。

在蒙特卡羅粒子的追蹤過程中，需抽樣給出碰撞距離。對于經(jīng)典IMC方法，若介質(zhì)的不透明度為χ0，則粒子飛行的軌跡長度ds為：

(37)

其中，ξ3為0～1之間的隨機(jī)數(shù)。即粒子飛行ds的距離后發(fā)生1次碰撞，在實驗室坐標(biāo)系下，介質(zhì)的不透明度與粒子運動方向相關(guān)，兩種坐標(biāo)系下的轉(zhuǎn)化關(guān)系為：

χ0ν0=χν

(38)

聯(lián)合式(29)、(38)可給出實驗室坐標(biāo)系下的不透明度，碰撞距離抽樣仍可仿照式(37)，將流體靜止系下的不透明度替換為實驗室坐標(biāo)系下的不透明度即可。在追蹤粒子的運動過程中，需統(tǒng)計粒子的能量和動量沉積。碰撞一般分為吸收和散射。對于吸收事件，粒子飛行距離ds后，其能量沉積密度ΔEm和動量沉積密度Δp分別為：

(39)

(40)

其中，上標(biāo)p、f分別表示碰撞前、后的物理量。

(41)

在此過程中，能量沉積密度和動量沉積密度分別為：

(42)

(43)

結(jié)合式(35)～(36)和式(39)～(43)，可給出1個時間步總的能量沉積密度ΔEm,tot和動量沉積密度Δptot。更新流體物理量，在此過程中流體的密度保持不變，則：

ρ(un+1-un)=Δptot

(44)

(45)

最后通過狀態(tài)方程更新流體的溫度和壓力。

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

3.1 輻射平衡問題

測試氣體在輻射作用下的加熱和冷卻過程，在此過程中可驗證輻射輸運部分蒙特卡羅程序的正確性。本文選擇Turner&Stone算例[17]，計算區(qū)域氣體始終處于靜止?fàn)顟B(tài)，即不考慮相對論效應(yīng)的影響，其熱容為9×10-7J·cm-3·K。區(qū)域內(nèi)存在一均勻的、各向同性的輻射場，輻射能密度為105J/cm3，對應(yīng)的輻射溫度為3.4×106K。忽略輻射對物質(zhì)動量的影響，這種情況下流體的能量方程為：

(46)

本文考慮兩種情況：1) 氣體被輻射加熱，取初始時刻的氣體溫度為11 K；2) 氣體冷卻過程，取初始時刻氣體的溫度為1.1×109K。兩種情況下輻射能密度均遠(yuǎn)大于氣體能量密度，這意味著能量交換的過程中輻射能密度是近似不變的，因此氣體溫度最終接近于輻射溫度。圖1為IMC方法計算結(jié)果與解析解結(jié)果的比較，其中時間步長為10-15s，蒙特卡羅粒子數(shù)取為10 000，可看出，兩者結(jié)果吻合良好。

圖1 IMC方法計算結(jié)果與解析解的比較Fig.1 Comparison of IMC calculation and analytic solution results

3.2 穩(wěn)態(tài)輻射波

Lowrie等[18]給出了經(jīng)典擴(kuò)散近似下一維輻射流體力學(xué)方程組關(guān)于穩(wěn)態(tài)波結(jié)構(gòu)的半解析解結(jié)果，在其分析中，方程組的穩(wěn)態(tài)解可由4個物理量參數(shù)決定：上游氣體的馬赫數(shù)M，光速和上游氣體的聲速的比率C，上游輻射壓力(乘以3)與上游氣體壓力的比率P0，光學(xué)厚度τ。

根據(jù)文獻(xiàn)[18-19]，取P0=10-4、C=1.732×103、τ=577。初始時刻上游氣體處于輻射平衡狀態(tài)，溫度和密度均為1，根據(jù)狀態(tài)方程可知此時的聲速為1，上游氣體速度取為馬赫數(shù)，下游氣體的狀態(tài)根據(jù)間斷跳躍條件計算得到，IMC方法基于式(1)～(4)進(jìn)行模擬。圖2為M=2時物質(zhì)密度、速度、溫度及輻射溫度的穩(wěn)態(tài)解，可看出，計算結(jié)果與半解析解結(jié)果基本一致，初步驗證了程序的正確性。

3.3 非穩(wěn)態(tài)輻射波

算例來自文獻(xiàn)[20]，文獻(xiàn)[15-16，21-22]先后對此算例進(jìn)行了測試，并對模型進(jìn)行了簡化，本文采用文獻(xiàn)[22]中的簡化模型。一維計算區(qū)域長度為8.7×108m，左邊界為反射邊界，右邊界為自由邊界，計算區(qū)域內(nèi)氣體為氫原子，其初始溫度為50 K，密度為7.78×10-5kg/m3，絕熱系數(shù)Γ=1.4，不透明度為3.1×10-8m-1，初始時刻氣體以-6×103m/s或-2×104m/s的速度向左面運動，兩種速度下計算得到的輻射波分別為亞臨界輻射波和超臨界輻射波。圖3分別為亞臨界輻射波和超臨界輻射波中物質(zhì)溫度和輻射溫度的分布[16]，可看出，兩者結(jié)果基本一致，但在輻射前驅(qū)波區(qū)域，本文計算結(jié)果偏大，這可能是由于輻射輸運模型的差異導(dǎo)致的。圖4為未考慮和考慮相對論效應(yīng)IMC方法的計算結(jié)果比較，可看出，在早期兩者的結(jié)果基本一致，但當(dāng)計算到1.3×104s時，兩者結(jié)果出現(xiàn)了較為明顯的差異，考慮相對論效應(yīng)的IMC方法物質(zhì)溫度升得更快，盡管物質(zhì)速度(-2×104m/s)僅約為光速的萬分之一，但對計算結(jié)果已有明顯的影響。

圖2 物質(zhì)密度、速度、溫度及輻射溫度的穩(wěn)態(tài)解Fig.2 Steady state solution of material density, velocity, temperature and radiation temperature

圖3 亞臨界(a)和超臨界(b)輻射波中物質(zhì)溫度和輻射溫度分布Fig.3 Material and radiation temperature distributions under subcritical shock (a) and supercritical shock (b)

圖4 兩種IMC方法計算結(jié)果的比較Fig.4 Comparison of result by two IMC methods

4 結(jié)論

本文研究了輻射流體力學(xué)數(shù)值模擬中的IMC方法。介紹了在不考慮流體運動過程的IMC方法，給出了該方法在輻射流體力學(xué)問題中詳細(xì)執(zhí)行過程。與熱輻射輸運問題相比，IMC方法需考慮輻射流體力學(xué)的相對論效應(yīng)。在模擬過程中給出了粒子的能量和動量沉積的統(tǒng)計方法，這種處理能保證守恒性。最后給出了幾個典型問題的數(shù)值模擬結(jié)果，并與相關(guān)的參考解進(jìn)行對比，對比結(jié)果驗證了算法和程序的正確性。