閆 凱，田 宙，郭永輝，董 楠

（西北核技術(shù)研究所，陜西 西安710024）

輻射輸運(yùn)的模型問題一直是強(qiáng)爆炸火球數(shù)值模擬中的一個主要難點［1］。在早期的火球計算中，對輻射輸運(yùn)的描述往往采用根據(jù)區(qū)域劃分為不同近似的方法。H.L.Brode等［2－3］對輻射輸運(yùn)過程的描述采用雙流近似和較簡單的P1近似方法。陳健華等［4］依據(jù)溫度的分布采用輻射熱傳導(dǎo)近似、輻射熱輸出近似描述輻射輸運(yùn)過程?；鹎虻陌l(fā)展過程中光學(xué)薄和光學(xué)厚區(qū)域同時存在，采用上述簡單劃分的模型不能準(zhǔn)確描述火球的輻射場，并且由于輻射場是連續(xù)的，因此對輻射場的描述也應(yīng)該采用連續(xù)的模型。王心正等［5］、田宙等［6］采用最大熵變Eddington因子P1近似，即 Minerbo近似，描述輻射輸運(yùn)過程，該模型可以連續(xù)地描述輻射場。

Eddington近似由于其具有計算簡單、適合大規(guī)模并行計算等優(yōu)點，在天體物理、強(qiáng)爆炸火球等領(lǐng)域的數(shù)值模擬中得到了廣泛的應(yīng)用［5－6，8－11］。如何構(gòu)造一個合適的Eddington因子一直是該方法的一個難點。C.D.Levermore等［12］和A.M.Anile等［13］依據(jù)輻射熵原理分別獨立地提出了一個 Eddington因子近似模型。該模型又被稱為M1近似，目前已被廣泛應(yīng)用在天體物理的數(shù)值模擬中［8－11，14－15］，但其在強(qiáng)爆炸火球數(shù)值計算中的應(yīng)用還未見報道。本文中，將M1近似引入到強(qiáng)爆炸火球的數(shù)值計算中，并將計算結(jié)果與現(xiàn)有輻射輸運(yùn)模型的計算結(jié)果進(jìn)行比較。

1 計算模型

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

本文中，采用Euler型一維球?qū)ΨQ輻射流體動力學(xué)方程組，其基本表達(dá)式由輻射輸運(yùn)方程：

和Euler方程：

耦合構(gòu)成。式中：r為空間坐標(biāo)，t為時間坐標(biāo)，c為光速，υ為光子運(yùn)動的頻率，Ω 為光子運(yùn)動的方向，i為輻射強(qiáng)度，j為發(fā)射率，κ為吸收系數(shù)；p、ρ、u、T 分別為氣體的壓力、密度、速度和溫度；Em為氣體體積內(nèi)能；I為單位矩陣；Er，F(xiàn)r和pr分別為輻射能密度、輻射能流和輻射壓力，其定義如下：

一般來說，由于光子頻率、光子運(yùn)動方向等引起的復(fù)雜性，直接求解輻射流體方程組是極其困難的.在強(qiáng)爆炸和天體物理的數(shù)值模擬中，一個常見的做法就是對輻射輸運(yùn)方程做角度積分，這樣可以得到輻射輸運(yùn)的零階矩方程和一階矩方程組成的方程組：

方程組（4）是輻射輸運(yùn)矩方程在實驗室坐標(biāo)系下的形式，如果考慮流體速度帶來的相對論效應(yīng)，方程組形式將更復(fù)雜［16－17］。如果忽略散射和相對論效應(yīng)，采用局部熱力學(xué)平衡、灰體近似條件，輻射輸運(yùn)矩方程組將簡化為：

式中：κ′為考慮到受激輻射后的吸收系數(shù)，a為Boltzmann常數(shù)。

1.2 Eddington方法

輻射輸運(yùn)矩方程組的一個基本問題是它不封閉，以方程組（5）為例，除去流體的物理量，矩方程組有3個未知量，卻只有2個方程。為了封閉方程組，一個流行的作法是引入Eddington張量f，令：

對于多維模擬，f是一個對稱張量。本文中只討論一維情形，此時f為標(biāo)量。Eddington因子可以定義為：

如果f＝1/3，那么該近似就是P1近似，又稱為常Eddington因子近似。對于輻射近似各向同性的輻射場，P1近似是一個很好的近似；同時，P1近似方程簡單，計算量小，在強(qiáng)爆炸火球數(shù)值模擬早期了得到廣泛應(yīng)用［2－4］。P1近似的缺點也很明顯，輻射波速度恒定為倍光速，它不適合描述光子在光學(xué)薄區(qū)域的傳播。

20世紀(jì)70年代，G.N.Minerbo［7］發(fā)現(xiàn)輻射強(qiáng)度角分布與輻射能流Fr定義的方向的方位角無關(guān)，其根據(jù)最大熵原則對Eddington因子f給出了一個有理逼近式：

式中：R為各向異性因子，R＝0～1。R＝0，意味著輻射場為各向同性；R＝1，即光子朝一個方向運(yùn)動。該模型構(gòu)建的f具有一般Eddington因子都具有的性質(zhì)，即f（0）＝1/3，f（1）＝1，f（R）隨R 的增大而增大。王心正等［5］和田宙等［6］利用該模型對強(qiáng)爆炸早期火球進(jìn)行了數(shù)值模擬，取得了較好的效果。

1.3 M1近似

A.M.Anile等［13］分析了一維輻射輸運(yùn)矩方程組，提出Eddington因子應(yīng)該滿足以下條件：

（1）由于f是輻射輸運(yùn)方程歸一化的二階矩，因此需要滿足R2≤f≤1；

（3）在光子朝一個方向運(yùn)動時，應(yīng)滿足f（1）＝1；

（4）f應(yīng)該是R 的一個凸函數(shù)，即f″（R）≥0；

（5）在自由流中，奇異點應(yīng)該以光速傳播；

（6）能流限制器應(yīng)隨各向異性因子的增大而減小，由此推出f′（R）－2R≤0。

最后，A.M.Anile等根據(jù)輻射熵原理構(gòu)造的Eddington因子表達(dá)式為：

該表達(dá)式最早被 C.D.Levermore［12］提出，目前已被廣泛應(yīng)用于天體物理的數(shù)值模擬中［8－11］，在強(qiáng)爆炸火球的數(shù)值計算中則未見報道。

2 一維輻射輸運(yùn)矩方程的特征值問題

假定輻射波在真空中傳播，則此時可忽略耦合項，輻射壓力采用Eddington因子的表達(dá)式，則一維輻射輸運(yùn)矩方程組可以寫成：

容易證明P1近似和M1近似下的一維輻射輸運(yùn)矩方程屬于雙曲型方程。對于雙曲方程，其特征值可以反映方程特征波的波速。針對不同的輻射壓力模型，可以計算得到其最大特征值隨各向異性因子R的變化關(guān)系。其法向通量的Jacobi矩陣為：

考慮Fr＞0的情況，此時：

通過化簡Jacobi矩陣的變量，可以得到：

對于 Minerbo近似［7］：

對于 M1近似［13］：

Jacobi矩陣的特征值可以表示為：

則矩陣的最大特征值應(yīng)該為：

圖1給出了不同近似模型中Eddington因子f隨各向異性因子R的變化關(guān)系。從圖中可以看出，變Eddington因子P1近似和M1近似所構(gòu)造的Eddington因子均符合一般Eddington因子的性質(zhì)，但不能區(qū)分兩者之間優(yōu)劣。

圖2給出了不同近似下輻射輸運(yùn)矩方程的最大特征值λmax與各向異性因子R的關(guān)系。從圖中可以看出，P1近似下的輻射波波速與R無關(guān)，其速度為恒定值；Minerbo近似在R比較大的時候出現(xiàn)了較明顯的超光速現(xiàn)象，R＝1時其輻射波波速接近1.2倍光速；M1近似下的輻射波波速在R＝0時，即各向同性時，與P1近似下的輻射波波速一致，R＝1時，即光子朝一個方向運(yùn)動時，輻射波波速和光速一致。由于輻射輸運(yùn)方程組的最大特征值實際反映了輻射波波速，因此可以看出，在M1近似下得到的輻射波波速更符合物理規(guī)律。對于Fr＜0的情況，可以得到類似的結(jié)果，這里不再贅述。

圖1 不同近似下Eddington因子隨各向異性因子的變化Fig.1Eddington factor varied with anisotropy factor for different approach models

圖2 不同近似下輻射輸運(yùn)矩方程的最大特征值隨各向異性因子變化Fig.2The largest eigenvalues of moment equations of radiative transfer varied with anisotropy factor for different approach models

3 強(qiáng)爆炸火球的數(shù)值模擬

本文中，利用M1近似構(gòu)造的輻射輸運(yùn)模型對一維強(qiáng)爆炸火球問題進(jìn)行數(shù)值模擬，并將計算結(jié)果與已有結(jié)果進(jìn)行比較。

計算中選取的強(qiáng)爆炸當(dāng)量為1kt TNT，等壓球半徑為0.75m。源區(qū)物質(zhì)等價于實際氣體，其密度為海平面空氣密度；初始時刻火球為等壓球，處于輻射平衡狀態(tài)；火球內(nèi)邊界采用對稱條件，外邊界采用透射條件。圖3給出了采用M1近似計算得到的1kt TNT當(dāng)量下強(qiáng)爆炸的火球陣面和沖擊波走時以及 H.L.Brode［2］和田宙等［6］的計算結(jié)果。

H.L.Brode［2］對輻射輸運(yùn)方程的描述采用了雙流近似，即火球內(nèi)部采用P1近似、外部采用輻射熱輸出近似。圖3中火球陣面的傳播過程分為輻射擴(kuò)張階段、過渡階段和沖擊波傳播階段，對應(yīng)火球擴(kuò)張的主要原因分別是輻射輸運(yùn)、輻射輸運(yùn)與流體力學(xué)過程共同作用和沖擊波擴(kuò)張。在火球擴(kuò)張的前0.2μs，采用M1近似得到的火球陣面剛開始超過了H.L.Brode［2］計算得到的火球陣面，這是由于在早期X射線火球處于非平衡輻射擴(kuò)散中，各向異性因子R接近1，因此采用M1近似計算得到的輻射熱波傳播速度高于采用雙流近似的情況。0.2μs后的一段時間內(nèi)，火球內(nèi)部開始接近輻射平衡，M1近似下的輻射熱波速度下降很快，而雙流近似下的輻射熱波傳播速度受輻射輸運(yùn)的影響較小，因此出現(xiàn)M1近似下的輻射波陣面落后雙流近似的情況。幾微秒后，火球進(jìn)入過渡階段，輻射輸運(yùn)對火球陣面的擴(kuò)張作用越來越小，因此采用這2種方法計算得到的火球陣面趨于一致。由于雙流近似中沒有考慮輻射壓力對熱空氣動量的影響，因此雙流近似下得到的彈殼沖擊波陣面明顯落后于M1近似下得到的彈殼沖擊波陣面。后期，輻射輸運(yùn)的影響越來越弱，因此采用這2種近似得到的結(jié)果又趨于一致。

在強(qiáng)爆炸火球數(shù)值模擬中，本文中采用的空氣狀態(tài)方程和吸收系數(shù)均與田宙等［6］采用的條件一致，因此兩者的結(jié)果最具有對比性。田宙等［6］在強(qiáng)爆炸火球數(shù)值模擬中對輻射輸運(yùn)方程的描述采用了變Eddington因子P1近似，即Minerbo近似。圖3中M1近似下的火球陣面落后于Minerbo近似下的計算結(jié)果［6］；而M1近似下的彈殼沖擊波的產(chǎn)生時間及走時均與Minerbo近似下的計算結(jié)果［6］一致。值得注意的是，在火球發(fā)展的中后期，即過渡階段及以后，M1近似下的輻射波和沖擊波走時與雙流近似下的計算結(jié)果［2］和Minerbo近似下的計算結(jié)果［6］均符合較好，這說明不同的輻射輸運(yùn)模型對輻射擴(kuò)張階段計算結(jié)果的影響較明顯，而對火球發(fā)展的過渡階段及沖擊波擴(kuò)張階段計算結(jié)果的影響開始減弱。

圖3 采用M1近似計算得到的強(qiáng)爆炸的火球陣面和沖擊波走時與已有結(jié)果［2，6］的比較Fig.3Fireball and shock wave fronts calculated by MI approach compared with existent results［2，6］

4 結(jié) 束 語

針對P1近似只適合于光學(xué)厚區(qū)域，Minerbo近似在光學(xué)薄區(qū)域出現(xiàn)了超光速現(xiàn)象等問題，本文中將M1近似模型引入到強(qiáng)爆炸火球數(shù)值計算中，并比較了3種不同輻射輸運(yùn)模型下火球陣面和沖擊波陣面走時的計算結(jié)果。數(shù)值實驗結(jié)果表明，在火球輻射的擴(kuò)張階段采用不同的輻射輸運(yùn)模型得到的計算結(jié)果有較明顯的區(qū)別，初步定性分析了本文中采用模型的合理性，下一步將進(jìn)一步定量驗證M1近似模型用于強(qiáng)爆炸火球數(shù)值模擬的合理性。

［1］ 喬登江.核爆炸物理概論［M］.北京：國防工業(yè)出版社，2003：169 －262.

［2］ Brode H L.Fireball phenomenology［R］.The RAND Corporation.AD0612197，1965.

［3］ Brode H L，Hillendahl R W，Landshoff R K.Thermal radiation phenomena.Volume V：Radiation hydrodynamics of high temperature air［R］.AD0672837，1967.

［4］ 陳健華，王心正，謝龍生，等.均勻空氣中的強(qiáng)爆炸一維輻射流體力學(xué)數(shù)值解［J］.爆炸與沖擊，1981（2）：37 －49.Chen Jian －h(huán)ua，Wang Xin －zheng，Xie Long －sheng，et al.A one －dimensional radiation hydrodynamic numerical solution for a strong explosion in uniform atmosphere［J］.Explosion and Shock Waves，1981（2）：37 －49.

［5］ 王心正，隋衛(wèi)星.高空核爆炸火球的二維輻射流體力學(xué)計算［J］.計算物理，1987，4（2）：159 －168.Wang Xin －zheng，Sui Wei －xing.Two －dimension radiation hydrodynamics calculation of the high －altitude fireball［J］.Chinese Journal of Computational Physics，1987，4（2）：159 －168.

［6］ 田宙，喬登江，郭永輝.強(qiáng)爆炸早期火球現(xiàn)象的一維數(shù)值研究［J］.計算物理，2010，27（1）：8 －14.Tian Zhou，Qiao Deng－jiang，Guo Yong －h(huán)ui.A one －dimensional numerical study on early fireball in strong explosion［J］.Chinese Journal of Computational Physics，2010，27（1）：8 －14.

［7］ Minerbo G N.Maximum entropy Eddington factors［J］.Journal of Quantitative Spectroscopy and Radiative Transfer，1978，20（6）：541 －545.

［8］ Kumholz M R，Klein R I，Mckee C F，et al.Equations and algorithms for mixed －frame flux －limited diffusion radiation hydrodymics［J］.The Astrophysical Journal，2007，667（1）：626 －643.

［9］ Seaid M，Klar A，Dubroca B.Flux limiters in the coupling of radiation and hydrodynamic models［J］.Journal of Computational and Applied Mathematics，2004，168（1/2）：425 －435.

［10］ Buer C，Despres B.Asymptotic preserving and positive schemes for radiation hydrodynamics［J］.Journal of Computational Physics，2006，215（2）：717 －740.

［11］ Swesty F D，Myra E S.A numerical algorithm for modeling multigroup neutrino －radiation hydrodynamics in two spatial dimensions［J］.The Astrophysical Journal Supplement Series，2009，181（1）：1 －52.

［12］ Levermore C D.Relating Eddington factors to flux limiters［J］.Journal of Quantitative Spectroscopy and Radiative Transfer，1984，31（2）：149 －160.

［13］ Anile A M，Pennisi S，Sammartino M.A thermodynamical approach to Eddington factors［J］.Journal of Mathematical Physics，1991，32（2）：544 －555.

［14］ Brunner T A，Holloway J P.One －dimensional Riemann solvers and the maximum entropy closure［J］.Journal of Quantitative Spectroscopy and Radiative Transfer，2001，69（5）：543 －566.

［15］ Buet C，Despres B.Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics［J］.Journal of Quantitative Spectroscopy and Radiative Transfer，2004，85（3/4）：385 －418.

［16］ Castor J I.Radiation hydrodynamics［M］.Cambridge，UK：Cambridge University Press，2004：213 －245.

［17］ Pomraning G C.The equations of radiation hydrodynamics［M］.Dover：Dover Publications Inc，2005：427 －505.