1 引言

随着火药或推进剂组分的复杂化，燃烧产物也日趋复杂，从十几种甚至上百种产物都有，需要计算机进行编程计算，目前计算方法有很多种，如内能法，迭代法、吉布斯自由能法等^[1-3]，但不管那种方法都存在收敛速度慢、或者对初值要求严格、或者收敛精度差、或者得到全局最优解困难等问题^[4]。

其中目前用得最多的吉布斯自由能法^[5-7]，其基本原理就是燃烧产物的吉布斯自由能最低并符合元素原子守恒，且吉布斯自由能最低也是方程反应的一个基本条件。但采用该方法进行计算机编程计算中，对初值要求严格，而且几十种或上百种的产物用手工给定初值是不现实的，更不要说是否接近实际解。因此本文提出采用随机方向法进行计算机赋初值，对最小吉布斯自由能法求燃烧产物的平衡组成进行改进。

2 数学模型 2.1 吉布斯自由能法

根据热力学原理，在高温条件下，火药的燃烧气体可视为理想气体，当体系达到化学平衡时，体系的自由能为最小。因此，在一定的温度和压力条件下，求出既能使体系的自由能最小，又能满足体系质量守恒的一组组分值，即为该条件下的平衡组分，系统总的自由能就等于组成该体系各组分的自由能之和，这就是最小自由能法计算燃烧产物平衡组成的基本原理^[6]。

为了讨论方便，不失一般性，设反应物含有C、H、O、N四种元素，则等效分子式可以写为：C_n1H_n2O_n3N_n4燃烧产物有：H₂，O₂，N₂，CO, CO₂，H₂O，等，其反应方程式为：

$ \begin{array}{l} {{\rm{C}}_{n1}}{{\rm{H}}_{n2}}{{\rm{O}}_{n3}}{{\rm{N}}_{n4}} = {x_1}{{\rm{H}}_2} + {x_2}{{\rm{O}}_2} + {x_3}{{\rm{N}}_2} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{x_4}{\rm{CO}} + {x_5}{\rm{C}}{{\rm{O}}_2} + {x_6}{{\rm{H}}_2}{\rm{O}} + \cdots \cdots \end{array} $

(1)

设第i种组分的自由能为：

$ {G_i} = {x_i}\left( {{C_i} + \ln {x_i} - \ln \overline x } \right) $

(2)

其中

$ \overline x = \sum\limits_{i = 1}^m {{x_i}} = {x_1} + {x_2} + \cdots + {x_i} + \cdots + {x_m} $

(3)

式中，G_i为所有组分的吉布斯自由能；C_i为第i种气体组分的吉布斯自由能函数；x_i为第i种气体组分物质的量；m为产物个数。

系统总的自由能为：

$ G = \sum\limits_{i = 1}^m {\left[{{x_i}\left( {{C_i} + \ln {x_i}-\ln \overline x } \right)} \right]} $

(4)

2.2 反应过程的物质守恒

由物质守恒，并设有j种元素(这里j为4)，则有：

$ \sum\limits_{i = 1}^m {{a_{ij}}{x_i}} = {n_j}\;\;\left( {j = 1, 2, \cdots, k} \right) $

(5)

其中, a_ij是第i种气体中的j元素的物质的量；n_j是第j种元素原子总的物质的量。

2.3 数学模型

$ {\rm{目标函数:}}\;G{\rm{ = }}\sum\limits_{i = 1}^m {\left[{{x_i}\left( {{C_i} + \ln {x_i}-\ln \overline x } \right)} \right]} $

(6)

$ {\rm{等式约束条件:}}\;\sum\limits_{i = 1}^m {{a_{ij}}{x_i}} = {n_j} $

(7)

$ {\rm{不等式约束条件:}}\;{x_m} \ge 0 $

(8)

这是一个带等式条件的极值问题，求解此问题一般采用成熟的优化方法。如果考虑到每组组份都必须大于零，则上式是一个带有等式约束和不等式约束优化的问题。采用拉格朗日因子求解虽然速度非常快，但存在需要手工计算初值，并且收敛值和初值选取关系较大，有时会收敛于局部最优值，而不是全局最优值，甚至初值选取不合适会迭代不收敛而导致求解失败。而这些恰恰是随机方向法可以解决的。

3 随机方向法求解 3.1 随机方向法算法

随机方向法完全不考虑复杂的极值条件，其思路是运用随机数方法在每一次符合条件的随机方向中选取最小向量作为搜索方向进行寻优。具体如下：

在可行域内选择一个初始点X₀，利用随机数的概率特性，产生若干个可行的随机方向，并从中选择一个能使目标函数值下降最快的随机方向作为可行搜索方向S，从初始点X₀出发，沿方向S按给定的初始步长α取试探点：

$ X = {X_0} + \alpha S $

(9)

检查X点的适用性和可行性，若满足则继续按上面的迭代式在S₁方向搜索，直至到达某迭代点不能满足条件为止；若不满足则再产生另一个随机方向S₂，重复以上过程，直至得到成功点为止。

3.2 初值个数的确定

对于有m种产物生成的体系，由于n种元素守恒，提供了n个等式约束，传统方法是用拉格朗日求解，需要引入n个拉格朗日参数。由极值条件可知：

$ \frac{{\partial G}}{{\partial {x_i}}} = 0\left( {i = 1, 2, 3 \cdots m} \right)\;\;\;\frac{{\partial G}}{{\partial {\lambda _j}}} = 0\left( {i = 1, 2, 3 \cdots n} \right) $

(10)

这样共有m+n个等式。把初值看成一个向量E，由m+n个分量组成，那么由E构成一个解的空间，我们需要的是符合极值条件的一组向量。

但是细致分析，其实守恒式子可以提供n个等式。那么这n个分量可由其它m个算出。如此一来，向量分量便由m+n降为m个，在可行的空间中随机给出m个初值，由m个算出其它n个，当都大于零时则赋值成功。

3.3 可行区间的确定

随机方向法可以产生较好的伪随机数，而且周期足够大，符合统计规律。生成随机数的方法有很多，这里采用的是均布随机数法^[8]。输入设计变量估计的上限和下限值，a_i≤x_i≤b_i(1, 2, ……n)，在区间(0，1)内产生n个伪随机数r_i(i=1, 2, ……n)，计算随机点X的各分量x_i=a_i+r_r(b_i-a_i)(1, 2, ……n)，可以映射到任意区间(a_i, b_i)。

由于可行区间直接影响初值大小的选取，选取不当则有可能使得选取失败。一种方法是手工试探选取区间，直到程序可以运行为止；另一种则是通过简单的分析来确定合适的区间。根据物质守恒，在反应式(1)中，与碳原子有关的组份永远不会超过n₁，那么理论上与碳相关产物的生产区间在(0，n₁)即可。但事实上以理论和经验来看，这样的区间可能过大，生成的随机数比较稀疏，很可能导致计算机赋值失败，实际应用中为了保证程序的顺利进行，可以根据产物个数来调整，如含碳原子的产物有j个，区间可缩为(0，n₁/j)，最后一种元素由其它几种算出，这样就大大缩小了区间。

3.4 计算程序

设第一次符合初值的随机数为m+n组，比较m+n组矢量，找出自由能最小的一组作为优化方向，并以一定优化步长寻优，直至该方向失效，则重新进行第二组随机数计算，并重复上一次计算，反复迭代，直至找到符合精度的一组向量值作为计算结果。本文只是对其随机方向法进行应用，此处不再赘述，具体可以参考优化设计相关书籍^[8-9]。

4 计算验证

针对式(6)~(9)表示的数学模型，运用C语言编写随机方向法改进吉布斯自由能算法计算程序，计算肼和氧的等物质的量混合物燃烧产物在温度3500 K和压力5.168 MPa时的平衡组成，主要产物有：H₂，O₂，N₂，NO，OH，H₂O，H，O，N，NH。由于反应物中只有H、O、N三种元素，则j=3，产物有10种组分，即m=10。

化学方程式可以写为：

$ \begin{array}{l} {{\rm{C}}_{n1}}{{\rm{H}}_{n2}}{{\rm{O}}_{n3}}{{\rm{N}}_{n4}} = {x_1}{{\rm{H}}_2} + {x_2}{{\rm{O}}_2} + {x_3}{{\rm{N}}_2} + {x_4}{\rm{NO}} + \\ \;\;\;\;\;\;{x_5}{\rm{OH}} + {x_6}{{\rm{H}}_2}{\rm{O}} + {x_7}{\rm{H}} + {x_8}{\rm{O}} + {x_9}{\rm{N}} + {x_{10}}{\rm{NH}} \end{array} $

(11)

系统总自由能为：

$ G = \sum\limits_{i = 1}^m {\left[{{x_i}\left( {{C_i} + \ln {x_i}-\ln \overline x } \right)} \right]} $

(12)

根据最小自由能法原理，在给定的温度和压力条件下，通过随机方向法赋初值，计算得到的结果见表 1。

从表 1可以看出，采用随机方向法赋初值，并采用随机方向法来计算得到的结果与文献[1]采用吉布斯自由能法计算结果一致，而且该方法对初值没有限制，使得计算更加便捷。对于中小型的优化计算，随机方向法可以完全替代吉布斯自由能法。

5 结论

采用随机方向法赋初值，对最小吉布斯自由能法求火药燃烧产物的平衡组成进行了改进，通过对初值个数和可行区间的分析，缩小了计算区间。采用改进后的算法对肼和氧燃烧产物的计算结果与文献[1]计算结果一致。

随机方向法不仅可以对给定优化迭代求解方程初值时提供比较好的初值，而且可以扩展到大型计算中，与吉布斯自由能法或者其他优化方法结合应用来提高计算速度。