強爆炸光輻射脈沖輻照特征與爆炸當量的相關(guān)性*

高銀軍，田 宙，閆 凱，劉 峰

(1.北京理工大學材料學院,北京 100081；2.西北核技術(shù)研究所，陜西 西安 710024)

強爆炸光輻射脈沖輻照特征與爆炸當量的相關(guān)性*

高銀軍1,2，田 宙2，閆 凱2，劉 峰2

(1.北京理工大學材料學院,北京 100081；2.西北核技術(shù)研究所，陜西 西安 710024)

為定量分析強爆炸光輻射輻照特征及其與爆炸當量的關(guān)系，建立了用于描述光輻射輸運過程的輻射流體力學模型。在算子分裂方法運用的基礎(chǔ)上，采用溫度梯度作為指示子進行并行區(qū)域的動態(tài)劃分，從而實現(xiàn)較高效率的并行求解。在此基礎(chǔ)上數(shù)值計算了千噸～兆噸當量下強爆炸光輻射的發(fā)展過程，分析表明：光輻射強度隨時間呈現(xiàn)“雙脈沖”變化，強度極小和強度第2極大時刻與當量的某次方成正比。光輻射總功率變化歷程與光輻射強度變化歷程相似，但受輻射源半徑隨當量變化的影響，其極值時刻會出現(xiàn)差異。

強爆炸;光輻射;脈沖輻照;爆炸當量

光輻射是強爆炸重要毀傷效應之一。受大氣作用的影響，低空大氣中的強爆炸光輻射，其發(fā)展過程具有脈沖輻照的特征[1-2]：一方面光輻射向外發(fā)展造成輻射源溫度降低，致使其強度減弱，另一方面爆炸中的沖擊波擴張使得波后空氣對輻射的吸收由強變?nèi)?，導致光輻射形成了具?次極大的脈沖輻照特征。光輻射的脈沖輻照特征與當量具有密切關(guān)系[3-4]?，F(xiàn)有光輻射特征與當量關(guān)系的研究，主要是基于早期實驗結(jié)果總結(jié)的經(jīng)驗關(guān)系[1,5]。這些經(jīng)驗公式在最初的強爆炸現(xiàn)象研究中發(fā)揮了積極作用，但由于缺少必要的理論計算和分析，導致對于光輻射特征及其物理規(guī)律的研究不夠深入?，F(xiàn)階段，數(shù)值模擬是光輻射輻照特征研究的唯一方法，其定量分析結(jié)果可以直接彌補原有經(jīng)驗規(guī)律研究的不足，為進一步的光輻射特征及其物理規(guī)律研究提供有效的技術(shù)手段。光輻射特征的計算有其自身的復雜性[6-7]：一是描述強爆炸過程的輻流方程組具有強非線性和剛性[8-11]，使得輻射能流的計算很不穩(wěn)定；二是光輻射發(fā)展過程是一個大空間尺度、長時間過程，相對于其求解步長而言計算量巨大，需要發(fā)展高效并行方法。對于輻射能流計算的不穩(wěn)定性問題，閆凱等[12]、高銀軍等[13]引入了算子分裂方法，在火球參量研究中開展了初步分析，但尚未對光輻射輻照特征進行計算。在大規(guī)模并行求解方面，由于算子分裂方法中源項求解在迭代過程中計算量差異很大，極易造成并行計算的負載不平衡。本文在算子分裂法運用的基礎(chǔ)上，通過對矩方程源項特征及其與流體狀態(tài)參量的對比分析，以溫度梯度作為指示子進行區(qū)域劃分，較好實現(xiàn)負載平衡。在此基礎(chǔ)上計算不同當量下強爆炸光輻射輻照特征，研究光輻射極值時刻隨當量的變化，并對光輻射相關(guān)物理參量的空間分布開展定量分析。

1 計算模型

灰體近似下強爆炸火球光輻射發(fā)展的一維球?qū)ΨQ輻射流體力學方程組(滿足局域熱動平衡假定)[3,12,14]如下：

(1)

式中：ρ為空氣密度，u為r方向空氣速度，p為空氣壓強，e為空氣體積總能量，κ為靜止坐標系中考慮了受激發(fā)射的空氣吸收系數(shù)，pr為光輻射壓強，c為光速;φ=aT4為溫度T下的黑體輻射能密度,a=7.56×10-16J/(m3·K4);Er為光輻射能密度，F(xiàn)r為光輻射能流，即單位時間內(nèi)通過垂直于傳播方向單位面積上光輻射的能量，是表征光輻射強度的主要參量。采用最大熵Eddington因子下的輻射壓強pr=χEr封閉方程組[3,14-15]， 空氣的狀態(tài)方程采用實際空氣狀態(tài)方程[5]。

對于計算中的初始條件，假定爆炸釋放的能量集中于一定半徑的等壓球內(nèi)(體積為V)，輻射能(ErV)與流體內(nèi)能(eIV)之和為爆炸總能量(相當于爆炸當量所對應能量的85%，1 kt爆炸當量所對應的能量約為3.55×1012J)。

2 求解方法

圖1 強爆炸過程中溫度、吸收系數(shù)及源項的空間分布特征Fig.1 Spatial distribution of temperature, absorption coefficient and source terms in intense explosion

運用算子分裂法對上述方程進行求解，其核心是將方程分裂為對流項和剛性源項2步求解[12,14]：第1步求解對流項，由于不含相互作用項，求解精度和穩(wěn)定性都較高；第2步化為常微分方程求解，根據(jù)方程特性可以進一步轉(zhuǎn)化，從而實現(xiàn)快速穩(wěn)定求解。

光輻射輻照特征的計算，需要發(fā)展高效并行算法。對于算子分裂方法求解來說，通常的區(qū)域劃分方法會由于各個節(jié)點計算量的不同導致求解效率大大降低。通過對分裂求解的進一步分析，可以發(fā)現(xiàn)：對流項求解幾乎不消耗計算資源，影響計算效率的主要因素在于源項求解中的迭代計算。以零階矩源項M0(單位時間內(nèi)輻射與流體能量交換項)為例，其表達式(忽略相對論修正項)為：

(2)

由于吸收系數(shù)κ也與溫度T有關(guān)，所以上述源項的大小與溫度直接相關(guān)，從而最終影響源項求解所需要的計算時間。對于光輻射發(fā)展過程而言，溫度在空間分布上存在明顯的非均勻性，且梯度差異明顯，與溫度有關(guān)的物理參量也呈現(xiàn)類似的分布特征，見圖1。

由圖1可見：溫度梯度較大的區(qū)域，吸收系數(shù)與源項均有大的變化。從物理過程上來說，也就是在溫度變化非常劇烈的區(qū)域(如波陣面附近)，輻射與流體強烈相互作用，導致源項極其巨大，最終造成該區(qū)域單個物理量的求解所需計算量要比其他區(qū)域大得多。針對這一特點，在并行求解中為較大程度的實現(xiàn)負載平衡，采用以溫度梯度作為區(qū)域劃分的標準：溫度梯度越大，所求解的區(qū)域較??；反之則越大。另外，考慮到光輻射發(fā)展過程中，物理參量隨時間不斷變化，因此需要對溫度進行動態(tài)監(jiān)測，利用上述標準實現(xiàn)區(qū)域的自動劃分。采用該方法對一維情況下光輻射過程進行求解，與均勻區(qū)域劃分相比，當計算節(jié)點達到10時并行效率顯著提升。

3 計算分析

爆炸高度取h=0.5 km，在kt～1 Mt當量范圍內(nèi)，求解相應的強爆炸光輻射輸運過程，空氣初始狀態(tài)參量如表1所示，cs為聲速，爆炸高度小于1 km時狀態(tài)參量取h=0 km時的數(shù)值。

表1 不同高度下空氣初始狀態(tài)參數(shù)

圖2給出了不同當量下光輻射強度(即火球半徑處光輻射能流)隨時間的變化曲線。從不同曲線的變化規(guī)律來看，幾種當量下光輻射強度隨時間變化趨勢基本一致：均是先增大然后減小，到某一極小值后又逐漸增大，最后減小，即雙脈沖波形。

圖2 不同爆炸當量下光輻射強度隨時間變化關(guān)系Fig.2 Relationship between intensity of thermal radiation and time at different yields

對于第1個脈沖來說，主要是由于輻射源不斷向外輻射能量導致溫度降低，從而引起光輻射強度下降，第2個脈沖的形成過程則與沖擊波的作用密切相關(guān)。圖3給出了1 kt爆炸情況下，第2個脈沖形成前、形成后典型時刻流體物理參量(密度、壓力和溫度)的空間分布特征。從圖中可以看出：在第2個脈沖形成前，沖擊波(密度、壓力突變處)位于輻射波(溫度突變處)之后，向外發(fā)出的光輻射不會受到?jīng)_擊波的影響；而當沖擊波趕上并超過輻射波之后，盡管內(nèi)部輻射溫度很高，但流體的吸收作用決定了光輻射強度的變化，最終導致了第2個脈沖的出現(xiàn)。

光輻射強度極值所對應的時刻與當量有關(guān)。爆炸當量越大，光輻射極值時刻不斷增加。圖4給出了不同當量下光輻射第1極小和第2極大所對應的時刻。在對數(shù)坐標下，極值時刻與當量有著近似線性的關(guān)系。通過擬合并經(jīng)過相應的轉(zhuǎn)化，得到：

(3)

圖3 爆炸當量為1 kt的情況下第2個脈沖形成前、后流體參量的空間分布特征Fig.3 Spatial distribution characteristics of fluid parameter before and after the formation of second pulse at 1 kt yield

文獻[1,6-7]中認為：光輻射極小和第2極大時刻與當量有著形如t=aQb的關(guān)系。不同文獻中參數(shù)a有所不同，但參數(shù)b相差不大。文獻[1]給出的光輻射極值時刻與當量的經(jīng)驗公式中，參數(shù)b分別為0.44和0.46，計算結(jié)果與文獻結(jié)果基本一致。

對于光輻射第1極大時刻(見圖5)，計算結(jié)果與經(jīng)驗公式相比出現(xiàn)較大不同。從計算的角度來講，主要是由于光輻射第1極大時刻出現(xiàn)較早，該時刻空氣實際溫度很高(約百萬K)，吸收系數(shù)的非線性影響很大，導致輻射能流的計算存在不穩(wěn)定性。另外，根據(jù)文獻[7]，光輻射第1極大時刻受到除當量外的多種因素影響(如彈體質(zhì)量等)，經(jīng)驗公式僅是針對當量的近似估計。

圖4 光輻射強度的極值時刻與爆炸當量的關(guān)系Fig.4 Thermal radiation extremum times at different explosion yields

圖5 不同當量下不同方法得到的光輻射強度第一極大時刻Fig.5 First maximum times of thermal radiation at different yields by different methods

光輻射總功率是光輻射輻照特征的一個關(guān)鍵參量。若光輻射源(即火球)半徑為rb，則光輻射向外輻射的總功率為：

(4)

可見，光輻射功率的計算需要首先確定火球半徑。圖6(a)給出了不同當量下火球半徑隨時間變化關(guān)系。進一步根據(jù)式(4)計算得到光輻射總功率的變化歷程，見圖6(b)。由圖6(b)和圖2可以看出，光輻射功率走時與其強度的變化趨勢相似，但其極值時刻有所不同。其根本原因是光輻射源(火球)半徑隨時間的變化與當量有關(guān)。當量越大，相同時刻的火球半徑越大。整體上來看，光輻射功率隨當量的增大而不斷增大。

圖6 不同爆炸當量下火球半徑及光輻射總功率隨時間的變化Fig.6 Variation of fireball radius and thermal radiant power with as time at different yields

4 小 結(jié)

(1)采用溫度梯度作為指示子進行并行區(qū)域的動態(tài)劃分，在算子分裂方法基礎(chǔ)上達到了較高效并行求解，實現(xiàn)了光輻射特征及其當量相關(guān)性的數(shù)值模擬。

(2)不同當量下光輻射均呈現(xiàn)雙脈沖，其數(shù)值大小與當量有關(guān)。曲線擬合的結(jié)果表明：光輻射極小時刻和第二極大時刻分別與當量的0.443和0.451次方成正比。

(3)光輻射總功率走時與其強度走時有著類似的變化趨勢，但隨當量的變化關(guān)系有所不同，這主要是由于輻射源半徑隨當量變化所引起的。

[1] 喬登江.核爆炸物理概論[M].北京:國防工業(yè)出版社,2003:225-262.

[2] Brode H L. Fireball phenomenology: AD612197[R]. 1964.

[3] 田宙,喬登江,郭永輝.不同當量強爆炸早期火球現(xiàn)象的數(shù)值模擬[J].爆炸與沖擊,2009,29(4):408-412. Tian Zhou, Qiao Dengjiang, Guo Yonghui. Numerical simulation on early fireball phenomenology of strong explosions for different yields[J]. Explosion and Shock Waves, 2009,29(4):408-412.

[4] 田宙,喬登江,郭永輝.強爆炸早期火球現(xiàn)象的一維數(shù)值研究[J].計算物理,2010,27(1):9-14. Tian Zhou, Qiao Dengjiang, Guo Yonghui. A one dimensional numerical study on early fireball in strong explosion[J]. Chinese Journal of Computational Physics, 2010,27(1):9-14.

[5] 屠琴芬.核爆炸火球的輻射流體力學計算中的幾個問題[R].西安:西北核技術(shù)研究所,1986.

[6] Pomraning G C. The equations of radiation hydrodynamics[M]. Mineola: Pergamon Press, 1973:309-330.

[7] Symbalisty E M D, Zinn J, Whitaker R W. RADFLO physics and algorithms: LA-12988-MS[R]. 1995.

[8] Crowley B K, Glenn H D, Marks R E. An analysis of marvel: A nuclear shock-tube experiment[J]. Journal of Geophysical Research, 1971,76(14):3356-3374.

[9] Marrs R E, Moss W C, Whitlock B. Thermal radiation from nuclear detonations in urban environments: UCRL-TR-231593[R]. 2007.

[10] Lowrie R B, Edwards J D. Radiative shock solutions with grey nonequilibrium diffusion[J]. Shock Waves, 2008,18:129-143.

[11] 陳健華,王心正,謝龍生.均勻大氣中的強爆炸一維輻射流體力學數(shù)值解[J].爆炸與沖擊,1981,1(2):37-49. Chen Jianhua, Wang Xinzheng, Xie Longsheng. An one-dimensional radiation hydrodynamic numerical solution for a strong explosion in uniform atmosphere[J]. Explosion and Shock Waves, 1981,1(2):37-49.

[12] 閆凱,李若,田宙,等.強爆炸火球數(shù)值模擬中的一種算子分裂方法[J].計算物理,2013,30(3):379-386. Yan Kai, Li Ruo, Tian Zhou, et al. An operator splitting method for numerical simulation of strong explosion fireball[J]. Chinese Journal of Computational Physics, 2013,30(3):379-386.

[13] 高銀軍,閆凱,田宙,等.強爆炸早期火球光輻射能譜的數(shù)值計算[J].爆炸與沖擊,2015,35(3):289-295. Gao Yinjun, Yan Kai, Tian Zhou, et al. Numerical calculation of early fireball radiation spectrum in strong explosion[J]. Explosion and Shock Waves, 2015,35(3):289-295.

[14] 閆凱,田宙,郭永輝,等.M1近似應用于強爆炸火球輻射輸運數(shù)值模擬[J].爆炸與沖擊,2014,34(2):241-246. Yan Kai, Tian Zhou, Guo Yonghui, et al. Application of M1 approach to numerical simulation of radiative transfer in strong explosion fireball[J]. Explosion and Shock Waves, 2014,34(2):241-246.

[15] Minerbo G N. Maximum entropy Eddington factors[J]. Journal of Quantitative Spectroscopy & Radiative Transfer, 1977,20(6):541-545.

(責任編輯 張凌云)

Correlation between pulse irradiation characteristic of thermal radiation in intense explosion and explosion yield

Gao Yinjun1,2, Tian Zhou2, Yan Kai2, Liu Feng2

(1.SchoolofMaterialScienceandEngineering,BeijingInstituteofTechnology,Beijing100081,China;2.NorthwestInstituteofNuclearTechnology,Xi’an710024,Shaanxi,China)

A radiation hydrodynamics model was established to analyse the pulse irradiation characteristic of thermal radiation in intense explosion and its dependence on explosion yield. Based on the splitting method, the temperature gradient was set as the indicator for the dynamic regional division to achieve highly efficient parallel calculation, on the basis of which the thermal radiation evolutions in intense explosion at different yields were calculated. The results show that the intensity of thermal radiation at different times exhibits a two-pulse pattern. The intensity extremums and extremum times vary with the change of the explosion yields. The minimum time and second maximum time of the thermal radiation are proportional to the power of the explosion yields. The radiant power history exhibits similar results with the thermal radiation intensity. However, the extremum times may differ due to the dependence of the fireball radius on the explosion yield.

intense explosion; thermal radiation; pulse irradiation; explosion yield

10.11883/1001-1455(2017)03-0549-05

2015-11-02；

2016-05-26

國家自然科學基金項目(91330205)

高銀軍(1983— )，男,博士研究生，gyj@mail.ustc.edu.cn。

O381 國標學科代碼： 13035

A