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

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

?

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

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

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

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

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

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

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

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

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

1.1 BKW狀態(tài)方程

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

(1)

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

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

1.2 求解爆轟產(chǎn)物組分的推導(dǎo)過(guò)程

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

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

(2)

(3)

(4)

利用拉格朗日乘數(shù)法把式(2)～(4)的條件極值轉(zhuǎn)化為非條件極值，根據(jù)拉烏爾定律和牛頓迭代法處理得到以下公式：

(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(ν)代表第ν次迭代的爆轟組分物質(zhì)量；

(8)

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

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

由能量守恒法則知爆炸前后系統(tǒng)能量守恒可得[9]:

(9)

式中：E0為炸藥內(nèi)能，E為爆炸產(chǎn)物內(nèi)能，Eg為爆炸產(chǎn)物氣體內(nèi)能，Es為爆炸產(chǎn)物固體內(nèi)能。

(10)

(11)

炸藥爆轟氣相產(chǎn)物分子體系必須滿足Hugniot關(guān)系[5-7]:

(12)

式中：E1為爆轟產(chǎn)物內(nèi)能，pH爆轟產(chǎn)物CJ條件的爆轟壓力，p0為初始?jí)毫Γ琕H為爆轟產(chǎn)物CJ條件下的體積，V0為初始體積。爆速表達(dá)式為:

(13)

由于CJ點(diǎn)的爆速是最小的，利用拋物線最小法由Hugoniot曲線上的3個(gè)點(diǎn)便可求得爆轟CJ點(diǎn)的爆速DCJ。

3 PETN的計(jì)算結(jié)果

為了驗(yàn)證本文計(jì)算方法的準(zhǔn)確性，采用具有實(shí)驗(yàn)數(shù)據(jù)的密度為1.77 g/cm3的PETN進(jìn)行驗(yàn)證。本文中依據(jù)純物質(zhì)熱化學(xué)數(shù)據(jù)手冊(cè)[12]，根據(jù)(1)～(13)式，自編程序求得PETN在CJ爆轟點(diǎn)的爆轟參數(shù)。密度為1.77 g/cm3、初始摩爾內(nèi)能為-399.572 kJ/mol的PETN炸藥的計(jì)算結(jié)果見(jiàn)表2和表3。本文中得到的壓力、爆速和溫度的誤差分別為0.01%、0.51%和0.21%。利用本文方法還能確定密度為1.77 g/cm3的PETN在CJ爆轟點(diǎn)的產(chǎn)物密度為2.43 g/cm3，與以往學(xué)者認(rèn)為的爆轟產(chǎn)物密度范圍2～3 g/cm3相符[6]。

表2 每摩爾PETN炸藥CJ點(diǎn)爆轟產(chǎn)物組成

表3 PETN(1.77 g/cm3)炸藥爆轟CJ點(diǎn)的爆轟參數(shù)

4 結(jié) 論

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

(2) 根據(jù)自行推導(dǎo)的計(jì)算方法，自編程序?qū)崿F(xiàn)了爆轟參數(shù)計(jì)算，用密度為1.77 g/cm3的PETN驗(yàn)證計(jì)算方法的準(zhǔn)確性，所得CJ點(diǎn)對(duì)應(yīng)爆壓、爆溫和爆速與實(shí)驗(yàn)值的誤差均小于1%，并得到密度為1.77 g/cm3的PETN在CJ爆轟點(diǎn)的產(chǎn)物密度為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] 趙艷紅,劉海風(fēng),張弓木.PETN炸藥爆轟產(chǎn)物狀態(tài)方程的理論研究[J].高壓物理學(xué)報(bào),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] 趙艷紅,劉海風(fēng),張弓木.基于統(tǒng)計(jì)物理的爆轟產(chǎn)物物態(tài)方程研究[J].物理學(xué)報(bào),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ó)外民爆器材法規(guī)和技術(shù)資料匯編[M].北京:中國(guó)爆破器材行業(yè)協(xié)會(huì),2010:290-298.

[10] 楊向東,謝文,武保劍.液氮的沖擊壓縮理論計(jì)算[J].高壓物理學(xué)報(bào),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].高壓物理學(xué)報(bào),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].四川師范大學(xué)學(xué)報(bào),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.

(責(zé)任編輯 王小飛)

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

國(guó)家自然科學(xué)基金項(xiàng)目(51134012)

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

O381 國(guó)標(biāo)學(xué)科代碼： 13035

A