球幾何中輻射源粒子抽樣方法的改進*

許育培 李樹

1) (北京應用物理與計算數學研究所， 北京 100094)2) (中國工程物理研究院研究生院， 北京 100088)(2020 年1 月5日收到; 2020 年4 月17日收到修改稿)

利用隱式蒙特卡羅方法模擬熱輻射輸運問題時， 輻射源粒子的抽樣對結果的正確性很重要. 在球幾何熱輻射輸運隱式蒙特卡羅模擬中， 通常的做法是假設單個網格(球殼)的溫度不隨空間變化， 即源粒子在球殼內均勻分布. 對于網格內部溫度分布梯度不大的情況， 這種處理不會造成太大誤差， 然而當物質的吸收系數很大或球殼較厚時， 那么單個網格內溫度的空間變化也可能較大， 這種處理將導致模擬的熱輻射傳播速度比實際的要快. 本文分析了導致這種偏差的原因， 并基于輻射能量密度分布， 推導了球幾何下輻射源粒子發射位置的概率密度函數， 提出了兩種輻射源粒子空間抽樣方法， 設計了抽樣流程， 計算了典型的熱輻射輸運問題.數值計算表明， 這兩種新的源粒子空間抽樣方法均能修正輻射傳播速度與參考值的偏差， 解決計算結果過分依賴于網格剖分的問題. 同時， 在網格數較少、模擬粒子數較少的情況下新方法仍能得到較準確的結果， 有效提高了計算效率.

1 引 言

激光間接驅動的慣性約束聚變(inertial confinement fusion， ICF)中， 黑腔的溫度可達上百萬度， 腔壁熱輻射(簡稱輻射)光子(X光)均勻地照射裝有DT聚變材料的靶丸， 驅動飛層壓縮靶丸， 實現聚變點火. 因此， 研究輻射在物質中的輸運及輻射與物質的相互作用是ICF的重要研究方向. 熱輻射的傳播行為可用輻射輸運方程來描述，理論上可通過求解輻射輸運方程及與之相耦合的物質能量方程來獲得輻射場及物質溫度場的時空演化規律. 由于輻射強度與物質溫度的強耦合性以及輻射輸運方程的非線性特性[1]， 即便數值方法求解也是非常困難的， 蒙特卡羅(Monte Carlo， MC)方法是求解輻射輸運問題的重要方法之一[2，3].

傳統的直接MC方法于20世紀六十年代由Fleck[4]提出， 但由于在很多情況下通過直接MC方法求得的結果存在過大的漲落和能量不均衡等缺點[4，5]， 因此很少被采用. 20世紀七十年代初， 加利福尼亞大學的Fleck和Cummings[6]提出了隱式蒙特卡羅(implicit Monte Carlo， IMC)方法， 該方法的主要思想是將物質的真實吸收截面由有效吸收截面和有效散射截面代替， 并引入與隱式迭代類似的手段， 故具有天然的穩定性. IMC方法與直接MC方法相比更穩定， 計算精度、效率更高，因此得到了廣泛的應用， 國外的ICF數值模擬程序[7]及其他輻射輸運程序[8-10]均用到了IMC方法解決其中的輻射輸運問題， 部分文獻[11-15]將IMC方法的模擬結果當作其他數值方法的參考結果. 國內已有學者研究并開發了IMC輻射輸運模擬程序[16]， 該程序可以進行一維或三維的輻射輸運數值模擬， 目前已經應用于ICF中黑腔物理的研究[17，18].

Fleck和Cummings[6]推導IMC輻射輸運方程時， 認為“某一時間步內的某網格中溫度是不隨空間變化的”， 即對單一網格做了“空間等溫假設”，因此網格內的輻射粒子的發射位置可由空間均勻抽樣得到， 例如在一維平板幾何中， 輻射粒子的發射位置是網格左右邊界之間的均勻分布函數， 本文將其稱為“等溫法”. 國外的三維輻射輸運模擬程序Milagro[7]、國內的半隨機模擬方法[19]以及其他基于IMC的輻射輸運模擬方法[6，11]都采用該抽樣方法. 在物質的吸收系數不大或空間網格較小的情況下， 這種處理方法不會給計算結果帶來明顯偏差， 然而對于強吸收、大網格問題， 繼續采用這種基于“等溫假設”的空間均勻抽樣法將導致較大誤差. 為此， 文獻[20]提出了基于輻射能量密度分布的新抽樣法， 推導了輻射能量空間分布密度函數和輻射源粒子空間抽樣公式， 解決了正交幾何下， 同一網格中溫差太大導致IMC模擬結果與實際結果不符的問題. 然而該方法并不適用于球幾何情況，因此本文將推導球幾何情況下輻射能量的空間分布密度函數， 并提出相應的輻射源粒子抽樣方法.

本文的主要內容如下， 第2節推導了球幾何情況下“等溫法”的輻射源粒子空間概率密度分布函數及抽樣方法， 探討了該方法的不足， 并通過數值手段獲得了一個典型球對稱熱輻射輸運問題的參考解. 第3節推導了基于能量密度分布的輻射源空間概率密度函數， 分析了IMC模擬中熱輻射波傳播偏快的原因， 提出了兩種新的輻射源粒子抽樣方法， 并設計了抽樣流程. 第4節通過典型算例測試了本文提出的兩種輻射源粒子抽樣方法的正確性，驗證了本項研究工作對球幾何情況下的IMC輻射源粒子抽樣精度及計算效率有顯著的改進效果.

2 等溫法的輻射源粒子抽樣及其缺點

IMC輻射輸運模擬中的輻射源粒子可分成四種[16，21]: 來自獨立外源的粒子; 從網格邊界進入的粒子; 從上一時間步存活下來的粒子; 從物質輻射出的粒子(以下稱為“輻射源粒子”). 前三種粒子均有確定的位置參數， 而輻射源粒子在當前時間步開始時需通過隨機抽樣確定其初始位置.

在一維球幾何情況下， IMC輸運方程與物質能量方程可寫為:

其中， 所有帶下標n的變量均表示第n個離散時間步中對應的變量值， I為輻射強度， c為真空中的光速， t為時間， r為粒子所在位置的半徑， μ為方向余弦， σ為物質的吸收截面， b為歸一化Planck函數， f為Fleck因子， Uγ為輻射能量密度， ζ為局域再發射系數， ν為光子頻率， Q為獨立外源， Cv為比熱容， T為物質溫度， σP為Planck平均吸收截面.

(1)式等號右邊第一項為相空間(r， μ， ν， t)處的輻射源粒子發射密度S(r， μ， ν， t)，

其中

式中， a為輻射常數，Tn(r) 是網格溫度.

在 r 處dr范圍內物質輻射出的能量 En(r) 為

在球幾何情況下，.

IMC方法中， 輻射源粒子數目與其輻射的能量 En(r) 呈正比， 因此， 某網格內的輻射源粒子關于 半徑 r 的空間分布概率密度函數可寫為

其中r1， r2分別是網格的內、外半徑.

因此在等溫假設下的輻射源粒子的位置變量 r可 用反函數法直接抽樣得到，

式中ξ為在[0， 1]上均勻分布的隨機數. 同一網格中輻射源粒子初始位置的半徑 r 均可由(9)式抽樣獲得， 這實際上等同于球殼內的空間均勻抽樣.

在物質吸收系數小或網格較小(球殼較薄)的情況下， 網格內部的溫度梯度不大， 即溫度隨空間(半徑r)的變化不顯著，近似引入的偏差不大， 故這種等溫法抽樣對計算結果影響較小.然而當物質的吸收系數很大或網格較大時， 網格內部溫度差可能會比較大， 若用網格平均溫度 Tn取代方程(7)中的 Tn(r) ， 則計算結果將產生較大誤差. 因此， 在等溫法中， 為了避免網格內部溫差大所導致的問題， 就必須盡量細分網格， 使網格內溫度差盡可能小， 令等溫假設盡可能接近成立. 那么，究竟要將網格剖分至多小才能使計算結果的誤差不至于太大呢？如果網格太小(網格數量多)導致的問題又是什么呢？

為此， 本文設計了一個熱輻射傳播的球對稱問題， 以分析等溫法抽樣對網格剖分的依賴情況，并找到盡可能接近真解的模擬結果. 本文所有問題的計算平臺為個人電腦(CPU型號Intel(R)Core(TM) i7-3770@3.4 GHz; 核數8; 內存4.00 GB;操作系統為Windows 7旗艦版; 編譯器為Intel Visual Fortran).

模型為半徑0.2 cm的一維球， 其芯部(中心網格)為一溫度恒為1 keV的輻射源， 一維球外表面為真空邊界， 材料比熱為Cv= 0.1 gJ/(cm3·keV)，吸 收 系 數 與 溫 度 的 三 次 方 呈 反 比， σ = σ0/T3，σ0= 200， 單位為keV3/cm， 溫度T的單位為keV.

系統的初始物質溫度與輻射溫度將分別為1和0 eV， 計算時間步長為0.01 ns， 總時長為10 ns，網格數分別設置為20， 25， 40， 80， 100， 160， 200，250， 320， 400， 不同網格數的計算結果如圖1和圖2所示.

從圖1和圖2可以看出， 輻射源粒子采用等溫法抽樣， 不同網格剖分情況下的計算結果差異是很顯著的， 這說明等溫法的計算結果對空間剖分的依賴性很強， 顯然這對數值模擬來說是非常不好的性質. 同時可以看出， 隨著網格數的增加， 模擬結果是逐漸收斂的， 對于本問題， 當網格數達到200之后， 模擬結果就基本趨于一致了， 因此， 為了避免因網格剖分問題引入誤差， 本問題的網格數至少為200.

圖 1 不同網格數的物質溫度空間分布(t = 10 ns)Fig. 1. Material temperature with different cell numbers (t =10 ns).

圖 2 物質溫度的收斂情況 (a)不同網格數情況下r =0.05 cm處的物質溫度隨時間的變化; (b) r = 0.05 cm處物質溫度隨網格數的變化(t = 5 ns)Fig. 2. The convergence of material temperature: (a) Material temperature change with time in r = 0.05 cm; (b) material temperature change with cell number in r = 0.05 cm (t =5 ns).

從圖1還能看出， 當網格數目較少(網格較大)時， 輻射傳播的速度是偏快的， 這里稍做分析.溫度相對較高的區域輻射出的能量(粒子)更多，其與附近溫度相對較低的物質相互作用并加熱低溫區， 被加熱的區域輻射出能量加熱更遠的區域，熱輻射得以向前傳播. 在IMC模擬中， 系統區域被剖分為離散的網格， 在等溫假設下， 對于某個網格， 抽樣出的輻射源粒子均勻地分布在網格各處.然而， 在輻射的傳播過程中， 即使是某單一網格其內部也應該是有溫差的， 那么輻射源粒子應該更多地分布在網格中溫度較高處. 對于本問題， 更多的輻射源粒子本應分布在網格中半徑更小的位置， 因此在等溫假設下， 粒子的半徑抽樣值事實上是較理論值偏大了， 即粒子的輻射位置總體上比理論值靠前， 最終使得熱輻射的傳播比實際更快. 增加網格數， 即減小網格厚度可使網格內溫差減小， 等溫法抽樣產生的誤差隨之減小， 理論上網格越小誤差越小. 因此， 本文將網格數400的計算結果作為問題解的參考解. 為更清晰地比較各個溫度曲線與參考解的偏差， 表1列出了不同網格數時溫度曲線相對參考解的標準差和最大誤差， 容易看出， 隨著網格數的增加， 溫度分布曲線相對參考解的偏差確實變小了.

表 1 不同網格數時的溫度曲線相對參考解的標準偏差和最大誤差Table 1. Relative to the reference solution， the standard deviation and the maximum error of temperature curves with different cell numbers.

既然網格尺度大了有問題， 網格小了計算結果有保證， 那么是不是應該一開始就將網格分得很小呢？這樣做理論上是行得通的， 但是實踐中將會帶來弊端: 網格數增加后， 若保持每個時間步模擬粒子數不變， 每個網格能夠分配到的粒子數就會變少， 模擬結果的統計漲落將變大. 因此為避免統計漲落過大， 單個網格分配到的粒子數必須足夠多.表2是保證單個網格分配粒子數足夠多的情況下，不同網格剖分數對應的總模擬粒子數和計算時間，從表中可以看出， 隨著網格數的增加， 模擬的總粒子數相應增加， 由此帶來的后果是更多的模擬粒子花費更多的計算時間以及計算機內存， 導致模擬效率降低.

表 2 不同網格數的計算時間Table 2. Computation time with different cell numbers.

從前面的模擬結果及分析得知， 網格尺度不能太大， 否則計算誤差偏大. 但是如果網格尺度太小，計算開銷又太大. 那么“等溫法”抽樣就存在這樣的問題: 究竟要剖分多少網格合適？有沒有辦法讓計算結果不太依賴于網格剖分， 或者說在網格尺度較大的情況下計算誤差仍然較小呢？下面將回答這個問題.

3 輻射源粒子抽樣的兩種新方法

假設同一網格內溫度與半徑r的關系是線性的(以下稱為“線性假設”)， 這種假設在絕大多數情況下比“等溫假設”的近似效果更好. 如圖3所示，某網格的內外邊界半徑分別為r1， r2， 內外邊界物質溫度分別為T1， T2， 則溫度與r的關系可表示為

其中k = (T2— T1)/(r2— r1)， b = T1— kr1. 將方程(10)代入方程(7)中得

圖 3 網格內溫度與空間的關系可近似為線性關系Fig. 3. The dependence of temperature on space is approximately linear.

為分析方程(11)與等溫假設下推導出的(8)式在描述輻射源粒子空間分布時的差異， 下面以第2節熱輻射輸運問題中網格數為40的計算結果為例， 選取了第9， 18， 22， 26個網格(球的中心網格為第1個網格)作為輻射波的代表點. 圖4為相應網格的輻射源粒子空間概率密度分布.

圖4中曲線與橫坐標所圍成的面積代表輻射源粒子初始位置在空間分布的概率， 容易看出， 在越靠近輻射波波頭、網格內溫差越大的地方(圖4(d))，等溫假設下輻射源粒子總體位置與線性假設偏離越遠; 而在遠離波頭、網格內溫差越小的地方(圖4(a))，兩者越接近. 因此等溫假設下IMC模擬中輻射波傳播速度偏快主要是輻射波波頭處網格輻射源粒子總體位置與實際偏差太大造成的， 印證了第2節的討論結果.

對于方程(11)的概率密度分布函數， 如果用反函數法來抽樣r值， 則需要求解一元七次方程，比較困難. 下文給出兩種新的抽樣法.

1) 乘抽樣法

方程(11)可寫成

其中

方程(14)滿足概率密度函數的條件， 且H(r) ≥ 0，令為H(r)的上界， 可以證明[22]: 從f1(r)得到的抽樣值rf若滿足則rf為f(r)的抽樣值. 實際上f1(r)是球殼內均勻分布函數， 容易得到其隨機抽樣值

因此， 乘抽樣法步驟為:

① 抽樣得到f1(r)的抽樣值rf;

② 將rf代入方程(13)中得到H(rf) ;

③ 取出一隨機數ξ， 與H(rf)/ 比較;

④ 若ξ ≤ H(rf)/， 則rf為f(r)的抽樣值，否則返回步驟①.

圖 4 輻射波不同位置的輻射源粒子空間分布概率密度 (a)網格9， 波后處; (b)網格18， 波后處; (c)網格22， 波頭處; (d)網格26， 波頭處Fig. 4. Spatial probability density distribution of radiation source particle in different positions of radiation wave: (a) Cell 9， in the behind of wave; (b) cell 18， in the behind of wave; (c) cell 22， in the head of wave; (d) cell 26， in the head of wave.

2) 階梯近似抽樣法

將區間[r1， r2]分成m個子區間， 則每個子區間長度δr = (r2— r1)/m， 第i個子區間為[r1+(i — 1)δr，r1+ iδr]， 對方程(11)在[r1， r2]范圍內積分，

因此階梯近似抽樣法的抽樣步驟為:

需要注意的是， 對于不同問題， 乘抽樣法和階梯近似抽樣法的抽樣效率可能不相同， 因此在實際計算時， 選擇哪種抽樣法可視情況而定.

4 數值結果

為檢驗乘抽樣法和階梯近似抽樣法能否有效減小空間均勻抽樣帶來的偏差， 本節仍以第2節中的熱輻射輸運問題作為算例.

模型參數及計算參數與第2節中描述的相同，進行熱輻射輸運模擬時， 分別采用乘抽樣法和階梯近似抽樣法抽取網格輻射源粒子的初始位置， 并將計算結果與參考解比較. 圖5為分別采用兩種新抽樣方法在不同網格數時計算得到的物質溫度空間分布. 圖6為兩種新抽樣方法在網格數為40時的計算結果與參考解的比較.

圖 5 不同網格數時的物質溫度空間分布 (t = 10 ns) (a)乘抽樣法; (b)階梯近似抽樣法Fig. 5. Material temperature with different cell numbers(t = 10 ns): (a) Multiplying sampling method; (b) stepped approximation sampling method.

圖 6 網格數為40時兩種新抽樣法計算得到的物質溫度空間分布(t = 10 ns)Fig. 6. Results of two new sampling methods with 40 cells(t = 10 ns).

從圖5可以看出， 除網格數為20時的計算結果外， 兩種新抽樣法的計算結果都沒有出現明顯的輻射波傳播速度過快的問題. 值得注意的是， 當網格數為40時， 兩種新抽樣法的計算結果已經與參考解基本一致(如圖6所示)， 僅僅在輻射波頭處與參考解有所偏差， 其原因與文獻[20]中的類似， 即在波頭處溫度隨半徑r的變化規律偏離了線性， 溫度線性變化假設不太適用(事實上網格數為20時的計算結果的誤差也是由該原因引起的). 波頭處的偏差可通過加密網格或引入更復雜的T-r關系及抽樣方法解決， 該偏差對大多數熱輻射輸運問題的影響不大. 僅從圖5和圖6還不能明顯地看出乘抽樣法和階梯近似抽樣法對模擬結果修正的差異，為了更加直觀地比較兩者的模擬結果以及兩者相對等溫抽樣法的修正效果， 表3列出了各溫度曲線相對基準解的偏差， 可以看出， 兩種新抽樣法得到的溫度曲線同樣隨著網格數的增加而減小， 另外在40網格時兩者的溫度曲線與參考解的標準偏差已經接近甚至小于等溫抽樣法200網格時溫度曲線的標準偏差(如表1)， 而最大誤差更是明顯小于后者， 說明兩種新抽樣方法在40網格時得到的溫度曲線已經和參考解相當. 至于乘抽樣法和階梯近似抽樣法哪個更優， 后者的標準偏差和最大誤差僅比前者略小， 因此并不能確定后者比前者更優， 具體采用哪種方法可根據實際情況確定. 總體上， 線性假設是球殼中溫度空間分布規律的一種較好的近似， 由此推導出的兩種新抽樣方法也比等溫法更符合源粒子的發射規律.

表 3 不同網格數時的溫度曲線相對基準解的標準偏差和最大誤差Table 3. Relative to the reference solution， the standard deviation， and the maximum error of temperature curves with difference cell numbers.

輻射源粒子抽樣方法改進后， 計算結果不太依賴于網格剖分， 即只要網格尺度不是特別大(例題的20個網格)， 各種網格尺度的計算結果均能基本保持一致， 因為線性假設下網格內輻射源粒子的空間分布更加符合實際， 受網格尺度的影響很小. 下面從節省計算時間角度來看新方法的改進效果，表4是采用新舊抽樣法模擬不同網格剖分模型的計算時間， 可以看出計算費用基本上均與模擬粒子數呈線性增長關系. 另一方面， 根據圖6的比較可以知道， 在網格相對比較大(40個網格)的情況下，計算結果的精度已經能夠得到保證， 故利用新方法來模擬時可采用相對較少的網格和較少的粒子， 由此可以節省大量的計算機時. 對于本文例題， 如果采用舊方法， 則網格至少為200， 對應的計算機時為2.28 × 104s， 采用新方法， 則網格只需要40即可， 對應的計算機時為2.56 × 103s， 這大致相當于效率提升了約9倍.

表 4 計算問題所花時間Table 4. Computation time of the problem.

5 結 論

本文研究了球幾何中輻射源粒子的空間分布，分析了IMC方法中傳統源粒子抽樣方法的不足，推導了球幾何中基于能量密度分布的源粒子空間概率密度函數及相關參數的計算公式， 提出了兩種新的源粒子空間抽樣方法及流程， 新的空間概率密度分布函數能更準確地描述源粒子輻射的物理規律， 而且兩種新抽樣方法也能避免傳統抽樣方法在“強吸收”或“大網格”情況下計算結果偏差太大的問題. 由于在球幾何中確定論方法難以得到精確的熱輻射輸運的解析解， 本文在傳統的源粒子空間均勻抽樣法下， 通過精細劃分網格使計算結果收斂于真解， 并將該收斂后的結果作為兩種新的源粒子抽樣方法計算結果的參考解. 結果表明: 兩種新抽樣方法的計算結果與參考解相符， 很好地解決了“強吸收”或“大網格”情況下IMC模擬的熱輻射傳播速度過快的問題以及計算結果依賴于網格剖分的困擾; 同時由于新抽樣法只需較少的網格、較少的粒子便可得到傳統抽樣法需要較多網格、較多粒子才能獲得的準確結果， 故新抽樣法還顯著提高了計算效率， 大幅節省了計算資源. 下一步將開展復雜T-r關系下輻射源粒子精確抽樣方法研究.