烟火药在定容条件下的爆燃p-t 曲线仿真

2020-06-15赵象润李朝振闫利伟赵变玲王德君

含能材料 2020年6期

赵象润，严楠，李朝振，闫利伟，赵变玲，王德君

（1.北京理工大学爆炸科学与技术国家重点实验室，北京100081；2.辽宁北方华丰特种化工有限公司研发中心，辽宁抚顺113003）

1 引言

随着数值仿真技术的发展，采用非线性显式动力学软件对烟火药作用过程的可视化建模及分析提供了必要手段［1］。目前常用的非线性显式动力学仿真软件（如LS-DYNA，Autodyn，ABAQUS 等）中内置的与药剂有关的状态方程诸如JWL 状态方程、γ 律状态方程（Gamma Law Equation of State）、点火增长模型（Ignition and Growth Explosive Material Model）等都是描述炸药高速爆轰的作功过程［2-4］。这些模型并不适合描述烟火药爆燃的作功过程，因此现阶段对于烟火药爆燃过程的模拟大多采用基于经典内弹道方程组的理论体系，利用MATLAB 进行数值计算，获得燃烧室内的压力时间历程（p-t 曲线）［5-7］。在进行以烟火药为驱动能源的作动装置（如分离螺母）作用过程的动力学仿真时，将此p-t 曲线作为动力学仿真的输入条件，代替药剂的作用过程［8-10］。这种方法涉及两个不同软件的联合仿真，过程较为繁琐，而且这种方法不能获得压力在燃烧室内的空间分布。

采用非线性显式动力学软件MSC.Dytran 新开发的爆燃状态方程（Deflagration equation of state）可直接用来模拟烟火药作功过程。本研究中采用燃速较慢热稳定性良好的炭黑/硝酸钾（C/KNO3）和燃速较快的铝/高氯酸钾（Al/KClO4）两种烟火药，利用爆燃状态方程模拟定容条件下的p-t 曲线［11-13］。并对比了仿真结果与密闭爆发器试验结果，验证仿真方法的准确性和适用性［14］。以期为烟火药的内弹道仿真以及烟火药驱动的分离装置作功过程的动力学仿真提供一种新的方法。

2 实验部分

2.1 试剂与仪器

炭黑：工业级，GB3778-2011，天津亿博瑞化工有限公司；硝酸钾：分析纯，GB/T 647-2011，北京化工厂；铝粉：GB/T 2085.2-2007，沈阳东轻铝材经销有限公司；高氯酸钾：HG 3-1395-1981，分析纯，天津市风船化学试剂科技有限公司。

装药量：200 mg；装药方式：压装。密闭爆发器：（10±0.5）mL。

压力传感器：KISTLER601A 电荷型压电式压力传感器；量程：0～25 MPa；线性误差：±0.27%；灵敏度：150 pC/MPa；工作温度：（-196～200）℃；固有频率：150 kHz；抗冲击10000 g；过载：50 MPa。

信号调理器：丹麦戴维创DEWE-PACK-16。动态分析仪：丹麦戴维创DEWE-5001-TR，8 个通道，每个通道最高采样速率为10 MS/s/CH，单通道最高采样速率为25 MS/s/CH，12bit A/D。

压力药筒p-t 曲线测试装置如图1 所示。

图1 p-t 曲线测试装置Fig.1 Test device of p-t curves

2.2 仿真模型的建立

MSC.Dytran 内置的爆燃状态方程是用于模拟烟火药或推进剂爆燃过程的数学模型，可模拟烟火药燃烧时产生燃气的压力和反应速率，还能获得燃烧过程的能量特征（如温度、能量值等）［15］。该爆燃状态方程基于经典内弹道学，不考虑热散失，并假定密闭爆发器不变形、不漏气。在该模型中，烟火药作用产生的燃气的压力由诺贝尔-阿贝尔状态方程定义［16-17］：

式中，p 为压力，Pa；R 为与气体组分有关的气体常数，J·（kg·K）-1；T 为爆温，K；v 为比容，m3·kg-1；b 为余容，表示与单位质量气体分子体积有关的修正量，m3·kg-1；γ 为比热容比；E 为单位质量药剂内能，J·kg-1。

为描述火药燃烧过程，假设：①装药的物理化学性质均一，几何形状相同；②所有药粒表面同时着火，且燃烧过程按药粒表面平行层或同心层逐层燃烧的；③所有药粒的燃烧环境和条件相同，各方向的燃速相同；④不论火药燃烧期间还是燃烧结束后，生成物的成分始终保持不变［18］。设mp为装药初始质量，kg；mb为一段时间内火药的燃去量，kg；Vp为装药初始体积，m3；V 为一段时间内装药已燃去体积m3；n 为装药单体数目；ρp为装药密度，kg·m-3。则可以得到式（3）［19］：

式中，ψ 为装药已燃去质量与装药初始质量之比，为也称为火药燃去百分比。将ψ 对时间微分可得式（4）：

为了导出dV/dt，设某时刻单体药粒的表面积为Sb，m2，在dt 时间内燃去的厚度为de。则得到式（5）：

式中，e 为燃去厚度，m。将e 对时间微分，则得到火药的线性燃速表达式（6）：

式中，u 为火药燃烧的线性燃速，m·s-1；联立式（5）～（6）则得到式（7）：

由于燃烧涉及多种复杂化学反应，实际上难以根据火药的化学和物理性质预测火药的线性燃速［20］。在内弹道学上，火药的线性燃速常用圣罗伯特方程（Saint-Robert equation）表示［21］：

式中，w 为燃速系数，m·s-1·Pa-β，与火药性能、燃烧压强范围、装药初温等有关；β 为燃速压强指数，为量纲为1 的量，对于特定的火药，在一定压强范围内是常数［22］。

火药在燃烧时快速产生气体的能力由活度（vivacity）表示为式（9）［23］：

式中，ξ 为药剂活度，s-1；SA是药剂初始表面积（Sp）与Vp的比率。在燃烧过程中，随着碎裂程度的增加，火药燃速更快，活度更大。SA可以基于颗粒几何形状计算初始表面积与体积比，每种几何形状都具有特定的SA值［24］。

定义相对燃烧表面σ=Sb/Sp，其物理意义是燃烧期间当前表面积与初始表面积之比，它提供火药颗粒几何形状以影响燃烧部分的燃烧速率［25］，得到式（10）：

因为在分析过程中当前的颗粒表面积是未知的，为了简化输入，使用两个参数X 和s 输入到式（11）来近似替代典型的颗粒几何燃烧表面［26］：

式中，X、s 均为与形状相关的参数，软件帮助文档提供了各种形状的X，s 值供选用。

密闭爆发器燃烧室简化为圆柱形，共有15180 个节点和13664 个网格单元。按照密闭爆发器的形状设置欧拉域，再使用EOSDEF 模型定义炭黑/硝酸钾的物理和化学属性。在燃烧期间，燃气流入欧拉单元。仿真模型所需药剂输入参数通过公开文献或试验获取，部分能量参数应用Real 热力学程序进行求解得到。软件中每个参数都有对应的域名称（field heading），仿真输入的参数如表1 所示。

3 结果与讨论

3.1 炭黑/硝酸钾装药压力药筒升压过程

炭黑/硝酸钾的压力云图、p-t 曲线如图2～图3 所示。由于仿真计算不考虑桥丝升温点火过程，因此燃烧室从0 s 开始就产生压力。

图2 表明密闭爆发器内的压力传递过程和空间分布状态，整体上压力是近似逐层分布的。图3 可以看出仿真与试验的p-t 曲线之间的差异。二者的量化对比主要考察峰值压力（pm）和对应的作用时间t 的接近程度。在压力快速上升阶段（从0 到0.5pm），仿真结果和实测结果的升压速率（曲线的斜率）几乎重合。压力上升到0.5pm时，仿真与试验结果的作用时间仅相差0.03 ms，相对误差为3.4%；压力到达pm的时刻，时间差扩大到0.3 ms，相对误差为4.8%。实测pm比仿真结果高约0.5 MPa，压力偏高的原因是由于实测的压力包括了始发装药（镁点火药）产生的压力。通过单独的压力测试得到镁点火药自身产生的压力约为1.1 MPa，可见仿真的结果实际上偏高约0.6 MPa。这主要是由于仿真分析没有考虑热散失，且仿真分析中输入的能量参数是由Real 软件按照化学反应完全进行的理想条件下计算得出的。而实际作用过程中，由于制药工艺、装药条件等因素的影响，总会有能量损失。扣除发火药压力贡献的因素，仿真得到的峰值压力相对误差为10%。可见MSC.Dytran 对于火药燃烧过程的模拟准确度较好。

表1 炭黑/硝酸钾仿真模型输入参数Table 1 Input parameters of the simulation model of the carbon black/potassium nitrate

图2 炭黑/硝酸钾装药压力药筒在密闭爆发器内的压力-时间云图Fig.2 Pressure-time fringe/deformation diagram of the pressure cartridge with carbon black/potassium nitrate charge in closed bomb

仿真和测试p-t 曲线的明显不同之处在于实测曲线在压力达到峰值后开始下降，而仿真曲线并没有下降趋势，但压力波动减小。这是由于仿真运算是把整个反应过程按照绝热过程进行的，没有考虑密闭爆发器器壁的热散失，而测试结果的压力降一般认为是热散失导致的［28-30］。仿真曲线在压力达到峰值后震荡减小是因为软件在运算过程中，达到峰值后反应已经完全结束，没有燃气产生，因而密闭爆发器内的气体流动趋于稳定。

图3 炭黑/硝酸钾装药压力药筒p-t 曲线Fig.3 p-t curves of the pressure cartridge with carbon black/potassium nitrate charge

3.2 铝/高氯酸钾装药压力药筒升压过程

为验证仿真方法的准确性和对不同药剂的适用性，在相同条件下对铝/高氯酸钾进行了试验和仿真（试验条件与炭黑/硝酸钾相同）并进行分析，结果如图4 所示。

图4 铝/高氯酸钾装药压力药筒p-t 曲线Fig.4 p-t curves of the pressure cartridge with aluminum/potassium perchlorate charge

图4 表明，铝/高氯酸钾到达压力峰值的时间约为炭黑/硝酸钾的1/5，压力峰值也要比炭黑/硝酸钾高4.3 MPa。压力上升到1/2pm时，仿真结果的作用时间与实测结果相差0.02 ms，相对误差为22%；压力达到pm时，时间差扩大到0.1 ms，相对误差为7.7%。铝/高氯酸钾相对误差偏大的原因是其反应速率比炭黑/硝酸钾快很多，因此计算时的基数小导致。扣除发火药压力贡献的因素，仿真所得峰值压力比实测峰值压力高0.4 MPa，相对误差为4%。由铝/高氯酸钾的仿真也可以看出，MSC.Dytran 对于升压阶段的仿真结果与实测结果一致性很好。

上述分析表明，利用MSC.Dytran 的爆燃状态方程来模拟定容条件下火药的作用过程（仅指升压过程）有较好的准确性。事实上，利用该方法对于变容过程也是适用的，因为软件本身也可以解算不同体积时的药剂作用过程并进行迭代。从航天系统的实际应用情况来看，使用MSC.Dytran 软件进行分离螺母作用过程的动力学仿真本身是有较好的准确性的，整个过程也无需进行多软件的联合仿真［31］。