杜明燃，汪旭光,2，郭子如，顏事龍

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

?

爆轟產物組成和爆轟參數(shù)計算方法的理論研究

杜明燃1，汪旭光1,2，郭子如3，顏事龍1

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

為實現(xiàn)爆轟產物組成和爆轟參數(shù)的計算，采用拉格朗日乘數(shù)法和牛頓迭代的方法預測爆轟產物組成，利用BKW狀態(tài)方程預測爆轟參數(shù)，在0～600 GPa和300～15 000 K壓力溫度范圍內選取金剛石作為碳的生成相；對爆轟產物系統(tǒng)采用最小自由能原理，結合牛頓迭代法求解爆轟產物的化學平衡方程組；對BKW狀態(tài)方程參數(shù)提出修訂，取α=0.5，β=0.298，θ=6 620，κ=9.50；采用自編程序實現(xiàn)計算過程。使用此方法和Hugoniot關系計算密度為1.77 g/cm3的PETN爆轟CJ點爆轟參數(shù)驗證計算精度，結果顯示計算與實驗結果的誤差均小于1%。利用此方法結合Hugoniot關系預測出爆轟CJ點的產物密度為2.43 g/cm3。

爆炸力學；爆轟參數(shù)；熱力學性質；狀態(tài)方程；吉布斯自由能

爆轟產物組成和爆轟參數(shù)是決定炸藥爆炸性能重要參數(shù)，也是設計和改進炸藥的重要依據。從20世紀40年代開始，很多國外學者在爆轟產物狀態(tài)方程研究領域作了大量的工作，提出了很多的數(shù)學物理模型[1]，對爆轟產物組成和爆轟參數(shù)進行研究，如BKW[2]、LJD、JCZ3及WCA模型，并發(fā)展了有關的計算方法及計算程序，如BKW、APPEGE、RUBY、LAMINEUR、TIGER、CHEQ[3]程序都能較成功地解決這一領域的一些問題，其中BKW、APPEGE、RUBY、LAMINEUR、TIGER等Fortran程序中應用了BKW狀態(tài)方程，而CHEQ程序中主要應用了液體狀態(tài)方程的微擾理論，并且還考慮了碳的石墨-金剛石-液碳三相組成。

從20世紀60年代至今，隨著計算機科學技術的發(fā)展，中國學者在爆轟產物狀態(tài)方程研究領域也作了大量的研究。X.Wu[4]利用BKW狀態(tài)方程模型，從熱力學著手，確定各種可能混合物的狀態(tài)方程，計算混合物的自由能，找出具有最小自由能的組成，最終求出爆轟CJ點爆速、爆壓及爆轟溫度，但計算爆轟溫度和實驗結果差別比較大，而爆速的誤差也達到500 m/s。李德華等[5]利用WAC狀態(tài)方程作為爆轟氣相產物的物態(tài)方程，考慮化學平衡條件和化學計算條件，基于Gibbs自由能最小原理，編程對化學平衡狀態(tài)的炸藥進行數(shù)值模擬計算，對幾種單質炸藥的爆轟參數(shù)作預言，得到的爆轟CJ點的爆速與實驗值得誤差小于3%、爆壓與實驗值的誤差小于4%、爆轟溫度與實驗值的誤差小于5%。趙艷紅等[6]采用van der Waals等效單組分流體模型和Ross硬球微擾理論軟球修正模型計算爆轟氣相產物的狀態(tài)方程；用石墨相、金剛石相、類石墨液相和類金剛石液相4種相態(tài)描述凝聚成分，由Gibbs自由能最小確定不同狀態(tài)下的凝聚產物相態(tài)，對爆轟產物混合系統(tǒng)采用Gibbs自由能最小原理，通過化學平衡方程組求解炸藥爆轟產物系統(tǒng)的平衡組分，計算結果與BKW和LJD的結果相近。

雖然針對爆轟產物組成的研究和編程很多，但很少有人公開自己的程序，且用基于最小自由能原理建立的非線性方程組研究爆轟產物和爆轟參數(shù)時，需要計算出各產物的最小自由能，求解過程復雜。本文中嚴格按照數(shù)學推導過程，自行設計更易于求解的理論計算模型求解多元非線性方程，選用BKW狀態(tài)方程，依據碳的相圖和文獻[8]選定碳的4種單質相態(tài)，將一般爆轟壓力(0≤p≤600 GPa)和溫度(300 K≤T≤15 000 K)下游離態(tài)的碳當作金剛石處理，采用自編程序對爆轟產物組成和爆轟參數(shù)進行計算。

1 數(shù)值分析求解最小自由能組分

1.1 BKW狀態(tài)方程

BKW狀態(tài)方程最初是由Becker提出的方程，后來經過Kistiakowsky和Wilson等的修正，最后確定為以下形式：

(1)

式中：p為壓力；Vm為爆轟氣體的摩爾體積；R為理想氣體常數(shù)；T為爆轟氣體的熱力學溫度；X=κ∑xiki/[Vm(T+θ)α]，κ、α和θ為經驗常數(shù)，ki為第i種物質的余容，xi為第i種物質的摩爾分數(shù)。BKW狀態(tài)方程4個常數(shù)α、β、κ和θ的值通常如表1所示，本文中對歐洲民用爆炸物品系列標準[9]取值進行修正，修正前α、β、θ和κ對應的取值分別是0.5、0.298、6 620和10.50。通過自編程序計算發(fā)現(xiàn)，若κ取9.50，計算的爆轟參數(shù)與實驗值更接近，本文中采用歐洲標準修正后的取值作為BKW狀態(tài)方程參數(shù)，如表1所示，然后計算爆轟參數(shù)關系。

表1 BKW狀態(tài)方程參數(shù)取值

1.2 求解爆轟產物組分的推導過程

以往的研究人員根據化學平衡方程組建立的非線性方程組，根據化學平衡條件推導出爆炸產物組成與自由能關系。在假定壓力和溫度下由狀態(tài)方程或者高壓狀態(tài)下分子作用模型確定假定爆轟產物組分總自用能[10-11]，根據求得假定組分自由能求出一組可能的爆轟產物組成，由熱力學關系確定假定爆轟組分的壓力和溫度，再利用求出的可能爆轟產物組成求出爆轟組分總自由能，如此利用程序反復迭代，直到求得爆轟產物組分總自由能、壓力和溫度在迭代程序中達到一定精度，最終得到爆轟產物組成、壓力和溫度。

本文中利用化學平衡條件和數(shù)學轉化推導出非線性方程組，根據質量守恒和物質自由能公式可得：

(2)

(3)

(4)

利用拉格朗日乘數(shù)法把式(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)

利用數(shù)學的方法設計了計算爆轟產物組成的方法，利用數(shù)學方法計算爆炸產物組成，不需要考慮單個組分的自由能。前期學者研究爆轟產物組成考慮單個組分的自由能，由于分子的作用勢和自由能轉換計算復雜，在一定程度上會給計算帶來誤差。

2 爆轟參數(shù)的確定

由能量守恒法則知爆炸前后系統(tǒng)能量守恒可得[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的計算結果

為了驗證本文計算方法的準確性，采用具有實驗數(shù)據的密度為1.77 g/cm3的PETN進行驗證。本文中依據純物質熱化學數(shù)據手冊[12]，根據(1)～(13)式，自編程序求得PETN在CJ爆轟點的爆轟參數(shù)。密度為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點的爆轟參數(shù)

4 結 論

(1) 本文中基于最小自由能原理，自行設計理論計算模型，求解爆轟產物和爆轟參數(shù)，選取BKW方程作為產物狀態(tài)方程，狀態(tài)方程參數(shù)采用歐洲標準取值，并對歐洲標準參數(shù)取值進行修正。利用數(shù)值分析的方法，推導求解爆轟產物的方程組，在計算爆轟產物化學平衡組成時不需要單獨計算每個產物組分的自由能，具有優(yōu)化計算減小誤差的優(yōu)點，且容易編程實現(xiàn)。

(2) 根據自行推導的計算方法，自編程序實現(xiàn)了爆轟參數(shù)計算，用密度為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炸藥爆轟參數(shù)的數(shù)值模擬[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炸藥爆轟產物狀態(tài)方程的理論研究[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] 趙艷紅,劉海風,張弓木.基于統(tǒng)計物理的爆轟產物物態(tài)方程研究[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] 楊祖一.國外民爆器材法規(guī)和技術資料匯編[M].北京:中國爆破器材行業(yè)協(xié)會,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] 劉福生,陳先猛,陳攀森,等.液態(tài)CO2高溫高密度狀態(tài)方程研究[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炸藥爆轟參數(shù)的數(shù)值模擬[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