基于多維空間的燃燒平衡產物組成的計算方法

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

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

?

基于多維空間的燃燒平衡產物組成的計算方法

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

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

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

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

引　言

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

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

1　求解方法

1.1數學模型的建立

燃燒平衡產物的求解模型滿足質量守恒、反應前后總焓相等以及產物吉布斯自由能最小等條件。其數學模型可表示如下：

(1)

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

1.2基本可行解的篩選

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

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

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

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

1.3搜索方法

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

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

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

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

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

2　計算驗證

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

采用本方法計算結果與文獻[13]中采用遺傳算法的計算結果進行對比，結果如表2所示，其中燃燒產物壓強為30MPa。

表1　基本可行解篩選結果

表2　燃燒產物計算結果

從表2可以看出，采用基于多維空間的燃燒平衡產物的計算方法，得到的結果與文獻[11]采用遺傳算法計算的結果基本一致，差異稍大的為PETN，其次為RDX，但燃燒溫度誤差均不超過3%。分析文獻[11]的產物數據可以發現，由于PETN及RDX存在較多的次要產物可能未被考慮(即表中標記為“-”的位置，其總數量分別為0.54mol、0.32mol)，故引起結果與本研究不一致。文獻[13]采用平衡常數法計算PETN的燃燒產物組成(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的燃燒產物組成(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。可見采用本方法的計算結果大體上介于文獻[11]和文獻[13]結果之間，因此具有較好的可信性。

搜索各平衡產物過程中，r和T的總更新次數分別為：TNT(377，14)、Tetryl(465，14)、TNBA(874，14)、RDX(2114，15)、PETN(278，15)、Nitromethane(454，14)，括號中的第1個數值為r的總更新次數，第2個數值為T的總更新次數。除RDX外，其他含能材料燃燒平衡產物的搜索，r的更新次數在200～900之間。在Intel core i5處理器下，程序計算TNT、Tetryl、TNBA、RDX、PETN、Nitromethane所需要的時間分別為469、1078、1500、3282、750、1343ms。可見利用本計算方法，程序在經過搜索步數不是很大的情形下，就能快速找出符合精度要求的最優解。

3　結　論

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

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

[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]. 火炸藥學報,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]潘功配，楊碩. 煙火學[M].北京:北京理工大學出版社,1997:42-52.

[5]田德余. 氣態物系平衡組成的計算[J]. 計算機與應用化學, 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]宋東明，潘功配，王乃巖. 基于最小自由能法的煙火藥燃燒產物預測模型[J]. 彈箭與制導學報,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]王雪，潘功配，李毅. 最小自由能法煙火劑平衡產物計算[J]. 火炸藥學報,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]許小平，張唯. 懲罰函數法計算燃燒產物的平衡組成[J]. 宇航學報,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]范磊，潘功配，歐陽的華，等. 基于支持向量機的Mg/PTFE煙火藥燃燒特性預測[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].含能材料,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] 謝中元，周霖，王浩，等. 遺傳算法在燃燒產物平衡組分計算中的應用[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] 周霖. 爆炸化學基礎[M].北京:國防工業出版社,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

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

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

TJ55;O643.2+1

A

1007-7812(2016)04-0046-05