基于BKW狀態方程的爆轟產物及參數的改進算法

何偉平，黃 菊，陳厚和，劉曉靜，王德堂,4

(1.徐州工業職業技術學院化學工程學院，江蘇 徐州 221140； 2.徐州工程學院化學化工學院，江蘇 徐州 221111； 3.南京理工大學化工學院，江蘇 南京 210094； 4.徐州工業職業技術學院江蘇省化工新材料工程技術研發開發中心，江蘇 徐州 221140)

基于BKW狀態方程的爆轟產物及參數的改進算法

何偉平1，黃 菊2，陳厚和3，劉曉靜1，王德堂1,4

(1.徐州工業職業技術學院化學工程學院，江蘇 徐州 221140； 2.徐州工程學院化學化工學院，江蘇 徐州 221111； 3.南京理工大學化工學院，江蘇 南京 210094； 4.徐州工業職業技術學院江蘇省化工新材料工程技術研發開發中心，江蘇 徐州 221140)

為了降低爆轟產物及爆轟參數的求解難度，通過對質量守恒方程的基本可行解進行線性組合，得到了爆轟產物的平衡組成，并在此基礎上進一步獲得了爆轟參數。其主要實現方法為：由最小自由能原理對基本可行解進行篩選，然后根據最大放熱原則確定初始解，并在最小自由能原則的引導下，由初始解和基本可行解的線性組合獲得爆轟產物的平衡組成，以上操作步驟均由自編程序完成。應用支持向量機(SVM)線性模型對BKW狀態方程參數進行了調整，并詳細介紹了其主要步驟。使用此方法預測了PETN、CL-20和含鋁炸藥的爆轟產物及爆轟參數，經參數調整后，發現預測結果與實驗值吻合良好；通過與單質炸藥爆轟實驗數據對比，發現調整BKW狀態方程參數時，應當盡可能使用爆轟產物中氣體含量相近的含能材料對SVM模型進行訓練；若預測含鋁炸藥，應當使用鋁氧比接近待測炸藥的樣品來訓練SVM模型。

爆炸力學；BKW狀態方程；基本可行解；支持向量機；SVM；吉布斯自由能

引 言

含能材料爆轟參數及產物的計算離不開爆轟產物的狀態方程[1]。從20世紀40年代起，各國研究者對爆轟產物狀態方程作了大量研究，提出了各種數學物理模型，如美國的BKW、JCZ、JWL等[2-4]，日本的KHT[5]，前蘇聯的BKWRR，以及我國的VLW等[6-8]，并根據狀態方程發展了各種計算程序[3,9-11]。其中，BKW狀態方程作為典型的爆轟產物狀態方程，在國際上取得了廣泛的應用和發展[12-14]。

由于BKW狀態方程對維里方程取一級近似，對工業炸藥和高能炸藥往往不能準確預測[1,15]，研究者通常采取調整系數的方法來進行改進[16-18]，但關于系數的調整方法卻鮮見報道。另一方面，通常采用最小自由能法來求解爆轟產物及參數問題，并引入拉格朗日因子將條件極值轉化為非條件極值，通過反復迭代進行計算，過程復雜而且會帶來誤差[13,18]，甚至會因為收斂困難而導致求解失敗。由于傳統的求解方式借助牛頓迭代法完成，其收斂速度與目標函數的偏導函數有關，故具有一定的不確定性，并且涉及到偏導函數的計算，其過程十分復雜。本文在前期研究[19]的基礎上，利用爆轟產物的基本可行解，通過簡單的線性組合來獲取爆轟產物的平衡組分，該計算方法原理可靠、過程簡潔，得到的爆轟產物組成嚴格遵守質量守恒原則；在此基礎上計算爆轟參數，可降低計算的復雜程度和避免迭代帶來的誤差；此外，采用支持向量機理論(SVM)[20]，結合自編程序調整了BKW狀態方程參數。

1 爆轟產物和參數的計算方法

1.1 問題描述

采用最小自由能法求解爆轟產物，其數學模型可簡述如下：

(1)

式中：下標g、s分別表示氣態產物和凝聚態產物；G為吉布斯自由能；x為產物的物質的量；μ為化學勢；aik、ajk分別為每摩爾第i、j種產物中含有第k種原子的物質的量；bk為體系第k種原子的總物質的量；E為爆轟產物比內能；E0為炸藥比內能；P、ν分別為爆轟狀態的壓強、產物比體積；P0、ν0分別為初始的壓強、炸藥比體積。其中，μgi、μsj、E均為關于溫度(T)、壓強(P)、產物組成(xg1，…，xgm，xs1，…，xgn)的函數，其具體表達形式見文獻[16]。

1.2 基本可行解的處理

計算爆轟參數的關鍵是確定爆轟產物的平衡組成。在以往研究中，通常采用拉格朗日乘數法將式(1)的質量守恒方程與吉布斯自由能組合為新的目標函數，通過牛頓迭代法，將Hugniot關系作為終止判據進行求解，但其計算過程極為繁瑣[18]。本研究直接求解質量守恒方程，利用線性組合搜索最優解，降低了計算的復雜程度；并通過優化搜索方法，保證了計算效率和精度。

該算法的基本原理為：爆轟產物的任意一種可能組成Xi=(xgli，…，xgmi，xsli，…，xgni)，均可視為式(1)中質量守恒方程的(m+n-l+1)個基本可行解{X1，X2，…，Xm+n-l+1}的線性組合。本研究基本可行解是指：從(m+n)種產物中選取l種，代入求解質量守恒方程，若解存在，則該解稱為基本可行解。

通常在假定溫度和壓強下求解平衡組分，此時吉布斯自由能為X=(xg1，…，xgm，xsl，…，xgn)的連續函數；另一方面，爆轟產物的平衡組分(Xopt)處于吉布斯自由能的最低點，故在其他條件未知的情況下，選取吉布斯自由能最小的(m+n-l+1)個基本可行解進行線性組合較為合理。

1.3 組分迭代

初始解(X0)的選取也很重要，若將X0設定在Xopt附近，則有利于求解。同時考慮程序編寫的簡便性，本研究采用“最大放熱原則”確定X0，計算結果表明可以很好地收斂到Xopt。

將篩選出的(m+n-l+1)個基本可行解按其吉布斯自由能升序排列，依次與X0線性組合：Xt1=(1-r)X0+rXi(0

若一輪比較結束后，X0仍未發生變化，則令r′=0.5r(r初始值設為0.1)，將r′作為新的r，進入下一輪比較和替換，直至r達到設定的精度則停止搜索，即認為得到該溫度和壓強下的最優解Xopt=X0。

1.4 溫度和壓強迭代

獲得某溫度和壓強下的最優解Xopt后，計算f(P,T)=(E-E0)-0.5(P+P0)(ν0-ν)。若f(P,T)>0，則降低溫度；若f(P,T)<0，則升高溫度。重新進行平衡組分的計算，直至f(P,T)達到預定的精度，該操作可采用二分法完成。

爆轟產物的比體積(ν)可以由BKW狀態方程和凝聚態產物的狀態方程(如Cowan方程)計算得到。爆轟產物平衡組分、爆轟溫度分別由前述方式確定的條件下，改變P，即可得到相應的ν。由于Hugoniot曲線上C-J點的爆速最小，利用最小拋物線法反復迭代，即可確定爆轟C-J點[18]。

2 BKW狀態方程的參數調整

2.1 獲取訓練數據

根據平衡組分的搜索方法，利用Visual Basic 6.0編制爆轟參數計算代碼，涉及產物的熱力學數據來自美國國家標準與技術局(NIST)數據庫。保持α=0.5，選取β=0.140～0.180，θ=300～500，κ=10～12，根據{β，θ，κ}的不同組合，計算若干個炸藥的爆壓(P)和爆速(D)并與實驗值比較，將結果記為“偏低”和“偏高”兩類。

2.2 構造SVM模型

支持向量機(SVM)線性規劃模型[21]可表示為式(2)，其對應最優超平面為A(ω+-ω-)-γ=0。

(2)

(3)

構造向量C、W(待求向量)、b、LB及系數矩陣A(上標M、H表示所屬分類)：

(4)

(5)

2.3 訓練結果

根據構造的SVM模型，利用Matlab求解并繪制最佳分割平面，其結果如圖1所示。

由圖1(a)～(f)可見，利用SVM線性模型可以完整地劃分兩類結果(“偏低”、“偏高”)，故其最佳分割平面即為BKW狀態方程參數的合理取值范圍。由圖1(g)～(h)可知，各最佳分割平面沿β、κ方

向升降較為明顯，而沿θ方向較為平坦；故在調整參數時，θ的調整幅度不宜過小。另外，對比圖1(g)～(h)可見各最佳分割平面的相互關系并無一致規律，其具體參數見表1。

表1 最佳分割平面的相關參數Table 1 Relevant parameters of the optimal classifying planes

2.4 數據處理

為了使調整后的參數{β，θ，κ}盡可能滿足同類含能材料爆轟參數的計算，選取點(β，θ，κ)使其與表1中6個最佳分割平面的距離di總和為最小，即：

(6)

保持α不變，可得β=0.184，θ=1080，κ=10.76，以此作為BKW狀態方程參數進行爆轟產物及參數的計算。

3 算法驗證及優化

3.1 PETN和CL-20炸藥爆轟參數的驗證及優化

為了驗證本算法的準確性，采用標準生成焓為-502.8kJ/mol、密度為1.77g/cm3的PETN[22]，以及標準生成焓為460.0kJ/mol、密度為2.04g/cm3的CL-20[22-24]進行驗證，計算得到爆轟參數見表2。

由表2可見，PETN的預測爆壓準確性雖高于BKW方法，但低于文獻方法；預測爆速的準確性雖高于文獻方法，但低于BKW方法。分析其原因，可能是BKW狀態方程本身存在缺陷[1]，使得調整的參數不具備普適性。故需要根據類似炸藥或自身實驗進一步調整BKW狀態方程參數。

表2 PETN和CL-20炸藥爆轟CJ點的爆轟參數Table 2 Detonation parameters of PETN and CL-20 at C-J point

參照爆轟產物中氣體所占的質量分數以選取BKW參數的方法，單位質量PETN、RDX、HMX、CL-20、TATB生成的氣體產物質量分別為1.206、1.172、1.172、1.077、0.971。CL-20的產氣質量位于RDX、HMX與TATB之間，故BKW參數不再進行調整，CL-20的預測值與文獻數據相比，爆壓、爆速的誤差分別為4.65%、1.14%；但PETN與TATB差異較大，與RDX、HMX則較為接近，故選取點(β，θ，κ)使其與RDX、HMX涉及的4個最佳分割平面的距離總和最短，可得β=0.176、θ=540、κ=10.81。重新計算PETN爆轟參數得P=34.0GPa，D=8388m/s，T=2963K，PETN爆轟產物物質的量見表3。經過參數調整后，本方法計算爆壓、爆速的誤差分別為1.49%、1.06%，可見與實驗值吻合良好。

表3 1 摩爾 PETN炸藥C-J點的爆轟產物物質的量Table 3 The amount of detonation product for 1 mol PETN at C-J point

實際上，若利用PETN的實驗值對不同BKW狀態方程參數組合下的爆壓、爆速和爆溫進行檢驗，可得SVM最佳分割平面如圖2所示，BKW狀態方程參數的較好取值為(β,θ,κ)=(0.149,1425,12.26)。經自編程序驗證，在該BKW狀態方程參數取值下，PETN的爆壓、爆速和爆溫預測值分別為P=34.0GPa、D=8326m/s、T=3432K，誤差分別為1.49%、0.31%、0.94%。

利用CL-20的文獻數據對不同BKW狀態方程參數組合下的爆壓和爆速進行檢驗，可得SVM最佳分割平面如圖3所示，BKW狀態方程參數的較

好取值為(β,θ,κ)=(0.148,15,10.76)。在該BKW狀態方程參數取值下，CL-20的爆壓、爆速預測值分別為P=44.0GPa、D=9600m/s，誤差分別為2.3%、0。

3.2 含鋁炸藥爆轟參數的驗證

與前述理想炸藥不同，由于含鋁炸藥中的鋁顆粒在爆轟過程中的化學反應滯后于C-J爆轟反應，因此一般認為鋁粉在爆轟過程中是惰性的，故求解過程中鋁粉按惰性物質處理。為了方便與實驗數據對比，選取文獻[25]編號為2～6的RDX基含鋁炸藥為研究對象，其配方及爆轟參數測定結果[25]如表4所示。

表4 RDX基含鋁炸藥配方及爆轟參數測定結果Table 4 The formulation and the testing results of detonation parameters of RDX-based aluminized explosives

保持α=0.5，β=0.160，選取θ=2500～3500，κ=10.00～12.00，根據{θ，κ}的不同組合，計算4號含鋁炸藥的爆壓和爆速，并與實驗值進行對比分類，結果見圖4。

在觀測范圍內，(θ,κ) 的一個較好取值為(3100, 10.75)。以此作為BKW狀態方程參數，計算得到2～6號含鋁炸藥的爆轟參數見表5。

表5 RDX基含鋁炸藥爆轟參數的計算結果Table 5 The calculation results of detonation parameters for RDX-based aluminized explosive

計算結果表明，當鋁氧比接近4號含鋁炸藥時，爆壓的預測值與實測值較為接近，反之預測值誤差較大；因此，使用調整BKW參數的方法預測含鋁炸藥的爆壓時，應當在待測炸藥的鋁氧比接近訓練樣本的情況下使用。爆速的預測值則較為準確，最大誤差為5.90%，而且預測值與實測值變化趨勢基本一致。

4 結 論

(1)介紹了利用SVM模型優化BKW狀態方程參數的方法。直接求解質量守恒方程，利用其基本可行解進行線性組合，并充分結合最小自由能原理，可以方便地求解爆轟產物的平衡組分。

(2)根據自編程序驗證，使用SVM線性模型可以有效地獲得BKW狀態方程參數的合理取值。

(3)通過對比發現，利用BKW狀態方程預測含能材料的爆轟參數時，應當使用爆轟產物中氣體所占質量分數相近的含能材料的爆轟參數訓練SVM模型(若預測含鋁炸藥，則應當使用鋁氧比接近待測炸藥的樣品訓練SVM模型)，獲取的優化參數才具有適用性。

[1] 吳雄，龍新平，何碧，等.VLW狀態方程的回顧與展望[J].高壓物理學報，1999，13(1)：55-58. WU Xiong, LONG Xin-ping, HE Bi, et al.Preview and look forward to the progress of VLW equation of state[J]. Chinese Journal of High Pressure Physics, 1999, 13(1): 55-58.

[2] Mader C L.Recent advances in numerical modeling of detonations[J].Propellants, Explosives, Pyrotechnics, 1986, 11(6): 163-166.

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

[4] 宋浦, 楊凱, 梁安定,等. 國內外TNT炸藥的JWL狀態方程及其能量釋放差異分析[J]. 火炸藥學報, 2013, 36(2):42-45. SONG Pu, YANG Kai, LIANG An-ding, et al.Difference analysis on JWL-EOS and energy release of different TNT charge [J].Chinese Journal of Explosives & Propellants(Huozhayao Xuebao), 2013, 36(2):42-45.

[5] Katsum, Tanaka.Detonation properties of high explosives calculated by revised Kihara-Hikita equation of state[C]∥The Eighth Symposium(International) on Detonation. Annapolis:Naval Surface Weapons Center,1985: 548-557.

[6] 高大元, 呂春緒, 董海山,等. 工業炸藥的爆轟性能研究[J]. 火炸藥學報, 2003, 26(1):26-29. GAO Da-yuan, Lü Chun-xu, DONG Hai-shan, et al.Study on detonation performance of industrial explosives [J].Chinese Journal of Explosives & Propellants(Huozhayao Xuebao),2003, 26(1):26-29.

[7] 龍新平，何碧，蔣小華，等.論VLW狀態方程[J].高壓物理學報，2003，17(4)：247-254. LONG Xin-ping, HE Bi, JIANG Xiao-hua, et al.Discussions on the VLW equation of state[J].Chinese Journal of High Pressure Physics, 2003, 17(4): 247-254.

[8] Wu Xiong, Long Xin-Ping, He Bi, et al.VLW equation of state of detonation products[J].Science in China Series B：Chemistry, 2009, 52(5): 605-608.

[9] Howard W M, Fried L E, Souers P K.Kinetic modeling of non-ideal explosives with CHEETAH, DE200114716[R]. Livemore:LLNL,2001.

[10] 魏賢鳳, 龍新平, 韓勇. VLWR程序計算CHNO炸藥爆轟性能碳相態的選擇[J]. 含能材料, 2013, 21(5)：604-608. WEI Xian-feng, LONG Xin-ping, HAN Yong. Selection of carbon phase in calculation of detonation performance by VLWR program for CHNO explosives[J]. Chinese Journal of Energetic Materials, 2013, 21(5)：604-608.

[11] 張為鵬, 畢福強, 王永順,等. 1,1′-二羥基-5,5′-聯四唑二羥胺鹽理論爆速的計算[J]. 火炸藥學報, 2015, 38(6):67-71. ZHANG Wei-peng, BI Fu-qiang, WANG Yong-shun, et al. Calculation of theory detonation velocity of dihydroxylammonium 5,5′-bistetrazole-1,1′-diolate [J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao), 2015, 38(6):67-71.

[12] 韋秉旭，陳壽如.用BKW規則和BKW程序計算乳化炸藥爆轟產物的比較[J].湘潭師范學院學報(自然科學版)，2003，25(4)：50-54. WEI Bing-xu, CHEN Shou-ru.Comparison of detonation products of emulsion explosives computed with B-W regulation and BKW program[J].Journal of Xiangtan Normal University(Natural Science Edition), 2003, 25(4): 50-54.

[13] 李德華.炸藥爆轟參數、生成熱及爆熱的理論研究[D].成都：四川大學，2005.

[14] Muhamed Suceska, How-Ghee Ang, Hay Yee Serene Chan.Study of the effect of covolumes in BKW equation of state on detonation properties of CHNO explosives[J].Propellants, Explosives,Pyrotechnics, 2013, 38: 103-112.

[15] 韓早，王伯良.混合炸藥爆速預報的新方法[J].爆炸與沖擊，2014，34(4)：421-426. HAN Zao, WANG Bo-liang.A new method for predicting detonation velocity of composite explosive[J].Explosion and Shock Waves, 2014, 34(4): 421-426.

[16] 王小紅，李曉杰，李瑞勇，等.特殊炸藥的爆轟參數計算[J].高壓物理學報，2013，27(1)：119-124. WANG Xiao-hong, LI Xiao-jie, LI Rui-yong, et al.Calculation of detonation parameters of special explosives[J].Chinese Journal of High Pressure Physics, 2013, 27(1): 119-124.

[17] 韋秉旭，吳雄，唐健軍.用BKW狀態方程計算乳化炸藥爆轟參數的研究[J].爆破器材，2001，30(6)：1-4. WEI Bing-xu, WU Xiong, TANG Jian-jun, et al.The calculation research of emulsion explosive detonation parameters with BKW equation[J].Explosive Materials, 2001, 30(6): 1-4.

[18] 杜明燃，汪旭光，郭子如，等.爆轟產物組成和爆轟參數計算方法的理論研究[J].爆炸與沖擊，2015，35(4)：449-453. DU MING-ran, WANG Xu-guang, GUO Zi-ru, et al.Theoretical studies for calculating the detonation products and properties of explosives[J].Explosion and Shock Waves, 2015, 35(4): 449-453.

[19] 何偉平, 黃菊, 劉曉靜, 等. 基于多維空間的燃燒平衡產物組成的計算方法[J]. 火炸藥學報, 2016, 39(4)：46-50. HE Wei-ping, HUANG Ju, LIU Xiao-jing, et al. Calculation method of combustion equilibrium product compositions based on multidimensional space [J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao),2016, 39(4)：46-50.

[20] Narwaria, Manish, Lin Wei-si.Objective image quality assessment based on support vector regression[J].Neurai Networks, 2012, 21(3): 515-519.

[21] Bosch R, Smith J A.Separating hyperplanes and the authorship of the disputed federalist papers[J].American Mathematical Monthly, 1998, 105(7): 601-608.

[22] 歐育湘, 王艷飛, 劉進全, 等. 近20年問世的5個新高能量密度化合物[J]. 化學通報, 2006, 69：1-7. OU Yu-xiang, WANG Yan-fei, LIU Jin-quan, et al. Five novel high energy density compounds synthesized in the latest 20 years[J]. Chemistry Bulletin, 2006, 69：1-7.

[23] 歐育湘, 孟征, 劉進全. 高能量密度化合物CL-20應用研究進展[J]. 化工進展, 2007, 26(12)：1690-1694. OU Yu-xiang, MENG Zheng, LIU Jin-quan. Review of the development of application technologies of CL-20[J]. Chemical Industry and Engineering Progress, 2007, 26(12)：1690-1694.

[24] 陳魯英, 楊培進, 張林軍, 等. CL-20炸藥性能研究[J]. 火炸藥學報, 2003, 26(3)：65-67. CHEN Lu-ying, YANG Pei-jin, ZHANG Lin-jun, et al. Study of the performance of explosive CL-20[J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao),2003, 26(3)：65-67.

[25] 項大林, 榮吉利, 李健, 等. 黑索今基含鋁炸藥的鋁氧比對爆轟性能及其水下爆炸性能的影響[J]. 兵工學報, 2013, 34(1): 45-50. XIANG Da-lin, RONG Ji-li, LI Jian, et al. Effect of Al/O ratio on detonation performance and underwater explosion of RDX-based alumiized explosive[J]. Acta Armamentarii, 2013, 34(1): 45-50.

Improved Algorithm of Detonation Products and Parameters Based on the BKW Equation of State

HE Wei-ping1, HUANG Ju2, CHEN Hou-he3, LIU Xiao-jing1, WANG De-tang1,4

(1. School of Chemical Engineering，Xuzhou College of Industrial Technology，Xuzhou Jiangsu 221140, China；2. School of Chemistry & Chemical Engineering, Xuzhou Institute of Technology, Xuzhou Jiangsu 221111,China; 3. School of Chemical Engineering, Nanjing University of Science & Technology, Nanjng 210094,China; 4.Jiangsu Province Engineering Technology Research and Development Center of New Chemical Materials, Xuzhou College of Industrial Technology, Xuzhou Jiangsu 221140, China)

To reduce the difficulty of predicting the detonation products and solving detonation parameters, the equilibrium compositions of detonation products were achieved by linear combination of the basic feasible solutions, which were obtained from the mass conservation equations; and the detonation parameters were further obtained based on equilibrium compositions. The major process was executed as follows: the basic feasible solutions were selected out by the principle of minimum free energy, and the initial solution was given by the principle of largest heat release. The equilibrium compositions of detonation products were linearly searched by uniting the initial solution with the basic feasible solutions, and the above-mentioned operation steps were completed by using self-made program. The parameters of the BKW equation of state were adjusted applying the linear support vector machine (SVM), and its main steps were introduced in detail. The detonation products and parameters of PETN, CL-20 and aluminized explosives were predicted with this method, and after parameter adjustment, it is found that the predicted results and the experimental ones are in better agreement. In comparison with the detonation experiment data of single compound, it is found out that when the BKW equation parameters are adjusted, the energetic materials with more similar percentage of gas fraction in detonation mass to the explosives predicted should be used as the training set of the LS-SVM model. If the detonation parameters of aluminized explosives are predicted, it should use the Al/O ratio close to measured explosive to train the SVM model.

explosion mechanics; BKW equation of state; basic feasible solution; support vector machine;SVM; Gibbs free energy

10.14077/j.issn.1007-7812.2017.03.009

2016-09-01；

2017-01-04

徐州市科技計劃社會發展項目(KC15SH064);徐州工業職業技術學院科技基金資助項目(XGY201607)

何偉平(1983-)，男，碩士，講師，從事含能材料理論計算研究。E-mail: 252927740@qq.com

TJ55;O381

A

1007-7812(2017)03-0053-07