1 引言

高超声速飞行时，超燃冲压发动机燃烧室壁面承受着极大的热负荷，如果没有有效的热防护措施，即使采用耐高温性能出色的复合材料也无法保证燃烧室的正常工作^[1]，当前普遍认为利用飞行器自身携带的碳氢燃料对燃烧室壁面进行主动冷却具有较高的可行性^[2]。在主动冷却系统中，低温的碳氢燃料在注入燃烧室中燃烧前，一部分将作为冷却剂流入分布在燃烧室内壁上的再生冷却通道，与高温壁面进行对流换热，在换热过程中，再生冷却通道的结构参数对碳氢燃料的换热效率具有较大的影响^{[2, 3]}，因此，国内外开展了诸多以强化换热为目的的对再生冷却通道的设计研究，Zhang等^[3]认为矩形通道的截面积和肋片厚度的减小会有效降低壁温，但前者会造成流动压力损失的增大；牛禄等^[4]认为薄壁厚、大高宽比的矩形通道能的冷却效果更好，增大人工粗糙度可以增强流体扰动从而提高对流换热能力；吴峰等^[5]通过两种优化方案改变冷却通道的高宽比，发现存在着一个最佳冷却通道个数，使得推力室壁面的冷却效果达到最佳，在相同质量流量下，通过降低通道高度能够强化推力室传热，同时会造成流动压力损失的增大。但以上的研究都只是对几种不同的模型进行对比，却不能得到再生冷却通道的准确的最优解。鲍文等^[6]基于等壁温设计准则对再生冷却通道的高度和人工粗糙度进行局部强化换热优化设计；王厚庆等^[7]基于遗传算法，以惰性质量最小、冷却液流量系数最低为目标，对再生冷却通道的内壁厚度、槽宽、肋宽和肋高进行了优化设计。

目前多数的对再生冷却结构的优化设计均是在基于对流换热经验关联式上建立的简化的一维传热模型，忽略温度在横截面的变化及壁面在轴向的热传导等，而且没有重视流动压力损失这一关键问题。正如文献[8]所提，为得到可靠的冷却通道参数结果，三维因素必须予以考虑。另外，通常为了避免煤油“拟沸腾”造成的局部换热恶化，同时保持煤油的喷射特性，要在不同工况下始终保持冷却通道中的煤油处于较高的压力（一般要求处于超临界压力状态）^[9]，过大的流动压力损失可能会导致管道内一部分碳氢燃料处于亚临界压力状态。

本文为改善设计工况下的再生冷却通道的综合流动换热性能，将通道燃气侧平均壁温最低和流动压力损失最小作为优化目标，以单根再生冷却通道的肋高、槽宽和肋厚为设计变量，利用CFD软件开展对热/流/固耦合场的三维数值计算，采用中心复合设计法（Central Composite Design，CCD）并结合完全二阶多项式响应面模型（Standard Response Surface-Full 2nd Order Polynomials）拟合出响应面，详细分析了各结构参数对优化目标的影响，然后运用多目标遗传优化算法（Multi-Objective Genetic Algorithm，MOGA）进行迭代计算，得到Pareto最优解集，进而获得冷却平板的高度、通道数量的设计方案及允许的进口质量流量范围。

2 计算模型和数值方法 2.1 再生冷却通道模型和边界条件

本文参照文献[10]，把超燃冲压发动机燃烧室的四个壁面均看作是矩形平板，含有多个再生冷却通道的平板模型如图 1所示，沿宽度方向均匀分布N个大小相同的冷却通道，一定质量流量的碳氢燃料流入冷却平板后，均匀流入每个冷却通道，单根冷却通道的截面几何参数和受热环境如图 2所示，δ为肋厚，H为肋高，W为槽宽，d为外壁厚度，e为内壁厚度，与燃气侧接触的壁面受热，外壁和侧壁均视为绝热壁面。由于这些通道具有相同的受热环境，因此可取单根再生冷却通道作为研究对象进行优化设计，进而可求得平板的几何参数和通道数量。由于燃烧室结构是一定的，即冷却平板的宽度W_N是定值，故在冷却平板进口质量流量m_N一定的情况下，冷却平板的高度H_N，通道数量N和单根通道的进口质量流量m_in会随肋高H、槽宽W和肋厚δ的改变而改变，满足以下约束

$\begin{array}{*{20}{l}} {{H_N} = e + H + d}\\ {N\left( {W + \delta } \right) = {W_N} = {\rm{Const}}}\\ {N{m_{{\rm{in}}}} = {m_N} = {\rm{Const}}} \end{array}$

(1)

若给出初始参数N₀，W₀，δ₀，则优化过程中可求得

$\begin{array}{*{20}{l}} {{m_{{\rm{in}}}} = \frac{{{m_N}}}{{{N_0}}} \cdot \frac{{W + \delta }}{{{W_0} + {\delta _0}}}}\\ {N = {N_0} \cdot \frac{{{W_0} + {\delta _0}}}{{W + \delta }}} \end{array}$

(2)

本文给定一组初始结构参数，其中W₀=2mm，δ₀=4mm，H₀=2mm，d₀=3mm，e₀=2mm，N₀=100（W_N=600mm），l=2m，并假定设计工况下的m_N=4kg/s。

根据上述条件，单根再生冷却通道的边界条件设置如下:

（1）入口边界条件为质量流量入口，m_in=m_N(W+δ)/600kg/s，入口温度T_in=300K。

（2）出口边界条件设置为压力出口，p_out=3MPa，以此保证管内煤油始终处于超临界压力状态（RP-3煤油临界压力为2.2MPa^[11]）。

（3）在Ma > 6时，超燃冲压发动机燃烧室壁面峰值热流密度可达3MW/m²以上^[12]，本文将燃气侧壁面受热条件简化成恒定的热流，取q_w=3MW/m²。

（4）其余壁面热流密度均设置为0。

2.2 材料热物性

燃烧室壁面采用同一种材料—不锈钢1Cr18Ni9Ti，耐温极限为1473K^[13]。密度和比热容按常数处理，分别取值为7.9×10³kg/m³和502J/（kg·K），热导率按温度的线性函数处理^[14]，即

${\lambda _{\rm{w}}} = 0.01525T + 10.6$

(3)

单位为W/（m·K）。

碳氢燃料选用RP-3航空煤油，由于RP-3航空煤油是由几百上千种碳氢化合物组成的复杂物质，难以直接获得热物性数据，因此采用文献[11]提出的RP-3航空煤油的十组分替代模型，在NIST SUPERTRAPP计算程序中获得3MPa压力下煤油的密度、定压比热、黏度、热导率随温度的变化规律，如图 3所示。通过分段多项式拟合，可获得煤油的各物性参数与温度的函数关系式。

2.3 网格划分和数值方法

计算域包括流体域和固体域，均采用结构化网格划分，如图 4所示为计算域的截面网格。计算中考虑到近壁区的黏性影响层，控制流体域第一层网格到壁面的距离，使y⁺≤1，并在黏性影响层（即y⁺≤11）内划分10层以上网格，以此保证对流体域的精确计算。为确保计算结果的准确性，进行了网格的无关性验证。以进口质量流量为40g/s下的初始通道为例，共比较了5组网格，网格划分数量与计算结果如表 1所示，T_wg为燃气侧平均壁温，T_{w, max}为最高壁温，Δp为流动压力损失。5组网格的计算结果中，网格1与网格2的T_wg和T_{w, max}误差分别达到了2.14%，2.46%，Δp误差达到了5.24%，而网格2，3，4，5之间T_wg，T_{w, max}和Δp的最大误差均不超过1%，已与网格无关。对于优化过程中不同通道网格的划分，固体域和流体域截面网格数分别不低于1388和1734，轴向网格划分数均为600。

湍流模型采用对模拟变物性流动具有较高精度的RNG k-ε两方程模型，近壁区采用适用于低雷诺数流动的增强壁面处理法。通过有限容积法离散三维的Navier-Stokes方程，运用SIMPLEC算法求解压力-速度耦合方程，动量和能量方程均采用二阶迎风格式。

2.4 数值模拟验证

仲峰泉等^[11]和李中洲等^[15]分别在二级加热圆管内和长为575mm，直径为1mm的微小圆管内进行的煤油传热特性实验和流阻特性实验具有典型性，为验证本文选用的煤油热物性、湍流模型及数值方法的可靠性，分别对两者实验条件进行数值模拟，并与实验结果进行比较验证，结果如图 5所示。温度的计算值与实验值最大误差为8.3%；流动压力损失计算值与实验值的吻合度较好，最大误差为10.8%，故认为数值模拟结果是可靠的。

3 冷却通道的优化方案 3.1 优化设计过程及原理

图 6给出了优化设计的具体过程。在生成响应面的过程中，试验样本点的选择是非常重要的，由于CCD方法具有精度高、试验次数少、预测性好等优点，尤其适合对影响因素比较敏感的指标进行优化，故本文采用CCD方法生成试验样本点。对于k个设计变量，试验样本点由2^k-ξ个析因点加上2k个坐标轴点（±α，0，…，0），（0，±α，…，0），…，（0，0，…，±α）和1个中心点所组成，如图 7所示，其中ξ为析因系数，当0 < n≤4时，ξ=0，4 < n≤7时，ξ=1。

采用对未知试验点预测准确度较高的完全二阶多项式响应面模型，对试验样本点进行拟合生成响应面。对于n个设计变量，函数模型可表示为

$y\left( X \right) = {\beta _0} + \sum\limits_{i = 1}^n {\beta _i}{x_i} + \sum\limits_{i = 1} ^n {\beta _{ii}}x_i^2 + \sum\limits_{i = 2} ^n \sum\limits_{j = 1} ^{i - 1} {\beta _{ij}}{x_i}{x_j}$

(4)

式中y(X)是拟合函数，X=（x₁，x2，x₃，…x_n）是设计变量，β₀，β_i，β_ii，β_ij是待定系数，通过最小二乘法拟合确定。

全局寻优采用的MOGA具有思想简单、易于实现、全局搜索等优点，是一种应用较为普遍的寻优算法。利用MOGA对再生冷却通道进行循环优化设计，可在较少的计算量下准确地找到最优结果。

3.2 优化目标及设计变量

多目标优化的数学模型可表示为

$\begin{array}{*{20}{l}} {{\rm{min}}F\left( X \right) = \left[ {{f_1}\left( X \right), {f_2}\left( X \right), \cdots , {f_m}\left( X \right)} \right]}\\ {\left( {m = 1, 2, \cdots , M} \right)} \end{array}$

(5)

${\rm{s}}.{\rm{t}}.\;\;\;{g_i}\left( X \right) \le 0, \;\;\left( {i = 1, 2, \cdots N} \right)$

(6)

$ {h_j}\left( X \right) = 0, \;\;\left( {j = 1, 2, \cdots r} \right) $

(7)

$ \mathit{\boldsymbol{X}} = {({x_1}, {x_2}, \cdots {x_k})^{\rm{T}}} $

(8)

再生冷却通道的主要功用是保证不同工况下超燃冲压发动机燃烧室的壁温始终低于材料的耐温极限，由于与燃气直接接触的壁面温度最高，经受的考验最大，且考虑到燃气侧平均壁温更能反映壁面整体冷却效果，因此把T_wg作为优化目标之一，同时将最高壁温T_{w, max}作为反映局部冷却效果的输出参数；另外，为使通道中的煤油始终处于超临界压力状态，要尽可能地减小流动压力损失，因此目标函数为

$ {\rm{min}}{T_{{\rm{wg}}}}\left( X \right) = \frac{1}{{{A_{\rm{g}}}}}\mathop \smallint \nolimits^ T{\rm{d}}A $

(9)

$ {\rm{min\Delta }}p\left( X \right) = {p_{{\rm{in}}}} - {p_{{\rm{out}}}} $

(10)

式中A_g为通道燃气侧受热面积，p_in和p_out分别为进口压力和出口压力。

多目标优化最终会得到一个优化解集，为避免一些优化解的T_{w, max}超过壁面材料的耐温极限，或Δp过大（文献[10]认为高于1MPa是严重的），为优化解集增加了一组筛选条件

$\begin{array}{*{20}{l}} {{T_{{\rm{w}}, {\rm{max}}}} < 1473{\rm{K}}}\\ {{\rm{\Delta }}p < 1{\rm{MPa}}} \end{array}$

(11)

通过以往的研究发现，外壁厚度d和内壁厚度e对流动压力损失的影响可忽略^[16]，对换热效果呈单调影响，因此本文取对通道传热和流动均有较大影响的肋高H，槽宽W和肋厚δ作为待优化设计变量。即

$\mathit{\boldsymbol{X}} = {[{x_1}, {x_2}, {x_3}]^{\rm{T}}} = {\left[ {H, W, \delta } \right]^{\rm{T}}}$

(12)

考虑到室壁的承载能力等因素，冷却通道结构参数的取值范围有一定的设计要求。因此，取如下约束条件

$\begin{array}{*{20}{l}} {2{\rm{mm}} \le H \le 8{\rm{mm}}}\\ {1{\rm{mm}} \le W \le 9{\rm{mm}}}\\ {1{\rm{mm}} \le \delta \le 7{\rm{mm}}} \end{array}$

(13)

4 结果分析 4.1 响应面的检验

响应面的优劣程度将直接影响优化结果的准确性，因此需要对生成的响应面检测评判，一般采用可决系数R²和调整可决系数R_adj²来评价响应面模型对试验样本点的拟合程度，其定义详见文献[17]。通常要求R²，R_a²的值在0.9以上^[18]（最优值为1）。对于本文建立的响应面模型，T_wg，T_{w, max}和Δp的R²分别为0.9977，0.9941和0.9812，说明响应面预测值和试验样本点有高度的相关性，模型拟合程度较好；T_wg，T_{w, max}和Δp的R_adj²分别为0.9957，0.9930和0.9647，说明响应面模型分别能解释99.57%的T_wg，99.30%的T_{w, max}和96.47%的Δp响应值的变化，有总变异0.43%的T_wg，0.70%的T_{w, max}和3.53%的Δp不能用该模型来解释。

4.2 各参数对目标函数的影响

研究不同进口质量流量下各参数对目标函数的影响，限于篇幅，以m_N=4kg/s和m_N=6kg/s为例进行分析，δ=2mm时H，W共同作用下与T_wg，Δp的三维响应面如图 8，9所示，三维响应面能直观地表现设计变量对目标函数的影响。可以看出，随着H，W的增大，T_wg呈升高趋势，Δp呈减小趋势，T_wg，Δp的峰值分别出现在H，W的最大值和最小值处，这是由于H，W的增大使流道截面积增大进而导致流速减小，换热能力下降，流动压力损失减小；H方向的响应曲面倾斜度比W方向的更高，即T_wg，Δp变化得更快，这是由于进口质量流量m_in与H的大小无关，但会随着W的增大以相对较小的斜率增大，故H比W对流速的影响大，T_wg，Δp随H的增大变化得更快；相比图 8，图 9的T_wg的变化范围减小，而Δp的变化范围扩大。

图 8（a），图 9（a）中，W较大的通道T_wg随H的增大升高得更快，这是由于影响对流换热的因素包括流速、流体的物理性质、换热表面的几何因素等^[19]，在H相等时，W较大的通道流速小，换热能力低，此时流速的变化对T_wg的影响起主导作用，W较小的通道流速大，流速的变化对T_wg的影响变小，此时其它因素对T_wg的影响起主导作用，故随H的增大，W较大的通道T_wg受流速的影响程度更大，升高得更快；同样的，对于H较大的通道，T_wg随W的增大升高得更快。图 8（b）、图 9（b）中，随着H或W的减小，Δp增大的趋势变快，与文献[20]的研究结论吻合，这是由于流速随着H或W的减小线性增大，根据流体力学公式（14）可知

${f_{\rm{l}}} = \frac{{2{\rm{\Delta }}p \cdot D}}{{l \cdot {\rho _{\rm{l}}} \cdot {v^2}}}$

(14)

Δp与流速的平方呈线性关系，故流速的均匀增大将导致流动压力损失增大的趋势变快；还可明显看出，H较小的通道Δp随W的减小增大得更快，这是由于其截面积更小，流速更大，截面积以同等斜率变化时Δp受到的影响更大；同样的，W较小的通道Δp随H的减小增大得更快。

图 10给出了在m_N=4kg/s和m_N=6kg/s条件下的T_wg，Δp与δ的二维响应曲线。图 10（b）中，在不同进口质量流量下，随着δ的增大，Δp增大，这是由于δ的增大导致m_in增大，进而导致流速及流阻的增大；但在图 10（a）中，δ对T_wg的影响与H和W的大小密切相关，当H=7mm，W=7mm时，T_wg主要受流速的影响，随δ的增大而单调下降；当H=4mm，W=4mm和H=2mm，W=2mm时，T_wg随δ的增大先略微下降，而后逐渐升高，说明此时流速已不是影响到换热能力的主导因素。根据文献[20]所述，受导热和冷却效率的限制，肋片太厚，就不能够有效冷却，肋厚只有近似地满足如下条件才能发挥强化传热的作用

$\delta \le {\delta _{{\rm{cr}}}} = \frac{{2{\lambda _{\rm{w}}}}}{h}$

(15)

式中h是冷却液与壁面之间的对流换热系数。可以看出，临界肋厚δ_cr与壁面热导率λ_w成正比，与h成反比。由于H和W较小的通道流速大，湍流强度和h较大，故δ_cr较小，随着δ的增大，热阻的增大对换热的负作用逐渐超过了流速的增大对换热的强化作用，导致H和W较小的通道会出现T_wg随δ的增大而升高的现象。以初始通道为例，当m_N=4kg/s时，在通道的l/2处截面，壁面平均温度为644K，流体平均温度为485K，热导率λ_w=0.01525×644+10.6=21.03W/（m·K），粗略估算δ_cr < 3mm，故初始通道的肋厚对换热起到了负作用。

通过上述三图还可看出，不同工况下各设计变量对优化目标的影响规律是相似的。

图 11给出了m_N=4kg/s时目标变量T_wg和Δp对设计变量H，W和δ的灵敏度。由图可知，T_wg与H，W呈负相关，与δ呈正相关，Δp与三者的关系相反；T_wg和Δp对设计变量的灵敏度绝对值大小次序是相同的，依次是H，W和δ，因此对冷却通道进行设计时可优先考虑改变H值。

4.3 优化结果及分析

MOGA的具体设置为种群规模1000，迭代样本数200，Pareto最大许用值为70%，最大迭代数20。图 12给出了m_N=4kg/s时Pareto最优解集，在红线以左的点T_wg和Δp均小于初始通道的。颜色显示了与既定设计目标的匹配程度，蓝的表示好，红的表示差。可以看出，两个优化目标之间是背道而驰的，T_wg的降低就意味着Δp的增大，反之亦然，不存在一个最优设计点使两者同时达到最优状态；1000K的T_wg对应的Δp为121kPa，但是要用50kPa的Δp实现同样的T_wg是不可能的；当T_wg较高时，随T_wg降低，Δp略微增大，当T_wg较低时，T_wg的小幅下降也会造成较大的Δp。根据实际工作需要，可从Pareto解集中甄选出合适的设计方案。

表 2给出了初始通道与三组优化通道的结构参数及计算结果，可以看到，三组优化通道的δ都非常小，接近取值范围的最小值。相比初始通道，Opt1通道的流动和换热性能均有一定程度改善。图 13给出了两者的燃气侧壁面温度分布，Opt1通道的燃气侧壁温明显下降，而且在x方向的温差变化更小。图 14给出了两者在流体域y=H/2截面上的速度-压力分布，沿z增大方向，由于油温的升高导致密度的下降，进而导致流速逐渐增大，Opt1通道内的流速和流动压力损失相比初始通道显著减小。根据响应面模型计算的结果，Opt1通道比初始通道的T_wg下降了21.8K，T_{w, max}下降了23.5K，由于Opt1通道比初始通道的当量直径仅减小0.1mm，而流速减小了约1倍，故Δp大幅减小1.0088MPa（63.7%）。尽管流速减小，但壁温没有升高反而降低，导致这一现象主要的原因就是初始通道的肋片过厚，对换热起到较大的负作用，而Opt1通道的肋片相对较薄，强化了换热。

Opt2通道的T_wg和Δp在Pareto解集中相对折衷，相比初始通道，T_wg升高20.9K，而Δp减小了1.4272MPa，因此也可以认为Opt2通道的综合流动换热性能得到一定改善；Opt3通道更关注于流动压力损失。将4组通道进行完整的CFD计算，与模型计算结果相比，T_wg，T_{w, max}的最大误差分别为1.4%，1.3%，Δp的最大误差为5.7%，由于Δp有一个很大的变化（从3.07kPa~6.34241MPa），因此认为Δp的误差可接受。由此说明了本文采用的响应面方法结合MOGA可以有效地优化冷却通道。

实际上受飞行高度、马赫数等多种因素影响，再生冷却平板并不是一直在设计工况（m_N=4kg/s）下工作的，在非设计工况下工作可能会导致最高壁温或流动压力损失超出限制，因此需要给出设计通道能够正常工作的进口质量流量范围。考虑到Opt2通道的两个目标值相对折衷，工作范围更广，故以Opt2通道为例进行分析。首先将Opt2通道进行圆整，之后得到冷却平板的设计方案如表 3所示。

采用同样的响应面方法，以m_in为设计变量，求得圆整后的Opt2通道在（10，100）g/s范围内的T_{w, max}和Δp随m_in的变化，如图 15所示。当m_in≤15.7g/s时，T_{w, max}超过了1473K，当m_in≥98g/s时，Δp超过了1MPa，说明圆整后的Opt2通道在m_in处于（15.7，98）g/s内可正常工作，则对于冷却平板能够正常工作的m_N范围是（1.539，9.604）kg/s。

5 结论

（1）对再生冷却平板结构参数及通道数量的优化设计，可以通过对单根再生冷却通道的结构参数优化计算来实现，从而增大了优化设计的可行性，降低了复杂度，节约了计算时间。

（2）采用CCD法结合完全二阶多项式响应面模型拟合生成的响应面来解决超燃冲压发动机再生冷却结构的优化设计问题，克服了传统优化方法对经验关联式的依赖，计算精度高。

（3）不同进口质量流量下的再生冷却通道设计变量对目标函数的影响规律是相似的，影响程度的大小顺序依次是肋高、槽宽、肋厚，设计工况下三者对燃气侧平均壁温的灵敏度分别是0.712，0.284，-0.072，对流动压力损失的灵敏度分别是-0.647，-0.504，0.270。

（4）得到Pareto最优解集，从中可选取出多组综合流动换热性能优于初始通道的结构。

（5）得到一组冷却平板的设计方案及正常工作允许的进口质量流量为1.539~9.604kg/s。