基于机器学习建模的液体火箭发动机喷管内型面优化设计

李晨沛，周晨初，高玉闪，胡海峰

(西安航天动力研究所，陕西 西安 710100)

0 引言

喷管是液体火箭发动机产生推力的重要部件。推进剂燃烧所产生的高温高压燃气在喷管中迅速膨胀，经由喷管出口排出，产生发动机所需推力。喷管型面的结构将直接影响燃烧所产生的燃气在喷管中的流动情况，进而对发动机的性能产生影响[1-2]。根据液体火箭发动机喷管的结构特点，可将喷管分为收缩段和扩张段两个部分，气流在喷管扩张段处于超声速流动，所以一般情况下，扩张段对喷管的性能和尺寸影响最大，如果扩张段型面设计得不合理，喷管效率就会受到明显的影响。

通常设计喷管型面是基于一维分析理论进行型面优化，并利用二维或三维CFD流场分析的方法来进行修正。马国宝等[3]利用一维等熵流的气动力公式、二次抛物线型面的几何关系等，优化得到了比冲最优的液体火箭发动机型面结构。方杰等[4]根据发动机计算公式，利用单级优化算法实现了某型发动机喷管的多学科设计优化。上述方法可有效实现发动机内型面设计，但一维方法不能实现对喷管性能的准确评估。

雷静等[5]指出多学科设计优化技术可较好地平衡计算精度与成本，在设计域内，采用高精度分析并行计算大量样本数据，利用这些数据构建代理模型，通过代理模型对喷管性能进行优化。代理模型是一种用于建立输入和输出之间关系的模型，可基于各种数据和算法，如机器学习技术，构造一个计算规模小、求解时间短，且计算精度与复杂物理场仿真模型一致的数学模型。利用这一技术众多学者考虑将CFD仿真与喷管优化设计相结合，基于CFD仿真平台获得大量仿真样本点数据，基于机器学习的方式从数据中发现规律，拟合出反映喷管规律的代理模型，并寻优得到优化后的型面结构[5-9]。为进一步提高样本点计算效率，周正[10]、苗萌等[11]、高太元等[12]在优化吸气式冲压发动机后体尾喷管结构时提出自动生成计算网格，极大程度地提高了计算效率。

目前针对液体火箭发动机喷管内型面优化设计的研究较少，本文采用B样条曲线对抛物面型线进行参数化，基于一体化优化设计平台，批量计算样本点CFD流场数据，通过机器学习，构建喷管内型面结构参数与发动机比冲之间的代理模型，以比冲为优化变量对喷管性能进行评估。利用网格变形技术提高计算效率，基于流程算法在设计平台上自动寻优，获得最佳的喷管内型面结构。

1 喷管内型面优化设计流程

本文基于代理模型代替CFD计算，采用全局优化算法完成喷管内型面的优化设计，具体流程如图 1所示。该流程图主要包括两个部分，分别是初始样本的生成及优化迭代更新。

如图1所示流程图的第一部分中，首先针对喷管内型面优化问题选取设计变量，并确定各个设置变量的取值范围。在已知室压、喉部直径、燃气物性、喷管长度等条件的前提下，对喷管的比冲进行优化，选取喷管内型面几个控制点作为设计参数，由三次B样条曲线生成喷管的内型面结构。设计变量的范围可以设置为初始设计点±10%～20%。其次，在完成设计变量的选取后，需要进行试验设计。试验设计需要确定样本点的个数及样本点在取样空间内的分布，本文采用拉丁超立方方法，在设计空间内生成了200个样本点数据。接着，在完成试验设计后，针对每个样本点的几何相关参数建立喷管的几何模型和计算网格。在生成新的喷管几何外形后，建立计算网格，即初始网格，如图 2所示。当样本点数目较多时，通常采用动网格技术生成各个样本点的网格，由此避免手动生成网格时导致的网格质量和数目的不可控。根据初始网格结构，利用径向基函数插值算法，在不改变网格原拓扑关系的基础上，根据变形的物面(面网格)重新生成新的计算网格，如图 3所示。之后，针对各个样本点的几何和网格，将开展CFD计算提取比冲和流场信息，由此构成初始样本库。

图1 基于代理模型的喷管内型面优化流程

图2 初始网格结构

图3 变形后网格

图4 B样条曲线示意图

在完成初始样本库的计算后，将进入图 1的第二部分，即优化计算。首先，基于当前样本库训练代理模型，本文应用Kriging模型、响应面特征模型(Polynomial Response Surface，PRS)和径向基函数特征模型(Radial Basis Function，RBF)训练样本点数据得到代理模型，并基于最大相对偏差、平均偏差等，对代理模型的误差进行初步分析；接着，采用遗传算法(Genetic Algorithm，GA)和粒子群优化算法(Particle Swarm Optimization，PSO)进行优化分析，输出最优解或Pareto解集；之后，采用CFD对最优解进行校核计算。若本次CFD计算得到的目标值与上一次优化的值满足收敛条件，则停止计算并输出最优解。若不满足收敛条件，则把本次CFD计算解添加到样本库中，重新构建代理模型及优化分析。

下面将对液体火箭发动机喷管内型面优化设计流程中几个重要过程进行详细阐述。

1.1 喷管内型面设计

针对喷管扩张段内型面结构，本文采用的参数化方法为B样条曲线方法。计算公式如下：

Pi，k(u)=[Bi-k+1，k(u)Bi-k+2，k(u)…Bi，k(u)]×[pipi+1…pi-k+1]T

(1)

式中，Bj，k(u)表示第j(j=i-k+1，…，i)条k阶B样条基函数有关段；pn表示第n(n=i，…，i-k+1)个顶点；u表示设计变量。

本文给定喷管喉部直径、出口直径以及喷管长度，采用四阶三次B样条曲线表征喷管内型面结构，如图 4所示。此外，本文将第一个控制点设置在喉部附近，并保证曲线与喉部最小部分相切，使喷管喉部光滑过渡。

1.2 网格变形

在喷管的优化设计过程中，需要对不同的设计参数进行多次CFD计算。针对该优化设计问题，最简单的方式是每次CFD计算时重新生成网格，但对于复杂几何问题或网格数量大等情况该方式通常非常耗时。因而，有学者提出使用网格变形方法，其特点是在不改变拓扑关系和节点总数的基础上，移动网格节点的位置，以满足计算域变化[13]。

本文提供径向基函数插值算法，可以根据变形的物面(面网格)重新生成新的计算体网格。径向基函数插值算法对每个坐标方向采用一个插值函数来计算任意网格节点在该方向上的偏移量。

S(X)=p(X)+∑λjΦ||X-Xj||

(2)

其中，

p(X)=c1x+c2y+c3z+c4

(3)

满足插值条件

S(Xj)=fj

(4)

以及正交条件

(5)

求解如下方程可得到系数λj和c1，c2，c3，c4。

(6)

针对数据量较大的情况，径向基函数方法处理难度较大，因此本文选用局部径向基函数插值方法，将全局空间划分为若干子区间，对每个子区间利用径向基函数进行插值运算，进而提高插值效率[14]。该方法网格变形质量高，适用于显著弯曲和扭转的几何非线性大变形情况，可用于结构和非结构网格变形。

1.3 代理模型参数训练及参数优化方法

在指定样本数量情况下，径向基函数模型和响应面模型的训练参数可以采用数值方法进行求解。响应面模型基于最小二乘法计算多项式的拟合系数。径向基函数模型求解一个维度等于样本数目的关联矩阵，采用全选主元高斯消去法求解。Kriging模型训练参数θ体现在关联函数中无解析解，是一个无约束问题，最终转换为求解极大值问题。采用梯度类优化算法和全局优化算法求解参数的最优解，具体包括：遗传算法、粒子群优化算法等。

1.3.1 遗传算法

遗传算法采用群体搜索技术，通过对当前群体实施选择、交叉和变异等一系列遗传操作，从而产生新一代群体，并使群体逐步进化到包含或接近最优解的状态[15]。

在遗传算法中，将n维设计变量用向量X=[x1，x2，…，xn]表示，每个X可以表示一个个体(又被称为染色体)，包含n个遗传基因(设计变量x)。遗传算法求解最优解的过程为模拟生物进化过程，通过群体按照优胜劣汰的进化，将适应度高的基因更多地遗传到下一代，在最终代的群体中达到或接近最优解。

遗传算法的流程如图5所示。

图5 遗传算法流程

(1)初始化。设置最大进化代数G，设置个体个数N，设置进化迭代器t=0，随机生成N个初始个体P(t)。

(2)个体评价。计算各个个体的适应度。

(3)选择运算。将选择算子作用于群体，根据个体适应度，按照一定规则选择父代和母代。常见的选择算子包括轮盘赌选择法、随机遍历抽样法、局部选择法及锦标赛选择法等。

(4)交叉运算。将交叉算子作用于父和母，以一定的概率交换它们之间的染色体，从而产生新的个体。

(5)变异运算。将变异算子用于交叉产生的新个体，以一定的概率改变个体的基因。经过交叉和变异操作后，最终形成下一代种群P(t=t+1)。

(6)终止条件判断。若t

采用实数编码，选择算子采用“君主方案”，即对群体按照适应度排序，用这一代中最优个体一直作为父，其他偶数位个体依次作为母，二者进行交叉后产生新的个体。新个体采用多点变异产生下一代子群。本模块采用了两种不同的变异操作，一种为恒定变异概率，另一种采用大变异遗传算法。大变异遗传算法中，当群体种群过于集中时改变变异概率，改变变异概率条件为：

am·fmax

(7)

式中am为用户输入参数，fmax为种群中最大的适应度，favg为种群中平均的适应度。满足该条件时采用一个更大的变异概率。

1.3.2 粒子群优化算法

粒子群优化算法(PSO)是一种进化计算技术，源于对鸟群捕食的行为研究。其基本思想是：通过个体间的协作和信息共享，在群体中找最优解[16]。PSO的优势在于其易操作性、易实现性，调整参数较少，不受函数约束条件限制。PSO可广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向[16-17]。在设置空间范围内，单独搜索每个粒子的最优解，标记为该粒子的极值，然后共享每个粒子的极值与粒子群的其他粒子进行比较，在整个粒子群中寻找当前全局最优解，然后让粒子群中的所有粒子根据当前自己的极值和粒子群的当前全局最优解来调整自己的速度和位置[17]。算法流程如图6所示。

图6 粒子群优化算法流程

图7 0.5万和7.1万网格计算结果对比

(1)初始化，设置最大进化代数G，设置粒子群体个数N，设置进化迭代器t=0，随机生成每个粒子在各个设计空间维度j的位置xj和速度vj；

(2)计算各个粒子的适应度；

(3)根据适应度更新各个粒子的历史最优pbest；

(4)根据适应度更新整个群体粒子的全局最优gbest；

(5)根据位置和速度公式对每个粒子的速度和位置进行更新；

(6)进行边界条件处理；

(7)若达到最大迭代次数，则输出最优解并退出算法，否则回到步骤(2)。

步骤(5)中粒子速度和位置的更新公式如下：

vi，j(t+1)=w·vi，j(t)+c1r1[pi，j(t)-xij(t)]+c2r2[pg，j(t)-xi，j(t)]

(8)

式中，i为粒子编号，j为设计变量的维度，w为惯性权重，c1和c2为学习因子，r1和r2为[0，1]范围内的均匀随机数，pi，j为粒子i历史最优的第j分量，pg，j为全局最优粒子的第j分量。在得到速度后，采用下式可以得到更新后粒子的位置分量：

xi，j(t+1)=xi，j(t)+vi，j(t+1)

(9)

惯性权重w对算法的收敛速度有很大影响，本模块采用了动态惯性权重值，w的计算公式为：

(10)

式中，wmax为最大惯性权重，wmin为最小惯性权重，Tmax为最大进化次数。在大多数应用中wmax=0.9，wmin=0.4。

2 优化结果分析

2.1 CFD样本计算

本研究利用Fluent的TUI命令，批量计算了200个样本点。由于燃气在喷管扩张段属于超音速流动，初始样本算例设置选取Standard k-ε模型，考虑能量方程。由于本文优化的目标参数为比冲，在燃烧室室压、喉部直径、喷管长度等参数确定的情况下，比冲与出口燃气流速成正比，因此本文选取出口流速作为优化参数。为验证不同网格数对计算结果的影响，本文测试了网格数目为0.5万、1.1万、2.1万、3.3万、4.2万、5.6万、7.1万的情况，在初始结构下，喷管出口流速如表1所示。在网格数目大于2万后，计算得到的质量加权速度偏出小于2 m/s。图 7对比了网格数0.5万和7.1万情况下的速度分布，因此本文选择7.1万网格数进行所有样本点计算。

表1 不同网格数对计算结果的影响

2.2 优化结果

本算例优化算法采用遗传算法(GA)和粒子群算法(PSO)，种群数为500，进化代数设置为1 500。由于比冲正比于出口的速度，本算例的优化目标选为质量加权的出口轴向速度，图8给出了该速度随着进化代的演化过程(仅显示前100代)。由图可知在进化大约30代之后速度值达到收敛。

图9展示了优化前后喷管的外形对比，图10给出了优化前后喷管内的速度分布对比。基于优化前后的设计参数进行CFD计算，给出了优化前后喷管气动性能的对比。优化后出口质量加权速度为3 452.22 m/s，比初始情况增加了1.88 %。针对优化结果进行CFD校核计算，得到的出口质量加权速度为3 453.52 m/s，与基于代理模型的结果对比偏差小于0.04%。图11对比了优化后输出的云图和CFD计算结果的对比，由图可知，二者基本一致。这一现象表明：基于代理模型的预测云图可有效预测CFD仿真结果，用于仿真结果的有效预测评估。

图9 优化前及优化后的喷管外形对比

图10 优化前及优化后速度分布对比

图11 优化结果的预测云图与CFD计算结果对比

另外将本文所用的基于CFD一体化优化流程的优化结果与特征线法优化设计的喷管型面结果对比，喷管内型面结构相当，如图12所示。表2将不同优化算法的出口流速与特征线法优化的出口流速进行比较，本文使用的优化算法与特征线模型计算结果相当，最大相对误差为0.28 %。该结果表明：基于CFD一体化优化流程的优化结果较传统特征线法计算结果具有一致性，可优化得到喷管内型面的最优结构，并获取最优解下的喷管流场信息。

表2 不同优化算法计算结果对比

图12 特征线法优化结果与本文优化流程结果对比

对比遗传算法与粒子群优化算法的计算结果，两种算法的优化时间相近，计算时间均小于5 min。对比表2中两种优化算法的计算结果发现，遗传算法与粒子群算法计算结果相近，粒子群算法略优于遗传算法。

3 结论

本文基于三次B样条法构建喷管内型面结构，通过控制点坐标设计样本点，同时结合网格变形技术，实现样本点喷管内型面网格的自动变形，完成样本点的CFD自动计算。然后基于样本点的CFD计算结果，构建代理模型，代替优化过程中大规模的CFD计算，通过机器学习，在较短时间内实现了发动机喷管内型面结构的自动优化设计。通过本文研究得出以下结论：

(1)采用三次B样条法可有效描绘液体火箭发动机内型面结构，同时采用径向基函数插值法实现样本点网格的自动生成，能高效地实现前处理的自动化，有利于工程应用。

(2)基于代理模型的优化方法计算结果与CFD仿真结果相当，误差为0.04%，证实基于代理模型的优化结果可有效代替大规模的CFD计算优化。本文所用方法在30代后，优化结果基本保持稳定，高效地完成了发动机喷管内型面优化设计。

(3)基于代理模型的优化方法与传统特征线法计算相比较，优化得到的喷管内型面结构基本保持一致，对比发动机比冲的相对误差为0.28%，进一步证实了本文所使用的基于代理模型的喷管内型面优化方法可较精确、高效地实现液体火箭发动机喷管内型面优化设计，具有较强的工程应用性。