基于中尺度天气模型的核爆炸烟云演化模型

The Nuclear Explosion Cloud Evolution Model Based on the Mesoscale Weather Model

摘要

核爆炸是在有限体积内瞬时间释放出大量能量的过程，释放能量很大，杀伤力很强。核爆炸过程中和核爆炸过程后，会产生多种物理效应向远离爆炸点区域传播。核爆炸产生的烟云既会污染空气，又会造成地面大范围的沾染。因此研究核爆炸烟云的运动规律是非常必要的。本文开展了对核爆炸烟云的模拟研究，在天气研究和预报模型中加入火球初始化，模拟了不同爆炸当量、不同爆炸高度、不同气象条件的核爆炸产生的烟云的演化过程，从垂直方向和水平方向全面演示，并模拟了流场随时间的变化。模拟实验结果表明，核爆炸产生的烟云演化的过程符合相关理论知识，基本达到实验预期。本研究为了解核爆炸烟云的运动规律提供了依据。

Abstract

A nuclear explosion is a process in which a large amount of energy is released instantly in a finite volume. Various physical effects can be produced during and after a nuclear explosion, spreading far away from the explosion area. The nuclear cloud generated by nuclear explosions can pollute the air and cause large-scale ground contamination. Therefore, it is indispensable to study the motion law of nuclear clouds. This paper mainly studies the formation of nuclear clouds from nuclear explosions. Fireball initialization is added to the weather research and forecast model to simulate the evolution process of nuclear explosion smoke clouds with different explosion yields, height and meteorological conditions. In this paper, the vertical and horizontal directions are fully demonstrated, and the change of flow field with time is simulated. The experimental results show that the process of nuclear cloud evolution of the nuclear explosion conforms to the relevant theoretical knowledge, and the experiments are basically desired. Overall, this study provides a basis for improving the capability of nuclear weapons capability and understanding the movement of nuclear cloud.

Graphical abstract

关键词

中尺度天气模型 / 天气研究和预报模型 / 核爆炸 / 烟云演化

Key words

the mesoscale weather model / weather research and forecasting model / nuclear explosion / nuclear cloud evolution

0　 引 言

自第二次世界大战后，研究核爆炸物理现象及其规律受到了人们广泛的重视，得到了相当充分的发展，取得了大量的成果，但是由于专业交流的限制，更主要的是由于核爆炸现象本身的复杂性，还有不少问题停留在经验公式上，需要深入讨论核爆炸现象，探讨杀伤破坏效应的规律，为研究核爆炸效应及防护提供理论依据^[1]。

核武器是通过能自持进行的核裂变或聚变反应，瞬间释放巨大能量，产生爆炸作用，并且具有大规模杀伤破坏效应的武器的总称^[2]。核爆炸是在有限体积内瞬时释放出大量能量的过程，释放的能量比常规炸药爆炸更大，能量密度也更高^[3]。核武器技术的发展会带动核试验技术的发展，核试验技术的发展又会反过来促进核武器技术的发展。继续研究核试验技术及其相关技术研究的理论与实践成果，不仅是禁核武器试验形势下培养核试验人才的需要，也是保留和发展核试验能力的需要，并且也会对新形势下研究和发展新的实验技术起到重要的指导作用^[1]。

大气层核爆炸产生的烟云承载了核裂变以后的绝大部分剩余放射性产物，其在大气的作用下向下风方向飘散，既污染了空气，又造成了地面大范围的沾染，因此核爆炸烟云的运动规律不可忽视^[4]。Morton等^[5]在1956年进行了一系列实验，完成了烟云上升的流体力学方面的开创性理论工作。Caulfield和Woods^[6]在1998年将Morton等人的理论推广为包含一个非线性分层的背景密度场，并且之后在理论中加入了化学反应。Onufriev^[7]提出了核爆炸烟云上升的涡环运动理论。Chang^[8]采用了类似的方法模拟大的污染云内部的流场。Kansa^[9]发展了一个爆炸源的烟云模式。Makhviladze等^[10]、Krispin等^[11]和Brown等^[12]采用多相流方法模拟了爆炸烟云的运动过程，但由于计算量过大，他们的模式只限于二维模拟。

以前的研究分别使用了分析模型、高分辨率计算流体动力学（Computational Fluid Dynamics,CFD）代码和历史核试验数据，以建立对烟云上升流体动力学的理解^[13]。例如，Kanarska等^[14]在2009年结合了可压缩GEODYN代码和基于Bell等^[15]的低马赫数代码模拟烟云上升过程的各个阶段。Song等^[16]与Won和Lee^[17]采用了类似的方法，模拟了几次历史核试验期间的烟云上升，模拟结果与观测结果吻合良好。

虽然上述理论模型揭示了烟云上升流体动力学的重要细节，但它们在核响应建模中的应用可能受到限制。因此，需要一种高保真度的方法来理解和预测在更复杂、更现实的大气条件下的核烟云上升过程。天气研究和预报模型（Weather Research and Forecasting，WRF）^[13]是一种常用的数值天气预报模型，也能够进行大涡模拟，非常适合这一目的。

本论文针对核爆炸场景，开展了对核爆炸多物理场效应的模拟研究，使用WRF模型模拟了理想大气中核爆炸烟云的上升过程。

1　 WRF烟云模型

1.1　核爆炸烟云产生过程

核爆炸烟云是核爆后冷却了的火球以及卷入的周围环境物质组成。当爆炸高度高，卷入的地面土壤或岩石少时，烟云卷入的主要是水汽，烟云颜色较浅；当爆炸高度低，火球接触地面，烟云颜色为灰黑色^[3]。

核爆炸瞬间在有限空间内释放出巨大能量，形成了高温、高压的火球。火球急剧膨胀、扩大，迅速上升，体积增大和温度降低，逐渐变成一团烟云。烟云内部发生如图1所示的激烈的内循环运动。

由于这种激烈的内循环运动，加上烟云的迅速上升，会在爆心投影点附近产生一股强烈的、向上抽吸的上升气流，使更多的空气从环底吸入，并将地面掀起的尘土、碎石卷进去，形成从地面升起的尘柱。尘柱带来了大量地面物质，熔融后与放射性裂变产物相混形成了颗粒较大的灰粒子，并较快地沉降在爆区和近云迹区的地面上。

烟云上升到一定高度后，就不再上升，此时的烟云称为稳定烟云。随后，烟云边缘起毛，轮底开始模糊，并逐渐被风吹散^[18]。稳定时烟云的高度与当量有关，表1展示了爆炸当量为0.3、0.5和1.0 Mt的核爆炸烟云的稳定顶高。

1.2　在WRF中包含火球的模型

1.2.1　火球近似

WRF模型可以生成具有理想初始化的模拟和利用真实数据的模拟。本文使用的是理想初始化模拟。

WRF版本4解决了非静液压可压缩的纳维⁃斯托克斯方程，但不能模拟高马赫数爆轰动力学。因此，初始爆炸近似为一个高温球体，假设爆轰中1/3的能量用于温度扰动。这对于表面相互作用最小的空气爆发来说是一个很好的近似值^[19，20]，但不适合普遍应用，特别是近地表爆炸。Kanarska等^[14]使用了类似的火球近似，温度扰动约为1 300 K，半径为400 m，用于与明确捕捉到早期爆炸动力学的模拟进行比较。

Molenkamp^[20]开发的迭代程序中，根据核装置产量估计了适当的火球温度扰动\mathrm{\Delta }\theta和半径R。由于\mathrm{\Delta }\theta和R都是自由参数，我们先验地选择\mathrm{\Delta }\theta，并使用下面描述的过程估计相应的R值。考虑火球的温度场，首先计算：

\theta ={\theta }_{\mathrm{b}\mathrm{g}}+\frac{\mathrm{\Delta }\theta }{\mathrm{2}}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }r}{R}\right)+\mathrm{1}\right]（1）

其中，{\theta }_{\mathrm{b}\mathrm{g}}为三维背景电位温度场(在爆炸时取自 WRF)，r为火球中心的径向距离，即距离火球中心的距离。然后，计算由于火球存在而增加的热能为：

∆E=c_{\mathrm{p}}∆TM（2）

\mathrm{\Delta }TM={\int }_V \left(\theta -{\theta }_{\mathrm{b}\mathrm{g}}\right){\left(\frac{p}{p_{\mathrm{0}}}\right)}^{\frac{r_{\mathrm{d}}}{c_{\mathrm{p}}}}\frac{p}{r_{\mathrm{d}}\theta {\left(\frac{p}{p_{\mathrm{0}}}\right)}^{\frac{r_{\mathrm{d}}}{c_{\mathrm{p}}}}}\mathrm{d}V（3）

其中，r_{\mathrm{d}}为干气普适气体常数，c_{\mathrm{p}}为恒定压强下空气的比热，p_{\mathrm{0}}是参考海平面压强（p_{\mathrm{0}}=\mathrm{1.01}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{P}\mathrm{a}），（3）式是对体积积分在WRF模型网格上的离散计算。火球半径R在方程（1）中迭代更新直到（2）式中\mathrm{\Delta }E\approx Y/\mathrm{3}，其中Y(\mathrm{J})为器件产量^[4]。

1.2.2　WRF中的火球初始化

由于WRF采用了一种基于地形跟随压强的混合垂直坐标，该坐标由一个流体静力平衡的参考状态和扰动场进行初始化，因此有必要修改WRF的初始化例程以包含较大的火球温度扰动。WRF可以使用用户定义的数据（例如输入探测）进行初始化，以实现理想化的仿真。

下面对WRF初始化程序进行概述，并介绍包含高温火球所需的修改。

在本文中，WRF的初始化例程被修改为包括一个高温火球。在理想WRF仿真中，大气状态通常由用户定义的一维输入测深计算，其中包括在指定高度的位势温度{\theta }_m（\mathrm{K}）、水汽混合比q_{\mathrm{v}}（\mathrm{g}/\mathrm{k}\mathrm{g}）、水平分速度u(\mathrm{m}/\mathrm{s})\mathrm{和垂}\mathrm{直分}\mathrm{速度}v（\mathrm{m}/\mathrm{s}），以及压强（\mathrm{h}\mathrm{P}\mathrm{a}）、温度（\mathrm{K}）和水蒸气混合比（\mathrm{g}/\mathrm{k}\mathrm{g}）。WRF将热力学变量分为参考状态和扰动状态。在理想情况下，参考状态是由干测深计算，而扰动状态是去除参考状态的湿状态。

由于WRF中的垂直坐标既跟随地形，又基于压强，因此当大气柱膨胀时，垂直坐标水平上升，当大气柱收缩时，垂直坐标水平下降，例如随着白天的加热和夜晚的冷却。垂直坐标由变量\eta给出，其范围从地形表面的1到模型顶部的0，并根据模型输入指定，对应的干空气静水压强为：

p_{\mathrm{d}}\left(x,y,\eta \right)=B\left(\eta \right)\left(p_{\mathrm{s}}\left(x,y\right)-p_{\mathrm{t}}\right)+

[\eta -B\left(\eta \right)]\left(p_{\mathrm{0}}-p_{\mathrm{t}}\right)+p_{\mathrm{t}}（4）

其中，p_{\mathrm{s}}\left(x,y\right)和p_{\mathrm{t}}表示地表和模型顶部的干静水压强。系数B是混合坐标的权重因子，因此当B\left(\eta \right)=\eta时，坐标是跟随地形的，当B\left(\eta \right)=\mathrm{0}时，坐标是基于静水压强。在WRF模型的控制方程中出现了垂直坐标变换的度量项{\mu }_{\mathrm{d}}，定义为：

{\mu }_{\mathrm{d}}\left(x,y,\eta \right)=\frac{\partial p_{\mathrm{d}}\left(x,y,\eta \right)}{\partial \eta }=\frac{\partial B\left(\eta \right)}{\partial \eta }{p}_{\mathrm{c}}\left(x,y\right)+

\left[\mathrm{1}-\frac{\partial B\left(\eta \right)}{\partial \eta }\right]\left(p_{\mathrm{0}}-p_{\mathrm{t}}\right)（5）

其中，p_{\mathrm{c}}\left(x,y\right)=p_{\mathrm{s}}\left(x,y\right)-p_{\mathrm{t}}是大气的柱压。注意，在WRF版本4中，垂直坐标度量从二维形式{\mu }_{\mathrm{d}}\left(x,y\right)变成了三维形式{\mu }_{\mathrm{d}}\left(x,y,\eta \right)。这个对应WRF中代码为：

{\mu }_{\mathrm{d}}\left(x,y,\eta \right)=C_{\mathrm{1}}\left(\eta \right){\mu }_{\mathrm{d}}\left(x,y\right)+C_{\mathrm{2}}\left(\eta \right)（6）

其中，C_{\mathrm{1}}(\eta )=\frac{\partial B\left(\eta \right)}{\partial \eta }和C_{\mathrm{2}}\left(\eta \right)=[\mathrm{1}-\partial B\left(\eta \right)/\partial \eta ](p_{\mathrm{0}}-p_{\mathrm{t}})。在接下来的讨论中，柱压强p_{\mathrm{c}}(x,y)对应于WRF代码中的二维变量{\mu }_{\mathrm{d}}(x,y)。

由于从模型输入中已知坐标水平\eta和加权项B，因此WRF初始化的第一步是使用（4）和（5）式确定干静水压强p_{\mathrm{d}}和柱压强p_{\mathrm{c}}。这就需要确定地表的静水压强p_{\mathrm{s}}和来自输入测深的模型顶部p_{\mathrm{t}}。第二步是将其用于模型方程（以及来自输入测深的信息）以初始化密度和位势。采用压强作为自变量，将测深变量({\theta }_m,q_{\mathrm{v}},u,v)插值到模型水平上。

为了确定模型顶部和地表的压强值，WRF集成了流体静力学方程

\frac{\partial p}{\partial z}=-{\rho }_{\mathrm{d}}g\left(\mathrm{1}+q_{\mathrm{v}}\right)（7）

通过探测的垂直大气柱（即这种积分不是在模型水平上进行的），在垂直积分过程中，利用状态方程迭代得到干密度{\rho }_{\mathrm{d}}，该状态方程还包含了输入测深的位温。首先以输入测深的海平面压强和温度作为底部边界条件，向上进行垂直积分。干大气的第二次积分(即q_{\mathrm{v}}=\mathrm{0})是向下进行垂直积分，使用第一次积分的压强作为顶部边界条件。如上所述，来自这些计算的压强值与来自输入测深的变量一起使用可以初始化WRF参考（基数）和扰动状态。

当火球存在时，火球内空气密度降低，输入测深所得的海平面压强不能再作为静水力方程（7）的底边界条件，之前的方法不再成立，因此，将火球温度扰动场叠加到输入测深的位温场上，建立了一种新的位温三维测深方法，如（1）式所示。首先完成未修改的干测深和湿测深的计算，并将参考状态定义为标准WRF初始化。然后，对新三维场的每一列进行独立的垂直积分，利用流体静力学方程（7）求解压强。在垂直积分中，积分方向是反向的，从输入测深顶部开始，以未修改的输入测深所得的湿压强作为顶边界条件。这是假设在探测顶部的大气状态不受火球的增加的影响。湿压和干压都是通过大气柱向下积分计算的。利用这个加入火球后的压强场，WRF模型就可以继续计算柱压、密度和位势的扰动值，并可以将输入测深的温度、湿度和速度值插值到模型水平。

理想WRF仿真中的水汽混合比q_v和水蒸气混合比是湿空气中的水汽质量与干空气质量之比，公式如下：

q_{\mathrm{v}}=\mathrm{622}\frac{e}{p-e}\left(\mathrm{g}/\mathrm{k}\mathrm{g}\right)（8）

式中，p/\mathrm{h}\mathrm{P}\mathrm{a}为气压，e/\mathrm{h}\mathrm{P}\mathrm{a}为水汽压，q_{\mathrm{v}}的单位为\mathrm{g}/\mathrm{k}\mathrm{g}。水汽压e是指湿空气中汽态水本身的压强，它是空气中水汽所产生的分压强。当空气饱和时便是饱和水汽压。水蒸气混合比和水汽比区别在于水蒸气混合比是参考状态，与高度无关，是地表初值，由比湿求得；水汽混合比q_{\mathrm{v}}是扰动状态，与高度有关，由相对湿度求得。比湿（specific humidity）是指湿空气中水汽的质量与该团空气总质量（水汽质量加上干空气质量）的比值，本文以\mathrm{g}/\mathrm{k}\mathrm{g}为单位，通常大气中比湿都小于40 \mathrm{g}/\mathrm{k}\mathrm{g}，比湿记录大气水汽状况的指标，它的计算公式如下：

q=m_{\mathrm{w}}/\left(m_{\mathrm{w}}+m_{\mathrm{d}}\right)（9）

式中，m_{\mathrm{w}}为该团湿空气中水汽的质量；m_{\mathrm{d}}为该团湿空气中干空气的质量，q的单位为\mathrm{g}/\mathrm{k}\mathrm{g}。

根据（9）式和气体状态方程可得：

q=\mathrm{622}\times \frac{e}{p}（10）

式中，p/\mathrm{h}\mathrm{P}\mathrm{a}为气压，e/\mathrm{h}\mathrm{P}\mathrm{a}为水汽压。

在WRF中可以得到地表的比湿q和气压p，因此可以通过比湿q和气压p，求得对应的水汽压，从而求得地表的水蒸气混合比。

相对湿度（relative humidity）是指空气中的实际水汽压与同温度下的饱和水汽压的比值:

f=\frac{e}{E}\times \mathrm{100}\mathrm{\%}（11）

式中，e为水汽压，单位h\mathrm{P}\mathrm{a}；E为饱和水汽压，是水汽达到饱和时的水汽压强，单位为h\mathrm{P}\mathrm{a}，是一个与温度有关的函数，其经验计算公式为：

\mathrm{l}\mathrm{n}E=\mathrm{53.679}\mathrm{57}-\frac{\mathrm{6}\mathrm{743.769}}{T}-\mathrm{4.845}\mathrm{1}\mathrm{l}\mathrm{n}T（12）

其中，T的单位是绝对温标\mathrm{K},与摄氏温度t（℃）的关系是：

T=\mathrm{273.15}+t（13）

在WRF中可以得到温度T和相对湿度f，因此可以先求饱和水气压E，再求水汽混合比q_v。

综上所述，将火球位势温度叠加在环境位势温度上，并进行一系列的计算来平衡流体静力参考状态和更新扰动场^[4]即可完成理想的WRF的初始化工作。

2　 实验系统及测量结果

2.1　实验数据

本文利用天气研究和预报模型模拟了核爆炸烟云上升的过程以及大气扰动引起空气对流的流场随时间的变化。气象数据是从国家环境预测中心/国家大气研究中心再分析项目的再分析数据下载的。此全球数据集从1948年起，每6 h可以在水平网格间距为209 km、垂直高度为28个的网格上获取，下载的数据格式为netcdf4，经度范围是0°~ 360°，纬度范围是-\mathrm{90}°~90°（北纬为正），可以在linux系统里用python读取。WRF模型中按气压分层，并会给出每一层对应的高度和地表的初始状态值作为输入数据。

2.2　实验结果与分析

本次模拟选择了2016年3月1日的数据，其地表的表面压强为985.72 hPa，位势温度为285.76 K,水汽混合比为10 \mathrm{g}/\mathrm{k}\mathrm{g}，高度范围为0~20 km，均保留6位小数。模拟的火球当量为0.3、0.5和1.0 Mt，对应的扰动温度分别为80、100、200 K，火球半径为200 m，爆高为10、500和1 000 m，风速分别为0、15和20 m/s，分别从垂直方向和水平方向展示烟云上升过程，并模拟了垂直方向和水平方向的流场随时间的变化。本节的实验中相对湿度设为0，即不考虑空气中的水汽对热源上升的影响。垂直方向的模拟图中，x轴表示水平方向距离，单位为m，左边y轴是用气压表示垂直方向距离，因为WRF模型是以气压分层，气压随着高度增加而减小，右边y轴表示高度，颜色表示温度，图中的黑色箭头表示流场；水平方向的模拟图是取对应时间的垂直方向模拟图中能量最大的一层展示，并标注了所在层的压强，图中的黑色箭头表示流场。

图2表示当量为0.3 Mt、高度为10 m、风速为0 m/s的烟云上升过程。风速为0 m/s时，仿真的烟云垂直上升。垂直方向上爆炸瞬间，火球急剧上升和膨胀，逐渐变成一团烟云，并且在约1 min时上升速度最大。随着时间往后，烟云上升，上升速度逐渐减小，烟云逐渐扩散，温度也越来越低，到达一定高度，烟云停止上升，之后颗粒慢慢沉降。当量为0.3 Mt的烟云稳定高度约15 000 m。在到达稳定高度之前烟云上方的流场是辐散的，烟云下方的流场是辐合的。在烟云上升过程中，烟云顶的流场更加密集；因此，温度主要集中在烟云顶。在到达稳定高度后，随着烟云沉降，流场密集程度差异不大，因此温度差异也不是很大。在水平方向上，烟云能量最大的一层中心保持不变，在到达稳定高度之前，烟云顶部逐渐扩大，但整体保持一个圆形，到达稳定高度之后，烟云开始沉降，逐渐扩散满整个区域。

图3表示当量为0.3 Mt、高度为10 m、风速为15 m/s的烟云上升过程。风速为15 m/s时，仿真的烟云在上升的同时受风速影响移动。垂直方向上，以900 m/min向右移动，超出右边界后从左边界出现继续显示。上升的趋势和风速为0 m/s时相同。在到达稳定高度之前，温度主要集中在烟云顶，在到达稳定高度后，温度差异不是很大。在水平方向上会受x轴和y轴两个方向的风速影响从而向斜对角移动，x轴和y轴的移动速度均为900 m/min，同样超出边界后从另一边继续显示。在到达稳定高度前，烟云顶部逐渐扩大，但整体保持一个圆形，到达稳定高度之后，烟云开始沉降，逐渐扩散满整个区域。

图S1（见附加材料）表示当量为0.3 Mt、高度为10 m、风速为20 m/s的烟云上升过程。风速为20 m/s时，仿真的烟云在上升的同时受风速影响移动。垂直方向上，以1 200 m/min向右移动，超出右边界后从左边界出现继续显示。上升的趋势和风速为0 m/s时相同。在水平方向上，受x轴和y轴两个方向的风速影响，向斜对角移动；在x轴和y轴的移动速度均为1 200 m/min，同样超出边界后从另一边继续显示，变化趋势和风速为15 m/s时相同。

图S2（见附加材料）表示当量为0.5 Mt、高度为10 m、风速为0 m/s的烟云上升过程。风速为0 m/s时，仿真的烟云垂直上升。垂直方向上，上升的趋势和当量为0.3 Mt相同。当量为0.5 Mt的烟云稳定高度约17 000 m。在水平方向上，变化趋势和当量为0.3 Mt相同。

图S3（见附加材料）表示当量为0.5 Mt、高度为500 m、风速为0 m/s的烟云上升过程。风速为0 m/s时，仿真的烟云垂直上升。垂直方向上，上升的趋势和爆炸高度为10 m时相同。在本文模拟的爆炸高度范围内，相同当量，不同爆炸高度的烟云的稳定高度相同，并且上升趋势相同，只是起始爆炸高度不同。在水平方向上，变化趋势和当量为0.3 Mt相同。

图S4（见附加材料）表示当量为0.5 Mt、高度为1 000 m、风速为0 m/s的烟云上升过程。风速为0 m/s时，仿真的烟云垂直上升。垂直方向上，上升的趋势和爆炸高度为10、500 m时相同。本文模拟的爆炸高度为10、500和1 000 m，在此范围内，相同当量，不同爆炸高度的烟云的稳定高度相同，并且上升趋势相同，只是起始爆炸高度不同。在水平方向上，趋势和当量为0.3 Mt相同。

图S5（见附加材料）表示当量为1.0 Mt、高度为10 m、风速为0 m/s的烟云上升过程。风速为0 m/s时，仿真的烟云垂直上升。垂直方向上，上升的趋势和当量为0.3、0.5 Mt时相同。当量为1.0 Mt的烟云稳定高度约18 700 m。爆炸当量对烟云的稳定高度有影响，并且爆炸当量越大，爆炸形成的烟云、尘柱越大；烟云上升高度越高，烟云消散的越慢。在水平方向上，趋势和当量为0.3、0.5 Mt相同。

3　 讨 论

图2、3、图S1~图S5展示了不同爆炸当量、不同爆炸高度、不同气象条件的核爆炸形成的烟云的上升过程，表2对其进行了总结对比。从中可以得出以下结论：爆炸当量、爆炸高度和气象条件对烟云有影响，并且WRF模型模拟的烟云上升过程与1.1节中描述的核爆炸烟云产生过程符合。

根据表1，爆炸当量对烟云的稳定顶高有影响，同时会使爆炸形成的烟云、尘柱更大，烟云上升高度更高，烟云消散的更慢。

本文模拟的爆炸高度不影响烟云稳定高度，爆炸高度分别为10、500和1 000 m时，烟云上升趋势是一致的，只是起始高度不同。

不考虑风速时，烟云垂直上升，水平方向中心点保持不变。考虑风速时，烟云在上升过程中受风速影响移动，当风速为15 m/s时，烟云在垂直方向模拟中以900 m/min的速度向右移动；当风速为20 m/s时，烟云在垂直方向模拟中以1 200 m/min的速度向右移动，超出右边界后从左边界出现继续显示，在水平方向受x轴和y轴两个方向的风速影响，向斜对角移动，同样超出边界后会从另一边继续显示。因此，加入风速会使烟云飘散范围更广，造成更大面积的沾染。

从流场看，流场可以通过密集程度来展示温度的变化。在烟云上升过程中，温度主要集中在烟云顶，在到达稳定高度后，温度差异不是很大。

4　 验 证

将本文的研究结果与郑毅等^[3]利用RAMS系统模拟的烟云运动趋势进行对比。RAMS的模拟结果如图4所示。模式的计算范围为20 km×20 km×20 km，网格为101×101×101，网格距为200 m，设火球高度为800 m，火球半径为500 m，扰动温度为100 K，分别模拟了垂直方向标量场及流场随时空的变化和水平方向涡度场随时间变化。本文用WRF模型模拟了火球半径为500 m、扰动温度为100 K、高度为800 m、风速为0 m/s的烟云上升过程，如图5。

从图4可以看出，当在大气中瞬时释放热量后，加热后的空气迅速上升，RAMS系统模拟的烟云在约3 min后上升速度达到最大，之后上升速度逐渐减小。在到达最大高度之前烟团上方流场辐散，烟团下方流场是辐合的。水平方向上，稳定层结的涡度场在平衡高度上下振荡，逐步波浪式沿水平方向向外传播。用WRF模型模拟的烟云扰动温度100 K，爆炸高度为800 m，在约1 min左右上升速度达到最大，之后上升速度逐渐减小。在到达最大高度之前烟团上方流场辐散，烟团下方流场是辐合的。在水平方向上，在到达稳定高度之前，烟云顶部逐渐扩大，但整体保持一个圆形。两者相比，模拟的烟云上升趋势一致，但本文不仅模拟了烟云垂直方向上的运动过程及流场时空分布和烟云水平方向上的运动过程和流场的时空分布，还模拟了烟云随温度和压强的变化，并且对比了不同爆炸当量、不同爆炸高度、不同气象条件的烟云运动过程，更加全面形象。此外，WRF模型可以选择爆炸时间、爆炸经纬度，大气环境可以考虑水汽混合比、风速等条件，模拟的理想大气中的爆炸背景更接近真实爆炸的背景，考虑的参数更详细全面。

AD报告《DNA 1251-1-EX_Compilation_Fallout_Data_CONUS》对1945年7月16日的数据进行了模拟，爆高为100英尺（约30.5 m），爆炸当量为19 kt，地理坐标为\mathrm{33}°\mathrm{40}\text{'}\mathrm{31}\text{'}\text{'}\mathrm{N}，\mathrm{106}°\mathrm{28}\text{'}\mathrm{29}\text{'}\text{'}\mathrm{W}，云顶高度为35 000英尺（约10 668 m）。报告给出了烟云的运动轨迹，如图6所示（图中MLS为平均海拔高度）。

图6中，烟云随着时间的增加而迅速上升，随着时间继续增加，烟云继续上升，并且上升速度不断减小，大约10 min后，烟云继续上升，上升速度又逐渐增大，直至达到稳定高度，烟云高度处于平衡状态。用与图6相同的条件，本文在WRF1模型中模拟烟云上升的过程，每一分钟保留一张图片，共计算了0~11 min的结果，部分结果如图7所示。

取每一分钟对应的烟云顶层高度，绘制出烟云顶高随时间的变化，如图8所示。图8中，烟云随着时间的增加而迅速上升，随着时间继续增加，烟云继续上升，但上升速度逐渐减小，11 min后，达到稳定高度，烟云高度处于平衡状态，烟云顶高最大值约33 500英尺，约10 210 m。因此，图8和图6趋势一致，量级一致。所以，本文使用的WRF模型模拟的结果真实性可靠。

5　 结 语

本文对核爆炸多物理场效应进行了模拟，用WRF模型模拟了不同爆炸当量、爆炸高度和气象条件的烟云上升过程。本文主要做了如下工作：

1） 对关于研究烟云上升过程的方法进行了总结。

2） 详细介绍了核爆炸中烟云上升的理论过程，并阐述了用WRF模型模拟烟云上升的基本原理和方法。在WRF模型模拟烟云上升原理中，介绍了火球近似模型以及在WRF中包含火球的模型开发。

3） 基于前人的研究，对核爆炸产生的烟云的上升过程进行了研究，并分析对比了不同爆炸当量、爆炸高度和气象条件产生的烟云，为提高核武器的使用能力和了解核爆炸烟云的运动规律奠定重要基础。

4） 将WRF模拟结果和RAMS系统模拟的结果以及真实核爆炸数据的模拟结果进行对比，证明了本文模拟结果的可靠性。

在烟云上升模拟研究中还存在可以改进的地方。本文使用WRF模型模拟核爆炸烟云上升的过程时忽略了水汽混合比，爆炸火球的扰动温度取值也较小，但实际大气中是存在水汽和风的。因此，为了更贴近实际大气中的核爆炸烟云上升过程，在之后的实验中可以考虑水汽混合比和逐渐提高火球扰动温度进行模拟，并将其与实际核爆炸产生的烟云上升过程进行对比研究，使模拟研究的结果更加准确。

参考文献

[1] 乔登江. 核爆炸物理概论［M］. 北京： 原子能出版社， 1988.

[2] QIAO J D. Introduction to Nuclear Explosion Physics ［M］. Beijing： Atomic Energy Press， 1988（Ch）.

[3] 吴卓明，石洪祥，侯明远，等. 核化生防护大辞典［M］. 上海： 上海辞书出版社， 2000.

[4] WU Z M， SHI H X， HOU Y M， et al. Dictionary of Nuclear， Chemical and Biological Protection［M］. Shanghai： Shanghai Lexicographical Publishing House， 2000（Ch）.

[5] 郑毅. 核爆炸监测技术概论［M］. 北京： 国防工业出版社， 2019.

[6] ZHENG Y. Introduction to Nuclear Explosion Monitoring Technology［M］. Beijing： National Defense Industry Press， 2019 （Ch）.

[7] 潘自强. 核与辐射恐怖事件管理［M］. 北京： 科学出版社， 2005.

[8] PAN Z Q. Nuclear and Radiation Terrorism Incident Management［M］. Beijing： Science Press， 2005（Ch）.

[9] MORTON B R， TAYLOR G， TURNER J S. Turbulent gravitational convection from maintained and instantaneous sources［J］. Proceedings of the Royal Society of London Series A， 1956， 234（1196）： 1-23. DOI： 10.1098/rspa.1956.0011 .

[10] CAULFIELD C P， WOODS A W. Turbulent gravitational convection from a point source in a non-uniformly stratified environment［J］. Journal of Fluid Mechanics， 1998， 360（1）： 229-248. DOI： 10.1017/S0022112098008623 .

[11] ONUFRIEV A T. Theory of the motion of a vortex ring under gravity. Rise of the cloud from a nuclear explosion［J］. Journal of Applied Mechanics and Technical Physics， 1967， 8（2）： 1-7. DOI： 10.1007/BF00918021 .

[12] CHANG S H. Motion of a large dusty buoyant thermal with a vortex ring［J］. Journal Applied Mechanics， 1978，45（1）：711-716. DOI: 10.1115/1.3424407 .

[13] KANSA E J. A time-dependent buoyant puff model for explosion sources［R/OL］.［2024-01-01］. DOI: 10.2172/302869 .

[14] MAKHVILADZE G M， ROBERTS J P， YAKUSH S E. Modelling of atmospheric pollution by explosions［J］. Environmental Software，1995，10（1）：117-127. DOI: 10.1016/0266-9838(94)00005-r .

[15] KRISPIN J， COLLINS J P. Simulations of dusty flows in an incompressible gas using the projection method［C/OL］.［2024-01-01］. DOI: 10.2514/6.1995-2164 .

[16] BROWN B P， KRISPIN J， POTTS M A. Projection methods for incompressible multiphase cloud-rise phenomena［C/OL］.［2024-01-01］. DOI: 10.2514/6.1999-3321 .

[17] ARTHUR R S， LUNDQUIST K A， MIROCHA J D， et al. Simulating nuclear cloud rise within a realistic atmosphere using the Weather Research and Forecasting model［J］. Atmospheric Environment， 2021， 254（1）：118-363. DOI: 10.1016/j.atmosenv.2021.118363 .

[18] KANARSKA Y， LOMOV I， GLENN L， et al. Numerical simulation of cloud rise phenomena associated with nuclear bursts［J］. Annals of Nuclear Energy， 2009， 36（10）： 1475-1483. DOI: 10.1016/j.anucene.2009.08.009 .

[19] BELL J B， DAY M S， RENDLEMAN C A， et al. Adaptive low Mach number simulations of nuclear flame microphysics［J］. Journal of Computational Physics， 2004， 195（2）： 677-694. DOI： 10.1016/j.jcp.2003.10.035 .

[20] SONG S， LI Y B， LEE C， et al. Effect of surface conditions on blast wave propagation［J］. Journal of Mechanical Science and Technology， 2016， 30（9）： 3907-3915. DOI： 10.1007/s12206-016-0802-5 .

[21] WON S， LEE C. Simulation of the mushroom cloud generated from a high-energy explosion using large-eddy simulation［J］. Journal of Mechanical Science and Technology， 2020， 34（6）： 2443-2453. DOI： 10.1007/s12206-020-0520-x .

[22] 王坚， 李路翔. 核武器效应及防护［M］. 北京： 北京理工大学出版社， 1993.

[23] WANG J， LI L X. Nuclear Weapons Effect and Protection［M］. Beijing： Beijing Institute of Technology Press， 1993（Ch）.

[24] GLASSTONE S， DOLAN P J. The effects of nuclear weapons［R/OL］.［2024-01-01］. DOI: 10.21236/ada087568 .

[25] MOLENKAMP C R. OCTET User’s Manual［R］. Lawrence Livermore National Laboratory ：Technical Report UCRL-SM-226802， 2006.