遺傳算法在燃燒產物平衡組分計算中的應用

謝中元， 周 霖, 王 浩, 李守殿

(1. 西安近代化學研究所, 陜西 西安 710065； 2. 北京理工大學爆炸科學與技術國家重點實驗室， 北京 100081； 3. 中國兵器工業(yè)規(guī)劃研究院， 北京 100053)

1 引 言

燃燒產物平衡組分的準確計算是火炸藥配方設計與性能預估的重要研究內容。目前，燃燒產物組分含量的計算方法主要有平衡常數(shù)法、最小自由能法等，其中最小自由能法由于不需要知道化學反應的中間過程且不涉及化學反應平衡常數(shù)的選取而應用最為廣泛[1]，如美國國家航空和宇宙航行局采用最小自由能法計算了包含20種元素的燃燒產物組成[2]； 國內焦清介[3]、劉晶如[4]、龐維強[5]等應用最小自由能法計算了煙火藥的燃燒產物組成，得到了與實驗相符的結果。然而，最小自由能法求解過程中由于采用泰勒展開法、Newton法、布羅依登法[1]等方法對非線性方程組進行求解，其初始值設置難度增加，且對于多組分多相態(tài)的復雜的化學反應，其收斂也比較困難。

遺傳算法是Holland J教授[6-7]提出的基于生物遺傳和進化機制的適合于復雜系統(tǒng)優(yōu)化的自適應概率優(yōu)化計算方法，包含變量編碼、復制、交叉、變異等過程，具有全局優(yōu)化、操作簡便以及較強的魯棒性等特點，在單目標尋優(yōu)、多目標尋優(yōu)、參數(shù)優(yōu)化、計劃調度等方面具有廣泛的應用。近年來，遺傳算法逐漸開始應用于化學平衡計算，安維中[8]、潘功配[9]、賈冬梅[10]等應用遺傳算法對含化學反應體系的相平衡進行求解，取得了良好效果。

基于遺傳算法的優(yōu)良特性以及火炸藥的燃燒特性，本研究擬在最小自由能法的基礎上引入遺傳算法建立燃燒產物平衡組分計算方法，并以文獻數(shù)據(jù)與實驗測試結果驗證計算結果的正確性。

2 計算模型

2.1 問題描述

熱力學平衡理論[1]表明，當燃燒產物體系達到動態(tài)平衡時，自由能取得最小值，反之，當體系自由能達到最小值時，燃燒產物體系即為系統(tǒng)平衡體系，這就是最小自由能法求解燃燒產物平衡組分的基本原理。為表述清楚，設氣態(tài)產物組分編號為1,2,…,N，固態(tài)產物的組分為N+1,N+2,…,M，則體系自由能G可表述為[9]：

(1)

(2)

(3)

式(1)中氣態(tài)燃燒產物與固態(tài)燃燒產物同時需滿足以下約束條件：

(4)

式中，aij表示第i種燃燒產物中第j種元素的個數(shù);Nj為火炸藥配方體系中第j種元素的物質的量，mol;l為火炸藥配方體系中元素的種類數(shù)。燃燒產物組分計算即在滿足式(4)條件下，使得式(1)達到最小值時所求得的組分ni。

2.2 計算方法

計算按以下步驟進行： (1)根據(jù)變量數(shù)NVAR、種群大小NIND、變量字節(jié)數(shù)PRECI以及變量變化范圍對基因編碼; (2)計算目標函數(shù)值并將其作為初始值，進入包含適應度、選擇、交叉、變異的循環(huán)計算，當條件滿足時跳出計算; (3)輸出計算結果。如圖1所示。

圖1 遺傳算法計算框圖

Fig.1 Calculation flow chart of genetic algorithm

下面主要闡述遺傳算法中自變量的選取與變化范圍的確定、適應度函數(shù)的確定以及燃燒產物平衡組分計算方法。

2.2.1 自變量的選取與變化范圍的確定

為敘述簡便，以火炸藥中應用廣泛的CHON類化合物為例闡述自變量選取與變化范圍確定的方法。

根據(jù)文獻[1]中的研究結果，將1 kg CHON類化合物分子式表示為CaHbOcNd(下標表示化合物中相應元素的物質的量)，燃燒產物組分(14種)表示為：nC(solid)，nCO2，nCO，nCH4，nH2O，nH2，nH，nHO，nNH3，nN2，nN，nNO，nNO2，nO。鑒于各組分含量均滿足約束條件，將式(4)展開為：

(5)

由于式(5)中包含4個線性化約束方程，因此可將4種組分表示為其它組分的線性關系式，進而使得獨立燃燒產物組分變量由14個減少為10個，最小自由能函數(shù)(1)變?yōu)?0個自變量的函數(shù)，且求解問題由線性化約束問題轉化為線性化非約束問題。

變量的選取與取值范圍確定，原則上可以從燃燒產物中任取10種組分，只要這10種組分涵蓋所有元素即可，但是由于各組分含量差異較大，任意選定的自變量不容易確定其變化范圍。基于此，根據(jù)文獻[1,3,11,12]，首先將上述14種燃燒產物組分分為兩組，第一組為主要組分，包含nC(solid)、nCO2、nCO、nH2O、nH2、nN2，第二組為次要成分，包含nCH4、nH、nHO、nNH3、nN、nNO、nNO2、nO，其中，主要組分含量總和占燃燒產物總含量的95%以上; 然后根據(jù)次要組分含量少、變化范圍小且容易確定等特點，選取包含不同元素的部分次要組分作為函數(shù)自變量，同時補充一定的主要成分作為自變量，本研究選取自變量為:nC(solid)，nCO2，nCH4，nH，nHO，nNH3，nN，nNO，nNO2，nO; 最后通過對文獻[11]中火炸藥燃燒產物組分含量的統(tǒng)計與整理，以較大上限確定各組分含量變化范圍，如表1所示。

表1中a、b、c、d與CHON類化合物分子式下標相同。經(jīng)計算表明，各變量的變化范圍除nC(solid)和nCO2需要根據(jù)不同物質稍作調整外，其他變量的取值范圍不需調節(jié)就可以滿足不同化合物燃燒產物計算需求。顯然，相對于泰勒展開法、Newton法以及布羅依登法[13]等方法對于空間點的初始值設置，本方法對少數(shù)自變量變化范圍的確定要簡便得多。

表1 變量取值范圍

Table 1 The range of input variables

variablesvariationrangevariablesvariationrangenC(solid)(0,0.4a)nCO2(0,0.2a)nHO(0,0.1b)nH(0,0.1b)nCH4(0,0.05a)nO(0,0.05c)nNH3(0,0.1b)nN(0,0.05d)nNO2(0,0.05d)nNO(0,0.05d)

2.2.2 適應度函數(shù)確定

適應度是與目標函數(shù)值(體系自由能G)接近程度相對應的評價個體優(yōu)劣的參數(shù)，越接近于目標函數(shù)值，則個體對應的適應度值越大，個體越容易繼承。在Matlab遺傳算法工具箱中是通過Ranking(ObjV)函數(shù)來計算個體適應度值的，其返回一個包含對應個體適應度的FitnV列向量。這個函數(shù)是從最小化方向對個體適應度進行排序的，即函數(shù)值越小，個體適應度越大，這正好與本研究的問題相吻合，因此可直接將目標函數(shù)(式(1))作為判別個體優(yōu)劣的適應度函數(shù)。然而，Ranking函數(shù)是對正值數(shù)列進行求解的，而自由能為負值，因此，需要在最小自由能值前加上一個較大的正數(shù)使得最小自由能為正值，考慮實際計算情況，本研究增加的正數(shù)為2000。另外，由于各變量在其變化范圍內是隨機抽取的，因此很有可能出現(xiàn)組分含量為負值而不滿足式(4)，為使其在迭代過程中被淘汰，需要設定一個大的數(shù)值使其適應度變小，本研究設定其為10000。適應度函數(shù)如式(6)所示，其中，“其它”是指燃燒產物組分含量為負值的情況。

(6)

2.2.3 燃燒產物平衡組分計算方法

圖2 總程序計算框圖

Fig.2 Block diagram of the total calculation program

表2 輸入變量

Table 2 Input variables

variablenamesymbolcompoundformulaCaHbOcNdenthalpyofformationΔrHΘmcombustionpressurepinitialtemperatureT0bytecountofVariablePRECIpopulationsizeNINDmaximumiterationMAXGENvariablenumberNVARGenerationgapGGAPvariationrangeofnC(solid)[0,X*a](0

表3 輸出參數(shù)

Table 3 Output parameters

parameterssymbolisobariccombustiontemperatureTpcombustionheatQvmolenumberoftotalcombustionproductsn0molenumberofcombustionprod-uctsnC(solid),nCO2,nCO,nCH4,nH2O,nH2,nH,nHO,nNH3,nN2,nN,nNO,nNO2,nO

3 實例計算與分析

3.1 CHON類化合物燃燒產物平衡組分計算

為驗證本方法的可行性及準確性，針對貧氧化合物(TNT，Tetryl，TNBA)與富氧化合物(RDX，PETN，Nitromethane)，采用本方法計算結果與文獻[11]中零近似法、平衡常數(shù)法計算結果進行對比，結果如表4、表5，其中，燃燒產物壓力均為30 MPa，變量數(shù)NVAR為10、種群大小NIND為200、變量字節(jié)數(shù)PRECI為20、最大迭代次數(shù)MAXGEN為3000、代溝GGAP為0.7、T0設置為3000 K、nC(solid)變化范圍設置為[0, 0.4a]、nCO2變化范圍設置為[0, 0.4a]。

從表4可以看出，針對貧氧化合物，零近似法計算、平衡常數(shù)法、遺傳算法三者計算結果非常接近，但仍然存在一定的差別。其中，相對于TNT，nH2O與nCH4的平衡常數(shù)法計算結果分別為0.18,0.25 mol，而零近似法計算結果均為零，這是因為零近似計算法由于其自身假設與平衡常數(shù)法在微量元素計算方面存在較大的誤差。遺傳算法與平衡常數(shù)法計算結果基本吻合，組分含量誤差在1%以內，溫度誤差在1.5%以內。

表5可以看出，針對富氧化合物，零近似法、平衡常數(shù)法、遺傳算法三者計算結果存在一定差異。相對于PETN，零近似法與平衡常數(shù)法計算所得高含量組分nCO之間誤差最高為20%，溫度誤差為12%; 平衡常數(shù)法與遺傳算法計算所得高含量組分nCO2之間誤差最大為7.7%，溫度誤差為4.9%。分析誤差產生的原因在于: 零近似法是基于最大放熱原則而提出的一種計算方法，當混合物體系為正氧平衡時，這種假設存在一定的不合理性從而導致計算結果誤差較大; 而平衡常數(shù)法與遺傳算法是基于化學動態(tài)平衡的兩種不同形式而提出的計算方法，致使兩者計算結果誤差相對較少，另外，由于本方法只涉及到初始狀態(tài)與最終狀態(tài)，不受化學反應過程以及平衡參數(shù)的影響，因此計算結果更具可信性。

表4 貧氧化合物計算結果

Table 4 Combustion products of negative-oxygen compounds obtained by caculation

molecularformula(1kg)ΔrHΘm/kJ·kg-1calculationmethodn0/molnH2/molnH2O/molnN2/molnCO/molnCO2/molnCH4/molnC(solid)/molTp/Kzeroapproximation[11]44.0311.0106.6026.42004.402034TNT(C30.84H22.03O26.43N13.22)-262[11]equilibriumconstant[11]44.9210.360.186.6126.550.090.254.152042geneticalgorithm43.4910.260.186.6026.060.100.284.402068zeroapproximation[11]41.796.412.308.7123.191.18002793Tetryl(C24.39H17.42O27.94N17.42) 111[11]equilibriumconstant[11]41.836.362.298.7123.131.19002784geneticalgorithm41.846.402.268.7123.171.220.000202791zeroapproximation[11]38.894.081.755.8325.092.14002137TNBA(C27.34H11.67O31.13N11.67)-1548[11]equilibriumconstant[11]38.944.071.765.8425.122.14002131geneticalgorithm38.914.101.735.8425.072.160.001502134

表5 富氧化合物計算結果

Table 5 Combustion products of oxygen-enriched compounds obtained by caculation

molecularformula(1kg)ΔrHΘm/kJ·kg-1calculationmethodn0/molnH2/molnH2O/molnN2/molnCO/molnCO2/molnHO/molnNO/molTp/KRDX(C13.51H27.03O27.03N27.03)322[11]zeroapproximation[11]40.523.519.9913.519.993.51003511equilibriumconstant[11]40.913.549.6213.4610.123.370.360.13375geneticalgorithm40.803.519.8313.459.973.5400.183439PETN(C15.82H25.32O37.97N12.66)-1681[11]zeroapproximation[11]34.790.7311.936.335.6010.22003897equilibriumconstant[11]36.061.0411.046.176.958.870.960.323473geneticalgorithm35.540.9211.616.006.219.6100.653646Nitromethane(C16.39H49.18O32.79N16.39)-1502[11]zeroapproximation[11]49.1510.8113.768.1913.762.62002628equilibriumconstant[11]49.2110.8013.918.1913.732.65002615geneticalgorithm49.2010.8713.688.1913.682.7100.0012634

3.2 CHONAlClCs類推進劑熱力學參數(shù)計算

爆熱是表征推進劑熱力學性能的重要參數(shù)，其數(shù)值大小與推進劑生成焓、燃燒產物組分含量等參數(shù)密切相關?；诖耍槍τ诤珻sNO3、鋁粉的磁流體發(fā)電機專用推進劑配方[12]，本研究首先采用氧彈量熱儀測試了不同CsNO3含量推進劑的爆熱值，然后應用本計算方法對上述配方的爆熱值進行計算，結果如表6所示。計算中各配方的生成焓由組分生成焓按含量加和得到，其它計算條件設置與2.1節(jié)相同。

表6可以看出，計算結果與實驗結果呈現(xiàn)相同的變化趨勢，隨著CsNO3含量的增加，推進劑爆熱逐漸降低; 計算結果與實驗結果誤差較少，當CsNO3含量為4%時，兩者誤差最大值為4.27%，表明本方法可用于含多種元素、多種固態(tài)物質的推進劑熱力學參數(shù)以及燃燒產物組分含量的計算。

表6 推進劑配方及其性能

Table 6 Performance of propellants contained different CsNO3content

hydroxyl-terminatedpo-lybutadiene/%Di-n-octylsebacate/%toluenediisocyanate/%content/%AlAPCsNO3density/g·cm-3heat/kJ·kg-1testcalculationerror/%13.810.61864.621.78162076004.73.2613.810.61862.641.80261675903.44.2713.310.61861.161.81460805831.24.0913.810.61858.681.83058305769.51.0313.810.61856.6101.85158185658.62.74

4 結 論

(1)在最小自由能法的基礎上，引入遺傳算法建立了火炸藥燃燒產物組分計算方法，并采用文獻數(shù)據(jù)與含銫鹽(CsNO3)推進劑實驗測試數(shù)據(jù)驗證了計算結果的正確性。

(2)相比于泰勒展開法、牛頓法等傳統(tǒng)燃燒產物求解方法中嚴格的初始值設置要求，本研究建立的方法僅需確定少量組分的變化范圍，計算過程難度降低。

(3)本研究建立的燃燒產物組分計算方法可適用于含C、H、O、N、Al、Cl、Cs等元素的復雜化學反應平衡體系的熱力學參數(shù)以及燃燒產物組分含量的計算，具有精度高、計算簡便等特點。

參考文獻:

[1] 潘功配, 楊碩. 煙火學[M]. 北京:北京理工大學出版社, 1997: 42-52.

[2] White W B, Johnson S M, Dantzig G B. Chemical equilibrium in complex mixtures[J].TheJournalofChemicalPhysics, 1958, 28(5): 751-755.

[3] 崔慶忠，焦清介. 基于最小自由能原理設計黑火藥組成[J]. 含能材料, 2004, 12(4): 214-217.

CUI Qing-zhong, JIAO Qing-jie. Design of the ingredients of black powder based on the least free-energy law[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2004, 12(4): 214-217.

[4] 劉晶如, 呂勇, 辛偉. 含能熱塑性聚氨酯推進劑的能量計算與分析[J]. 固體火箭技術, 2012, 35(3): 376-381.

LIU Jing-ru, Lü Yong, XIN Wei. Calculation and analysis on energy characteristics of energetic thermoplastic polyurethane propellant[J].JournalofSolidRocketTechnology, 2012, 35(3): 376-381.

[5] 龐維強, 張曉宏, 樊學忠. 含Cs鹽的HTPB/AP/Al復合推進劑特性研究[J]. 固體火箭技術, 2011, 34(5): 619-622.

PANG Wei-qiang, ZHANG Xiao-hong, FAN Xue-zhong. Study on properties of HTPB-based composite propellant with cesium salt[J].JournalofSolidRocketTechnology, 2011, 34(5): 619-622.

[6] Holland J H. Adaptation in natural and artificial systems [M]. Ann Arbor: University of Michigan Press, 1975.

[7] Montastruc L, Azzaro P C, Pibouleau L, et al. Use of genetic algorithms and gradient based optimization techniques for calcium phosphate precipitation[J].ChemicalEngineeringandProcessing, 2004, 43(10): 1289-1298.

[8] 安維中，胡仰棟，袁希鋼．多相多組分化學平衡和相平衡計算的遺傳算法[J]．化工學報，2003，54(5): 691-694．

AN Wei-zhong, HU Yang-dong, YUAN Xi-gang. Calculation of multiphase and multicomponent chemical equilibrium using genetic algorithm[J].JournalofChemicalIndustryandEngineering, 2003, 54(5): 691-694．

[9]范磊, 潘功配, 歐陽的華，等. 基于遺傳算法結合支持向量機的Mg/PTFE貧氧推進劑配方優(yōu)化[J]. 推進技術, 2012, 33(4): 620-624.

FAN Lei, PAN Gong-pei, OUYANG De-hua, et al. Application of genetic algorithm-support vector machine in formula optimization of Mg/PTFE fuel rich propellant [J].JournalofPropulsionTechnology, 2012, 33(4): 620-624.

[10] Jia D M, Deng W S. Phase and chemical equilibrium calculations by using of a hybrid method of the genetic algorithm[J].ComputersandAppliedChemistry, 2011, 28(9): 1175-1178.

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

[12] 周霖, 謝中元. 含Cs鹽推進劑燃燒產物導電特性研究[J]. 含能材料, 2009. 17 (1): 99-102.

ZHOU Lin, XIE Zhong-yuan. Research on electrical conductivity of combustion production of composite propellant containing Cs salt[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2009, 17 (1): 99-102.

[13] 趙慕愚, 徐寶琨. 復雜化學平衡計算 [M]. 長春: 吉林大學出版社, 1990: 1-10.