固体火箭推进剂配方性能最优化的直接搜索法

译自AIAA 82-1190.直接搜索法(又称直接法)是一种有效的最优化方法。它不用求目标函数和约束函数的偏导数。这种方法已广泛应用于各种工程问题中。

摘 要

本文研究的数值计算方法，用以确定固体火箭推进剂系统的最佳配比，其比冲最高。

根据对最优化理论的分析，本文采用直接搜索法结合计算机程序CEC 71(一种热力计算程序)[6],以获得固体推进剂在冻结流或平衡流条件下的最优配比。比较了计算结果和文献报导的数据，可看出直接搜索法用于确定固体推进剂系统的最优配比是有效的。

符 号

A 推进剂组分的总数

a 热力学常数

c 单位质量推进剂所含某组分的莫尔数，mole/kg

f 函数

Gf Gibbs自由能，cal/kgg 重力加速度

g 9.8m/sec²

H 推进剂组分的热焓，cal/mole

Isp 比冲，sec

J 热功当量

L J/calL 原子总数

M 分子量 kg/mole

m气相产物总数

N产物总数

n 莫尔数，mole

P 总压，atm

R 通用气体常数，8.314J/(mole.°K)

r 直接搜索法的步长

s 熵，cal/(mole.°K)

(Sr°) 标准状态下的熵，cal/(mole.°K)

T 温度，°K

X 质量分数

α,β每个元素在某种产物或某种组分中的原子数

λ 拉格朗日乘子

μ 化学势，eal/mole

δ 最小”步长

ρ 步长缩小因子，ρ<1

下标

c 燃烧室

e 出口

i 产物

j 推进剂组分

k 原子

一、引 言

在研制固体火箭推进剂的时候，推进剂化学方面的科学家或工程师们一直在寻求一些新的组分，能产生较高的比冲。这样的研究工作在试验室内既费时间又费钱，因而必须发展一种数学方法来预测一种固体推进剂系统的最佳配方，并使理论性能最高。

评价推进剂性能的第一步需要计算推进剂燃烧产物的平衡组成，在文献中已有许多计算方法1[2]。1958年以前，所有的平衡组成计算都是使用平衡常数公式作为控制方程。例如Brinkley[3]和Huff等[4]。但是，这个方法仅适用于简单的化学反应。1958年White等[5]建议用最小自由能法计算平衡组成，这个方法分析复杂燃料和推进剂特别适用。1971年Gordon等[6]采用最小自由能法计算了平衡流和冻结流条件下的化学平衡组成和理论的火箭发动机性能。编制的计算程序CEC71已广泛地用来预测火箭发动机性能[7-10]。

1978年Swaminaihan等曾采用平衡常数公式和梯度投影法确定冻结流条件下含金属的固体推进剂的最佳配比。但是，在应用这个方法时发现有以下一些不方便之处：

a)对不同推进剂系统的气相产物组成和其它性能参数必须分别作出初始估计，并满足原子守恒方程组。

b)使用梯度投影法时，偏导数的解析计算困难。

c)将这个方法加以发展，用于平衡流条件下的最优化问题时又非常费机时。

因此，本文的目的就在于提出一个有效的数值方法，且能在平衡流和冻结流两种条件下均能容易又方便地确定对于给定推进剂系统组分的最优配比，其比冲最高。根据对最优化理论的分析，12-131采用直接法结合计算机程序CEC71实现最优。所得数值结果与文献报导的数据吻合较好。

二、数学分析

1.目标函数

比冲是表征固体推进剂性能的主要参数之一。在这篇文章里它是最优化对象，其定义为

(a)在喷管膨胀过程，燃气总是化学平衡的(平衡流),有：

\({I}_{sp}={\left[\frac{2J}{{g}^{2}} \left(\sum\limits_{i=1}^{N}{{n}_{ic}{h}_{ic}}-\sum\limits_{i=1}^{N}{{n}_{ie}{h}_{ie}}\right)\right]}^{\frac{1}{2}}\)

(b)对喷管膨胀过程，燃气具有不变化学组成(冻结流),有：

\({I}_{sp}={\left[\frac{2J}{{g}^{2}}\sum\limits_{i=1}^{N}{{n}_{ci} \left( {h}_{ic}-{h}_{ie} \right) }\right]}^{\frac{1}{2}}\)

式中

\({h}_{i}={a}_{1i}{T}+\frac{{a}_{2i}}{2}{T}^{2}+\frac{{a}_{3i}}{3}{T}^{3}+\frac{{a}_{4i}}{4}{T}^{4}+\frac{{a}_{5i}}{5}{T}^{4}+{a}_{6i}\)

其中：\({a}_{1i},{a}_{2i},{a}_{3i},{a}_{4i},{a}_{5i},{a}_{6i}\)是常数。

2.约束条件

上面定义的目标函数受以下条件的约束不等式约束(·在实际应用中，由于力学性能、然烧性能和工艺性能等的要求，对某些组分还可加上一些不等式约束，例如\({a}_{j}≤{n}_{ij}≤{b}_{j}\)):

\({n}_{ic} \ge 0,{n}_{ie} \ge 0,{T}_{c} \ge 0,{T}_{e} \ge 0,\)

等式约束：

a)原子守恒方程

\(\sum\limits_{i=1}^{N}{{n}_{i}{\alpha}_{ki}}-\sum\limits_{j=1}^{A}{{C}_{j}{\beta}_{kj}}=0 \),K=1,2,…,L         (4)

(产物的某原子总数等于推进剂中该原子的总数)。

b) 系统的能量守恒

\(\sum\limits_{i=1}^{N}{{n}_{ic}{h}_{ic}}-\sum\limits_{j=1}^{A}{{C}_{j}{H}_{j}}=0 \)          (5)

(燃烧室内产物的总焓等于推进剂的总焓)。

c)喷管膨胀过程中平衡流条件下的等熵方程

\(\sum\limits_{i=1}^{N}{{n}_{ie}{s}_{ie}}-\sum\limits_{i=1}^{N}{{n}_{ic}{s}_{ic}}=0 \)          (6)

冻结流下则是

\(\sum\limits_{i=1}^{N}{n}_{ic}({{s}_{ie}-{s}_{ic}})=0 \)                     (7)

\({s}_{i}={({s}_{T}^{0})}_{i}-R\ln {\frac{{n}_{i}}{n}}-R \ln{p}\),   i=1,2,…,m

\({s}_{i}={({s}_{T}^{0})}_{i}\),   i=m+1,…,N                   (8)

方程(8)中，在考虑的可能的N个产物中，气相由1到m表示，凝聚相用m+1到N表示，式中

\(n=\sum\limits{i=1}^{N}{{n}_{i}}\)                           (9)

{({s}_{T}^{0})}_{i}是组分i标准状态下的熵。

\({s}_{T}^{0}={a}_{1i}{\ln{T}}+\frac{{a}_{2i}}{T}+\frac{{a}_{3i}}{2}{T}^{2}+\frac{{a}_{4i}}{3}{T}^{3}+\frac{{a}_{5i}}{4}{T}^{4}+{a}_{7i}\)

(d)化学平衡采用Gibbs自由能最小作为化学平衡的准则。Gibbs自由能为

\({G}_{j}=\sum\limits{i}^{N}{{\mu}_{i}{n}_{i}}\)

而这个自由能最小应受质量守恒方程(4)的约束。

我们定义G为，

\(G={G}_{j}+\sum\limits{k=1}^{N}{{{\lambda}_{k}} \left(\sum\limits_{i=1}^{N}{{n}_{i}{\alpha}_{ki}}-\sum\limits_{j=1}^{A}{{C}_{j}{\beta}_{kj}} \right)}=0\)        (11)

式中\({\lambda}_{k}\)是拉格朗日乘子。平衡条件就变成

\({\delta}{G}=\sum\limits_{i=1}^{N}{\left({\mu}_{i}+\sum\limits_{k=1}^{L}{{\lambda}_{k}{\alpha}_{ki}} \right)}\delta {n}_{i}+\sum\limits_{k=1}^{L}{\left( {n}_{i}{\alpha}_{ki}-\sum\limits_{j=1}^{A}{{c}_{j}{\beta}_{kj}}\right)}\delta {\lambda}_{k}=0\)               (12)

因为变分\(\delta {n}_{i}\);和\(\delta {\lambda}_{k}\)是独立的，所以方程(12)除得到原子守恒方程(4)外，还给出

\({\mu}_{i}+\sum\limits_{k=1}^{L}{{\lambda}_{k}{\alpha}_{ki}}=0\)         ,i=1,2,……,n

Gordon等编制了一个计算机程序[6]CEC71,用牛顿——拉福松叠代法解上述方程(1)

—(13),计算产物的平衡组成和比冲。这个程序是很成功的，我们在这里也予采用。

3.最优化方法

本文企图获得给定推进剂系统组分的最佳配比，使比冲最大。这一最优化问题可以再用数学公式描述如下：

把方程(5)代入(1)或(2),目标函数可写为

\({I}_{sp}=f({X}_{j})\),           j=1,2……,4                   (14)

这里组分j的质量分数是\({X}_{j}={C}_{j}{M}_{j}\),受如下的约束:

\(\sum\limits_{j=1}^{A}{{X}_{j}}=1\)            (15)

方程(14)和(15)可用一个三角形配比图(如图1)表示。图1是过氯酸氨/铝粉/聚乙烯推进剂的理论比冲分布图。由图可见，等比冲曲线是连续的，且彼此不相交。根据最优化理论，实现这一最优化的数值方法可分为两类：

——非直接法(或梯度法):最优解主要依据所含的函数的可微特性确定的。 Swamina-than等人曾采用梯度投影法来解决前面所述的推进剂性能最优化问题。

——等比冲曲线(压力比68/1)

—–直接搜索法的路径

图1过氯酸铵/铝粉/聚乙烯推进剂的性能分布和直接搜索的可行方向

—–直接法：其解需要直接用到性能计算和给定问题的约束方程，并采用系统的循环方法得到一个最优解。Swaminathan等[17]采用网格搜索法[13],这是一种最简单的直接搜索法，通过最小自由能法计算推进剂燃烧产物的平衡组成。但是，这个方法对具有许多独立变量的目标函数的最优化不太有效。另一方面，模式搜索法(Pattern search)是一种最广泛使用的直接法，通过一系列有序的外插和模式移动进行最优化。Eason等[18]比较了许多最优化方法，得到结论——模式搜索法总的来说比梯度法好。这使得模式搜索用在推进剂性能最优化问题上具有吸引力。

为了不违背约束方程(15),本问题要对模式搜索法做些修改。其可行方向取图1所示方向。计算框图见图2,其方法主要步骤可描述如下：

图2推进剂性能最优化所用直接搜索法的方框图

1.选择步长r、及步长缩小因子(ρ<1)和检验是否收敛使用的“最小”步长δ。

2.选择推进剂系统的初始配方(X₁^*,X₂^*,……,X_A^*)并使用计算程序CEC71计算其比冲I_sp^*。

3. 分别设定叠代计数器i=1和j=1。

4.计算X_i=X_i^*+r和组分变动后新的百分组成

\({X}_{k}={X}_{k}^{*}\left( \frac{1-{X}_{j}}{1-{X}_{j}^{*}} \right),k=1,2,…,A\)

\(k \ne i\)

5.使用计算程序CEC71计算这个试验配方的比冲I_sp。

6.如果\({I}_{sp} \ge {I}_{sp}^{*}\),设\({X}_{k}^{*}=Xk{}_{}\),k=1,2,……,A,

j=j+1,L, \({I}_{sp}^{*} = {I}_{sp}\),转第4步。

如果 \({I}_{sp} < {I}_{sp}^{*}\),

(1)如果j=1,设r=-r,j=j+1,转第4步。

(2)如果i>1转第7步。

7.如果i<A,设i=i+1,转第4步。 如果i>A转第8步。

8.如果r>δ,计算r=ρr(缩小步长)并设i=1,j=1,转第4步。

如果r<δ,停止叠代得到最优配比(X₁^*,X₂^*,……,X_A^*)。

三、计算结果及其比较(略)

四、结 论

本研究结合直接搜索法和CEC71计算程序，求解平衡流或冻结流条件下固体推进剂的配比优化问题。这个方法具有简单、精度较高和计算时间较少等优点。

主要参考文献

Gordon, s, and MoBride, B. J.”Computer program for Caloulation ofComplex Chemical Eguilibrium Compositions, Rooket Performance, Incidentand Reflected Shocks and Chapman-Jouguet Detonations” NASA Rep. SP-273,1971.