爆轰产物组成和爆轰参数是决定炸药爆炸性能重要参数, 也是设计和改进炸药的重要依据。从20世纪40年代开始, 很多国外学者在爆轰产物状态方程研究领域作了大量的工作, 提出了很多的数学物理模型^[1], 对爆轰产物组成和爆轰参数进行研究, 如BKW^[2]、LJD、JCZ3及WCA模型, 并发展了有关的计算方法及计算程序, 如BKW、APPEGE、RUBY、LAMINEUR、TIGER、CHEQ^[3]程序都能较成功地解决这一领域的一些问题, 其中BKW、APPEGE、RUBY、LAMINEUR、TIGER等Fortran程序中应用了BKW状态方程, 而CHEQ程序中主要应用了液体状态方程的微扰理论, 并且还考虑了碳的石墨-金刚石-液碳三相组成。

从20世纪60年代至今, 随着计算机科学技术的发展, 中国学者在爆轰产物状态方程研究领域也作了大量的研究。X.Wu^[4]利用BKW状态方程模型, 从热力学着手, 确定各种可能混合物的状态方程, 计算混合物的自由能, 找出具有最小自由能的组成, 最终求出爆轰CJ点爆速、爆压及爆轰温度, 但计算爆轰温度和实验结果差别比较大, 而爆速的误差也达到500 m/s。李德华等^[5]利用WAC状态方程作为爆轰气相产物的物态方程, 考虑化学平衡条件和化学计算条件, 基于Gibbs自由能最小原理, 编程对化学平衡状态的炸药进行数值模拟计算, 对几种单质炸药的爆轰参数作预言, 得到的爆轰CJ点的爆速与实验值得误差小于3%、爆压与实验值的误差小于4%、爆轰温度与实验值的误差小于5%。赵艳红等^[6]采用van der Waals等效单组分流体模型和Ross硬球微扰理论软球修正模型计算爆轰气相产物的状态方程; 用石墨相、金刚石相、类石墨液相和类金刚石液相4种相态描述凝聚成分, 由Gibbs自由能最小确定不同状态下的凝聚产物相态, 对爆轰产物混合系统采用Gibbs自由能最小原理, 通过化学平衡方程组求解炸药爆轰产物系统的平衡组分, 计算结果与BKW和LJD的结果相近。

虽然针对爆轰产物组成的研究和编程很多, 但很少有人公开自己的程序, 且用基于最小自由能原理建立的非线性方程组研究爆轰产物和爆轰参数时, 需要计算出各产物的最小自由能, 求解过程复杂。本文中严格按照数学推导过程, 自行设计更易于求解的理论计算模型求解多元非线性方程, 选用BKW状态方程, 依据碳的相图和文献[8]选定碳的4种单质相态, 将一般爆轰压力(0≤p≤600 GPa)和温度(300 K≤T≤15 000 K)下游离态的碳当作金刚石处理, 采用自编程序对爆轰产物组成和爆轰参数进行计算。

1.   数值分析求解最小自由能组分

1.1   BKW状态方程

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

式中:p为压力; V_m为爆轰气体的摩尔体积; R为理想气体常数; T为爆轰气体的热力学温度; X=κ∑x_ik_i/[V_m(T+θ)^α], κ、α和θ为经验常数, k_i为第i种物质的余容, x_i为第i种物质的摩尔分数。BKW状态方程4个常数α、β、κ和θ的值通常如表 1所示, 本文中对欧洲民用爆炸物品系列标准^[9]取值进行修正, 修正前α、β、θ和κ对应的取值分别是0.5、0.298、6 620和10.50。通过自编程序计算发现, 若κ取9.50, 计算的爆轰参数与实验值更接近, 本文中采用欧洲标准修正后的取值作为BKW状态方程参数, 如表 1所示, 然后计算爆轰参数关系。

表  1  BKW状态方程参数取值

Table  1.  The parameter values of BKW equition of state

炸药	α	β	θ	κ
RDX	0.5	0.16	400	10.91
工业炸药	0.5	0.16	2 890	10.91
本文中	0.5	0.298	6 620	9.50

下载: 导出CSV 

| 显示表格

1.2   求解爆轰产物组分的推导过程

以往的研究人员根据化学平衡方程组建立的非线性方程组, 根据化学平衡条件推导出爆炸产物组成与自由能关系。在假定压力和温度下由状态方程或者高压状态下分子作用模型确定假定爆轰产物组分总自用能^[10-11], 根据求得假定组分自由能求出一组可能的爆轰产物组成, 由热力学关系确定假定爆轰组分的压力和温度, 再利用求出的可能爆轰产物组成求出爆轰组分总自由能, 如此利用程序反复迭代, 直到求得爆轰产物组分总自由能、压力和温度在迭代程序中达到一定精度, 最终得到爆轰产物组成、压力和温度。

本文中利用化学平衡条件和数学转化推导出非线性方程组, 根据质量守恒和物质自由能公式可得:

式中:N_k⁰为每千克炸药配方含有第k种元素的总物质的量, N_k为每千克爆炸产物含有第k种元素的总物质的量, a_ik为每摩尔第i种气相组分含有k元素的物质的量, n_ig为每千克爆炸产物含有第i种气相组分的物质的量, a_jk为每摩尔第j种凝聚相组分含有k元素的物质的量, n_js为每千克爆炸产物含有第j种凝聚相组分的物质的量, G为爆炸产物总自由能, g_i为每摩尔i气相组分的吉布斯自由能, g_j为每摩尔j凝聚相组分的吉布斯自由能。

利用拉格朗日乘数法把式(2)~(4)的条件极值转化为非条件极值, 根据拉乌尔定律和牛顿迭代法处理得到以下公式:

式中:x=(n_1g, …, n_ig, …, n_mg, n_1s, …, n_js, …, n_ns, λ₁, …, λ_k, …, λ_l)^T, λ_k为拉格朗日因子; f_i(n, λ_k)=∂G(n, λ_k)/∂n_ig, i=1, …, m; f_j+m(n, λ_k)=∂G(n, λ_k)/∂n_js, j=1, …, n; f_k+m+n(n, λ_k)=∂G(n, λ_k)/∂λ_k, k=1, …, l; x^(ν)代表第ν次迭代的爆轰组分物质量;

g_i⁰和g_j⁰分别为第i种气体和第j种凝聚相组分的吉布斯自由能。

利用数学的方法设计了计算爆轰产物组成的方法, 利用数学方法计算爆炸产物组成, 不需要考虑单个组分的自由能。前期学者研究爆轰产物组成考虑单个组分的自由能, 由于分子的作用势和自由能转换计算复杂, 在一定程度上会给计算带来误差。

2.   爆轰参数的确定

由能量守恒法则知爆炸前后系统能量守恒可得^[9]:

式中:E₀为炸药内能, E为爆炸产物内能, E_g为爆炸产物气体内能, E_s为爆炸产物固体内能。

式中: E_i⁰为第i种气体的内能, E_imp为BKW状态方程对应的第i种气体的非理想条件项, E_j⁰为第j种固体的内能。根据BKW状态方程, 可以推导出E_imp表达式^[9]:

炸药爆轰气相产物分子体系必须满足Hugniot关系^[5-7]:

式中:E₁为爆轰产物内能, p_H爆轰产物CJ条件的爆轰压力, p₀为初始压力, V_H为爆轰产物CJ条件下的体积, V₀为初始体积。爆速表达式为:

由于CJ点的爆速是最小的, 利用抛物线最小法由Hugoniot曲线上的3个点便可求得爆轰CJ点的爆速D_CJ。

3.   PETN的计算结果

为了验证本文计算方法的准确性, 采用具有实验数据的密度为1.77 g/cm³的PETN进行验证。本文中依据纯物质热化学数据手册^[12], 根据(1)~(13)式, 自编程序求得PETN在CJ爆轰点的爆轰参数。密度为1.77 g/cm³、初始摩尔内能为-399.572 kJ/mol的PETN炸药的计算结果见表 2和表 3。本文中得到的压力、爆速和温度的误差分别为0.01%、0.51%和0.21%。利用本文方法还能确定密度为1.77 g/cm³的PETN在CJ爆轰点的产物密度为2.43 g/cm³, 与以往学者认为的爆轰产物密度范围2~3 g/cm³相符^[6]。

表  2  每摩尔PETN炸药CJ点爆轰产物组成

Table  2.  Compositions of detonation products per mole of PETN explosive at CJ point

方法	n(CO2)/mol	n(CO)/mol	n(N2)/mol	n(O2)/mol	n(H2)/mol
BKW[1, 2]	3.950 9	0.096 2	1.999 2	1.359 1×10-5	0.568 9×10-4
LJD[1, 2]	3.779 0	0.446 9	1.969 5	0	0.480 4×10-1
WCA[13]	3.934 9	0.154 6	1.991 7	4.494 3×10-5	3.989 5×10-4
本文	3.999 5	0.001 1	1.999 9	1.190 4×10-5	2.077 7×10-4
方法	n(H2O)/mol	n(NH3)/mol	n(NO)/mol	n(CH4)/mol	n(C)/mol
BKW[1, 2]	3.998 6	1.327 5×10-4	2.244 1×10-4	0.537 3×10-5	0.951 6
LJD[1, 2]	3.929 1	0	0.619 5×10-1	1.197 9×10-2	0.759 8
WCA[13]	3.975 2	1.629 0×10-2	3.141 7×10-4	0	0.910 5
本文	3.999 7	7.354 8×10-8	4.530 7×10-6	5.482 0×10-6	0.999 1

下载: 导出CSV 

| 显示表格

表  3  PETN(1.77 g/cm³)炸药爆轰CJ点的爆轰参数

Table  3.  Detonation parameter of CJ for PETN explosive(1.77 g/cm³)

方法	p/GPa	D/(m·s-1)	T/K	方法	p/GPa	D/(m·s-1)	T/K
WCA[2]	31.5	8 416	4 349	LJD[6]	28.5	8 087	4 378
BKW[5]	31.8	8 421	2 833	WCA[7]	33.49	8 540	3 237
实验[2]	33.5	8 300	3 400	本文	33.53	8 342	3 473

下载: 导出CSV 

| 显示表格

4.   结论

(1) 本文中基于最小自由能原理, 自行设计理论计算模型, 求解爆轰产物和爆轰参数, 选取BKW方程作为产物状态方程, 状态方程参数采用欧洲标准取值, 并对欧洲标准参数取值进行修正。利用数值分析的方法, 推导求解爆轰产物的方程组, 在计算爆轰产物化学平衡组成时不需要单独计算每个产物组分的自由能, 具有优化计算减小误差的优点, 且容易编程实现。

(2) 根据自行推导的计算方法, 自编程序实现了爆轰参数计算, 用密度为1.77 g/cm³的PETN验证计算方法的准确性, 所得CJ点对应爆压、爆温和爆速与实验值的误差均小于1%, 并得到密度为1.77 g/cm³的PETN在CJ爆轰点的产物密度为2.43 g/cm³。