杜明燃，汪旭光,2，郭子如，颜事龙

(1.安徽理工大学土木建筑学院，安徽 淮南 232001； 2.北京矿冶研究总院，北京 100044； 3.安徽理工大学化学工程学院，安徽 淮南 232001)



爆轰产物组成和爆轰参数计算方法的理论研究

杜明燃1，汪旭光1,2，郭子如3，颜事龙1

(1.安徽理工大学土木建筑学院，安徽 淮南 232001； 2.北京矿冶研究总院，北京 100044； 3.安徽理工大学化学工程学院，安徽 淮南 232001)

为实现爆轰产物组成和爆轰参数的计算，采用拉格朗日乘数法和牛顿迭代的方法预测爆轰产物组成，利用BKW状态方程预测爆轰参数，在0～600 GPa和300～15 000 K压力温度范围内选取金刚石作为碳的生成相；对爆轰产物系统采用最小自由能原理，结合牛顿迭代法求解爆轰产物的化学平衡方程组；对BKW状态方程参数提出修订，取α=0.5，β=0.298，θ=6 620，κ=9.50；采用自编程序实现计算过程。使用此方法和Hugoniot关系计算密度为1.77 g/cm3的PETN爆轰CJ点爆轰参数验证计算精度，结果显示计算与实验结果的误差均小于1%。利用此方法结合Hugoniot关系预测出爆轰CJ点的产物密度为2.43 g/cm3。

爆炸力学；爆轰参数；热力学性质；状态方程；吉布斯自由能

爆轰产物组成和爆轰参数是决定炸药爆炸性能重要参数，也是设计和改进炸药的重要依据。从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等的修正，最后确定为以下形式：

(1)

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

表1 BKW状态方程参数取值

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

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

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

(2)

(3)

(4)

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

(5)

(6)

(7)

式中：x=(n1g,…,nig,…,nmg,n1s,…,njs,…,nns,λ1,…,λk,…,λl)T，λk为拉格朗日因子；fi(n,λk)=∂G(n,λk)/∂nig,i=1,…,m；fj+m(n,λk)=∂G(n,λk)/∂njs,j=1,…,n；fk+m+n(n,λk)=∂G(n,λk)/∂λk,k=1,…,l；x(ν)代表第ν次迭代的爆轰组分物质量；

(8)

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

2 爆轰参数的确定

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

(9)

式中：E0为炸药内能，E为爆炸产物内能，Eg为爆炸产物气体内能，Es为爆炸产物固体内能。

(10)

(11)

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

(12)

式中：E1为爆轰产物内能，pH爆轰产物CJ条件的爆轰压力，p0为初始压力，VH为爆轰产物CJ条件下的体积，V0为初始体积。爆速表达式为:

(13)

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

3 PETN的计算结果

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

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

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

4 结 论

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

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

[1] Chirat R, Pittion-Rossillon G. A new equation of state for detonation products[J]. Journal of Chemical Physics, 1981,74(8):4634-4645.

[2] Mader C L. Numerical modeling of explosives and propellants[M]. 3rd Edtion. Boca Raton, FL, USA: CRC Press, 2008:31-63.

[3] Thiel M V, Ree F H. Nonequilibrium effects of slow diffusion controlled reactions on the properties of explosives [C]∥Proceedings of 9th International Symposium on Detonation. Portland, USA: Oregon, 1991:743-750.

[4] Wu X. BKW equation of state for detonation products[C]∥Proceedings of 8th International Symposium on Detonation. Portland, USA: Oregon, 1985:438.

[5] 李德华,程新路,杨向东,等.PETN、RDX和HMX炸药爆轰参数的数值模拟[J].爆炸与冲击,2005,25(4):325-328. Li De-hua, Cheng Xin-lu, Yang Xiang-dong, et al. Numerical simulation of detonation parameters for PETN, RDX and HMX explosives[J]. Explosion and Shock Waves, 2005,25(4):325-328.

[6] 赵艳红,刘海风,张弓木.PETN炸药爆轰产物状态方程的理论研究[J].高压物理学报,2009,23(2):143-149. Zhao Yan-hong, Liu Hai-feng, Zhang Gong-mu. Equation of state of detonation products for PETN explosive[J]. Journal of High Pressure Physics, 2009,23(2):143-149.

[7] Fried L E, Howard W M. Explicit Gibbs free energy equation of state applied to the carbon phase diagram[J]. Physical Review, 2000,B61(13):8734-8743.

[8] 赵艳红,刘海风,张弓木.基于统计物理的爆轰产物物态方程研究[J].物理学报,2007,56(8):4791-4797. Zhao Yan-hong, Liu Hai-feng, Zhang Gong-mu. Equation of state of detonation products based on statistical mechanical theory[J]. Acta Physica Sinica, 2007,56(8):4791-4797.

[9] 杨祖一.国外民爆器材法规和技术资料汇编[M].北京:中国爆破器材行业协会,2010:290-298.

[10] 杨向东,谢文,武保剑.液氮的冲击压缩理论计算[J].高压物理学报,1998,12(1):1-5. Yang Xiang-dong, Xie Wen, Wu Bao-jian. Theoretical calculation for the Hugoniot curves of liquid nitrogen[J]. Journal of High Pressure Physics, 1998,12(1):1-5.

[11] 刘福生,陈先猛,陈攀森,等.液态CO2高温高密度状态方程研究[J].高压物理学报,1998,12(1):28-33.Liu Fu-sheng, Chen Xian-meng, Chen Pan-sen, et al. Equation of state of liquid CO2at high temperatures and high densities[J]. Journal of High Pressure Physics, 1998,12(1):28-33.

[12] Ihsan B, Gregor P. Thermochemical Data of Pure Substances: Volumn 2 [M].3rd Edtion. New York: VCH, 1995.

[13] 李德华,程新路,杨向东,等.PETN炸药爆轰参数的数值模拟[J].四川师范大学学报,2005,28(4):448-451. Li De-hua, Cheng Xin-lu, Yang Xiang-dong, et al. Numerical simulation of detonation parameters for PETN explosive[J]. Journal of Sichuan Normal University, 2005,28(4):448-451.

(责任编辑 王小飞)

Theoretical studies for calculating the detonation products and properties of explosives

Du Ming-ran1, Wang Xu-guang1,2, Guo Zi-ru3, Yan Shi-long1

(1.SchoolofCivilEngineeringandArchitecture,AnhuiUniversityofScienceandTechnology,Huainan232001,Anhui,China; 2.BeijingGeneralResearchInstituteofMiningandMetallurgy,Beijing100044,China; 3.SchoolofChemicalEngineering,AnhuiUniversityofScienceandTechnology,Huainan232001,Anhui,China)

In order to calculate the detonation products and parameters, Lagrange multiplier and Newton iterative method were used to predict detonation products. The state equation of BKW was used to predict detonation parameters. In a range of pressure from 0 to 600 GPa and temperature from 300 to 15 000 K, diamond was intended as the elemental carbon product. Based on the principle of minimum free energy, the equilibrium compositions of detonation products were calculated by using Newton iterative method, which need not calculate the free energy of each composition. The parameters of the state equation of BKW were modified.α=0.5;β=0.298;θ=6 620;κ=9.50. Using self-made program, the detonation properties at CJ point of PETN, whose density is 1.77 g/cm3, were calculated with the theory in this paper and the equation of Hugoniot. The results show satisfactory agreement with the experimental data, with the error less than 1%. The density of detonation products is also predicted easily. When the density of PETN is 1.77 g/cm3, the density of detonation products is 2.43 g/cm3.

mechanics of explosion; detonation parameters; thermodynamic property; equation of state; Gibbs free energy

10.11883/1001-1455(2015)04-0449-05

2013-04-12；

2014-02-27

国家自然科学基金项目(51134012)

杜明燃(1987- )，男，博士研究生，dumingranaust@163.com。

O381 国标学科代码： 13035

A