輻射流體力學數值模擬中的隱式蒙特卡羅方法

閆 凱

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

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

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

1 輻射流體力學方程組

不考慮黏性和體積力，實驗室坐標系下的流體力學方程可寫為：

(1)

(2)

(3)

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

(4)

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

Fν=∮nIνdn

(5)

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

(6)

(7)

(8)

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

(9)

四元組表達式可簡化為：

(10)

(11)

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

2 數值方法

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

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

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

(12)

(13)

其中：I為輻射強度；μ為光子方向與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)

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

(20)

化簡后得到：

(21)

定義Fleck因子f為：

(22)

則式(21)可寫為：

(23)

t∈[tn,tn+1]

(24)

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

(25)

(26)

利用內能和平衡輻射能的轉化關系，并對能量平衡方程關于時間積分，最終得到內能的時間離散形式為：

(27)

(28)

2.2 輻射流體力學問題中的IMC方法

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

(29)

(30)

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

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

(31)

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

ε0≈fχ0PcurΔVΔt

(32)

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

(33)

每個蒙特卡羅粒子代表一束具有相同位置、發射時間和頻率的光子，因此粒子的能量權在兩種坐標系下的轉化關系為：

(34)

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

(35)

(36)

其中，下標i表示單元上的第i個粒子。

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

(37)

其中，ξ3為0～1之間的隨機數。即粒子飛行ds的距離后發生1次碰撞，在實驗室坐標系下，介質的不透明度與粒子運動方向相關，兩種坐標系下的轉化關系為：

χ0ν0=χν

(38)

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

(39)

(40)

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

(41)

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

(42)

(43)

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

ρ(un+1-un)=Δptot

(44)

(45)

最后通過狀態方程更新流體的溫度和壓力。

3 數值結果

3.1 輻射平衡問題

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

(46)

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

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

3.2 穩態輻射波

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

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

3.3 非穩態輻射波

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

圖2 物質密度、速度、溫度及輻射溫度的穩態解Fig.2 Steady state solution of material density, velocity, temperature and radiation temperature

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

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

4 結論

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