基于欧拉法的核爆烟尘颗粒群运动模拟

郭思禹，郑 伟，张 伟，郭 俊

（1.西安交通大学，西安，710025；2.北京宇航系统工程研究所，北京，100076）

0 引言

核爆烟尘是核爆炸火球熄灭后形成的一种蘑菇状、具有放射性的烟云和尘柱，持续时间长，覆盖范围广，烟尘颗粒物形成的热力学环境会对广域空间中的人员以及设备造成损伤［1］。

烟尘颗粒物的产生与分布涉及等离子体物理、辐射流体力学、凝聚态物理等多个学科。目前学术界尚难精确模拟核爆炸爆源初场［2］，主要研究核爆后烟尘颗粒物运动规律，以及沉降到地面后产生的辐射效应［3］，极少关注烟尘颗粒群在大气中的传输扩散。

刘朝辉等［4］模拟了烟尘中单个颗粒的大气运动，建立了一个拉格朗日模型，将颗粒分为大颗粒与小颗粒，分别给出其在烟尘中的轨迹。但是，无论是大、小颗粒的划分方法，还是颗粒初始条件的设置，都缺乏理论依据，文献中也没有涉及对颗粒物质量的模拟研究。卓俊等［5］建立了一个颗粒群拉格朗日模型，模拟颗粒速度、质量，但该模型以稳定烟尘为起始时刻，不能对烟尘上升阶段进行表征。郑毅等［6］采用中尺度气象软件RAMS模拟了上升阶段的烟尘颗粒流动，但仿真中将固定粒径颗粒的速度设为定值，且在模型中只考虑了大气中的放射性活度浓度，没有研究质量浓度。

本文在二维圆柱坐标系中基于欧拉网格方法建立了烟尘颗粒群运动模型，考虑了历史试验数据与真实大气条件，对颗粒速度、质量进行时空表征，研究了核爆烟尘颗粒群的传输扩散规律。

1 理论建模

地面爆和近地面爆的物理模型如图1所示，爆后蘑菇状烟尘的发展由大气中各种固体颗粒群的运动来体现，而地面沙粒、土壤、岩石等物质的成分是构成这些颗粒的主体。这些原本在地面的颗粒由于受到爆炸后气体的抽吸带动进入大气中，在烟尘稳定后颗粒群又开始沉降，沉降过程中伴随着扩散，其中，小粒径颗粒相比大粒径颗粒的迁移距离更远。

图1 地面和近地核爆炸颗粒运动物理模型Fig.1 Physical model of particle motion in ground and nearground nuclear explosions

同气体相一样，在气固两相流中，固体相同样满足质量守恒定律，固体相颗粒质量浓度c（单位为kg·m-3）满足［7］：

式中u，v，w分别为x，y，z方向上的颗粒平均移动速度。

如果核烟尘中只有均匀的气相流速，即假设层流条件下，颗粒群在气流中的运动如图2a 所示，事实上，大气中存在着剧烈的湍流运动，使颗粒群与空气之间强烈地混合和交换，使得颗粒不仅沿气流流动方向传输，而且还会向各个方向扩散，如图2b所示。

图2 层流和湍流情况下颗粒群的扩散Fig.2 Diffusion of particle population in laminar and turbulent flow conditions

考虑由湍流引起的速度脉动和浓度涨落，可将速度和浓度写成平均值与脉动值之和：

将式（2）代入式（1），并按流体力学中的雷诺平均法则取平均，经整理可得：

运用梯度输送理论［7］，任意物理量的脉动输送值与该特征量的平均值的梯度成比例关系。若该物理量为扩散物理的质量浓度c，则有：

式中Kx，Ky，Kz分别表示x，y，z方向的湍流扩散系数。

将式（4）代入式（3），有：

不考虑风速，认为核烟尘颗粒群在三维空间呈轴对称分布，即可在二维圆柱坐标系中建立颗粒群运动模型。

将颗粒都近似成球形，颗粒群按粒径划分成多个区间，以某一中值粒径下的颗粒浓度等效代替这一区间颗粒群的浓度，有：

式中Ci=Ci(z，r，t)为第i种粒径颗粒时刻t时在坐标点(r，z)的质量浓度；w(z，r，t)为颗粒群层流速度；Kz(z，r，t)，Kr(z，r，t)分别为颗粒湍流扩散系数的垂直分量和水平分量。

2 数值计算方法

在明确烟尘颗粒群运动物理模型的基础上，基于欧拉视角的数学求解方法按照如图3所示的流程进行。

图3 颗粒群运动模型算法流程Fig.3 Algorithm flow of particle population motion model

在二维圆柱坐标系中，将左边界取为核烟尘的对称轴，如图4所示，则左边界满足第二类边界条件：

图4 颗粒群运动模型示意Fig.4 Schematic of particle population movement model

下边界取为地面，地面与上方存在质量交换。上边界与右边界应取适当大的值，确保边界处浓度为0即可。

被带到大气中的土壤总质量为［3］

式中kH=0.077 41；Y（kt，TNT）为核爆当量；Hˉ为爆炸比高，0 ≤Hˉ≤7.19 m/t1/3.4。

对核烟尘放射性沾染尺寸分布的测量［8］表明，烟尘粒子的数量按粒径近似服从对数正态分布：

将粒子看成等密度球体，粒子表面积与粒径成二次方，体积与粒径成三次方，则粒子表面积与体积也均近似服从对数正态分布：

设定1 000 kt 内华达地面爆场景，烟尘颗粒平均密 度ρS=2 600 kg/m3，μN=ln0.407 μm，μS=2.94 μm，μV=5.61 μm，σ=ln 4 μm。

将烟尘颗粒群按粒径分为16 个区间，各区间的初始颗粒数如表1所示。

表1 1 000 kt内华达地面爆颗粒群分谱Tab.1 1 000 kt Nevada ground burst particle population spectra

假设初始烟尘是一个以爆心为球心、与地面相接的球状物，初始半径为

假定烟尘颗粒群在初始烟尘中均匀分布，初始时间为

本研究采用均匀正方形进行空间的网格划分，如图4 所示。应用时域有限差分法对式（6）进行离散求解，时变项采用显示格式，对流项采用一阶迎风格式，扩散项采用二阶中心差分格式［9］，如下：

简化得到离散的时间步进公式：

离散的左边界条件为

离散的初始质量为

在轴向上，颗粒会受到重力、浮力和流体的阻力等，当颗粒速度小于周围流体速度时，阻力表现为曳力。

忽略浮力及其他力，将w取为重力与阻力达到平衡时的值［10］：

式中wg为核烟尘中气体的发展速度与密度，根据典型核试验数据取值［8］；CD=0.44为阻力系数；ρg为大气背景密度［11］。

针对wg和w，取向上为正方向，经推导得到：

水平方向湍流扩散系数取为［12］

垂直方向湍流扩散系数取为

设置烟尘稳定前时间步长0.05 s，稳定后时间步长2 s，空间步长Δz=Δr=200 m，时间尺度1 h，空间尺度为z方向20.5 km、r方向12 km。

在物理上判断解是否收敛的方法是在每次时间步进计算后，即每对式（14）进行一次求解后，判断总颗粒质量是否不变，即质量是否守恒。

3 结果分析

基于上述数学物理方法，最终计算得到16 个粒径颗粒群的时空分布，粒径中值13.3 μm 颗粒群质量浓度在4个时刻的空间分布如图5所示。

图5 粒径中值13.3μm颗粒群数量时空分布Fig.5 Spatial distribution of the number of particle population with a median particle size of 13.3μm

在图5 中，标注了质量浓度大于1×10-3g/m3的颗粒群的轴向迁移速度，可以看出，由于上方大气空气稀薄，距地面更高处的速度大于低处的速度。

颗粒在升降过程中还会不断扩散，在空间近似满足高斯分布，由于在柱坐标系下建模，明显可见颗粒群向r轴正方向迁移。

在爆后60 s，中值粒径13.3 μm 颗粒群的质量浓度中心大概在地面上方约5 km处，爆后5 min升高到约13 km，之后缓慢降低，爆后1 h浓度中心为地面上方约8 km。此外，由于颗粒群的扩散，浓度中心处的质量浓度一直随时间减少，爆后60 s 为5～6 g/m3，爆后1 h为0.15 g/m3。

将模拟出的爆后1～5 min的5个时刻下2.37 μm颗粒群在空间的近似高斯分布的3σ值作为烟尘高度与宽度，并将此结果与已有文献中核试验数据［8］作对比，结果如表2 所示。由表2 可知，本模型得到的颗粒群运动仿真结果与已有核试验数据一致性良好。

表2 烟尘颗粒群模拟结果对比Tab.2 Comparison of simulation results of fallout particle population

4 结束语

本文基于欧拉方法构建了核爆烟尘颗粒运动模型，得到了与真实核爆炸运动过程相似的现象，模拟结果与已有核试验数据基本一致。因此，采用该方法模拟烟尘颗粒群空间传输是可行的，且可节省大量时间和费用。但由于没有采用计算流体力学方法，涡流现象难以模拟，也没有考虑风速的影响，这两个方面的研究有待于今后继续开展。