2. 北京理工大学爆炸科学与技术国家重点实验室，北京 100081;
3. 中国兵器工业规划研究院，北京 100053
2. State Key Laboratory of Explosion Science and Echnology, Beijing Institute of Technology, Beijing 100081, China;
3. Norincogroup Planning and Research Institute, Beijing 100053, China
燃烧产物平衡组分的准确计算是火炸药配方设计与性能预估的重要研究内容。目前，燃烧产物组分含量的计算方法主要有平衡常数法、最小自由能法等，其中最小自由能法由于不需要知道化学反应的中间过程且不涉及化学反应平衡常数的选取而应用最为广泛[1]，如美国国家航空和宇宙航行局采用最小自由能法计算了包含20种元素的燃烧产物组成[2]; 国内焦清介[3]、刘晶如[4]、庞维强[5]等应用最小自由能法计算了烟火药的燃烧产物组成，得到了与实验相符的结果。然而，最小自由能法求解过程中由于采用泰勒展开法、Newton法、布罗依登法[1]等方法对非线性方程组进行求解，其初始值设置难度增加，且对于多组分多相态的复杂的化学反应，其收敛也比较困难。
遗传算法是Holland J教授[6-7]提出的基于生物遗传和进化机制的适合于复杂系统优化的自适应概率优化计算方法，包含变量编码、复制、交叉、变异等过程，具有全局优化、操作简便以及较强的鲁棒性等特点，在单目标寻优、多目标寻优、参数优化、计划调度等方面具有广泛的应用。近年来，遗传算法逐渐开始应用于化学平衡计算，安维中[8]、潘功配[9]、贾冬梅[10]等应用遗传算法对含化学反应体系的相平衡进行求解，取得了良好效果。
基于遗传算法的优良特性以及火炸药的燃烧特性，本研究拟在最小自由能法的基础上引入遗传算法建立燃烧产物平衡组分计算方法，并以文献数据与实验测试结果验证计算结果的正确性。
2 计算模型 2.1 问题描述热力学平衡理论[1]表明，当燃烧产物体系达到动态平衡时，自由能取得最小值，反之，当体系自由能达到最小值时，燃烧产物体系即为系统平衡体系，这就是最小自由能法求解燃烧产物平衡组分的基本原理。为表述清楚，设气态产物组分编号为1, 2, …, N，固态产物的组分为N+1, N+2, …, M，则体系自由能G可表述为[9]：
| $ G = \sum\limits_1^N {{n_i}u_i^{\rm{g}}} + \sum\limits_{N + 1}^M {{n_i}u_i^{\rm{s}}} $ | (1) |
式中，气态产物组分化学势uig(J·mol-1)，表示为：
| $ u_i^{\rm{g}} = \left( {\frac{{{G_{\rm{m}}}}}{{RT}}} \right)_i^{\rm{g}} + RT\left[{\ln \left( {\frac{{n_i^{\rm{g}}}}{n}} \right) + \ln p} \right] $ | (2) |
固态产物组分化学势uis(J·mol-1)表示为：
| $ u_i^{\rm{s}} = \left( {\frac{{{G_{\rm{m}}}}}{{RT}}} \right)_i^{\rm{s}} $ | (3) |
式中，上标g和s分别表示为气态和固态; Gm为燃烧产物标准状态下的自由能，kJ·mol-1; R为气体常数，8.314 J·K-1·mol-1; T为体系温度，K; ni为燃烧产物组分的物质的量，mol; nig为燃烧产物组分中气态组分的物质的量，mol; n为气态产物物质的量的总和，mol; p为燃烧产物压强，MPa。
式(1)中气态燃烧产物与固态燃烧产物同时需满足以下约束条件：
| $ \left\{ {\begin{array}{*{20}{l}} {{n_i} \ge 0\;\;i = 1,2, \cdots ,M}\\ {\sum\limits_{i = 1}^N {{a_{ij}}{n_i}} + \sum\limits_{i = N + 1}^M {{a_{ij}}{n_i}} = {N_j}\;\;\;j = 1,2, \cdots ,l} \end{array}} \right. $ | (4) |
式中，aij表示第i种燃烧产物中第j种元素的个数; Nj为火炸药配方体系中第j种元素的物质的量，mol; l为火炸药配方体系中元素的种类数。燃烧产物组分计算即在满足式(4)条件下，使得式(1)达到最小值时所求得的组分ni。
2.2 计算方法计算按以下步骤进行：(1)根据变量数NVAR、种群大小NIND、变量字节数PRECI以及变量变化范围对基因编码; (2)计算目标函数值并将其作为初始值，进入包含适应度、选择、交叉、变异的循环计算，当条件满足时跳出计算; (3)输出计算结果。如图 1所示。
|
图 1 遗传算法计算框图 Fig.1 Calculation flow chart of genetic algorithm |
下面主要阐述遗传算法中自变量的选取与变化范围的确定、适应度函数的确定以及燃烧产物平衡组分计算方法。
2.2.1 自变量的选取与变化范围的确定为叙述简便，以火炸药中应用广泛的CHON类化合物为例阐述自变量选取与变化范围确定的方法。
根据文献[1]中的研究结果，将1 kg CHON类化合物分子式表示为CaHbOcNd(下标表示化合物中相应元素的物质的量)，燃烧产物组分(14种)表示为：nC(solid)，nCO2，nCO，nCH4，nH2O，nH2，nH，nHO，nNH3，nN2，nN，nNO，nNO2，nO。鉴于各组分含量均满足约束条件，将式(4)展开为：
| $ \left\{ \begin{align} &a={{n}_{{{\rm{C}}_{\left( \rm{solid} \right)}}}}+{{n}_{\rm{C}{{\rm{O}}_{2}}}}+{{n}_{\rm{CO}}}+{{n}_{\rm{C}{{\rm{H}}_{4}}}} \\ &b=2{{n}_{{{\rm{H}}_{2}}\rm{O}}}+2{{n}_{{{\rm{H}}_{\rm{2}}}}}+{{n}_{\rm{H}}}+{{n}_{\rm{HO}}}+3{{n}_{\rm{N}{{\rm{H}}_{3}}}}+4{{n}_{\rm{C}{{\rm{H}}_{4}}}} \\ &c={{n}_{{{\rm{H}}_{2}}\rm{O}}}+2{{n}_{\rm{C}{{\rm{O}}_{2}}}}+{{n}_{\rm{CO}}}+{{n}_{\rm{HO}}}+{{n}_{\rm{NO}}}+2{{n}_{\rm{N}{{\rm{O}}_{2}}}}+{{n}_{\rm{O}}} \\ &d={{n}_{\rm{N}{{\rm{H}}_{\rm{3}}}}}+{{n}_{\rm{NO}}}+{{n}_{\rm{N}{{\rm{O}}_{2}}}}+{{n}_{\rm{N}}}+2{{n}_{{{\rm{N}}_{2}}}} \\ \end{align} \right. $ | (5) |
由于式(5)中包含4个线性化约束方程，因此可将4种组分表示为其它组分的线性关系式，进而使得独立燃烧产物组分变量由14个减少为10个，最小自由能函数(1)变为10个自变量的函数，且求解问题由线性化约束问题转化为线性化非约束问题。
变量的选取与取值范围确定，原则上可以从燃烧产物中任取10种组分，只要这10种组分涵盖所有元素即可，但是由于各组分含量差异较大，任意选定的自变量不容易确定其变化范围。基于此，根据文献[1, 3, 11, 12]，首先将上述14种燃烧产物组分分为两组，第一组为主要组分，包含nC(solid)、nCO2、nCO、nH2O、nH2、nN2，第二组为次要成分，包含nCH4、nH、nHO、nNH3、nN、nNO、nNO2、nO，其中，主要组分含量总和占燃烧产物总含量的95%以上; 然后根据次要组分含量少、变化范围小且容易确定等特点，选取包含不同元素的部分次要组分作为函数自变量，同时补充一定的主要成分作为自变量，本研究选取自变量为: nC(solid)，nCO2，nCH4，nH，nHO，nNH3，nN，nNO，nNO2，nO; 最后通过对文献[11]中火炸药燃烧产物组分含量的统计与整理，以较大上限确定各组分含量变化范围，如表 1所示。
|
表 1 变量取值范围 Tab.1 The range of input variables |
表 1中a、b、c、d与CHON类化合物分子式下标相同。经计算表明，各变量的变化范围除nC(solid)和nCO2需要根据不同物质稍作调整外，其他变量的取值范围不需调节就可以满足不同化合物燃烧产物计算需求。显然，相对于泰勒展开法、Newton法以及布罗依登法[13]等方法对于空间点的初始值设置，本方法对少数自变量变化范围的确定要简便得多。
2.2.2 适应度函数确定适应度是与目标函数值(体系自由能G)接近程度相对应的评价个体优劣的参数，越接近于目标函数值，则个体对应的适应度值越大，个体越容易继承。在Matlab遗传算法工具箱中是通过Ranking(ObjV)函数来计算个体适应度值的，其返回一个包含对应个体适应度的FitnV列向量。这个函数是从最小化方向对个体适应度进行排序的，即函数值越小，个体适应度越大，这正好与本研究的问题相吻合，因此可直接将目标函数(式(1))作为判别个体优劣的适应度函数。然而，Ranking函数是对正值数列进行求解的，而自由能为负值，因此，需要在最小自由能值前加上一个较大的正数使得最小自由能为正值，考虑实际计算情况，本研究增加的正数为2000。另外，由于各变量在其变化范围内是随机抽取的，因此很有可能出现组分含量为负值而不满足式(4)，为使其在迭代过程中被淘汰，需要设定一个大的数值使其适应度变小，本研究设定其为10000。适应度函数如式(6)所示，其中，“其它”是指燃烧产物组分含量为负值的情况。
| $ \left\{ \begin{array}{l} {G_{\rm{s}}} = 2000 + G\;\;\;{n_{\rm{C}}}, {n_{{\rm{C}}{{\rm{O}}_{\rm{2}}}}}, \cdots, {n_{{{\rm{H}}_{\rm{2}}}{\rm{O}}}} \ge 0\\ {G_{\rm{s}}}{\rm{ = 10000}}\;\;\;\;\;\;\;\;\;\;{\rm{其他}} \end{array} \right. $ | (6) |
燃烧产物平衡组分计算方法总体计算框图如图 2所示，其中，ΔrHmΘ表示为生成焓，kJ·mol-1; ΔH为燃烧产物焓变，kJ·mol-1; T0为初始设置温度, K; Hi为燃烧产物在某一温度下的焓，Hi0为燃烧产物在标准状态下的焓, kJ·mol-1。
|
图 2 总程序计算框图 Fig.2 Block diagram of the total calculation program |
具体计算按以下步骤进行: (1)输入初始参数，具体如表 2所示，其中，初始温度T0设置为3000 K即可满足大部分火炸药燃烧温度的计算; 种群大小NIND即数据量可设置为200左右，以减少计算时间并保证计算精度; 最大迭代次数MAXGEN设置为3000左右即可保证计算精度; 变量字节数PRECI设置为20即可满足计算精度要求; 代沟GGAP可固定设置为0.7;变量变化范围除nC(solid)和nCO2需要根据不同物质中碳含量不同稍加修改上限和下限外，气体均按表 1取值。(2)利用燃烧产物焓变ΔH[11]与混合物生成焓ΔrHmΘ的比较作为判断条件，由遗传算法确定温度变化的上下界[T0, T1]。(3)利T0与T1之间差值作为判断条件，采用二分法[11]与遗传算法确定燃烧温度T以及燃烧产物各组分的物质的量ni; (4)输出计算结果并按要求作图，具体如表 3所示。
|
|
表 2 输入变量 Tab.2 Input variables |
|
|
表 3 输出参数 Tab.3 Output parameters |
为验证本方法的可行性及准确性，针对贫氧化合物(TNT，Tetryl，TNBA)与富氧化合物(RDX，PETN，Nitromethane)，采用本方法计算结果与文献[11]中零近似法、平衡常数法计算结果进行对比，结果如表 4、表 5，其中，燃烧产物压力均为30 MPa，变量数NVAR为10、种群大小NIND为200、变量字节数PRECI为20、最大迭代次数MAXGEN为3000、代沟GGAP为0.7、T0设置为3000 K、nC(solid)变化范围设置为[0, 0.4a]、nCO2变化范围设置为[0, 0.4a]。
|
|
表 4 贫氧化合物计算结果 Tab.4 Combustion products of negative-oxygen compounds obtained by caculation |
|
|
表 5 富氧化合物计算结果 Tab.5 Combustion products of oxygen-enriched compounds obtained by caculation |
从表 4可以看出，针对贫氧化合物，零近似法计算、平衡常数法、遗传算法三者计算结果非常接近，但仍然存在一定的差别。其中，相对于TNT，nH2O与nCH4的平衡常数法计算结果分别为0.18, 0.25 mol，而零近似法计算结果均为零，这是因为零近似计算法由于其自身假设与平衡常数法在微量元素计算方面存在较大的误差。遗传算法与平衡常数法计算结果基本吻合，组分含量误差在1%以内，温度误差在1.5%以内。
表 5可以看出，针对富氧化合物，零近似法、平衡常数法、遗传算法三者计算结果存在一定差异。相对于PETN，零近似法与平衡常数法计算所得高含量组分nCO之间误差最高为20%，温度误差为12%;平衡常数法与遗传算法计算所得高含量组分nCO2之间误差最大为7.7%，温度误差为4.9%。分析误差产生的原因在于:零近似法是基于最大放热原则而提出的一种计算方法，当混合物体系为正氧平衡时，这种假设存在一定的不合理性从而导致计算结果误差较大; 而平衡常数法与遗传算法是基于化学动态平衡的两种不同形式而提出的计算方法，致使两者计算结果误差相对较少，另外，由于本方法只涉及到初始状态与最终状态，不受化学反应过程以及平衡参数的影响，因此计算结果更具可信性。
3.2 CHONAlClCs类推进剂热力学参数计算爆热是表征推进剂热力学性能的重要参数，其数值大小与推进剂生成焓、燃烧产物组分含量等参数密切相关。基于此，针对于含CsNO3、铝粉的磁流体发电机专用推进剂配方[12]，本研究首先采用氧弹量热仪测试了不同CsNO3含量推进剂的爆热值，然后应用本计算方法对上述配方的爆热值进行计算，结果如表 6所示。计算中各配方的生成焓由组分生成焓按含量加和得到，其它计算条件设置与2.1节相同。
|
|
表 6 推进剂配方及其性能 Tab.6 Performance of propellants contained different CsNO3 content |
表 6可以看出，计算结果与实验结果呈现相同的变化趋势，随着CsNO3含量的增加，推进剂爆热逐渐降低; 计算结果与实验结果误差较少，当CsNO3含量为4%时，两者误差最大值为4.27%，表明本方法可用于含多种元素、多种固态物质的推进剂热力学参数以及燃烧产物组分含量的计算。
4 结论(1) 在最小自由能法的基础上，引入遗传算法建立了火炸药燃烧产物组分计算方法，并采用文献数据与含铯盐(CsNO3)推进剂实验测试数据验证了计算结果的正确性。
(2) 相比于泰勒展开法、牛顿法等传统燃烧产物求解方法中严格的初始值设置要求，本研究建立的方法仅需确定少量组分的变化范围，计算过程难度降低。
(3) 本研究建立的燃烧产物组分计算方法可适用于含C、H、O、N、Al、Cl、Cs等元素的复杂化学反应平衡体系的热力学参数以及燃烧产物组分含量的计算，具有精度高、计算简便等特点。
| [1] |
潘功配, 杨硕.
烟火学[M]. 北京: 北京理工大学出版社, 1997: 42-52.
|
| [2] |
White W B, Johnson S M, Dantzig G B. Chemical equilibrium in complex mixtures[J].
The Journal of Chemical Physics, 1958, 28(5): 751-755. DOI:10.1063/1.1744264 |
| [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]. Chinese Journal of Energetic Materials(Hanneng Cailiao), 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]. Journal of Solid Rocket Technology, 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]. Journal of Solid Rocket Technology, 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].
Chemical Engineering and Processing, 2004, 43(10): 1289-1298. DOI:10.1016/j.cep.2003.12.002 |
| [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]. Journal of Chemical Industry and Engineering, 2003, 54(5): 691-694. |
| [9] |
范磊, 潘功配, 欧阳的华, 等. 基于遗传算法结合支持向量机的Mg/PTFE贫氧推进剂配方优化[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]. Journal of Propulsion Technology, 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].
Computers and Applied Chemistry, 2011, 28(9): 1175-1178. |
| [11] |
周霖.
爆炸化学基础[M]. 北京: 国防工业出版, 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]. Chinese Journal of Energetic Materials(Hanneng Cailiao), 2009, 17(1): 99-102. |
| [13] |
赵慕愚, 徐宝琨.
复杂化学平衡计算[M]. 长春: 吉林大学出版社, 1990: 1-10.
|
Genetic algorithm for the calculation of combustion equilibrium composition was proposed based on the principles of chemical equilibrium and the method of minimum free energy. By determining the independent variables and its range, the fitness function was confirmed.