何偉平，黃　菊，劉曉靜，陳厚和

(1.徐州工業(yè)職業(yè)技術(shù)學(xué)院化學(xué)工程技術(shù)學(xué)院，江蘇 徐州 221140; 2.徐州工程學(xué)院化學(xué)化工學(xué)院，江蘇 徐州 221111;3.南京理工大學(xué)化工學(xué)院，江蘇 南京 210094)

?

基于多維空間的燃燒平衡產(chǎn)物組成的計算方法

何偉平1，黃菊2，劉曉靜1，陳厚和3

(1.徐州工業(yè)職業(yè)技術(shù)學(xué)院化學(xué)工程技術(shù)學(xué)院，江蘇 徐州 221140; 2.徐州工程學(xué)院化學(xué)化工學(xué)院，江蘇 徐州 221111;3.南京理工大學(xué)化工學(xué)院，江蘇 南京 210094)

為計算燃燒平衡產(chǎn)物的組成，在質(zhì)量守恒關(guān)系的基礎(chǔ)上，得到一系列基本可行解。根據(jù)吉布斯自由能，對基本可行解進(jìn)行優(yōu)化篩選，將篩選得到的基本可行解視為多維解空間的頂點，以吉布斯自由能降低為判據(jù)，利用基本可行解逐個與初始解線性組合，從而搜索得到最佳點，即為最終平衡產(chǎn)物組成。應(yīng)用文獻(xiàn)數(shù)據(jù)對計算結(jié)果的準(zhǔn)確性進(jìn)行驗證。結(jié)果表明，該方法計算平衡產(chǎn)物的組成與文獻(xiàn)數(shù)據(jù)基本一致，燃燒溫度最大誤差不超過3%。

燃燒；平衡產(chǎn)物；多維空間；搜索算法；最小自由能；吉布斯自由能

引　言

含能材料燃燒熱力學(xué)計算的關(guān)鍵在于平衡產(chǎn)物的確定，確定平衡產(chǎn)物的基本方法有平衡常數(shù)法和最小自由能法。目前應(yīng)用最多的是最小自由能法[1-3]，其基本原理是：在給定的溫度和壓強下，體系的吉布斯自由能達(dá)到最小值時，混合體系的組成即為平衡組成。應(yīng)用該原理求算平衡產(chǎn)物，其實質(zhì)是在滿足質(zhì)量守恒的條件下，解決自由能最小的極值問題，一般通過引入拉格朗日因子來處理該問題[4]，傳統(tǒng)的計算機(jī)求解方式有泰勒展開法、Newton法、布羅依登法[5-7]等，但這些求解方式存在初始解設(shè)置難度大，計算復(fù)雜或收斂困難的問題。

針對傳統(tǒng)求解方式存在的弊端，出現(xiàn)了一些新的解決辦法，如借助懲罰函數(shù)法[8]、支持向量機(jī)方法[9]、隨機(jī)方向法[10]、遺傳算法[11]等進(jìn)行不同程度的改進(jìn)，這些方法較好地解決了傳統(tǒng)求解方式的弊端。但也存在一些不足，如懲罰函數(shù)法的罰因子非常重要，然而較難設(shè)置合適的罰因子，并且罰因子每次迭代后，均需要重新求出最優(yōu)解[8]；支持向量機(jī)方法則由于涉及矩陣運算，當(dāng)訓(xùn)練樣本較多時會顯著增加計算復(fù)雜度，并且預(yù)測精度也嚴(yán)重依賴于訓(xùn)練樣本；隨機(jī)方向法與遺傳算法雖然具有良好的全局搜索能力，但其局部搜索能力較差，尤其在后期搜索效率較低，導(dǎo)致計算非常費時。本研究提出了一種基于多維空間的燃燒平衡產(chǎn)物的計算方法。其基本思路為：燃燒產(chǎn)物的所有可能解構(gòu)成多維解空間(為凸空間，其維度即為產(chǎn)物種類數(shù))，平衡組分的最優(yōu)解必然位于該解空間內(nèi)的某點，而質(zhì)量守恒方程組的基本可行解則可視為該解空間的頂點。由初始解出發(fā)，以吉布斯自由能下降為依據(jù)，向不同頂點所在的方位進(jìn)行搜索，則將收斂于吉布斯自由能最低點，即平衡組分的最優(yōu)解，為降低計算復(fù)雜程度、提高計算效率并保證計算精度等提供一種新方法。

1　求解方法

1.1數(shù)學(xué)模型的建立

燃燒平衡產(chǎn)物的求解模型滿足質(zhì)量守恒、反應(yīng)前后總焓相等以及產(chǎn)物吉布斯自由能最小等條件。其數(shù)學(xué)模型可表示如下：

(1)

式中：T為燃燒溫度；G為體系的吉布斯自由能;Hi為第i種產(chǎn)物的摩爾焓;Si為第i種產(chǎn)物的摩爾熵;ng為氣態(tài)產(chǎn)物種數(shù);ns為凝聚態(tài)產(chǎn)物種數(shù);aij為第i種產(chǎn)物分子式中的第j種元素原子的個數(shù);xi為第i種產(chǎn)物的摩爾數(shù);bj為第j種元素的總摩爾數(shù);m為體系的元素種數(shù);H0為含能材料的生成焓。

1.2基本可行解的篩選

平衡產(chǎn)物的任意一組解Xi=(x1,x2,…,xng,…,xng+ns)T均可表示為(ng+ns)維空間中的一點，并應(yīng)滿足xi≥0 (i=1,2,…,ng,…,ng+ns)。

實際上滿足式(2)的基本可行解的組合情形往往很多，故應(yīng)篩選出對后面搜索最為有利的組合情形。因為吉布斯自由能是連續(xù)函數(shù)，故最有利情形是指：式(1)的最優(yōu)解與式(2)所表示的向量的吉布斯自由能最為接近，這樣初始解經(jīng)過若干次搜索并修正后最容易收斂到最優(yōu)解上?？紤]到通常在假定溫度(如3000K)并在定壓條件(如30MPa)下求解平衡產(chǎn)物，故采取以下方式進(jìn)行篩選：

圖1　基本可行解的篩選流程Fig.1　Screening process of basic feasible solutions

計算表明，由于基本可行解之間的差異較大，在溫度、壓強變化不是特別顯著的情形下，篩選出的X″1～X″ng+ns不會發(fā)生明顯變化。故定壓條件下，分別計算溫度區(qū)間[T1，T2]兩個端點T1、T2的X″1～X″ng+ns，若無明顯差異(例如差異的個數(shù)少于n)，則認(rèn)為以其中一組X″1～X″ng+ns代替區(qū)間上各點的篩選結(jié)果，不會影響燃燒產(chǎn)物的搜索結(jié)果。

1.3搜索方法

(1)給定溫度和壓強，計算XT1=(1-r)X0+rX″1，若吉布斯自由能G(XT1)

(2)利用余下的X″2～X″ng+ns，依次對X0重復(fù)以上操作;

(3)重復(fù)(1)、(2)操作，直至X0不再被更新，轉(zhuǎn)向(4);

(4)令r=0.5r，繼續(xù)(1)～(3)的搜索操作，直至r

圖2　平衡組分搜索過程示意圖Fig.2　Schematic diagram of search process for equilibrium compositions

2　計算驗證

根據(jù)以上求解方法，利用Visual Basic 6.0語言編寫燃燒平衡產(chǎn)物搜索的程序代碼。為了驗證本方法的可行性和準(zhǔn)確性，選取文獻(xiàn)[11]討論的14種燃燒產(chǎn)物進(jìn)行研究，結(jié)合不同含能材料的元素組成，進(jìn)行基本可行解的篩選，并在此基礎(chǔ)上進(jìn)行計算。以PETN為例，在T=3000K，p=30MPa時篩選的基本可行解如表1所示(基本可行解總數(shù)N=163)，產(chǎn)物統(tǒng)計以1kg含能材料參與計算為基準(zhǔn)。

采用本方法計算結(jié)果與文獻(xiàn)[13]中采用遺傳算法的計算結(jié)果進(jìn)行對比，結(jié)果如表2所示，其中燃燒產(chǎn)物壓強為30MPa。

表1　基本可行解篩選結(jié)果

表2　燃燒產(chǎn)物計算結(jié)果

從表2可以看出，采用基于多維空間的燃燒平衡產(chǎn)物的計算方法，得到的結(jié)果與文獻(xiàn)[11]采用遺傳算法計算的結(jié)果基本一致，差異稍大的為PETN，其次為RDX，但燃燒溫度誤差均不超過3%。分析文獻(xiàn)[11]的產(chǎn)物數(shù)據(jù)可以發(fā)現(xiàn)，由于PETN及RDX存在較多的次要產(chǎn)物可能未被考慮(即表中標(biāo)記為“-”的位置，其總數(shù)量分別為0.54mol、0.32mol)，故引起結(jié)果與本研究不一致。文獻(xiàn)[13]采用平衡常數(shù)法計算PETN的燃燒產(chǎn)物組成(mol)：nO=36.06，nCO2=8.87，nCO=6.95，nH2O=11.04，nH2=1.04，nHO=0.96，nN2=6.17，nNO=0.32，燃燒溫度Tr=3473K；RDX的燃燒產(chǎn)物組成(mol)：n0=40.91，nCO2=3.37，nCO=10.12，nH2O=9.62，nH2=3.54，nHO=0.36，nN2=13.46，nNO=0.1，燃燒溫度Tr=3375K?？梢姴捎帽痉椒ǖ挠嬎憬Y(jié)果大體上介于文獻(xiàn)[11]和文獻(xiàn)[13]結(jié)果之間，因此具有較好的可信性。

搜索各平衡產(chǎn)物過程中，r和T的總更新次數(shù)分別為：TNT(377，14)、Tetryl(465，14)、TNBA(874，14)、RDX(2114，15)、PETN(278，15)、Nitromethane(454，14)，括號中的第1個數(shù)值為r的總更新次數(shù)，第2個數(shù)值為T的總更新次數(shù)。除RDX外，其他含能材料燃燒平衡產(chǎn)物的搜索，r的更新次數(shù)在200～900之間。在Intel core i5處理器下，程序計算TNT、Tetryl、TNBA、RDX、PETN、Nitromethane所需要的時間分別為469、1078、1500、3282、750、1343ms?？梢娎帽居嬎惴椒?，程序在經(jīng)過搜索步數(shù)不是很大的情形下，就能快速找出符合精度要求的最優(yōu)解。

3　結(jié)　論

(2)在給定溫度和壓強下，對基本可行解進(jìn)行篩選，可有效提高程序搜索效率和預(yù)測精度?？梢灶A(yù)見，該處理方式同樣適用于采取其他方法求解燃燒平衡產(chǎn)物的過程。

(3)不同于傳統(tǒng)方法將質(zhì)量守恒約束方程與目標(biāo)函數(shù)通過拉格朗日因子組合到一起，本方法將最優(yōu)解表示為約束方程可行解的線性組合，簡化了求解過程，同時保證了產(chǎn)物的計算精度；并通過優(yōu)化求解過程，保證了計算效率。

[1]Gordon S，McBride B J. Computer program for calculation of chemical equilibrium compositions and applications：I Analysis，NASA RP-1311[R].Washington D C：NASA，1994.

[2]McBride B J，Gordon S. Computer program for calculation of chemical equilibrium compositions and applications：Ⅱ Users manual and program description，NASA RP-1311[R]. Washington D C：NASA，1996.

[3]李猛，趙鳳起，徐司雨，等. 三種能量計算程序在推進(jìn)劑配方設(shè)計中的比較[J]. 火炸藥學(xué)報,2013,36(3): 73-77.

LI Meng，ZHAO Feng-qi，XU Si-yu，et al. Comparison of three kinds of energy calculation programs in formulation design of solid propellants[J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao)，2013,36(3): 73-77.

[4]潘功配，楊碩. 煙火學(xué)[M].北京:北京理工大學(xué)出版社,1997:42-52.

[5]田德余. 氣態(tài)物系平衡組成的計算[J]. 計算機(jī)與應(yīng)用化學(xué), 1986,3(4):310-313.

Tian De-yu. The computation of chemical equilibrium composition in gaseous systems[J]. Computers and Applied Chemistry, 1986,3(4):310-313.

[6]宋東明，潘功配，王乃巖. 基于最小自由能法的煙火藥燃燒產(chǎn)物預(yù)測模型[J]. 彈箭與制導(dǎo)學(xué)報,2006,26(1):120-122.

SONG Dong-ming，PAN Gong-pei，WANG Nai-yan. Calculation model of pyrotechnical combustion products based on Minimum free-energy[J]. Journal of Projectiles，Rockets，Missiles and Guidance，2006,26(1):120-122.

[7]王雪，潘功配，李毅. 最小自由能法煙火劑平衡產(chǎn)物計算[J]. 火炸藥學(xué)報,1994,22(4):61-63.

WANG Xue，PAN Gong-pei，LI Yi. Calculation of pyrotechnic composition equilibrium products based on principle of free-energy minimization[J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao)，1994,22(4):61-63.

[8]許小平，張唯. 懲罰函數(shù)法計算燃燒產(chǎn)物的平衡組成[J]. 宇航學(xué)報,1994(3)：90-95.

XU Xiao-ping，ZHANG Wei. Calculation of equilibrium composition of combustion products by penalty method[J]. Journal of Astronautics，1994(3)：90-95.

[9]范磊，潘功配，歐陽的華，等. 基于支持向量機(jī)的Mg/PTFE煙火藥燃燒特性預(yù)測[J].含能材料,2012,20(4):414-417.

FAN Lei，PAN Gong-pei，OUYANG De-hua，et al. Combustion characteristics forecast of Mg/PTFE pyrotechnic composition with support vector machine[J]. Chinese Journal of Energetic Materials，2012,20(4):414-417.

[10] 菅曉霞，劉慶上，鄒偉偉，等. 隨機(jī)方向法賦初值求火藥燃燒產(chǎn)物的平衡組成[J].含能材料,2012,20(2):202-204.

JIAN Xiao-xia，LIU Qing-shang，ZHOU Wei-wei，et al. Calculation of equilibriun composition of combustion products of gunpowder based on random direction method initialization[J]. Chinese Journal of Energetic Materials，2012,20(2):202-204.

[11] 謝中元，周霖，王浩，等. 遺傳算法在燃燒產(chǎn)物平衡組分計算中的應(yīng)用[J].含能材料,2015,23(4):340-345.

XIE Zhong-yuan，ZHOU Lin，WANG Hao，et al. Application of genetic algorithm in calculation of combustion equilibriunm composition[J]. Chinese Journal of Energetic Materials，2015,23(4):340-345.

[12] Kubota, Naminosuke. Propellants and Explosives: Thermochemical Aspects of Combustion[M].Hoboken:John Wiley & Sons, 2015: 23-29, 64-65.

[13] 周霖. 爆炸化學(xué)基礎(chǔ)[M].北京:國防工業(yè)出版社,2005:33-37.

Calculation Method of Combustion Equilibrium Product Compositions Based on Multidimensional Space

HE Wei-ping1， HUANG Ju2， LIU Xiao-jing1， CHEN Hou-he3

(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, Nanjing 210094, China)

To calculate the combustion equilibrium product compositions, a series of basic feasible solutions were obtained on the basis of conservation of mass. According to Gibbs free energy, the basic feasible solutions were optimized and selected. The basic feasible solutions obtained by selecting were considered as the vertex of the multidimensional space, the decrease of Gibbs free energy was considered as the criterion, the best point was obtained via search using linear combination one by one of basic feasible solution and initial solution, namely the final equilibrium product compositions were obtained. The accuracy of the calculation results was verified by the literature data. Results show that the compositions of equilibrium product calculated by this method are in agreement with the literature data, and the maximum error of combustion temperature is no more than 3%.

combustion；equilibrium product；multidimensional space; search algorithm；minimum free energy;Gibbs free energy

10.14077/j.issn.1007-7812.2016.04.009

2016-03-18;

2016-04-28

徐州市科技計劃社會發(fā)展項目(KC15SH064);徐州工業(yè)職業(yè)技術(shù)學(xué)院科技基金資助項目(XGY201409)

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

TJ55;O643.2+1

A

1007-7812(2016)04-0046-05