闫 凯，刘 钰，田 宙（西北核技术研究所，陕西西安 710024）

低空强爆炸中火球的一维数值模拟研究

闫 凯，刘 钰，田 宙
（西北核技术研究所，陕西西安 710024）

摘要：对海平面环境下当量为1Mt TNT的强爆炸进行一维辐射流体力学数值计算，得到了较为完整的火球扩张过程数值模拟结果。数值给出了火球不同发展阶段物质参数和辐射参数的空间分布，该结果与实际强爆炸中一些经验公式的计算结果吻合较好。给出了火球阵面和辐射功率随时间的变化关系，较合理地解释了低空强爆炸火球二次极大现象。

关键词：火球；辐射流体力学；数值模拟；二次极大

Numerical Simulation of Fireball in Strong Explosion at Low Altitude

YAN Kai，LIU Yu，TIAN Zhou
（Northwest Institute of Nuclear Technology，Xi’an710024，China）

火球现象是强爆炸早期最显著的特征之一。根据维恩位移定律，强爆炸瞬间辐射波以高频波为主，约75%～85%的爆炸能量为X射线［1］。在低空环境下，这种X射线的平均自由程只有m量级。X射线在传播过程中加热周围空气，形成一团高温高压的气体，即火球。火球迅速向外膨胀，压缩周围空气形成冲击波，同时向外发出光辐射。

低空强爆炸火球的重要特征是二次极大现象，即火球光辐射会出现两个脉冲。该现象的形成和空气在被加热过程中引起的光学性质变化有关［2－4］，文献［2］指出，随着冲击波卷入的空气层增厚，空气的不透明度增加，火球的亮度逐渐下降，达到极小值。但当火球阵面温度降至2 000K以下，波后空气对辐射的吸收越来越弱，逐渐变为透明，里面的火球一方面逐渐冷却，另一方面却越来越看得清楚。目前国内尚未见到对该现象的数值模拟结果。

数值模拟低空强爆炸火球现象的主要困难有两点。首先是辐射输运模型问题，由于低空环境下温度对空气不透明度有重要的影响，爆炸初期火球边缘附近的空气不透明度变化极为剧烈，因此辐射输运模型需能同时描述光学厚和光学薄介质中的光子输运过程。其次是刚性问题，辐射流体方程组涉及3种重要的动力学尺度：物质波波速（声速）、辐射波波速（光速）及源项［5］。当辐射与物质相互作用较强时，源项代表最快的尺度；当辐射与物质相互作用较弱时，源项代表最慢的尺度。因此如何处理源项问题是提高程序计算效率的关键。针对上述两点困难，本文采用一种最大熵变Eddington因子的M1近似来描述辐射场［6－7］，并设计一种分裂格式来处理源项问题［8－9］。通过对当量为1Mt TNT的强爆炸火球在海平面大气中的数值模拟，计算火球参数在不同发展阶段的时空分布，给出火球二次极大现象的数值解释。

1　辐射流体方程

本文采用文献［10］给出的一维球对称Euler型辐射流体动力学方程组。其中，流体力学方程忽略重力和黏性，辐射输运方程采用灰体近似和局部热力学平衡近似。由于爆炸初期火球的扩张速度很快，因此必须考虑相对论效应的影响。常用的做法是方程组采用互动坐标系或混合坐标系。本文考虑采用混合坐标系的形式，将相对论修正项当成源项的一部分，这样辐射流体方程组可写成双曲守恒律形式，即：

其中

此时辐射输运方程是不封闭的，本文引入

Levermore等［6－7］构造的变Eddington因子χ，令：

其中：r为空间坐标；t为时间坐标；p、ρ、u、T分别为气体的压力、密度、速度和温度；e、eI分别为单位体积气体的总能量和内能为单位体积气体的辐射能密度；为辐射能流；为辐射压力；c为光速；φ为平衡辐射能密度

2　计算方法

辐射流体方程组中源项的尺度对方程组的计算效率影响很大，本文应用以下算子分裂方法将耦合项从方程组中分裂开来：

其中

式（4）为双曲守恒律方程组，式（5）为常微分方程组。对于双曲守恒律方程组，本文采用有限体积5阶WENO格式进行数值求解；由于常微分方程组具有强刚性的特点，一般的显式迭代法时间步长必须非常苛刻才能得到收敛解。本文设计一种高效的温度迭代方法，数值计算的具体过程可参考文献［9－10］。

3　数值计算结果及分析

计算选取当量为1Mt TNT的强爆炸，火球半径取为0．75m。采用文献［11－12］的做法，数值计算中将源区物质等价为空气，源区总质量为1 000kg。强爆炸源区设为等压球，强爆炸能量均匀分布在等压球内，等压球外为实际空气。

初始时刻，等压球处于热力学平衡状态，等压球内辐射能与气体能量之和等于爆炸总能量。等压球内外的密度为大气密度，速度为0，温度、压力、内能利用状态方程给出。辐射能密度利用温度求出，辐射能流为0，吸收系数的计算采用灰体近似法。

本文计算采用一维球坐标系，其中内边界采用对称边界条件。对流体部分，压力、密度、能量对称相同，速度对称相反；对辐射输运部分，辐射能、辐射压力对称相同，辐射能流对称相反。外边界则采用透射边界条件。

从时间上来说，一个完整的火球发展过程应包括X射线火球阶段、辐射扩张阶段、过渡阶段及冲击波扩张阶段（也有文献将过渡阶段划分到辐射扩张阶段）。这4个阶段是以火球扩张的主要作用力的不同而划分的。在火球发展的早期阶段，约在前1～2μs，称为X射线火球阶段。X射线火球阶段尚难以用灰体近似法来描述，必须采用多群方法。由于本文仅采用灰体近似法，基于这种考虑，本文未计算X射线火球阶段，而是将X射线加热空气形成的火球作为初始条件。本文数值模拟的初始条件可看作辐射扩张阶段的开始。由于强爆炸火球参数的实际测量较困难，国内、外公开发表的文献中给出的结果大多数是定性的。本文将火球发展不同阶段的计算结果与文献［2］进行定性对比。

3．1 辐射扩张阶段

图1　 辐射扩张阶段火球的物质参数Fig．1 Fireball’s material parameters in radiation transport period

图1为辐射扩张阶段火球的物质参数。在该阶段，火球扩张主要是由辐射输运造成的，其特点是火球以辐射热波的形式向外传播。文献［2］指出该阶段波阵面上温度、压力和速度均有跃变，且波后温度大体是常数，即等温球。图1中各物理参数在5μs以前的空间分布基本反映了这一特点。图2为火球辐射参数的时空分布。可看出，火球阵面附近辐射能密度急剧降低，这是由于爆炸初期源区的高频射线不断被周围空气吸收并加热空气的结果。辐射能流在火球阵面附近有突变，这是由火球吸收系数的变化引起的。虽然火球内部的温度很高，但是由于源区内的密度较大，火球内部的吸收系数（灰体近似条件下）仍较大，因此火球内部并不透明；在火球阵面附近，气体密度较低，火球在该区域吸收系数较小；但当火球温度处于20 000～100 000K之间时，火球的吸收系数又会突然增大。随着火球的扩张，火球温度和辐射能密度迅速降低，其中辐射能密度下降的更快。该物理量实际反映火球的辐射扩张能力，它的减小预示着流体力学过程的作用开始增加，火球开始进入过渡阶段。

3．2 过渡阶段

随着时间的推移，流体力学过程的作用不能再忽略。表现在火球阵面的温度、压力不断下降，火球阵面逐渐向强冲击波过渡，即所谓过渡阶段。该阶段火球的增长是辐射扩展和流体力学过程的共同作用。文献［2］将火球发展过程分为辐射扩张阶段和冲击波扩张阶段，本文采用文献［12］的观点，认为在辐射扩张阶段和冲击波扩张阶段之间存在一流体和辐射共同作用的阶段，即过渡阶段，该阶段对应于文献［2］中辐射扩张阶段后期。

过渡阶段流体力学过程不能忽略，表现在火球阵面附近的密度逐渐有跃变，且密度压缩比逐渐趋向经典激波的极限，同时弹壳冲击波（又称等温激波）开始形成。图3为过渡阶段火球的物质参数。从图3中的温度和压力曲线上可明显看出弹壳冲击波的形成过程。在弹壳冲击波的阵面上各物理量均有突变，波后的温度、压力基本是均匀的，而且在火球内部会产生多个二次激波，图3中30、100、220μs时火球中心温度和压力均有突变，二次激波对火球中心物质参数的影响较大，但对辐射参数的影响并不大，如图4所示。随着时间的推移，二次激波对火球阵面的影响越来越小，在220μs时对火球阵面的影响已非常微弱。根据文献［2］，形成冲击波时相应的波阵面速度约为62km／s，阵面温度约为180 000K，大致可认为220μs左右时火球开始进入冲击波扩张阶段。

图2　 辐射扩张阶段火球的辐射参数Fig．2 Fireball’s radiation parameters in radiation transport period

图3　 过渡阶段火球的物质参数Fig．3 Fireball’s material parameters in transition period

图4　 过渡阶段火球的辐射参数Fig．4 Fireball’s radiation parameters in transition period

3．3 冲击波扩张阶段

图5为冲击波扩张阶段火球的物质参数。当等温激波赶上辐射波后，即进入火球发展的强冲击波扩张阶段。当冲击波超过辐射波后，在温度空间剖面上出现由冲击波阵面温度Ts到辐射波阵面温度TI的阶梯，形成冲击波扩张阶段火球特有的结构，即火球内部在辐射作用下为等温球，火球边缘基本对应冲击波波峰位置。随着火球的进一步扩张，冲击波阵面的温度压力不断下降，随着冲击波波速不断下降，约在10ms时冲击波阵面密度压缩比达到顶峰，随后开始下降。

图6为冲击波扩张阶段火球的辐射参数。对比图5、6可发现，此阶段辐射能密度的阵面已落后于冲击波阵面，其阵面大致与等温球阵面一致。虽然此时辐射对冲击波传播的影响已减弱（表现在辐射能流迅速减小），但对火球现象仍要用辐射流体力学来描述。这是因为火球内部温度仍较高，吸收系数仍很大，因此辐射仍是传热的主要方式，火球内部基本处于辐射平衡状态。在等温球内部，密度、速度均是均匀分布的，且随着时间单调下降。火球表面与等温球表面壳层间，压力和密度明显上升，温度则明显下降。

图5　 冲击波扩张阶段火球的物质参数Fig．5 Fireball’s material parameters in shock transport period

图6　 冲击波扩张阶段火球的辐射参数Fig．6 Fireball’s radiation parameters in shock transport period

3．4 火球冲击波阵面与文献［2］结果的对比

根据实测冲击波参数的分析，文献［2］提出应用以下经验公式描述强爆炸冲击波的传播规律。

式中：Q为爆炸总当量，kt；Rd为冲击波波峰距爆心的距离，m；p0为压力的峰值，105Pa。

图7为冲击波扩张阶段后（此时冲击波波峰距爆心约50m）的压力峰值与经验公式的对比，可看出，两者差异较小，说明本文的算法能较好地描述强爆炸冲击波。

图7　 冲击波压力峰值随冲击波半径的变化Fig．7 Peak pressure of shock wave versus radius

同时文献［2］提到：在冲击波扩张阶段，如果对模拟解不同时刻用不同γ的理想气体来逼近真实大气，那么冲击波阵面R（t）可写为：

图8为本文计算得到的冲击波阵面与经验公式计算结果的对比。可看出，在爆炸早期，本文计算结果较经验公式计算的结果略高，这是由于本文采用的空气吸收系数与实际有差异引起的；在冲击波阵面扩张至400m后，本文计算结果较经验公式计算的结果略低，这是由于随着计算区域的扩大网格尺寸变大（导致冲击波波峰不够锋锐，冲击波速度变慢）引起的。总体而言，两者计算结果差异不大，说明本文采用的物理模型可描述火球冲击波的扩张过程。

图8　 冲击波阵面随时间的变化Fig．8 Shock wave front versus time

3．5 火球的光辐射参数

光辐射的源是火球，因此火球的辐射功率与远方各点接收到的光辐射的能量有关。辐射功率可表示为：

式中：4πR2f为火球的面积；B为全波亮度，在灰体近似的条件下，可将辐射能流近似认为全波亮度，即：

图9　 辐射功率随时间的变化Fig．9 Radiation power versus time

图9为本文计算的辐射功率与文献［11］计算结果的对比。文献［13］主要介绍了RADFLO程序对强爆炸火球的一些数值模拟结果，该程序采用多群方法，辐射方面只考虑辐射能流对物质内能的影响。本文计算得到的第一极大，第一极小和第二极大时间分别为0．002、0．095和0．84s。文献［2］给出λ≈0．65μm的红光亮度的第一极小和第二极大到来时间的经验公式为：

对于当量为1Mt的强爆炸，计算得到红光亮度的第一极小和第二极大时间约为0．08s 和0．62s。此结果较本文的计算结果略小，这是由于本文采用灰体近似法，导致吸收系数的计算不准确引起的。总体来看，本文的计算结果与文献［2］中的经验公式及RADFLO的计算结果差异不大，验证了本文方法的合理性。

由于物理模型、空气状态方程、吸收系数等处理的不同，本文与RADFLO程序得到的辐射功率有差异。此外，在火球第一极大前，两者的结果均有震荡，这是由于此时火球温度的变化极为剧烈，而温度又影响网格的光子平均自由程造成的，原则上网格尺寸要小于光子平均自由程才能得到合理的数值解，但由于计算效率的原因网格尺寸不可能满足该条件。在火球第一极大后，随着火球阵面温度的降低，辐射功率的变化较为光滑。

图10为火球、冲击波及等温球半径随时间的变化。可看出，冲击波形成后，火球表面以冲击波速度向外扩张，而等温球界面的扩张速度明显减慢，约在2s等温球消失。在0．13s左右，冲击波脱离火球，此时对应的冲击波半径为500m左右，文献［2］给出的冲击波脱离火球的半径和时间的估算公式为：

图10　 冲击波阵面、火球阵面及等温球阵面随时间的变化Fig．10 Shock wave front，fireball front and isothermal sphere front versus time

对于当量为1Mt的强爆炸，可计算出脱离半径约为524m，脱离时间约为0．12s，此结果与本文的计算结果基本符合。该时间基本对应火球辐射功率极小值时间。

4　结论

采用Euler坐标系下的能描述光学厚和光学薄介质中光子输运的辐射输运方程，连同流体力学方程构成辐射流体力学方程组。数值计算给出了强爆炸火球现象中火球阵面、弹壳冲击波的形成发展过程及其物理参数的分布规律。计算结果与文献结果符合较好，证明了本文所采用物理模型、数值方法及计算程序的合理性。最后给出了低空强爆炸火球光辐射参数的长时间计算结果，计算结果表明：强爆炸初期，由于源区物质不透明度较大，X射线基本被吸收，辐射功率较小；随着强冲击波的形成，火球阵面附近的不透明度逐渐减小，辐射功率达到第一极大点；随着冲击波卷入空气的增厚，由于空气的不透明度（灰体近似）在20 000～100 000K时达到最大，因此空气又变得不透明，辐射功率达到了第一极小点；冲击波脱离火球后，火球阵面本身发光已微弱，且对波后辐射的吸收也逐渐微弱，火球外面有透明层出现，此时里面的火球逐渐露出，辐射功率达到第二极大点。

参考文献：

［1］ 欧阳建明，马燕云，邵福球，等．高空核爆炸X射线电离的时空分布数值模拟［J］．物理学报，2012，61（24）：242801．OUYANG Jianming，MA Yanyun，SHAO Fuqiu，et al．Numerical simulation of X－ray ionization and atmospheric temporal evolutions with high－altitude nuclear explosions［J］．Acta Phys Sin，2012，61（24）：242801（in Chinese）．

［2］ 乔登江．核爆炸物理概论［M］．北京：国防工业出版，2003：169－262．

［3］ 吴健辉．核爆炸光辐射特性及探测技术的理论与实验研究［D］．武汉：华中科技大学，2009．

［4］ GOLDBLAT J，COX D．Nuclear weapon tests：Prohibition or limitation［M］．ENG：Oxford University Press，1988：10－18．

［5］ SEKORA M，STONE J．A higher order Go－dunov method for radiation hydrodynamics：Ra－Computational Physics，2013，30（3）：379－386 diation subsystem［J］．Communications in Ap－（in Chinese）．plied Mathematics and Computational Science，

［10］SEKORA M．Algorithms for hyperbolic balance 2009，4（1）：135－152．laws with multiscale behavior：Applications in

［6］ LEVERMORE C D．Relating Eddington factorsradiation hydrodynamics［D］．United States：to flux limiters［J］．Journal of Quantitative Spec－Princeton University，2010．troscopy and Radiative Transfer，1984，31（2）：

［11］田宙，乔登江，郭永辉．不同高度强爆炸早期火149－160．球数值研究［J］．兵工学报，2009，30（8）：1 078－

［7］ ANILE A M，PENNISI S，SAMMARTINO M．1 084．A thermodynamical approach to Eddington fac－TIAN Zhou，QIAO Dengjiang，GUO Yonghui．tors［J］．Journal of Mathematical Physics，1991，Numerical investigation of early fireball of strong 32（2）：544－555．explosion for different altitudes［J］．Acta Arma－

［8］ LOWRIE R B，MOREL J E．Issues with high－mentarii，2009，30（8）：1 078－1 084（in Chiresolution Godunov methods for radiation hydro－nese）．dynamics［J］．Journal of Quantitative Spectrosco－

［12］BRODE H L．Thermal radiation phenomena，py and Radiative Transfer，2001，69（4）：475－VolumeⅤ：Radiation hydrodynamics of high 489．temperature air，AD672837［R］．United States：

［9］ 闫凯，李若，田宙，等．强爆炸火球数值模拟中的Lockheed Missiles and Space Company Sunny－算子分裂方法［J］．计算物理，2013，30（3）：379－vale，1967．386．

［13］EUGENE M D，JOHN Z，RODNEY W．RADYAN Kai，LI Ruo，TIAN Zhou，et al．An oper－FLO physics and algorithms，LA－12988－MS［R］．ator splitting method for numerical simulationofUnited States：Los Alamos National Laboratory，strong explosion fireball［J］．Chinese Journal of1995．

作者简介：闫 凯（1984—），男，河南民权人，助理研究员，硕士，从事辐射流体力学及磁流体力学数值模拟研究Abstract：A one－dimensional radiation hydrodynamics numerical computation was presented for a strong explosion at sea level，of which the equivalent is 1 Mt TNT．The spatial distributions of material and radiation parameters in strong explosion，which reflect a integrate process of fireball evolution，were given in terms of the simulation．The calculation results agree well with those of some empirical formulas．The results of optical observables－radiation power versus time and fireball radius versus time were described，and the second maximum phenomenon of the fireball at low altitude was given．Key words：fireball；radiation hydrodynamics；numerical simulation；second maximum

基金项目：国家自然科学基金资助项目（91330205）

收稿日期：2014－04－18；修回日期：2014－07－24

doi：10．7538／yzk．2015．49．08．1345

文章编号：1000－6931（2015）08－1345－09

文献标志码：A

中图分类号：O38；O242