炸药慢烤热点火起爆数值模拟及试验研究

王成，方钧宇，杨同会

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

大多数凝聚相炸药在持续热源的作用下会发生热分解，可能导致燃烧、爆燃甚至爆轰的发生，造成严重的损失.一般用烤燃试验[1-2]来测试并研究炸药在这种情况下的宏观响应规律，以减少意外事件的发生.但烤燃试验的缺点在于成本较高、周期较长，部分物理量测量难度较大，因此结合数值模拟可以对烤燃现象进行更全面的分析.

现阶段对烤燃过程的计算大多依赖商业软件.MAIENSCHEIN 等[3]使用ALE3D 对B 炸药和PBXN-109 炸药的烤燃进行了计算.常双君等[4]采用了Fluent 进行了LLM-105 基战斗部装药烤燃试验的计算.孙培培等[5]开展了烤燃试验，结果认为壳体厚度对烤燃的点火临界温度无显著影响；陈科全等[6]基于ABAQUS 软件针对不同壳体厚度的烤燃弹进行了数值模拟研究，得出了和孙培培等类似的结论.常见的凝聚相炸药多为低熔点炸药，比如TNT 炸药熔点为352～355 K，远低于其热分解温度，因此在烤燃过程中会发生从固相到液相的相变[7].陈朗等[8]、寇永锋等[9]参考了MCGUIRE 和TARVER[10]在1981 提出的用吸热描述炸药相变的简化模型，并结合观察烤燃试验中的升温曲线，发现相变过程不可忽略，并使用流体力学计算软件Fluent 进行了计算.商业软件为烤燃过程的计算提供了诸多便利，但大多商业软件无法在算法层面优化模拟过程，且这些软件在计算炸药点火后起爆过程时稍有欠缺，比如LS-DYNA 在计算时需要将点火区域换算成等当量点火药，与实际物理过程有所偏差.

为了能完整且精确地模拟出炸药慢烤热点起爆过程，形成完善且符合实际物理过程的慢烤热点起爆计算软件，本文针对柱形TNT 烤燃弹算例，基于其轴对称性质，发展了柱坐标下关于热传导方程的交替方向隐格式，对慢烤升温过程进行模拟.在炸药热点火后，采用化学反应流体动力学Euler 方程组描述起爆过程，使用5 阶WENO 格式，编写了炸药起爆与爆轰仿真程序.并且针对算例工况进行了慢烤试验，对计算得出的升温曲线、点火温度、点火时间进行了验证.

1 慢烤热点火起爆控制方程与数值格式

在慢烤过程中，凝聚相炸药主要经历两大阶段：升温阶段以及热点火后起爆、爆轰阶段.划分两个阶段的关键物理量为炸药反应速率.在慢烤升温阶段，炸药的热分解十分微弱[11]，反应速率几乎为0，炸药内部为低压零速的简单温度场，采用热传导方程即可描述这个过程；当温度到达一定阈值时，炸药发生热点火，反应速率就会急剧增长[12], 强烈的热分解反应导致炸药内部出现高温高压，粒子高速移动，热传导方程无法描述热点火后炸药内部演化过程，因此采用化学反应流体动力学Euler 方程组描述热点火后起爆、爆轰过程.

1.1 慢烤升温过程

1.1.1 慢烤升温过程计算模型

针对柱形装药，凝聚相炸药慢烤升温过程遵循的热传导方程为

dλ/dt实际上为Arrhenius 公式，表示化学反应速率，m为反应级数，本文采用一级反应模型进行计算：

1.1.2 慢烤升温过程计算方法

为实现慢烤升温过程与热点火后响应过程计算的耦合，保证数据交换的准确性，在两过程的计算中空间网格密度应当保持一致.慢烤升温过程持续时间较长，通常为几十分钟至数小时，若采用显式格式进行计算，计算过程会受限于CFL 条件，在高密度空间网格下，会产生极多的无用计算，浪费计算时间与储存空间.部分无条件稳定的隐式格式，在求解过程中需要求解大型方程组，计算效率低.

因此综合考虑计算精度、算法稳定性和求解效率，发展了如下针对慢烤升温过程的柱坐标下的交替方向隐格式：

以步长h划分空间，取时间步长 τ，定义导温系数D=λT/ρCv，将i、j、n分别作为划分z、r、t的标记，其中 0 ≤i≤imax，0 ≤j≤jmax，0 ≤n≤nmax，定义记号

该求解格式关于时间空间都为2 阶精度，无条件稳定，可任意设置时间步长，且在求解过程中产生的方程组全为三对角方程组，能有效地减少计算量，提高计算效率.

TNT 炸药等低熔点炸药在慢烤过程中会发生相变，在计算这类炸药的慢烤升温过程时，导温系数D会随着温度而改变，为此，定义液相率 β[8]，并根据焓值与温度、相变潜热的关系进行相变处理.

式中：TS为 炸药固-液相变起始温度；TL为炸药固-液相 变 终 止 温 度；λTL为炸药液态时的热导率；λTS为 炸药固态时的热导率；Qβ为炸药相变潜热.

1.2 热点火后起爆和爆轰过程

1.2.1 热点火后起爆和爆轰过程计算模型

针对柱形装药，凝聚相炸药热点火后起爆、爆轰过程遵循柱坐标系下Euler 方程组，描述为

式中：vr、vz分 别为径向和轴向的速度；p为压力；E为单位质量的总能；e为比内能.

对于凝聚相炸药与爆轰产物，通常采用JWL 状态方程来描述压强、密度、温度等物理量之间的关系：

式中：A、B、R1、R2、ω为 常数；V=ν/ν0表示当前比容与初始比容的比值.

在炸药热点火后的一段时间内，由于炸药内部压力与粒子速度较小，反应仍由温度主导，因此仍可采用Arrhenius 公式描述反应速率.当炸药内部压力成长到一定值时，压力将替代温度成为主导化学反应的因素，此时采用LEE 和TARVER 提出的点火增长模型描述反应速率：

式中I、a、b、c、d、e、g、x、y、z、G1、G2、λ1、λ2、λ3均为常数.3 个开关函数分别控制三项式点火增长模型中的热点成核项、慢速反应项、快速反应项.TNT 炸药等低熔点炸药在慢烤升温过程中会发生固-液相变，导致炸药内部无孔穴等结构缺陷，热点成核项不起作用，因此在计算这类问题时直接关闭热点成核项，取 Φ1=0.

因此采用如下修正后的反应模型计算凝聚相炸药烤燃热起爆和爆轰过程：

在该模型中，当温度主导化学反应时，采用Arrhenius 公式计算反应速率；当压力主导化学反应时，通过采用点火增长方程计算反应速率.

1.2.2 热点火后起爆和爆轰过程计算方法

在对热点火后起爆与爆轰过程的计算中，空间离散采用5 阶WENO 格式[13-14]，通过引入自适应权重，在强间断附近自动选择最光滑的模版进行通量重构，有效抑制了高精度计算中的非物理振荡，时间离散采用总变差递减的3 阶Runge-Kutta 法，进一步增强数值格式的稳定性.

热传导方程主要对温度进行了计算，而Euler方程组涉及到压力、速度、能量等物理量，因此在计算热点火后起爆与爆轰过程前，应当对升温过程的计算结果数据进行预处理.为此，设计如下数据接口：

式中：下标Euler 与heat 分别表示该变量在Euler 方程组与热传导方程中的取值；ε为小量，一般取0.01；αg为爆轰产物与炸药比热容的比值.

表1[8-9,15]给出了TNT 炸药计算过程中涉及到的的所有相关参数.

表1 TNT 炸药相关参数Tab.1 Related parameters of TNT

2 TNT 炸药烤燃试验

为验证算法可行性，设计并开展TNT 炸药烤燃试验.试验采用K 型WRNK-191 温度传感器进行测温，HIOKI LR8431-30 数据记录仪用于温度显示与采集.在壳体表面均匀缠绕电热线圈，外层包裹石棉保温层，图1 为TNT 烤燃弹加工前后图.

图1 TNT 烤燃弹弹体Fig.1 TNT cook-off projectile

图2 为TNT 烤燃弹弹体柱截面示意图，药柱直径为60 mm，高为240 mm，内部炸药为TNT 炸药，装药密度为1 624 kg/m3，外部壳体厚度为3 mm，材料为45 号钢.

图2 烤燃弹算例截面尺寸图Fig.2 Example section size diagram

在测温点A、B、C处放置温度传感器，其中A处温度传感器测得的温度数据用于分析温度变化规律，B、C处温度传感器测得的温度数据与温控系统进行实时交互，以控制升温速率.

烤燃试验起始环境温度为284 K，升温速率为1 K/min.位于测温点A处的温度传感器记录的温度数据于236 min 终止，记录的最高温度为509 K.此时从监控画面观测，TNT 烤燃弹已开始冒烟，温控系统显示壳体温度为520 K，随后TNT 烤燃弹发生剧烈反应，如图3 所示.

3 计算结果与分析

3.1 慢烤升温过程计算结果与对比分析

针对TNT 炸药慢烤升温过程进行数值模拟计算，并与试验测得数据进行对比，图4 为慢烤试验数据与计算结果对比.从图4(a)所示的计算结果曲线可以看出，测温点A处炸药升温曲线在112 ～122 min出现相变平台，升温速率减小至0.46 K/min，温度由352 K 升至355 K；随后升温速率跃升至2.7 K/min；直至146 min 时，升温速率逐渐减小至1.0 K/min.测温点A处试验结果曲线中，108 min 时，柱心传感器处的升温速率减小至0.5 K/min 左右；试验进行至125 min时，柱心传感器处的升温速率升高至6.0 K/min 左右；试验进行至133 min 时，柱心传感器处的升温速率恢复至1.0 K/min 左右直至纪录终止.可以看出测温点A处的计算结果在升温趋势、相变平台上与试验高度一致，在相变段的最大相对误差为5%左右，在非相变段的最大相对误差为3.5%左右.

图4(c)所示的计算结果中，壳体的升温速率保持在1 K/min 直至233 min，此时壳体温度为525 K；随后，壳体升温速率不断增大，在236 min 左右，壳体温度直线上升，说明炸药在236 min 左右发生热点火并起爆.在试验中，壳体升温速率在235 min 时由520 K 开始快速升高，与计算结果吻合度较高.

图4 TNT 烤燃试验数据与计算结果对比图Fig.4 Comparison of experimental data and calculated results

图5 为236 min 时TNT 烤燃弹温度分布图，此时由于TNT 炸药热分解产生的热量累计，弹体内部2个圆环区域的温度高于壳体，TNT 炸药发生热点火，点火时弹体内最高温度大约为550 K.点火区域是z=10 mm 与z=230 mm 位 置 处 半 径 为20 mm 的2 个点火圆环，圆环内径大约为3 mm.

图5 236 min 时TNT 烤燃弹内部温度云图Fig.5 Cloud diagram of temperature in TNT projectile at 236 min

3.2 热点火后起爆和爆轰过程计算结果与分析

以图5 对应的慢烤模拟结果为初始温度场，通过1.2.2 中的数据接口设置初始条件，同时设置刚性壁面条件，开展热点火后起爆、爆轰过程数值模拟.图6、图7 分别为TNT 炸药热点火后典型时刻内部温度、压力云图，高温区域逐渐扩大，且温度由600 K持续上升至1 800 K.高压区域伴随高温区域产生，在2.83 μs 时，高温区域产生了2.6 GPa 的高压，随后产生压力波向柱心传播，在3.80 μs 时，压力波达到6.5 GPa，并接近柱心.图8 为在4.28 μs 时弹体内部压力云图，此时冲击波在柱心汇集，产生26 GPa 的高压.

图6 热点火后典型时刻烤燃弹体内温度云图Fig.6 Cloud diagram of temperature in projectile at typical time

图7 热点火后典型时刻烤燃弹体内压力云图Fig.7 Cloud diagram of pressure in projectile at typical time

图8 4.28 μs 时烤燃弹体内压力云图Fig.8 Cloud diagram of pressure in projectile at 4.28 μs

图9 为TNT 烤燃弹热点火后压力、反应率随时间变化曲线.曲线对应的点处于柱心对称轴上z=10、30、50、70、90、110 mm 位置处.从计算结果可以看出，柱心对称轴上z=10 mm 位置处的TNT 炸药在3.85 μs 时达到压力峰值，峰值大约为15 GPa，此时该位置反应率接近100%，能量完全释放.结合点火圆环位置，推测出由点火圆环产生的压力波在柱心完成汇聚的时间点在热点火后3.85 μs 左右.柱心对称轴上z=30 mm位置处的TNT 炸药在4.94 μs 时达到压力峰值，峰值大约为33.7 GPa；柱心对称轴上z=50、70、90、110 mm 位置处的TNT 炸药几乎同时在5.74 μs 时达到压力峰值，且峰值几乎一致，大约为33.7 GPa.这是因为高压波阵面后稀疏波与反应流达到平衡，峰值保持稳定，说明TNT 炸药已形成稳定爆轰，起爆时间大约为4.94 μs.

图9 热点火后典型位置压力、反应率变化曲线Fig.9 Pressure and reaction rate curves at typical locations

4 结 论

结合试验与模拟的结果，得出以下主要结论：

①本文推导的交替方向隐格式在计算TNT 炸药慢烤时，具有较高精度且实现了相变，柱心传感器处的试验与模拟结果相对比，在非相变段的最大相对误差为3.5%左右，在相变段的最大相对误差为5.0%左右.计算得出点火区域为z=10 mm 与z=230 mm位置处半径为20 mm 的2 个点火圆环，圆环内径大约为3 mm.

②TNT 烤燃弹热点火起爆前壳体临界温度为520 K，柱心传感器处温度为509 K.

③热点火后，点火区域产生压力波向柱心传播，于3.85 μs 左右在柱心汇聚，产生15 GPa 的高压；TNT烤燃弹起爆时间大约为4.94 μs，形成稳定爆轰后，中心轴上峰值压力大约为33.7 GPa.