强爆炸光辐射脉冲辐照特征与爆炸当量的相关性*

高银军，田 宙，闫 凯，刘 峰

(1.北京理工大学材料学院,北京 100081；2.西北核技术研究所，陕西 西安 710024)

强爆炸光辐射脉冲辐照特征与爆炸当量的相关性*

高银军1,2，田 宙2，闫 凯2，刘 峰2

(1.北京理工大学材料学院,北京 100081；2.西北核技术研究所，陕西 西安 710024)

为定量分析强爆炸光辐射辐照特征及其与爆炸当量的关系，建立了用于描述光辐射输运过程的辐射流体力学模型。在算子分裂方法运用的基础上，采用温度梯度作为指示子进行并行区域的动态划分，从而实现较高效率的并行求解。在此基础上数值计算了千吨～兆吨当量下强爆炸光辐射的发展过程，分析表明：光辐射强度随时间呈现“双脉冲”变化，强度极小和强度第2极大时刻与当量的某次方成正比。光辐射总功率变化历程与光辐射强度变化历程相似，但受辐射源半径随当量变化的影响，其极值时刻会出现差异。

强爆炸;光辐射;脉冲辐照;爆炸当量

光辐射是强爆炸重要毁伤效应之一。受大气作用的影响，低空大气中的强爆炸光辐射，其发展过程具有脉冲辐照的特征[1-2]：一方面光辐射向外发展造成辐射源温度降低，致使其强度减弱，另一方面爆炸中的冲击波扩张使得波后空气对辐射的吸收由强变弱，导致光辐射形成了具有2次极大的脉冲辐照特征。光辐射的脉冲辐照特征与当量具有密切关系[3-4]。现有光辐射特征与当量关系的研究，主要是基于早期实验结果总结的经验关系[1,5]。这些经验公式在最初的强爆炸现象研究中发挥了积极作用，但由于缺少必要的理论计算和分析，导致对于光辐射特征及其物理规律的研究不够深入。现阶段，数值模拟是光辐射辐照特征研究的唯一方法，其定量分析结果可以直接弥补原有经验规律研究的不足，为进一步的光辐射特征及其物理规律研究提供有效的技术手段。光辐射特征的计算有其自身的复杂性[6-7]：一是描述强爆炸过程的辐流方程组具有强非线性和刚性[8-11]，使得辐射能流的计算很不稳定；二是光辐射发展过程是一个大空间尺度、长时间过程，相对于其求解步长而言计算量巨大，需要发展高效并行方法。对于辐射能流计算的不稳定性问题，闫凯等[12]、高银军等[13]引入了算子分裂方法，在火球参量研究中开展了初步分析，但尚未对光辐射辐照特征进行计算。在大规模并行求解方面，由于算子分裂方法中源项求解在迭代过程中计算量差异很大，极易造成并行计算的负载不平衡。本文在算子分裂法运用的基础上，通过对矩方程源项特征及其与流体状态参量的对比分析，以温度梯度作为指示子进行区域划分，较好实现负载平衡。在此基础上计算不同当量下强爆炸光辐射辐照特征，研究光辐射极值时刻随当量的变化，并对光辐射相关物理参量的空间分布开展定量分析。

1 计算模型

灰体近似下强爆炸火球光辐射发展的一维球对称辐射流体力学方程组(满足局域热动平衡假定)[3,12,14]如下：

(1)

式中：ρ为空气密度，u为r方向空气速度，p为空气压强，e为空气体积总能量，κ为静止坐标系中考虑了受激发射的空气吸收系数，pr为光辐射压强，c为光速;φ=aT4为温度T下的黑体辐射能密度,a=7.56×10-16J/(m3·K4);Er为光辐射能密度，Fr为光辐射能流，即单位时间内通过垂直于传播方向单位面积上光辐射的能量，是表征光辐射强度的主要参量。采用最大熵Eddington因子下的辐射压强pr=χEr封闭方程组[3,14-15]， 空气的状态方程采用实际空气状态方程[5]。

对于计算中的初始条件，假定爆炸释放的能量集中于一定半径的等压球内(体积为V)，辐射能(ErV)与流体内能(eIV)之和为爆炸总能量(相当于爆炸当量所对应能量的85%，1 kt爆炸当量所对应的能量约为3.55×1012J)。

2 求解方法

图1 强爆炸过程中温度、吸收系数及源项的空间分布特征Fig.1 Spatial distribution of temperature, absorption coefficient and source terms in intense explosion

运用算子分裂法对上述方程进行求解，其核心是将方程分裂为对流项和刚性源项2步求解[12,14]：第1步求解对流项，由于不含相互作用项，求解精度和稳定性都较高；第2步化为常微分方程求解，根据方程特性可以进一步转化，从而实现快速稳定求解。

光辐射辐照特征的计算，需要发展高效并行算法。对于算子分裂方法求解来说，通常的区域划分方法会由于各个节点计算量的不同导致求解效率大大降低。通过对分裂求解的进一步分析，可以发现：对流项求解几乎不消耗计算资源，影响计算效率的主要因素在于源项求解中的迭代计算。以零阶矩源项M0(单位时间内辐射与流体能量交换项)为例，其表达式(忽略相对论修正项)为：

(2)

由于吸收系数κ也与温度T有关，所以上述源项的大小与温度直接相关，从而最终影响源项求解所需要的计算时间。对于光辐射发展过程而言，温度在空间分布上存在明显的非均匀性，且梯度差异明显，与温度有关的物理参量也呈现类似的分布特征，见图1。

由图1可见：温度梯度较大的区域，吸收系数与源项均有大的变化。从物理过程上来说，也就是在温度变化非常剧烈的区域(如波阵面附近)，辐射与流体强烈相互作用，导致源项极其巨大，最终造成该区域单个物理量的求解所需计算量要比其他区域大得多。针对这一特点，在并行求解中为较大程度的实现负载平衡，采用以温度梯度作为区域划分的标准：温度梯度越大，所求解的区域较小；反之则越大。另外，考虑到光辐射发展过程中，物理参量随时间不断变化，因此需要对温度进行动态监测，利用上述标准实现区域的自动划分。采用该方法对一维情况下光辐射过程进行求解，与均匀区域划分相比，当计算节点达到10时并行效率显著提升。

3 计算分析

爆炸高度取h=0.5 km，在kt～1 Mt当量范围内，求解相应的强爆炸光辐射输运过程，空气初始状态参量如表1所示，cs为声速，爆炸高度小于1 km时状态参量取h=0 km时的数值。

表1 不同高度下空气初始状态参数

图2给出了不同当量下光辐射强度(即火球半径处光辐射能流)随时间的变化曲线。从不同曲线的变化规律来看，几种当量下光辐射强度随时间变化趋势基本一致：均是先增大然后减小，到某一极小值后又逐渐增大，最后减小，即双脉冲波形。

图2 不同爆炸当量下光辐射强度随时间变化关系Fig.2 Relationship between intensity of thermal radiation and time at different yields

对于第1个脉冲来说，主要是由于辐射源不断向外辐射能量导致温度降低，从而引起光辐射强度下降，第2个脉冲的形成过程则与冲击波的作用密切相关。图3给出了1 kt爆炸情况下，第2个脉冲形成前、形成后典型时刻流体物理参量(密度、压力和温度)的空间分布特征。从图中可以看出：在第2个脉冲形成前，冲击波(密度、压力突变处)位于辐射波(温度突变处)之后，向外发出的光辐射不会受到冲击波的影响；而当冲击波赶上并超过辐射波之后，尽管内部辐射温度很高，但流体的吸收作用决定了光辐射强度的变化，最终导致了第2个脉冲的出现。

光辐射强度极值所对应的时刻与当量有关。爆炸当量越大，光辐射极值时刻不断增加。图4给出了不同当量下光辐射第1极小和第2极大所对应的时刻。在对数坐标下，极值时刻与当量有着近似线性的关系。通过拟合并经过相应的转化，得到：

(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的关系。不同文献中参数a有所不同，但参数b相差不大。文献[1]给出的光辐射极值时刻与当量的经验公式中，参数b分别为0.44和0.46，计算结果与文献结果基本一致。

对于光辐射第1极大时刻(见图5)，计算结果与经验公式相比出现较大不同。从计算的角度来讲，主要是由于光辐射第1极大时刻出现较早，该时刻空气实际温度很高(约百万K)，吸收系数的非线性影响很大，导致辐射能流的计算存在不稳定性。另外，根据文献[7]，光辐射第1极大时刻受到除当量外的多种因素影响(如弹体质量等)，经验公式仅是针对当量的近似估计。

图4 光辐射强度的极值时刻与爆炸当量的关系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

光辐射总功率是光辐射辐照特征的一个关键参量。若光辐射源(即火球)半径为rb，则光辐射向外辐射的总功率为：

(4)

可见，光辐射功率的计算需要首先确定火球半径。图6(a)给出了不同当量下火球半径随时间变化关系。进一步根据式(4)计算得到光辐射总功率的变化历程，见图6(b)。由图6(b)和图2可以看出，光辐射功率走时与其强度的变化趋势相似，但其极值时刻有所不同。其根本原因是光辐射源(火球)半径随时间的变化与当量有关。当量越大，相同时刻的火球半径越大。整体上来看，光辐射功率随当量的增大而不断增大。

图6 不同爆炸当量下火球半径及光辐射总功率随时间的变化Fig.6 Variation of fireball radius and thermal radiant power with as time at different yields

4 小 结

(1)采用温度梯度作为指示子进行并行区域的动态划分，在算子分裂方法基础上达到了较高效并行求解，实现了光辐射特征及其当量相关性的数值模拟。

(2)不同当量下光辐射均呈现双脉冲，其数值大小与当量有关。曲线拟合的结果表明：光辐射极小时刻和第二极大时刻分别与当量的0.443和0.451次方成正比。

(3)光辐射总功率走时与其强度走时有着类似的变化趋势，但随当量的变化关系有所不同，这主要是由于辐射源半径随当量变化所引起的。

[1] 乔登江.核爆炸物理概论[M].北京:国防工业出版社,2003:225-262.

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

[3] 田宙,乔登江,郭永辉.不同当量强爆炸早期火球现象的数值模拟[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] 田宙,乔登江,郭永辉.强爆炸早期火球现象的一维数值研究[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].西安:西北核技术研究所,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] 陈健华,王心正,谢龙生.均匀大气中的强爆炸一维辐射流体力学数值解[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] 闫凯,李若,田宙,等.强爆炸火球数值模拟中的一种算子分裂方法[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] 高银军,闫凯,田宙,等.强爆炸早期火球光辐射能谱的数值计算[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近似应用于强爆炸火球辐射输运数值模拟[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