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

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

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

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

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

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

強(qiáng)爆炸;光輻射;脈沖輻照;爆炸當(dāng)量

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

1 計(jì)算模型

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

(1)

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

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

2 求解方法

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

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

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

(2)

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

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

3 計(jì)算分析

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

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

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

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

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

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

(3)

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

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

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

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

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

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

(4)

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

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

4 小 結(jié)

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

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

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

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

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

[3] 田宙,喬登江,郭永輝.不同當(dāng)量強(qiáng)爆炸早期火球現(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] 田宙,喬登江,郭永輝.強(qiáng)爆炸早期火球現(xiàn)象的一維數(shù)值研究[J].計(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] 屠琴芬.核爆炸火球的輻射流體力學(xué)計(jì)算中的幾個(gè)問(wèn)題[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] 陳健華,王心正,謝龍生.均勻大氣中的強(qiáng)爆炸一維輻射流體力學(xué)數(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] 閆凱,李若,田宙,等.強(qiáng)爆炸火球數(shù)值模擬中的一種算子分裂方法[J].計(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] 高銀軍,閆凱,田宙,等.強(qiáng)爆炸早期火球光輻射能譜的數(shù)值計(jì)算[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近似應(yīng)用于強(qiáng)爆炸火球輻射輸運(yùn)數(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.

(責(zé)任編輯 張凌云)

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

國(guó)家自然科學(xué)基金項(xiàng)目(91330205)

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

O381 國(guó)標(biāo)學(xué)科代碼： 13035

A