M1近似應用于強爆炸火球輻射輸運數值模擬＊

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

（西北核技術研究所，陜西 西安710024）

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

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

1 計算模型

1.1 輻射流體力學方程組

本文中，采用Euler型一維球對稱輻射流體動力學方程組，其基本表達式由輻射輸運方程：

和Euler方程：

耦合構成。式中：r為空間坐標，t為時間坐標，c為光速，υ為光子運動的頻率，Ω 為光子運動的方向，i為輻射強度，j為發射率，κ為吸收系數；p、ρ、u、T 分別為氣體的壓力、密度、速度和溫度；Em為氣體體積內能；I為單位矩陣；Er，Fr和pr分別為輻射能密度、輻射能流和輻射壓力，其定義如下：

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

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

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

1.2 Eddington方法

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

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

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

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

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

1.3 M1近似

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

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

（3）在光子朝一個方向運動時，應滿足f（1）＝1；

（4）f應該是R 的一個凸函數，即f″（R）≥0；

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

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

最后，A.M.Anile等根據輻射熵原理構造的Eddington因子表達式為：

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

2 一維輻射輸運矩方程的特征值問題

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

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

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

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

對于 Minerbo近似［7］：

對于 M1近似［13］：

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

則矩陣的最大特征值應該為：

圖1給出了不同近似模型中Eddington因子f隨各向異性因子R的變化關系。從圖中可以看出，變Eddington因子P1近似和M1近似所構造的Eddington因子均符合一般Eddington因子的性質，但不能區分兩者之間優劣。

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

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

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

3 強爆炸火球的數值模擬

本文中，利用M1近似構造的輻射輸運模型對一維強爆炸火球問題進行數值模擬，并將計算結果與已有結果進行比較。

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

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

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

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

4 結 束 語

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

［1］ 喬登江.核爆炸物理概論［M］.北京：國防工業出版社，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］ 陳健華，王心正，謝龍生，等.均勻空氣中的強爆炸一維輻射流體力學數值解［J］.爆炸與沖擊，1981（2）：37 －49.Chen Jian －hua，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］ 王心正，隋衛星.高空核爆炸火球的二維輻射流體力學計算［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］ 田宙，喬登江，郭永輝.強爆炸早期火球現象的一維數值研究［J］.計算物理，2010，27（1）：8 －14.Tian Zhou，Qiao Deng－jiang，Guo Yong －hui.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.