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

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

1 引 言

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

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

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

2 計(jì)算模型

2.1 問(wèn)題描述

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

(1)

(2)

(3)

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

(4)

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

2.2 計(jì)算方法

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

圖1 遺傳算法計(jì)算框圖

Fig.1 Calculation flow chart of genetic algorithm

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

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

為敘述簡(jiǎn)便，以火炸藥中應(yīng)用廣泛的CHON類(lèi)化合物為例闡述自變量選取與變化范圍確定的方法。

根據(jù)文獻(xiàn)[1]中的研究結(jié)果，將1 kg CHON類(lèi)化合物分子式表示為CaHbOcNd(下標(biāo)表示化合物中相應(yīng)元素的物質(zhì)的量)，燃燒產(chǎn)物組分(14種)表示為：nC(solid)，nCO2，nCO，nCH4，nH2O，nH2，nH，nHO，nNH3，nN2，nN，nNO，nNO2，nO。鑒于各組分含量均滿(mǎn)足約束條件，將式(4)展開(kāi)為：

(5)

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

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

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

表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 適應(yīng)度函數(shù)確定

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

(6)

2.2.3 燃燒產(chǎn)物平衡組分計(jì)算方法

圖2 總程序計(jì)算框圖

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 實(shí)例計(jì)算與分析

3.1 CHON類(lèi)化合物燃燒產(chǎn)物平衡組分計(jì)算

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

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

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

表4 貧氧化合物計(jì)算結(jié)果

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 富氧化合物計(jì)算結(jié)果

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類(lèi)推進(jìn)劑熱力學(xué)參數(shù)計(jì)算

爆熱是表征推進(jìn)劑熱力學(xué)性能的重要參數(shù)，其數(shù)值大小與推進(jìn)劑生成焓、燃燒產(chǎn)物組分含量等參數(shù)密切相關(guān)。基于此，針對(duì)于含CsNO3、鋁粉的磁流體發(fā)電機(jī)專(zhuān)用推進(jìn)劑配方[12]，本研究首先采用氧彈量熱儀測(cè)試了不同CsNO3含量推進(jìn)劑的爆熱值，然后應(yīng)用本計(jì)算方法對(duì)上述配方的爆熱值進(jìn)行計(jì)算，結(jié)果如表6所示。計(jì)算中各配方的生成焓由組分生成焓按含量加和得到，其它計(jì)算條件設(shè)置與2.1節(jié)相同。

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

表6 推進(jìn)劑配方及其性能

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 結(jié) 論

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

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

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

參考文獻(xiàn):

[1] 潘功配, 楊碩. 煙火學(xué)[M]. 北京:北京理工大學(xué)出版社, 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] 崔慶忠，焦清介. 基于最小自由能原理設(shè)計(jì)黑火藥組成[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ìn)劑的能量計(jì)算與分析[J]. 固體火箭技術(shù), 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] 龐維強(qiáng), 張曉宏, 樊學(xué)忠. 含Cs鹽的HTPB/AP/Al復(fù)合推進(jìn)劑特性研究[J]. 固體火箭技術(shù), 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] 安維中，胡仰棟，袁希鋼．多相多組分化學(xué)平衡和相平衡計(jì)算的遺傳算法[J]．化工學(xué)報(bào)，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]范磊, 潘功配, 歐陽(yáng)的華，等. 基于遺傳算法結(jié)合支持向量機(jī)的Mg/PTFE貧氧推進(jìn)劑配方優(yōu)化[J]. 推進(jìn)技術(shù), 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] 周霖. 爆炸化學(xué)基礎(chǔ)[M]. 北京:國(guó)防工業(yè)出版, 2005: 33-37.

[12] 周霖, 謝中元. 含Cs鹽推進(jìn)劑燃燒產(chǎn)物導(dǎo)電特性研究[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] 趙慕愚, 徐寶琨. 復(fù)雜化學(xué)平衡計(jì)算 [M]. 長(zhǎng)春: 吉林大學(xué)出版社, 1990: 1-10.