云爆燃料分散过程窜火机理的数值模拟

贾承志，张 奇

（北京理工大学爆炸科学与技术国家重点实验室，北京 100081）

1 引言

燃料空气炸药（Fuel Air Explosives，FAE）是一种新型爆炸能源，它利用云爆燃料分散与周围空气混合形成可燃云雾，再使可燃云雾发生爆轰达到毁伤目的，云雾范围是此类弹药威力的基础。燃料抛撒形成云雾的过程中，有可能提前发生爆燃现象，称为窜火。窜火的发生将大幅降低燃料空气炸药的爆炸性能，影响弹药威力，建立窜火预测模型，解决窜火问题对云爆装置的发展至关重要。Singh等［1］模拟了FAE装置在壳体破碎之后的近场阶段燃料抛撒模型，对燃料抛撒范围进行了研究。Sallam等［2］对燃料液滴在爆炸作用下的破碎机理进行了研究。肖邵清等［3-4］采用分散T型装药及不耦合装填中心分散药，以及在中心爆管内壁和分散药之间填充多孔惰性材料等新型装药方法，基本上使燃料分散过程中的窜火现象得到有效控制，有效遏制了窜火现象的发生。张奇等［5］认为抛撒过程中窜火现象的发生与中心抛撒药量有关，中心抛撒药量越大，越容易窜火。杜海文等［6］认为爆轰气体产物量是影响抛撒初始阶段窜火发生的主要因素，针对液固态燃料配方，认为最优的比药量范围为0.95%～1.70%。云爆燃料窜火机理复杂，国内外对窜火机理的研究公开报道很少，因此，本研究利用数值模拟软件Ansys-Fluent，针对以环氧丙烷为抛撒装药的云爆装置燃料抛撒过程进行了数值模拟，对燃料分散过程的窜火现象进行了分析，得到窜火规律，为小质量云爆装置燃料抛撒过程中窜火现象的预防提供了新的预测方法。

2 计算模型

2.1 二维模型的建立

对云爆装置进行燃料扩散及窜火研究时，为便于与已有实验数据进行对照，以2 kg云爆装置为研究对象，装置为高0.168 m，半径0.056 m的圆柱，云爆装置底部距离地面2 m，装置竖直放置，即云爆装置轴线与地面垂直。在中心装药作用下，研究液体燃料在开敞空间内抛撒的过程。由于云爆装置具有高度结构对称性，以云爆装置中心轴为旋转对称轴，建立轴对称旋转计算模型以减小计算量，并对二维计算网格进行划分，最小网格尺寸5 mm，最大网格尺寸200 mm，网格尺寸由中心装药位置至压力出口以1.04（104%）的增长率进行发散式非结构化网格划分，网格总数为14379个。云爆装置中心距离地面2 m，计算域为高6 m，长5 m的扇形区域。

云爆装置结构简化如图1所示。

图1 云爆装置简化结构图Fig.1 Diagram of simplified FAE structure

2.2 计算模型假设

研究中，假设中心装药起爆瞬间壳体已经破碎，燃料在中心装药作用下进行抛撒运动，液体燃料在抛撒过程中考虑液滴的破碎蒸发。

以中心装药药柱高度作为变量，为描述简便，以下使用图示高度h与云爆装置高度H的比值作为等效变量，在保持环境温度、环境压力、抛撒药装药量、中心装药药柱半径等条件不变的情况下，改变图1所示上下端燃料量高度h，即改变中心装药药柱高度，研究中心装药量对于窜火的影响。

云雾窜火的判据由窜火必须满足的两个必要条件［7］确定，即：燃料空气浓度处于爆炸上下限区间；温度不低于燃料云雾的自燃点。当云雾在抛撒过程中，满足上述两个条件时，即认为发生窜火。

2.3 控制方程

流体在流动过程中，需要遵循三个基本守恒定律，即：质量守恒定律、动量守恒定律、能量守恒定律。本研究通过Navier-Stokes方程（式（1）～（4）），求解环氧丙烷液雾抛撒过程中的气相流动［8］。

由质量守恒定律得到质量守恒方程：

式中，ρ为密度，kg·m－3；t为时间，s；Xi是空间坐标在 i方向的分量，m；Ui是速度矢量在 i方向的分量，m·s－1；Sm为源项，是来自液滴的质量通量。

由动量守恒定律得到动量守恒方程：

式中，p为静压，Pa；τij是应力张量；gi和 Fi分别是 i方向上的重力体积力和外部体积力，N。其中：

式中，μ为动力粘度，N·s·m－2；λ 为第二粘度，一般取-2/3。

式中，K是有效导热系数，无量纲；E=，h

eff为可压缩气体的显焓，kJ·mol－1；Jj是物质́的扩散通量，Sh为化学反应能量源项。

抛撒装药受力后很快破碎为液滴，在求解液滴的扩散时采用离散相模型，使用斯托克斯追踪（随机轨迹）轨迹模型，颗粒受作用力的平衡方程在笛卡尔坐标系下的形式为：

式中，u为气相速度，m·s－1；up为颗粒速度，m·s－1；ρp为颗粒密度，kg·m－3；FD（u-up）为颗粒的单位质量拖曳力，N；Fx为附加质量力，N。其中：

式中，dp为颗粒直径，m；Re为相对雷诺数，CD为拖曳力系数，且：

对于液滴的传热蒸发，当颗粒温度小于其蒸发温度（Tp＜Tvap）时，使用式（9）将液滴温度与液滴表面的热对流及辐射传热联系起来［9］。液滴与流体之间没有质量交换，当液滴穿过流体单元（计算网格）时，液滴吸收（释放）的热量作为源相作用到连续相的能量方程中：

式中，mp为液滴质量，kg；cp为液滴比热，J·（kg·K）－1；Ap为液滴表面积，m2；T∞为连续相的当地温度，K；Tp为液滴温度，K；h为对流传热系数，W·（m－2·K－1）；εp为颗粒黑度（辐射率），无量纲；σ为斯蒂芬孙-玻尔兹曼常数，为5.67×10－8W·（m－2·K－4）；θR为辐射温度，K；式（9）中假设液滴内部的热阻为零，即液滴内部温度处处一致。

当 Tvap＜Tp＜Tbp时，其中 Tbp为液滴的沸腾温度，液滴温度通过自身的热平衡得出，液滴的焓变与两相间的对流传热、汽化潜热通过热平衡的计算进行联系。

式中，（dmp）/dt为蒸发速率，kg·s－1；hfg为汽化潜热，J·kg－1。液滴与气相流体之间传质的计算中，液滴的蒸发量由梯度扩散确定，即从液滴向气相的扩散率与液滴与气流主流之间的蒸汽浓度梯度相关联：

式中，Ni为蒸汽的摩尔流率，kmol·（m－2·s－1）；ki为传质系数，m·s－1；C（i，s）为液滴表面的蒸汽浓度，kmol·m－3；C（i，∞）为气相主流的蒸汽浓度，kmol·m－3，对于第 i个组分，C（i，∞）由组分输运方程求解得到。液滴表面的蒸汽分压假定等于液滴温度Tp所对应的饱和压力psat，则此时有：

式中，R 为普适气体常数，8.314 J·（mol·K）－1。

气相流动涉及湍流现象，湍流流动是一种高度非线性的复杂流动。对湍流目前有多种计算方法，本文采用标准模型k-ε半经验公式模拟气相组分湍流运动过程。

式中，k表示湍流动能，m2·s－2；ε表示湍流耗散率，m2·s－3。

求解气态云雾浓度场与温度场的关系，在得到云团浓度场以及压力场数据的基础上，使用多方状态方程（式15）从压力场数据中取得对应的温度数据，确定温度场的分布，式中γ为多方指数。

考虑气体与液滴之间的辐射传热，采用P-1辐射模型，数值计算方法采用SIMPLE算法。

环氧丙烷爆炸极限采用式（16）、（17）、（18）确定［10］：

式中，L0下、L0上分别为环氧丙烷在常温常压下的爆炸上下限，由经验公式［11］确定。高温高压下的爆炸下限由式（16）确定，上限取（17）、（18）计算所得最大值。

2.4 参数设置

将2 kg环氧丙烷燃料进行抛撒，温度与大气压为常温常压，热容 1950 J·（kg·k）-1，密度 859 kg·m-3由文献［12］查得。

3 结果与讨论

计算中控制环境温度、环境压力、云爆装置总质量等其他参数不变，通过调整弹体高度H与中心抛撒装药距端板距离h的比值（i=H/h），改变中心装药量，得到不同比例下的燃料抛撒数据。

3.1 云雾抛撒半径

对2 kg质量的云爆装置燃料抛撒过程进行数值模拟计算，得到数值模拟和实验的云雾扩散半径与时间关系曲线。由图2可见，在保持其他条件不变，只变更中心装药量的情况下，云雾扩散半径并无明显变化，中心装药量对云雾扩散范围影响不大。数值模拟结果显示，扩散80 ms之后，云雾最大半径为3.397 m，最小为3.294 m，实验结果为 3.321 m［13］，误差小于 5%，在可接受范围内。

图2 云雾扩散半径与时间关系曲线Fig.2 Relationship between dispersion radius of propyleneoxide and corresponding time

由图2可以看出，数值模拟云爆燃料最终扩散半径约为3.3 m，燃料扩散速度在0～10 ms增长迅速，进一步计算表明，在10 ms左右达到速度最大值122 m·s-1，之后速度逐渐下降，但在35 ms时速度又有回升，这可能是中心大粒径液滴的影响，大粒径液滴受初始压力作用后加速较小液滴小，但小液滴在之后减速较快。

图3为相关文献中2 kg云爆装置燃料抛撒60 ms时刻的实验摄像图，图4为2 kg云爆装置燃料抛撒60 ms时刻的数值模拟云雾扩散浓度云图，从图像上可以看出数值模拟云雾扩散范围总体上与实验结果是吻合的。

图3 实验过程中2 kg云爆装置的扩散范围Fig.3 Dispersion range of a 2 kg FAE device in the experiment

图4 数值模拟2 kg云爆装置扩散浓度云图Fig.4 Simulated FAE fuel dispersion concentration contour

3.2 云雾场温度研究

由于温度只有在达到环氧丙烷自燃点（722.15 K）以上时，云雾才有可能发生窜火，为减少数据分析量，主要对燃料分散过程中云雾场的最高温度进行研究。由式（15）可知云雾场温度与压力正相关，为得到云雾场分散过程中的最高温度，对分散压力云图进行分析，以5 ms时刻的分散云团为例，图5所示为5 ms时刻燃料分散压力云图，以云爆装置初始时刻底部所在水平面建立横轴，以云爆装置对称轴为纵轴建立坐标，由式（15）计算得到此时刻云雾场的温度分布，可知此时温度最高点位于距离云爆装置底部约0.32 m，距中心轴约0.11 m处。

据此对数值模拟得到的不同时刻燃料分散云图进行分析，可得到不同中心装药量情况下燃料分散温度场最高温度与分散时间关系曲线，如图6所示。

从图6中可以看出，气相云雾温度场的温度在0～5 ms阶段先是从常温（298.15 K）骤增，在约2 ms时刻温度达到最大值，最高温度可达到3296.37 K，之后又急剧下降，在5 ms之后温度呈现缓慢下降趋势，最终趋于常温。在假定环氧丙烷自燃点不改变的情况下，可以看出10 ms之后云雾场温度一直小于自燃点温度（722.15 K［14］），即云雾在扩散 10 ms之后无窜火危险，这里重点研究云雾扩散前10 ms的特性。

图5 云爆装置燃料分散压力云图（t=5 ms）Fig.5 Pressure contour of FAE fuel dispersion（t=5 ms）

图6 云雾场最高温度与时间关系曲线Fig.6 Relationship between maximum temperature and time of FAE fuel dispersion field

3.3 窜火机理

一般云爆装置的窜火点往往与温度场温度最高点对应，针对不同中心装药量的云爆装置，对其燃料扩散过程中的温度场最高温度及附近气相浓度进行综合分析。

图7是i=4，中心药柱高为84 mm时燃料分散情况分析图，左侧Y轴为环氧丙烷燃料气相浓度、环氧丙烷爆炸上限以及爆炸下限，右侧Y轴为气态云雾场最高温度以及环氧丙烷自燃点。从图7可以看出，云雾场温度在约6.8 ms之前均高于环氧丙烷自燃点，但气相环氧丙烷浓度始终高于爆炸上限，故此种情况下不发生窜火。

图7 燃料分散情况分析曲线（i=4）Fig.7 Distribution of volume ratio and temperature of FAE fuel dispersion（i=4）

图8是i=5，中心药柱高度为100.8 mm时的燃料分散情况。从图8可以看出，云雾场内部最高温度在约5.8 ms之前均高于环氧丙烷自燃点（722.15 K），但云雾场内部气相环氧丙烷浓度在5.8 ms之前始终高于爆炸上限，此种情况下不存在同时满足窜火发生两个条件的时间段，故不发生窜火。

图 8 燃料分散情况分析曲线（i=5）Fig.8 Distribution of volume ratio and temperature of FAE fuel dispersion（i=5）

图9是i=6，中心药柱高度为112 mm时的燃料分散情况。从图9可以看出，云雾场内部最高温度在约5.6 ms之前均高于环氧丙烷自燃点，但云雾场内气相环氧丙烷浓度在此之前始终高于爆炸上限，窜火发生的两个条件无法同时被满足，不发生窜火。

图 9 燃料分散情况分析曲线（i=6）Fig.9 Distribution of volume ratio and temperature of FAE fuel dispersion（i=6）

图10是i=7，中心装药药柱高度为120 mm时的燃料分散情况。从图10中可以看出，虽然环氧丙烷气相浓度在4.8 ms～6.2 ms时处于爆炸上下限之间，气相占比在5 ms达到最小值27.77%，该时刻爆炸下限为1.6%，爆炸上限为33.5%，满足窜火发生的浓度条件，但云雾场内最高温度只有在约4.5 ms之前均高于环氧丙烷自燃点，在4.5 ms之后均低于环氧丙烷自燃点，此种情况下依然不存在同时满足窜火发生的温度与浓度条件的时间段，不发生窜火。

图10 燃料分散情况分析曲线（i=7）Fig.10 Distribution of volume ratio and temperature of FAE fuel dispersion（i=7）

当i=8，中心装药药柱高度为126 mm时的窜火情况分析曲线如图11所示。从图11可以看出，在约4 ms之前，温度场最高温度高于环氧丙烷自燃点，而气相浓度在约2.8～5.4 ms期间处于爆炸极限范围内，符合窜火发生的条件。在t=3 ms时刻，云雾场内部最高温度达到1244.7 K，环氧丙烷气相浓度为35.84%，爆炸上下限分别为1.58%和41.02%，此时云雾场内最高温度高于燃料环氧丙烷自燃点（722.15 K），气相燃料浓度在爆炸极限范围内，判定此种情况下发生窜火，t=3 ms时刻的云雾场压力云图如图12所示，图中红色线框所示为云爆装置初始时刻位置，由压力与温度的正相关关系可知，此时云雾场最高温度出现于约距中心轴0.01 m，距云爆装置底部0.42 m处，即为窜火点。

图11 窜火情况分析曲线（i=8）Fig.11 Distribution of volume ratio and temperature of FAE fuel dispersion（i=8）

图13为i=9，中心装药药柱高度为130.67 mm时的窜火情况分析曲线。由图13可以看出，在约4.3 ms之前，温度场最高温度位于环氧丙烷自燃点上方，气相浓度存在符合爆炸极限要求的阶段，在t=4 ms时刻，云雾场内最高温度达到1300.84 K，环氧丙烷气相浓度为39.24%，此刻的爆炸上下限分别为41.61%和1.17%，即云雾场温度高于燃料环氧丙烷自燃点（722.15 K），且气相燃料浓度在爆炸极限范围内，判定此种情况下也发生窜火。t=4 ms时刻的云雾场压力云图如图14所示，由图14可知，此时最高温度出现于约距中心轴0.03 m，距云爆装置底部0.26 m处，即为窜火点。

图 12 云爆装置燃料分散压力云图（i=8，t=3 ms）Fig.12 Pressure distribution contour of FAE fuel dispersion（i=8，t=3 ms）

图 13 窜火情况分析曲线（i=9）Fig.13 Distribution of volume ratio and temperature of premature-combustion（i=9）

图 14 云爆装置燃料分散压力云图（i=9，t=4 ms）Fig.14 Pressure distribution contour of FAE fuel dispersion（i=9，t=4 ms）

图15为i=10，中心装药药柱高度为134.4 mm时的窜火情况。，由图15可以看出，在约5 ms之前，温度场最高温度大于环氧丙烷自燃点，气相浓度在4.1 ms到4.5 ms阶段符合爆炸极限，在t=4.2 ms时刻，云雾场温度达到745.37 K，环氧丙烷气相浓度为29.97%，在此时刻的爆炸上下限分别为32.85%和1.54%，即云雾场温度高于燃料环氧丙烷自燃点（722.15 K），且气相燃料浓度在爆炸极限范围内，判定此种情况下同样发生窜火。t=4.2 ms时刻的云雾场压力云图如图16所示，由图16可知，窜火点出现于约距中心轴0.02 m，距云爆装置底部0.18 m处。

将上述数值模拟计算结果汇总，得到2 kg云爆装置窜火情况在表2中给出。

图15 窜火情况分析曲线（i=10）Fig.15 Distribution of volume ratio and temperature of premature-combustion（i=10）

图 16 云爆装置燃料分散压力云图（i=10，t=4.2 ms）Fig.16 Pressure distribution contour of FAE fuel dispersion（i=10，t=4.2 ms）

表2 数值模拟云爆装置窜火情况Table 2 Summary of numerical simulation about prematurecombustion

4 结论

（1）采用数值仿真模拟云爆装置燃料扩散过程，通过燃料分散压力云图导出温度场分布，再综合气相浓度确定窜火条件，分析了温度场与浓度场随时间的变化规律，通过监测不同时间点云雾场最高温度处燃料浓度，可对云爆装置燃料分散过程中的窜火现象进行预测。

（2）针对2 kg的云爆装置，中心装药药柱高度为126，130.67，134.4 mm时均发生窜火，发生窜火的位置位于云爆装置顶部附近，窜火时间集中在燃料分散开始后的3～4.2 ms。