炸药爆轰热力学数据非线性拟合方法研究❋

杜明燃 张金龙 梁琳娜 黄亮亮 韩体飞 刘 锋

安徽理工大学化学工程学院(安徽淮南，232001)

引言

研究炸药性能时，爆轰参数一直都是备受关注的课题，它是预测炸药性能和理论设计炸药的基础[1]。为计算炸药爆轰参数，国内外学者提出了BKW、JWL、LJD、JCZ3和 WCA 等物理模型[2-7]，并发展了相关计算程序，对理论研究爆轰参数作出了重大贡献。研究发现，利用BKW[3]、WCA和LJD等状态方程计算爆轰参数时，爆轰参数的求解过程是一个反复迭代的过程，求解爆温的热力学参数若满足一种函数关系，对设计程序自循环求解爆轰参数具有重要意义。

以往研究者在研究爆温时多采用Kast平均热容法[8-9]，这种计算方法基于平均热容建立，在精确计算中存在不足，导致爆轰参数计算得不准确。从理论角度分析，能量守恒法[9]计算爆温具有较高的精度，然而，由于理论数据的缺失，能量守恒最终确定的爆轰温度一般根据理论数据估算得出；同时，在编程计算中发现，爆轰产物一般十余种，若不对热力学参数进行处理，能量守恒法需要反复迭代大量数据，程序繁琐[10]且容易出错。因此，利用拟合的方法确定一种函数关系变得非常必要，在此基础上再建立求解爆轰参数的数学模型，对求解爆轰参数及爆炸物品的应用具有一定参考意义。

本文中，采用非线性拟合的方法处理爆轰产物热力学数据，比较非线性拟合法与Kast平均热容法计算物质的焓变(HT-H298)与理论值的差别。在保证非线性拟合方程精度的基础上，依据此非线性拟合方程建立求解爆轰参数的数学模型，模型选取BKW状态方程作为爆轰产物的状态方程，利用最小自由能原理和Hugoniot关系，最终编程实现此过程。利用此程序，求解RDX和硝基甲烷(CH3NO2)的爆温、爆压、爆速和爆轰产物组成，验证此方法的可行性。

1 炸药爆轰温度的确定方法

1.1 Kast平均热容法

计算炸药爆轰温度常用的方法是爆轰产物平均热容法[8]。爆轰产物平均热容法是利用在一定温度范围内物质的平均热容计算爆轰温度，物质的平均热容由经验公式确定。爆轰产物平均热容法的计算公式如下：

式中：QV为炸药定容爆热，kJ/mol；为爆轰产物的平均热容；T为爆温，K；T0=298.15 K；t为温差；A、B为平均热容的线性系数。

1.2 能量守恒法

能量守恒法[9]利用物质能量守恒确定爆温，能量守恒法满足如下计算公式：

1.3 非线性拟合法

能量守恒法具有严格的理论依据，但能量守恒法与Kast平均热容法相比也具有一个重要缺陷：理论数据不足，导致部分计算只能停留在估算阶段，给爆温的确定带来误差，导致计算其他爆轰参数也出现误差。因此，依据能量守恒法建立模拟物质的热力学关系的方程变得十分必要。依据式(4)、式(5)和式(6)，可以得出以下假设公式：

式中：A1、A2、A3、A4为待拟合系数，可以根据最小二乘法编程计算得到，编程依据热化学参数可由文献[11]获得。

式(7)、式(8)确定后，可以依据式(3)建立确定爆温的能量守恒公式。表1给出了爆轰产物对应系数的拟合结果，爆轰产物中的碳单质视为金刚石[12-18]。表1适应范围为298～5 000 K。为比较线性拟合法与平均热容法计算方法的准确度，分别利用Kast平均热容法和非线性拟合法计算CO2、H2O和N2的焓变(HT-H298)，共7组数据。表2列出两种方法计算所得的3种主要产物的焓变(HT-H298)与理论值[11]。

表2显示，对于CO2、H2O和N2，Kast平均热容法的计算误差远远大于本工作所采用的非线性拟合法。表2中，非线性拟合法的最大误差为0.5%；Kast平均热容法的计算最小误差为3.1%，最大误差为12.5%。

2 求解爆轰参数过程

2.1 BKW状态方程

BKW状态方程[10]最初是由Becker提出的，后来经过Kistiakowsky和Wilson等的修正，最后确定为以下形式：

表1 式(7)和式(8)系数的拟合值Tab.1 Fitting values of coefficient for Eq-7 and Eq-8

表2 焓变(HT-H298)的计算与理论比较Tab.2 Calculation results and theoretical values of(HT-H298) kJ/mol

式中：p为爆轰压力，Pa；β为常数；X=κ∑xiki/V(T＋θ)α；V为爆轰气体摩尔体积，κ、α、θ为经验常数，ki为第i种物质的余容，xi为这种物质的摩尔分数。

2.2 求解爆轰参数的数学模型及模拟算法

利用最小自由能原理建立非线性方程组。根据质量守恒和物质自由能公式，再运用拉格朗日乘数法把条件极值转化为非条件极值，可以得出与爆轰产物系统具有相同极值条件的公式：

式中：G、g、g0分别为爆轰体系自由能、物质自由能和标准状况下自由能，kJ/mol；Nk为1 kg爆炸产物含有第k种元素的物质的量；λ为拉格朗日因子；k为第k种元素序号；λk为第k种元素的拉格朗日因子；i为第i种气体产物序号；j为第j种固体产物序号；aik为每摩尔第i种气相组分含有k元素的物质的量，mol；ajk为每摩尔第j种凝聚相组分含有k元素的物质的量，mol；ng为气体产物总物质的量；nig为第i种气体的物质的量；njs为第j种固体的物质的量；l为炸药所含元素的种类数；m和n分别为气体和固体的种类。

式(10)对应极值条件的方程组为非线性方程组，为求解此方程组，根据本小组研究的牛顿迭代法求解非线性方程组的思路[19]，利用以下模型求解爆轰参数和爆轰产物组成。

首先，需要把 əG(n，λk)/ənig、əG(n，λk)/ənjs和əG(n，λk)/əλk看作是关于nig、njs和λk的方程组，令

把式(11)、式(12)、式(13)和式(14)代入式(16)，可求得F′(x)，进一步求得F′(x)-1，并记解为x(k＋1)，对应牛顿迭代法公式：x(k＋1)=xk-F′(xk)-1F(xk)，求解上述公式要先假设一组爆温T′和爆压p′，此温度和压力可以根据一般经验计算值确定，再取经验计算所得爆炸产物组分量并作一定变换后，作为初始近似根x1，设置计算精度ε=0.000 001，即xk和x(k＋1)对应元素的误差绝对值。

炸药爆轰气相产物分子体系必须满足Hugoniot关系[10，20]：

爆速表达式为

根据BKW状态方程、以上求解模型和本文中设计的爆轰产物非线性拟合方程，计算爆温T和爆压p，与假设爆温T′和爆压p′比较，反复修正，直到满足假设值与计算所得值相同。

以上求解爆轰参数的数学模型利用自编程序实现，计算过程详见本课题组前期研究工作[19]，其程序算法框图见图1。

图1 计算模拟程序图Fig.1 Computer simulation program chart

3 爆轰参数计算结果与分析

利用非线性拟合方法求得密度为1.80 g/cm3的RDX和密度为1.14 g/cm3的CH3NO2的爆轰产物组成，如表3所示，爆轰参数如表4所示。

计算过程BKW状态方程参数[10]α取0.5，β取0.16，θ取400，κ取10.91。碳的单质生成物为金刚石[12-18]。RDX和CH3NO2的生成热分别是-70.66 kJ/mol和 113.09 kJ/mol。

为对比本文中计算方法的精确性，表4还分别给出爆轰参数的实验结果和依据Kast计算爆温所得结果。

表4显示，采用非线性拟合法作为计算爆温的方法计算爆轰参数，比使用Kast平均热容法计算的数值更加接近实验值。

对于RDX，非线性拟合法计算爆温、爆压和爆速与实验值的误差分别为3.3%、0.4%和0.3%；Kast平均热容法计算爆温、爆压和爆速与实验值的误差分别为6.4%、3.4%和2.0%。

表3 每摩尔炸药CJ点爆轰产物组成Tab.3 Detonation products of explosives at CJ point mol

表4 RDX和CH3NO2的爆温、爆压和爆速Tab.4 Detonation temperature，pressure and velocity of RDX and CH3NO2

对于CH3NO2，非线性拟合计算的爆温值也比Kast平均热容法更接近实验值。

4 结论

1)采用最小二乘法法得到爆轰产物热力学参数的非线性拟合方程，与Kast平均热容法比较，非线性拟合法在计算物质焓变(HT-H298)时显示出更好的计算精确性。

2)在此基础上，利用自编程序实现爆轰产物组成、爆压、爆温和爆速的计算。密度为1.80 g/cm3的RDX和密度为1.14 g/cm3的CH3NO2的爆轰参数计算结果显示，采取非线性拟合法计算的爆轰参数比采用Kast平均热容法计算的结果更加接近实验值，显示出更好的应用前景，对炸药爆轰性能预测具有一定参考价值。