球幾何中輻射源粒子抽樣方法的改進(jìn)*

許育培 李樹(shù)

1) (北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所， 北京 100094)2) (中國(guó)工程物理研究院研究生院， 北京 100088)(2020 年1 月5日收到; 2020 年4 月17日收到修改稿)

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

1 引 言

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

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

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

本文的主要內(nèi)容如下， 第2節(jié)推導(dǎo)了球幾何情況下“等溫法”的輻射源粒子空間概率密度分布函數(shù)及抽樣方法， 探討了該方法的不足， 并通過(guò)數(shù)值手段獲得了一個(gè)典型球?qū)ΨQ(chēng)熱輻射輸運(yùn)問(wèn)題的參考解. 第3節(jié)推導(dǎo)了基于能量密度分布的輻射源空間概率密度函數(shù)， 分析了IMC模擬中熱輻射波傳播偏快的原因， 提出了兩種新的輻射源粒子抽樣方法， 并設(shè)計(jì)了抽樣流程. 第4節(jié)通過(guò)典型算例測(cè)試了本文提出的兩種輻射源粒子抽樣方法的正確性，驗(yàn)證了本項(xiàng)研究工作對(duì)球幾何情況下的IMC輻射源粒子抽樣精度及計(jì)算效率有顯著的改進(jìn)效果.

2 等溫法的輻射源粒子抽樣及其缺點(diǎn)

IMC輻射輸運(yùn)模擬中的輻射源粒子可分成四種[16，21]: 來(lái)自獨(dú)立外源的粒子; 從網(wǎng)格邊界進(jìn)入的粒子; 從上一時(shí)間步存活下來(lái)的粒子; 從物質(zhì)輻射出的粒子(以下稱(chēng)為“輻射源粒子”). 前三種粒子均有確定的位置參數(shù)， 而輻射源粒子在當(dāng)前時(shí)間步開(kāi)始時(shí)需通過(guò)隨機(jī)抽樣確定其初始位置.

在一維球幾何情況下， IMC輸運(yùn)方程與物質(zhì)能量方程可寫(xiě)為:

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

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

其中

式中， a為輻射常數(shù)，Tn(r) 是網(wǎng)格溫度.

在 r 處dr范圍內(nèi)物質(zhì)輻射出的能量 En(r) 為

在球幾何情況下，.

IMC方法中， 輻射源粒子數(shù)目與其輻射的能量 En(r) 呈正比， 因此， 某網(wǎng)格內(nèi)的輻射源粒子關(guān)于 半徑 r 的空間分布概率密度函數(shù)可寫(xiě)為

其中r1， r2分別是網(wǎng)格的內(nèi)、外半徑.

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

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

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

為此， 本文設(shè)計(jì)了一個(gè)熱輻射傳播的球?qū)ΨQ(chēng)問(wèn)題， 以分析等溫法抽樣對(duì)網(wǎng)格剖分的依賴(lài)情況，并找到盡可能接近真解的模擬結(jié)果. 本文所有問(wèn)題的計(jì)算平臺(tái)為個(gè)人電腦(CPU型號(hào)Intel(R)Core(TM) i7-3770@3.4 GHz; 核數(shù)8; 內(nèi)存4.00 GB;操作系統(tǒng)為Windows 7旗艦版; 編譯器為Intel Visual Fortran).

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

系統(tǒng)的初始物質(zhì)溫度與輻射溫度將分別為1和0 eV， 計(jì)算時(shí)間步長(zhǎng)為0.01 ns， 總時(shí)長(zhǎng)為10 ns，網(wǎng)格數(shù)分別設(shè)置為20， 25， 40， 80， 100， 160， 200，250， 320， 400， 不同網(wǎng)格數(shù)的計(jì)算結(jié)果如圖1和圖2所示.

從圖1和圖2可以看出， 輻射源粒子采用等溫法抽樣， 不同網(wǎng)格剖分情況下的計(jì)算結(jié)果差異是很顯著的， 這說(shuō)明等溫法的計(jì)算結(jié)果對(duì)空間剖分的依賴(lài)性很強(qiáng)， 顯然這對(duì)數(shù)值模擬來(lái)說(shuō)是非常不好的性質(zhì). 同時(shí)可以看出， 隨著網(wǎng)格數(shù)的增加， 模擬結(jié)果是逐漸收斂的， 對(duì)于本問(wèn)題， 當(dāng)網(wǎng)格數(shù)達(dá)到200之后， 模擬結(jié)果就基本趨于一致了， 因此， 為了避免因網(wǎng)格剖分問(wèn)題引入誤差， 本問(wèn)題的網(wǎng)格數(shù)至少為200.

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

圖 2 物質(zhì)溫度的收斂情況 (a)不同網(wǎng)格數(shù)情況下r =0.05 cm處的物質(zhì)溫度隨時(shí)間的變化; (b) r = 0.05 cm處物質(zhì)溫度隨網(wǎng)格數(shù)的變化(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還能看出， 當(dāng)網(wǎng)格數(shù)目較少(網(wǎng)格較大)時(shí)， 輻射傳播的速度是偏快的， 這里稍做分析.溫度相對(duì)較高的區(qū)域輻射出的能量(粒子)更多，其與附近溫度相對(duì)較低的物質(zhì)相互作用并加熱低溫區(qū)， 被加熱的區(qū)域輻射出能量加熱更遠(yuǎn)的區(qū)域，熱輻射得以向前傳播. 在IMC模擬中， 系統(tǒng)區(qū)域被剖分為離散的網(wǎng)格， 在等溫假設(shè)下， 對(duì)于某個(gè)網(wǎng)格， 抽樣出的輻射源粒子均勻地分布在網(wǎng)格各處.然而， 在輻射的傳播過(guò)程中， 即使是某單一網(wǎng)格其內(nèi)部也應(yīng)該是有溫差的， 那么輻射源粒子應(yīng)該更多地分布在網(wǎng)格中溫度較高處. 對(duì)于本問(wèn)題， 更多的輻射源粒子本應(yīng)分布在網(wǎng)格中半徑更小的位置， 因此在等溫假設(shè)下， 粒子的半徑抽樣值事實(shí)上是較理論值偏大了， 即粒子的輻射位置總體上比理論值靠前， 最終使得熱輻射的傳播比實(shí)際更快. 增加網(wǎng)格數(shù)， 即減小網(wǎng)格厚度可使網(wǎng)格內(nèi)溫差減小， 等溫法抽樣產(chǎn)生的誤差隨之減小， 理論上網(wǎng)格越小誤差越小. 因此， 本文將網(wǎng)格數(shù)400的計(jì)算結(jié)果作為問(wèn)題解的參考解. 為更清晰地比較各個(gè)溫度曲線與參考解的偏差， 表1列出了不同網(wǎng)格數(shù)時(shí)溫度曲線相對(duì)參考解的標(biāo)準(zhǔn)差和最大誤差， 容易看出， 隨著網(wǎng)格數(shù)的增加， 溫度分布曲線相對(duì)參考解的偏差確實(shí)變小了.

表 1 不同網(wǎng)格數(shù)時(shí)的溫度曲線相對(duì)參考解的標(biāo)準(zhǔn)偏差和最大誤差Table 1. Relative to the reference solution， the standard deviation and the maximum error of temperature curves with different cell numbers.

既然網(wǎng)格尺度大了有問(wèn)題， 網(wǎng)格小了計(jì)算結(jié)果有保證， 那么是不是應(yīng)該一開(kāi)始就將網(wǎng)格分得很小呢？這樣做理論上是行得通的， 但是實(shí)踐中將會(huì)帶來(lái)弊端: 網(wǎng)格數(shù)增加后， 若保持每個(gè)時(shí)間步模擬粒子數(shù)不變， 每個(gè)網(wǎng)格能夠分配到的粒子數(shù)就會(huì)變少， 模擬結(jié)果的統(tǒng)計(jì)漲落將變大. 因此為避免統(tǒng)計(jì)漲落過(guò)大， 單個(gè)網(wǎng)格分配到的粒子數(shù)必須足夠多.表2是保證單個(gè)網(wǎng)格分配粒子數(shù)足夠多的情況下，不同網(wǎng)格剖分?jǐn)?shù)對(duì)應(yīng)的總模擬粒子數(shù)和計(jì)算時(shí)間，從表中可以看出， 隨著網(wǎng)格數(shù)的增加， 模擬的總粒子數(shù)相應(yīng)增加， 由此帶來(lái)的后果是更多的模擬粒子花費(fèi)更多的計(jì)算時(shí)間以及計(jì)算機(jī)內(nèi)存， 導(dǎo)致模擬效率降低.

表 2 不同網(wǎng)格數(shù)的計(jì)算時(shí)間Table 2. Computation time with different cell numbers.

從前面的模擬結(jié)果及分析得知， 網(wǎng)格尺度不能太大， 否則計(jì)算誤差偏大. 但是如果網(wǎng)格尺度太小，計(jì)算開(kāi)銷(xiāo)又太大. 那么“等溫法”抽樣就存在這樣的問(wèn)題: 究竟要剖分多少網(wǎng)格合適？有沒(méi)有辦法讓計(jì)算結(jié)果不太依賴(lài)于網(wǎng)格剖分， 或者說(shuō)在網(wǎng)格尺度較大的情況下計(jì)算誤差仍然較小呢？下面將回答這個(gè)問(wèn)題.

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

假設(shè)同一網(wǎng)格內(nèi)溫度與半徑r的關(guān)系是線性的(以下稱(chēng)為“線性假設(shè)”)， 這種假設(shè)在絕大多數(shù)情況下比“等溫假設(shè)”的近似效果更好. 如圖3所示，某網(wǎng)格的內(nèi)外邊界半徑分別為r1， r2， 內(nèi)外邊界物質(zhì)溫度分別為T(mén)1， T2， 則溫度與r的關(guān)系可表示為

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

圖 3 網(wǎng)格內(nèi)溫度與空間的關(guān)系可近似為線性關(guān)系Fig. 3. The dependence of temperature on space is approximately linear.

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

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

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

1) 乘抽樣法

方程(11)可寫(xiě)成

其中

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

因此， 乘抽樣法步驟為:

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

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

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

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

圖 4 輻射波不同位置的輻射源粒子空間分布概率密度 (a)網(wǎng)格9， 波后處; (b)網(wǎng)格18， 波后處; (c)網(wǎng)格22， 波頭處; (d)網(wǎng)格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) 階梯近似抽樣法

將區(qū)間[r1， r2]分成m個(gè)子區(qū)間， 則每個(gè)子區(qū)間長(zhǎng)度δr = (r2— r1)/m， 第i個(gè)子區(qū)間為[r1+(i — 1)δr，r1+ iδr]， 對(duì)方程(11)在[r1， r2]范圍內(nèi)積分，

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

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

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

為檢驗(yàn)乘抽樣法和階梯近似抽樣法能否有效減小空間均勻抽樣帶來(lái)的偏差， 本節(jié)仍以第2節(jié)中的熱輻射輸運(yùn)問(wèn)題作為算例.

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

圖 5 不同網(wǎng)格數(shù)時(shí)的物質(zhì)溫度空間分布 (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 網(wǎng)格數(shù)為40時(shí)兩種新抽樣法計(jì)算得到的物質(zhì)溫度空間分布(t = 10 ns)Fig. 6. Results of two new sampling methods with 40 cells(t = 10 ns).

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

表 3 不同網(wǎng)格數(shù)時(shí)的溫度曲線相對(duì)基準(zhǔn)解的標(biāo)準(zhǔn)偏差和最大誤差Table 3. Relative to the reference solution， the standard deviation， and the maximum error of temperature curves with difference cell numbers.

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

表 4 計(jì)算問(wèn)題所花時(shí)間Table 4. Computation time of the problem.

5 結(jié) 論

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