2. 翼型、叶栅空气动力学重点实验室，西北工业大学航空学院，陕西西安 710072
2. National Key Laboratory of Science and Technology on Aerodynamic Design and Research, School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, China
高超声速技术是未来航空航天领域的战略制高点.随着X-43, X-51等高超声速飞行器成功进行飞行试验, 标志着高超声速飞行所需的基础技术已取得重要突破, 未来高超声速飞行器将向着速域更宽、空域更广的方向发展, 并最终实现飞行器单级入轨的能力.单级入轨空天飞行器实际飞行中必然要经历从亚声速、跨声速、超声速到高超声速巡航的多个飞行阶段.因此, 除了保证高超声速巡航设计点的飞行性能外, 这类飞行器的设计必须兼顾满足工程需要的亚、跨和超声速气动特性以确保良好的起降和加速性能.但是, 各个速域绕流的空气动力学特性有很大不同, 保证良好的气动性能所要求的外形/构型也存在矛盾, 同时还必须考虑到高超声速气动热防护, 这些给宽速域飞行器的设计带来了巨大的挑战.
高超声速飞行时, 翼型的作用相对减弱, 只要求尖前缘、薄机翼即能较好地满足高超声速飞行的需要.然而, 在亚声速及跨声速状态下, 翼型的设计对飞行器的升阻、操稳及控制特性具有重要的影响, 好的翼型设计是提高亚、跨声速气动性能的关键.另一方面, 尽管高超声速飞行器通常具有大后掠角、小展弦比的特征, 机翼的三维效应明显, 但从翼型入手研究宽速域飞行器流动机理仍然存在一定的价值和意义.随着高超声速飞行器对宽速域气动性能的需求, 近年来有学者开展了兼顾不同速域气动性能的宽速域翼型优化设计研究.国外, Ueno等[1]对高超声速翼型进行了优化设计, 在保证高超声速高升阻比的同时, 兼顾了翼型的跨声速气动性能.国内, 西北工业大学曹长强等[2]开展了超声速多边形翼型和双弧形翼型气动优化设计研究, 取得了一定的成果.西北工业大学孙祥程等[3]也开展了高超声速宽速域翼型优化设计研究, 并发明出一种能够兼顾跨声速与高超声速气动性能的宽速域翼型.研究发现, 相比普通翼型, 该翼型配置到三维机翼上仍能保持更好的宽速域性能.宽速域翼型优化设计是典型的多目标优化设计过程, 现有的宽速域翼型优化设计工作多采用线性加权法处理多目标优化问题, 尚未有宽速域翼型多目标优化问题的Pareto最优解集[4]的研究报道.但线性加权法处理多目标优化问题时, 优化效果强烈依赖于设计者的经验, 而且优化得到的结果比较单一.为了更好地理解宽速域翼型优化中亚、跨、超和高超声速状态下气动性能协调的空气动力学原理, 最大化地提升优化翼型的宽速域性能, 通过优化手段获取宽速域翼型多目标优化设计的Pareto前沿是十分有必要的.
本文首先对薄翼型在不同速域下的流动机理进行了分析, 总结了翼型在各速域下增升减阻的设计准则.然后利用课题组自主开发的RANS方程求解器程序PMNS2D[5-7], 结合基于代理模型的多目标多约束高效通用优化软件“SurroOpt”[8], 分别采用线性加权法与Pareto最优解集法开展了兼顾亚、超和高超声速气动性能并考虑了气动热防护的宽速域翼型优化设计研究.采用Pareto解集法的宽速域翼型优化设计研究中得到了包含一系列优化翼型的翼型簇.通过对Pareto最优化解集中的翼型进行研究, 分析了宽速域翼型协调各速域气动性能的空气动力学原理.预期本文发展的宽速域翼型优化设计流程对于多目标气动优化设计问题具有一定的工程参考价值.
1 流动数值模拟方法采用课题组自主开发的RANS方程流场求解器PMNS2D进行数值模拟.在连续介质假设下, 惯性坐标系下忽略彻体力和热源的二维非定常可压缩Reynolds平均Navier-Stokes方程可写成如下积分形式
| $ \iiint\limits_\mathit{\Omega } {\frac{{\partial \mathit{\boldsymbol{W}}}}{{\partial t}}{\text{d}}V} + \iint\limits_{\partial \mathit{\Omega }} {\mathit{\boldsymbol{H}} \cdot \mathit{\boldsymbol{n}}{\text{d}}S} - \iint\limits_{\partial \mathit{\Omega }} {{\mathit{\boldsymbol{H}}_{\text{v}}} \cdot \mathit{\boldsymbol{n}}{\text{d}}S} = 0 $ | (1) |
式中, W为流动守恒变量, H为无黏通量, Hv为黏性通量, Ω为控制体,
本文采用格心有限体积法在结构化网格上求解方程组(1), 时间推进采用隐式LU-SGS格式, 模型采用一方程Spalart-Allmaras[9]湍流模型实现控制方程的封闭.湍流模型实现控制方程的封闭.跨声速状态下, 空间离散采用中心格式, 高超声速状态下则采用高阶迎风AUSM+_up格式.
2 优化设计方法 2.1 CST几何参数化方法Kulfan等[10-11]提出了一种基于型函数/类函数变换的参数化方法(class function/shape function transformation, CST), 参数具有明确的几何意义, 控制参数少, 适应性强, 精度好.
采用CST函数方法对翼型进行参数化的表达式如下:
上表面
| $ {y_{\rm{u}}} = C\left( x \right) \cdot {S_{\rm{u}}}\left( x \right) + x \cdot {y_{{\rm{TEu}}}} $ |
下表面
| $ {y_{\rm{l}}} = C(x) \cdot {S_{\rm{l}}}(x) + x \cdot {y_{{\rm{TEl}}}} $ |
式中, C(x)为类函数, S(x)为型函数, yTEu, yTEl分别为上下表面后缘的坐标.
本文采用的是8阶CST参数化方法, 共18个设计变量.
2.2 约束翼型前缘半径的方法本文宽速域翼型优化设计中, 通过约束翼型前缘半径来限制前缘驻点处的气动加热.具体过程如图 1(a)对于由CST参数化生成的原始翼型, 在翼型前缘处按指定半径倒圆; 图 1(b)倒圆处理使弦长有一定的损失, 对翼型的升阻特性产生影响, 为了对比的公平性, 将倒圆后的翼型弦长重新变为“1”.具体的处理方法是:保持前缘圆形部分不变, 将翼型后段沿弦向均匀拉长使翼型弦长恢复到“1”.
|
| 图 1 约束翼型前缘半径的处理过程 Fig.1 Method to constrain leading edge radius |
Kriging模型是一种插值模型, 其插值结果定义为已知样本函数响应值的线性加权.该方法用相关函数描述不同空间位置间样本点的关联程度, 并寻找最优加权系数使预测函数与真实响应值的均差最小, 从而确定代理模型. Kriging模型预测值和均方误差公式分别为
| $ \hat y(\mathit{\boldsymbol{x}}) = {\beta _0} + {\mathit{\boldsymbol{r}}^{\rm{T}}}(\mathit{\boldsymbol{x}}){\mathit{\boldsymbol{R}}^{ - 1}}\left( {{\mathit{\boldsymbol{y}}_{\rm{s}}} - {\beta _0}\mathit{\boldsymbol{F}}} \right) $ |
| $ \begin{array}{l} {\rm{MSE}}\left\{ {\hat y(\mathit{\boldsymbol{x}})} \right\} = {s^2}(x) = \\ \;\;\;\;{\sigma ^2}\left\{ {1.0 - {\mathit{\boldsymbol{r}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{r}} + {{\left( {1 - {\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{r}}} \right)}^2}/{\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{F}}} \right\} \end{array} $ |
式中, R为相关矩阵, 由所有已知样本点之间的相关函数值组成; r为相关矢量, 由未知点与所有已知样本点之间的相关函数值组成; F为单位向量; 而β0=(FTR-1F)-1FTR-1ys, ys为已知样本点处的响应值.Kriging模型的具体细节可参考文献[12].
2.4 多目标优化方法在优化设计中同时要求几项设计指标达到最优的优化问题, 称为多目标优化问题.多目标优化也称为多准则优化或矢量优化, 用于找到满足约束值并且所有目标函数达到最优的决策向量解集.它可以定义为
| $ \left\{ \begin{array}{l} \min \;\;\;\;\;\mathit{\boldsymbol{F}}\left( X \right) = \left[ {{f_1}\left( X \right),{f_2}\left( X \right), \cdots ,{f_m}\left( X \right)} \right]\\ {\rm{s}}.{\rm{t}}.\;\;\;\;\;\;\mathit{\boldsymbol{X}} = {\left[ {{x_1},{x_2}, \cdots ,{x_n}} \right]^{\rm{T}}}\\ \;\;\;\;\;\;\;\;\;{g_i}\left( X \right) \le 0,\quad i = 1,2, \cdots ,p\\ \;\;\;\;\;\;\;\;\;{h_j}\left( X \right) = 0,j = 1,2, \cdots ,q \end{array} \right. $ |
式中, X⊂Rn为决策向量; F(X)⊂Rm为目标向量; gi(X)为第i个不等式约束, hj(X)为第j个等式约束, m为优化目标的个数, n为设计变量的个数.
2.4.1 线性加权法[4]线性加权法是将各个分目标函数加权到一个总的目标函数, fi(x)将多目标优化问题转化为单目标优化问题F(x), 于是优化过程中的步骤与单目标优化问题基本相同.在对F(x)进行极小化的过程中, 先将各向设计指标都转化为统一的无量纲值, 并将量级限定于0~1之间, 使目标规格化, 然后根据各分目标的重要性分别赋予加权因子.
本文使用线性插值的方法将各分目标函数规格化.对每项设计指标fi(x), 根据经验, 我们总能大致地预计其变化范围
| $ {\alpha _i} \le {f_i}\left( \mathit{\boldsymbol{x}} \right) \le {\beta _i}\;\;\;\;i = 1,2, \cdots ,q $ |
对fi(x)进行转换, 得到规格化的设计目标
| $ {{f'}_i}\left( \mathit{\boldsymbol{x}} \right) = \frac{{{f_i}\left( \mathit{\boldsymbol{x}} \right) - {\alpha _i}}}{{{\beta _i} - {\alpha _i}}}\;\;\;\;i = 1,2, \cdots ,q $ |
则多目标优化问题的数学模型变成以下形式
| $ \begin{array}{l} \min :\;\;\;\;F\left( \mathit{\boldsymbol{x}} \right) = \sum\limits_{i = 1}^q {{\omega _i}{{f'}_i}\left( \mathit{\boldsymbol{x}} \right)} \\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;\;{g_j}\left( \mathit{\boldsymbol{x}} \right) \ge 0,\;\;\;\;j = 1,2, \cdots ,{N_{\rm{G}}}\\ \;\;\;\;\;\;\;\;\;x_i^1 \le {x_i} \le x_i^{\rm{u}},\quad i = 1,2, \cdots ,n \end{array} $ |
式中，ωi是加权因子, 它是根据第i项设计指标在优化中的重要程度给定的, 加权因子的选择对优化结果具有决定性的影响; NG代表优化模型中约束函数的个数，nv是优化模型中设计变量的个数.
2.4.2 Pareto解集法[4]Pareto解集法是求解一次优化问题即可得到不同目标之间存在权衡的多个最优解(Pareto解集)的直接方法.一些关于Pareto解集的基本概念如下:
i. Pareto占优[24]. P=[p1, p2, …, pm]∈Rm和Q=[q1, q2, …, qm]∈Rm, 当且仅当∀i∈{1, 2, …, m}: pi≤qi∧∃j∈{1, 2, …, m}: pj＜qj满足时, 称Q比P占优或Q支配P, 记为P
ii. Pareto最优解.设Ω是Rm中满足约束条件的可行区域, 对于点X*∈Ω, 当且仅当∀i∈{1, 2, …, m}, ∀X∈Ω-{X*}: fi(X*)≤fi(X)∧∃j∈{1, 2, …, m}: fi(X*)＜fi(X)满足时, 称X*为该问题的Pareto最优解.
iii. Pareto最优解集.对于给定的多目标优化问题, Pareto解集是由所有的Pareto最优解组成的集合, 即PS*={X∈Ω|∃X′∈Ω:F(X′)
iv. Pareto前沿. Pareto前沿是Pareto最优解集PS*中使用决策变量获得的目标函数集合, 即: PF*={F(X)=(f1(X), f2(x), …, fm(X)):X∈PS*}.
(1) 多目标遗传算法: NSGA-Ⅱ[4]
多目标遗传算法是直接求解多目标最优解集(Pareto解集)的有效方法.由Deb等[25]于2000年提出来的NSGA-Ⅱ被公认为是最有效的多目标遗传算法之一, 从提出至今被广泛应用于各个领域.本文采用NSGA-Ⅱ求解多目标子优化问题.设种群规模为N, 进化代数G. NSGA-Ⅱ的实现步骤如下:
① 将第k代种群Pk进行“选择”, “交叉”, “变异”, 操作, 产生中间种群Qk;
② 将Pk与Qk合并成规模为2N的种群;
③ 对两倍的种群进行快速非支配排序, 将其按从优到劣分为若干等级;
④ 选择排序位于1~m的个体, 其中m是满足前m个等级个体总数Nm大于等于N的最小整数, 若Nm等于N, 则将这Nm个个体作为第k+1代种群Pk+1; 若Nm大于N, 则前m-1等级的直接进入第k+1代种群Pk+1, 将第m级的个体进行拥挤距离排序, 前(Nm-N)个个体进入第k+1代种群Pk+1.
⑤ 若k+1＜G, 转①, 否则, 停止.
图 2给出了从第k代进化到第k+1代的流程图.
|
| 图 2 多目标优化算法NSGA-Ⅱ流程图 Fig.2 Flow chart of multi-objective optimization algorithm NSGA-Ⅱ |
(2) 基于Kriging代理模型的多目标遗传算法
为了开展宽速域翼型的多目标优化问题, 本文采用了基于Kriging模型的NSGA-Ⅱ多目标优化方法, 具体的优化流程如图 3所示.
|
| 图 3 基于Kriging代理模型的NSGA-Ⅱ多目标优化流程图 Fig.3 Multi-objective optimization flow chart of NSGA-Ⅱ based on Kriging surrogate model |
要进行高超声速飞行器宽速域翼型优化设计研究, 首先须要对薄翼型在不同速域下的流动机理进行分析, 以便对宽速域翼型优化设计结果的正确性有定性的判断.选取4%厚度的NPU-Hyper-04翼型, NACA64A-204翼型和四边形翼型进行不同速域下薄翼型的流动机理分析, 其中NPU-Hyper-04[3]翼型是课题组自主开发的一种兼顾跨声速与高超声速气动特性的宽速域翼型.
3.1 数值模拟程序有效性验证对RAE2822翼型在跨声速状态下进行翼型绕流计算, 计算网格如图 4所示, 计算状态为Ma=0.73, α=2.79°, Re=6.5×106.根据图 5可以看出计算的压力分布与实验结果吻合较好, 验证了流场求解程序PMNS2D的可靠性.对比了实验值与PMNS2D求解器计算所得的力系数, 可以看出, 数值模拟得出的阻力系数稍微偏大, 而升力系数则十分接近实验值.
|
| 图 4 RAE2822亚临界翼型计算网格示意图(417×129) Fig.4 Computational grids of RAE2822 airfoil(417×129) |
|
| 图 5 RAE2822翼型计算压力分布与实验值对比 Fig.5 Comparison of pressure coefficient distribution of RAE2822 airfoil |
|
下载CSV 表 1 RAE2822翼型力系数计算结果与实验值对比 Tab.1 Comparisons of force coefficients of RAE2822 airfoil |
图 6为亚声速状态下NPU-Hyper-04翼型(专利号: ZL105752315.A)在不同迎角下的压力云图.随着迎角的增大, 翼型表面的流动是一个薄翼分离过程, 且当迎角达到10°以上时, 翼型上表面流动会完全分离. 图 7为NPU-Hyper-04翼型在不同迎角下的压力分布对比.随着迎角的增大, 翼型下表面的压力变化不大, 升力提升主要来源于上表面压力减小, 因此在亚声速状态下, 随迎角增大, 薄翼型的升力增量主要是由前缘分离涡使翼型上表面压力减小产生的.
|
| 图 6 NPU-Hyper-04翼型不同迎角下的压力云图(Ma=0.3, Re=5.59×107) Fig.6 Pressure contours of NPU-Hyper-04 airfoil at different angles of attack |
|
| 图 7 NPU-Hyper-04翼型不同迎角下的压力分布对比 Fig.7 Comparisons of pressure coefficient distributions of NPU-Hyper-04 airfoil at different angles of attack |
图 8为NPU-Hyper-04翼型与同等厚度的四边形翼型的压力云图对比, 由于NPU-Hyper-04翼型有一定的弯度, 因此在小迎角下, 其上表面前缘分离涡形成的低压区的面积更大, 因此亚声速小迎角下NPU-Hyper-04翼型具有比四边形翼型更大的升力系数.观察NPU-Hyper-04翼型在各迎角下的压力分布, 在其后缘处有一个局部高压区提供了额外的升力, 这是由该翼型在后缘处下弯的“S”外形而导致的后加载提供的附加升力.
|
| 图 8 四边形翼型与NPU-Hyper-04翼型压力分布对比 Fig.8 Comparsions of pressure coefficient distributions of NPU-Hyper-04 airfoil at different angles |
结合上述分析, 尽管上表面的分离涡是升力的主要来源, 但是在大迎角下(10°以上), 薄翼型上表面会完全发生分离, 上表面压力分布趋于一致(如图 8所示).因此大迎角下若要提升薄翼型的升力系数是很困难的, 且改进的余地非常有限, 因为此时只能通过改变翼型下表面的外形来寻求升力小量的提升.
3.3 超声速流动机理分析为了研究最大厚度位置对薄翼型超声速阻力特性的影响, 对具有不同最大厚度位置的9个4%厚度的四边形翼型进行计算分析. 图 9为9个翼型的外形示意图, 它们的最大厚度位置与弦长的比分别是0.1，0.2，0.3，0.4，0.5，0.6，0.7，0.8，0.9. 图 10是描述各翼型在0°迎角下的阻力系数柱状图, 可以发现翼型的阻力系数随着最大厚度位置的后移先减小后增大, 最大厚度位置在0.5c~0.6c(c为翼型的弦长, 下同)左右时翼型具有最小的阻力系数.
|
| 图 9 不同最大厚度位置的四边形翼型示意图 Fig.9 Quadrilateral airfoils with different maximum thickness locations |
|
| 图 10 不同最大厚度位置的四边形翼型的阻力系数对比(Ma=1.2, Re=8.16×107, α=0°) Fig.10 Comparisons of drag coefficients of quadrilateral airfoils with different maximum thickness positions (Ma=1.2, Re=8.16×107, α=0°) |
图 11对比了最大厚度位置在0.2c，0.5c和0.8c的四边形翼型在超声速0°迎角下的压力云图.最大厚度位置位于0.2c的翼型由于翼型头部的钝度较大, 超声速情况下在头部形成很强的激波, 使其激波阻力大大强于其他两个翼型; 最大厚度位置在0.8c的翼型的头部激波理论上是3个翼型中最弱的, 其激波阻力也最小.但由于其最大厚度位置过于靠后, 导致翼型最大厚度位置之后部分的膨胀波较强, 压力更小, 使翼型表面压力积分在阻力方向的分量较大.最大厚度位置在0.5c的翼型, 由于最大厚度位置比较适中, 其头部激波较弱, 激波阻力较小, 而且压力分布积分在阻力方向上的分量也较小.因此在三者摩擦阻力相差不大的情况下, 弦长位于0.5c的翼型具有最小的压差阻力, 它的阻力是3个对比翼型中最小的.上述分析结果表明最大厚度位置适当的后移可以减小翼型超声速时的阻力.
|
| 图 11 不同最大厚度位置的四边形翼型超声速压力云图对比(Ma=1.2, Re=8.16×107, α=0°) Fig.11 Comparison of pressure contours at supersonic state (Ma=1.2, Re=8.16×107, α=0°) |
高超声速状态下NACA64A-204翼型、四边形翼型与NPU-Hyper-04翼型的压力云图与压力分布对比如图 12所示. NACA64A-204翼型不是为高超声速飞行而设计的, 与其他两种翼型相比, 其最大厚度位置比较靠前, 前缘钝度比较大, 使翼型头部的激波非常强, 激波阻力远大于NPU-Hyper-04翼型及四边形翼型.不仅如此, 前缘钝度大使激波脱体, 会使翼型下表面前缘的压力减小, 却会增大上表面前缘的压力, 使翼型的升力有较大的损失.在该计算状态下, 3种翼型中NACA64A-204翼型不仅是阻力系数最大的翼型, 也是升力系数最小的翼型, 其高超声速升阻比性能极差.以上分析说明涉及高超声速飞行的翼型设计时应注意控制翼型前缘的钝度, 减小前缘的钝度有利于提高翼型的升阻特性.
|
| 图 12 不同翼型高超声速状态下压力云图与压力系数分布对比 Fig.12 Comparisons of pressure contours and pressure coefficients for different airfoils at supersonic state |
观察图 12中NPU-Hyper-04翼型的压力分布, 翼型下表面前、后缘各有一个局部压力增量, 大大提高了翼型的升力.这是由于NPU-Hyper-04翼型下表面双“S”外形使流动在下表面前缘等熵压缩, 起到增升减阻的作用, 在下表面后缘形成二次压缩增加翼型的升力. NPU-Hyper-04翼型下表面外形设计有利于提高翼型的高超声速升阻比.
4 宽速域翼型优化设计考虑到高超声速飞行器起飞、降落、加速、爬升以及高超声速巡航性能, 选择翼型优化的3个设计工况如下:
(1) Ma=0.3, Re=5.59×107, α=12°, 要求升力系数最大;
(2) Ma=1.2, Re=8.16×107, α=0°, 要求阻力系数最小;
(3) Ma=6, Re=2.47×107, α=4°, 要求升阻比最大.
其中, 工况(1)对应的是低速典型起飞状态, 工况(2)对应的是超声速加速状态, 工况(3)对应的是高超声速巡航状态.
4.1 线性加权法处理多目标优化问题以NACA64A-204翼型为基准翼型, 采用线性加权法开展高超声速飞行器宽速域翼型优化设计研究.
优化问题的数学模型为
| $ \begin{array}{l} \max \;\;\;\;{W_1} \cdot {A_1} \cdot {C_{1,Ma = 0.3}} - {W_2} \cdot {A_2} \cdot {C_{{\rm{d}},Ma = 1.2}} + \\ \;\;\;\;\;\;\;\;\;{W_3} \cdot {A_3} \cdot {\left( {{C_1}/{C_{\rm{d}}}} \right)_{Ma = 6}}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;t - 0.99{t_0} \ge 0\\ \;\;\;\;\;\;\;\;{r_{{\rm{le }}}} = 0.001\;25c \end{array} $ |
式中, Cl, Ma=0.3和Cd, Ma=1.2分别为优化翼型的亚声速升力系数与跨声速阻力系数, Cl, Ma=6为优化翼型在高超声速下的升力系数, Cl, Ma=0.8为优化翼型在跨声速下的升力系数, t为优化翼型的厚度. W1, W2, W3为各目标的加权系数, A1=1/Cl0, Ma=0.3, A2=1/Cd0, Ma=1.2, A3=1/(Cl/Cd)Ma=6.0为归一化系数, 下角标有“0”的为基准翼型的力系数和翼型厚度.
亚声速状态下的力矩参考点取在距离翼型前缘0.25c弦长处, 超声速状态与高超声速状态下的力矩参考点取在距离翼型前缘0.5c弦长处; 控制翼型的前缘半径为0.001 25c, 以限制高超声速时前缘驻点的热流.翼型厚度在(4±0.01)%c范围内.翼型外形的变形采用8阶直接CST参数化方法.在开始优化之前, 采用Latin超立方抽样生成初始样本点20个, 优化过程中同时采用EI+MSP+LCB+MSE+PI[13]5种加点准则每次增加5个点.
利用课题组自主开发的网格生成程序, 网格如图 13, 图 14和15所示.采用C型网格, 亚声速状态下网格远场距离为100c, 第1层网格高度为1×10-5, 网格量5.8×104; 超声速状态下远场距离为150c, 第1层网格高度为1×10-5, 网格量为5.8×104; 高超声速时采用C型网格, 对激波位置进行加密处理, 其远场距离变为20c, 第1层网格高度为1×10-5, 网格量为5.8×104.
|
| 图 13 亚声速状态计算网格示意图(361×161) Fig.13 Computational grids at subsonic state(361×161) |
|
| 图 14 跨声速状态计算网格示意图(361×161) Fig.14 Computational grids at transonic state(361×161) |
|
| 图 15 高超声速状态计算网格示意图(361×161) Fig.15 Computational grids at hypersonic state(361×161) |
采用上述优化模型, 总共进行了4轮优化.如图 16为优化翼型的几何外形对比图.相比基准翼型, 优化后翼型的弯度有显著的减小, 最大厚度后移到0.55c左右, 翼型上表面变化平缓, 以保证其超声速下具有良好的阻力特性.下表面前、后缘处略微向内凹陷, 以保证在阻力增加不多的情况下提高翼型在亚声速与高超声速状态下的升力, 这是与NPU-Hyper-04翼型相似的特征, 有助于提升翼型的亚声速升力与高超声速升阻比.
|
| 图 16 优化翼型与基准翼型及四边形翼型几何外形对比 Fig.16 Comparisons of geometric shape |
表 2列出了优化翼型与基准翼型各设计指标的对比.相比基准翼型, 优化翼型亚声速升力降低了2.91%, 但亚声速升力少量的损失使得其超声速阻力与高超声速升阻比性能得到巨大的改善.其中优化翼型的跨声速阻力减小了36.06%, 而高超声速升阻比提升了84.12%, 改善非常明显.这体现了宽速域翼性设计是在不同速域下的气动性能互相协调与权衡的过程. 表 3列出优化翼型与常规的高超声速四边形翼型各项设计指标的对比:优化翼型比四边形翼型的亚声速升力系数提高了4.46%, 但跨声速阻力系数增加了4.24%, 高超声速升阻比降低了2.86%.
|
|
下载CSV 表 2 优化翼型与基准翼型气动特性对比 Tab.2 Comparison of aerodynamic performance between baseline and optimized airfoil |
|
|
下载CSV 表 3 优化翼型与四边形翼型气动特性对比 Tab.3 Comparison of aerodynamic performance between quadrilateral airfoil and optimized airfoil |
图 17与18分别为优化翼型与基准翼型及四边形翼型在亚声速状态与高超声速状态下的压力分布对比.可以发现, 在亚声速状态下, 各翼型升力的差异主要来源于翼型下表面压力分布的不同, 优化翼型比四边形翼型亚声速升力更大的原因在于其下表面后缘处下弯的“S”外形提供了一个额外的压力增量.高超声速状态下, 优化翼型下表面也表现出了前、后加载效应, 有利提高翼型在高超声速状态下的升力, 进而提升其升阻比.
|
| 图 17 优化翼型与基准翼型及四边形翼型亚声速压力分布对比 Fig.17 Comparison of pressure coefficient distributions at subsonic state |
|
| 图 18 优化翼型与基准翼型及四边形翼型高超声速压力分布对比 Fig.18 Comparison of pressure coefficient distributions at hypersonic state |
采用线性加权法处理多目标优化问题时权重系数的选取严重依赖于设计者的经验, 极有可能因为权值系数选取不合适而无法找到综合性能最优的翼型.而且4.1节中通过线性加权法开展宽速域翼型优化设计仅得到一个优化翼型, 优化结果单一, 缺乏足够的信息来判断权重系数取得是否合适, 也难于判断优化翼型是否还有进一步改进的潜力.
现在采用Pareto解集法处理多目标优化问题, 将翼型的亚声速升力与超声速阻力作为优化目标, 而将高超声速升阻比作为约束(大于具有相同厚度的四边形翼型高超声速升阻比的90%).以4.1节中得到的优化翼型为基准翼型, 开展高超声速宽速域翼型多目标气动优化设计.优化问题的数学模型为
| $ \begin{array}{l} \max \;\;\;\;{C_{1,Ma = 0.3}}\\ \min \;\;\;\;\;{C_{{\rm{d}},Ma = 1.2}}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;\;\;{\left( {{C_1}/{C_{\rm{d}}}} \right)_{Ma = 6.0}} \ge 4.75\\ \;\;\;\;\;\;\;\;\;\;\;{r_{{\rm{le}}}} = 0.00125c\\ \;\;\;\;\;\;\;\;\;\;\;t - 0.99{t_0} > 0 \end{array} $ |
优化数学模型中各项的意义以及不同速域下所使用的计算网格与4.1节相同.翼型的变形采用8阶CST方法.在开始优化前, 使用Latin超立方抽样生成40个初始样本点. 图 19为优化过程中加点300次的优化历程, 优化得到的Pareto解集上有37个翼型, 它们是对翼型的亚声速升力指标与超声速阻力指标的重要性作不同的权衡时得到的优化翼型的集合.
|
| 图 19 两目标优化的加点历程 Fig.19 Process of adding new sample points by two-objective optimization |
图 20为经过优化设计后找到的Pareto前沿.图中还标注出了NACA64A-204翼型、四边形翼型和基准翼型(Opt-4)在目标空间中的位置.相比NACA64A-204翼型, Pareto前沿已经推进了很多(箭头方向是Pareto前沿推进的方向), 而四边形翼型与Opt-4翼型位于Pareto前沿上.
|
| 图 20 两目标优化得到的Pareto前沿 Fig.20 Pareto front obtained by two-objective optimization |
从Pareto解集中选出4个翼型(图 20中用1，2，3，4标注出的翼型)进行分析.其中1是Pareto解集中超声速阻力系数最小的翼型, 4是亚声速升力系数最大的翼型. 表 4给出4个翼型的外形示意图(表中从上到下依次是翼型1, 2, 3和4)及其气动力系数. 4个翼型均是最优翼型, 因为对各优化目标的侧重有所不同而导致外形与性能均有较大差异.从翼型1到翼型4, 翼型的超声速阻力增加, 亚声速升力得到改善.此过程中翼型外形的变化主要体现在前、后缘下弯以及下表面前缘附近向内凹陷, 即随着翼型弯度的变大以及下表面前缘附近向内凹陷量的增加, 翼型的亚声速升力性能得到改善, 这与流动机理分析中指的亚声速大迎角下要增大薄翼型升力应加强下表面设计的结论一致, 但却损失了超声速阻力性能.
|
|
下载CSV 表 4 两目标Pareto解集上的典型翼型及其气动性能对比 Tab.4 Typical airfoils on Pareto front and their aerodynamic performance |
图 21为翼型1与翼型4的亚声速压力分布与超声速的压力云图对比.由图 21(a)可见, 两个翼型亚声速升力的差异主要来自下表面压力分布的不同. 图 21(b)为翼型1(左)与翼型4(右)在超声速下的压力云图对比.翼型1的弯度很小, 上、下表面平滑, 除了头部有一道脱体激波之外, 翼面上没有强激波产生, 因此其超声速阻力很小.而翼型4是Pareto解集中亚声速升力系数最大的翼型, 结合上文分析, 翼型4通过前、后缘下弯, 下表面前段向内凹陷增加亚声速升力, 使翼型表面局部的变化较大.在超声速状态下翼型上表面前缘附近形成一道强激波, 极大地增加了翼型的超声速阻力.
|
| 图 21 翼型1与翼型4的亚声速压力分布与跨声速压力云图对比 Fig.21 Comparisons of pressure coefficient distributions and pressure contours of airfoil 1 and 4 |
四边形翼型是严格对称的翼型, 表面由简单的直线相连, 因此四边形翼型实际上是Pareto前沿上超声速阻力最小的翼型.对比表 3中由线性加权法得到的优化翼型(Opt-4)与四边形翼型的各设计指标的数据, 优化翼型的性能似乎比不上四边形翼型.观察图 20中Pareto前沿, 四边形翼型是Pareto前沿上超声速阻力最小的翼型, 相比于四边形翼型, Opt-4翼型处于Pareto前沿左下角近似水平段斜率即将快速增大的位置, Opt-4翼型在超声速阻力增加不大的情况下, 其亚声速阻力得到较大提升. Opt-4翼型在超声速阻力与亚声速升力之间做了非常理想的权衡.以Opt-4翼型为基准翼型, 采用Pareto解集法进行宽速域翼型优化设计, 优化结果表明, Pareto前沿并没有推进, 只是由Opt-4翼型展开为Pareto解集中的众多最优翼型, 说明采用线性加权法进一步改进Opt-4翼型的空间不会太大.
4.3 基于Pareto解集的三目标优化问题在4.2节中采用Pareto解集法处理两目标优化问题, 分别将翼型的亚声速升力与超声速阻力作为优化目标进行宽速域翼型优化设计.优化得到的Pareto解集中包含37个翼型, 实际工程应用中可以根据需要选择合适的翼型.通过分析Pareto解集上翼型的外形与对应的性能, 对亚声速升力指标与超声速阻力指标博弈的空气动力学原理有了更深的理解.
为了进一步明确3个设计指标间的博弈关系, 并充分挖掘翼型优化的潜能, 本节采用Pareto解集法, 将亚声速升力、超声速阻力、高超声速升阻比3个指标均作为优化目标, 进行高超声速宽速域翼型优化设计研究.优化问题的数学模型为
| $ \begin{array}{l} \max \;\;\;\;\;{C_{1,{M_a} = 0,3}}\\ \min \;\;\;\;\;{C_{d,{M_a} = 1.2}}\\ \max \;\;\;\;\;{\left( {{C_1}/{C_{\rm{d}}}} \right)_{{M_a} = 6.0}}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;\;\;{r_{{\rm{le}}}} = 0.00125c\\ \;\;\;\;\;\;\;\;\;\;\;t > 0.99{t_0} \end{array} $ |
翼型外形的变形采用8阶CST方法.在开始优化前, 使用Latin超立方抽样生成40个初始样本点. 图 22为优化过程中加点400次的优化历程.
|
| 图 22 三目标优化的加点历程 Fig.22 Process of adding new sample points by three-objective optimization |
图 23为经过440次CFD计算后得到的Pareto前沿面, Pareto解集中共有181个翼型.图中由D指向B的蓝色箭头的是Pareto前沿面的推进方向.在图 23中蓝色的四面体块代表的是四边形翼型, 紫色的小球体代表Opt-4翼型.从总体趋势来看:由A到C, 亚声速升力指标的大幅提升必然导致超声速阻力与高超声速升阻比的恶化, 而跨声速升力指标与高超声速升阻比指标的改善在一定程度上是可以统一的.这说明翼型的亚声速升力与超声速阻力以及高超声速升阻比之间的矛盾是主要的.在决策使用哪个翼型时, 设计者应首先确定翼型所需的亚声速升力, 之后再在跨声速阻力与高超声速升阻比间作出权衡.
|
| 图 23 三目标优化得到的Pareto前沿 Fig.23 Pareto front obtained by three-objective optimization |
图 24是从三目标Pareto解集中选取的3个翼型, 分别是:亚声速升力系数最大的翼型、超声速阻力系数最小的翼型、高超声速升阻比最大的翼型.亚声速升力最大的翼型最大厚度之前的部分几乎弯成弓形, 而且它与跨声速阻力最小翼型与高超声速升阻比翼型的外形差异较大, 后两者的外形差异较小, 这也说明了亚声速升力指标与其他两个指标存在的矛盾是主要的.
|
| 图 24 三目标Pareto解集上的典型翼型 Fig.24 Typical airfoils on Pareto front |
亚声速升力最大的翼型由于弯度太大, 明显会使超声速及高超声速阻力过大.
超声速阻力最小的翼型与高超声速升阻比最大的翼型的外形差异在于:超声速阻力最小的翼型很对称, 表面近似成直线, 外形与四边形翼型类似; 而高超声速升阻比最大的翼下表面前段微向内凹.这是由于高超声速状态下, 翼型主要靠下表面压缩产生升力, 并且翼型的乘波特性逐渐体现出来, 通过下表面前段微向内凹, 在翼型阻力增加不大的情况下使附体激波更好地与下表面前段贴合, 从而增加了翼型的升力, 相应地其升阻比也增大了.
图 25是三目标Pareto解集中高超声速升阻比最大的翼型、跨声速阻力最小的翼型和亚声速升力最大的翼型在亚声速与高超声速状态下压力分布对比.与之前的分析一致, 亚声速时升力差异主要来源于下表面压力分布的不同; 而在高超声速情况下, 下表面压缩是升力的主要来源, 高超声速升阻比最大的翼型利用了乘波效应.对比3个翼型在高超声速状态下的压力分布, 高超声速升阻比最大的翼型通过下表面前段的乘波外形设计获得了更大的升力. 图 26与图 27分别是跨声速阻力最小翼型与高超声速升阻比最大的翼型在超声速与高超声速状态下的压力云图对比.在超声速情况下, 高超声速升阻比最大的翼型由于上、下表面前段微向内凹, 导致其压差阻力大于超声速阻力最小的翼型.
|
| 图 25 三目标Pareto解集上3个典型翼型亚声速与高超声速下的压力分布对比 Fig.25 Comparisons of pressure coefficient distributions |
|
| 图 26 Pareto前沿面上翼型的超声速压力云图对比(Ma=1.2, Re=8.16×107, α=0°) Fig.26 Comparison of pressure contours at supersonic state |
|
| 图 27 Pareto前沿面上翼型的高超声速压力云图对比(Ma=6, Re=2.47×107, α=4°) Fig.27 Comparison of pressure contours at hypersonic state |
本文针对宽速域翼型优化问题, 进行了薄翼型在不同速域下的流动机理分析.然后分别采用线性加权法与基于Kriging模型的Pareto解集法开展了宽速域翼型优化设计研究, 得到一系列优化翼型.本文的一些研究结论如下:
(1) 设计了一种能很好权衡各速域气动性能的新翼型(Opt-4).相比常规的四边形翼型, 新翼型在超声速阻力与高超声速升阻比性能下降不多的前提下, 使亚声速升力得到了较大的提升, 具有更好的宽速域特性;
(2) 通过Pareto解集法优化设计得到两个翼型簇.这些翼型是根据各设计目标的重要程度作不同权衡时所得优化翼型的集合, 在工程应用可根据实际情况选择最合乎要求的翼型;
(3) 相比四边形翼型, 宽速域翼型优化应寻求在超声速阻力与高超声速升阻比增加不大的前提下尽可能增加翼型的亚声速升力, 优化翼型要全面超越四边形翼型是非常困难的;
(4) 采用Pareto解集法开展宽速域翼型优化设计研究, 不仅增加了解的多样性, 还揭示了优化问题的空气动力学的原理.实际处理多目标优化问题可以先用线性加权法找到Pareto前沿上的一个解, 再利用Pareto解集法扩展优化解的多样性, 并获得多目标优化问题更多的信息.
本文针对高超声速飞行器宽速域翼型进行了研究.揭示了部分宽速域流动机理, 然而由于高超声速飞行器本身的气动特征, 导致三维效应十分明显, 只从翼型进行研究是远远不够的.未来的研究中将进一步开展三维机翼宽速域流动机理研究, 并直接进行三维机翼优化设计研究.
| [1] |
Ueno A, Suzuki K. CFD-based shape optimization of hypersonic vehicles considering transonic aerodynamic performance[R]. AIAA 2008-288, 2008.
|
| [2] |
曹长强, 蔡晋生, 段焰辉. 超声速翼型气动优化设计[J]. 航空学报, 2015, 36(12): 3774-3784. Cao C Q, Cai J S, Duan Y H. Aerodynamic design optimization of supersonic airfoils[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(12): 3774-3784. (in Chinese) |
| [3] |
孙祥程, 韩忠华, 柳斐, 等. 高超声速飞行器宽速域翼型/机翼设计与分析[J]. 航空学报, 2018, 39(6): 26-37. Sun X C, Han Z H, Liu F, et al. Design and analysis of hypersonic vehicle airfoil/wing at wide-range Mach numbers[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(6): 26-37. (in Chinese) |
| [4] |
刘俊.基于代理模型的高效气动优化设计方法及应用[D].西安: 西北工业大学, 2015. Liu J. Efficient surrogate-based optimization method and its application in aerodynamic design[D]. Xi'an: Northwestern Polytechnical University, 2015(in Chinese). |
| [5] |
Han Z H, He F, Song W P, et al. A preconditioned multigrid method for efficient simulation of three-dimensional compressible and incompressible flows[J]. Chinese Journal of Aeronautics, 2007, 20(4): 289-296. DOI:10.1016/S1000-9361(07)60046-6 |
| [6] |
Xie F T, Song W P, Han Z H. Numerical study of high-resolution scheme based on preconditioning method[J]. Journal of Aircraft, 2009, 46(2): 520-525. DOI:10.2514/1.37976 |
| [7] |
Wu M M, Han Z H, Nie H, et al. A transition prediction method for flow over airfoils based on high-order dynamic mode decomposition[J]. Chinese Journal of Aeronautics, 2019. DOI:10.1016/j.cja.2019.03.020 |
| [8] |
Han Z H. SurroOpt: A generic surrogate-based optimiza-tion code for aerodynamic and multidisciplinary design[C]. 30th Congress of the International Council of the Aeronautical Sciences, Daejeon, 2016.
|
| [9] |
Spalart P, Allmaras S. A one-equation turbulence model for aerodynamic flows[J]. La Recherche Aérospatiale, 1994, 439(1): 5-21. |
| [10] |
Kulfan B M, Bussoletti J E. ″Fundamental″ parameteric geometry representations for aircraft component shapes[R]. AIAA 2006-6948, 2006.
|
| [11] |
Kulfan B M. A universal parametric geometry representation method-″CST″[R]. AIAA 2007-62, 2007.
|
| [12] |
韩忠华. Kriging模型及代理优化算法研究进展[J]. 航空学报, 2016, 37(11): 3197-3225. Han Z H. Kriging surrogate model and its application to design optimization:A review of recent progress[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3197-3225. (in Chinese) |
| [13] |
Han Z H, Zhang Y, Song C X, et al. Weighted gradient-enhanced Kriging for high-dimensional surrogate mode-ling and design optimization[J]. AIAA Journal, 2017, 55(12): 4330-4346. DOI:10.2514/1.J055842 |
| [14] |
Liu J, Song W P, Han Z H, et al. Efficient aerodynamic shape optimization of transonic wings using a parallel infilling strategy and surrogate models[J]. Structural and Multidisciplinary Optimization, 2017, 55(3): 925-943. DOI:10.1007/s00158-016-1546-7 |
| [15] |
Han Z H, Göertz S, Zimmermann R. Improving variable-fidelity surrogate modeling via gradient-enhanced Kriging and a generalized hybrid bridge function[J]. Aerospace Science and Technology, 2013, 25(1): 177-189. DOI:10.1016/j.ast.2012.01.006 |
| [16] |
Han Z H, Göertz S. Hierarchical Kriging model for variable-fidelity surrogate modeling[J]. AIAA Journal, 2012, 50(9): 1885-1896. DOI:10.2514/1.J051354 |
| [17] |
Han Z H, Zimmermann R, Göertz S. Alternative cokri-ging method for variable-fidelity surrogate modeling[J]. AIAA Journal, 2012, 50(5): 1205-1210. DOI:10.2514/1.J051243 |
| [18] |
Han Z H, Zhang K S. Surrogate based optimization, Intec book, Real-World Application of Genetic Algorithm[Z]. 2012: 343-362.
|
| [19] |
Han Z H, Xu C Z, Zhang L, et al. Efficient aerodynamic shape optimization using variable-fidelity surrogate models and multilevel computational grids[J]. Chinese Journal of Aeronautics, 2019. DOI:10.1016/j.cja.2019.05.001 |
| [20] |
Zhang K S, Han Z H, Gao Z J, et al. Constraint aggregation for large number of constraints in wing surrogate-based optimization[J]. Structural and Multidisciplinary Optimization, 2019, 29(2): 421-438. |
| [21] |
Han Z H, Chen J, Zhang K S, et al. Aerodynamic shape optimization of natural-laminar-flow wing using surrogate-based approach[J]. AIAA Journal, 2018, 56(7): 2579-2593. DOI:10.2514/1.J056661 |
| [22] |
Zhang Y, Han Z H, Zhang K S. Variable-fidelity expected improvement method for efficient global optimization of expensive functions[J]. Structural and Multidisciplinary Optimization, 2018, 58(4): 1431-1451. DOI:10.1007/s00158-018-1971-x |
| [23] |
韩忠华, 张瑜, 许晨舟, 等. 基于代理模型的大型民机机翼气动优化设计[J]. 航空学报, 2019, 40(1): 155-170. Han Z H, Zhang Y, Xu C Z, et al. Aerodynamic optimization design of large civil aircraft wings using surrogate-based model[J]. Acta Aeronauticaet Astronautica Sinica, 2019, 40(1): 155-170. (in Chinese) |
| [24] |
何衍儒, 宋保维, 曹永辉. 基于Pareto最优的翼身融合水下滑翔机结构优化设计[J]. 水下无人系统学报, 2017, 25(3): 243-249. He Y R, Song B W, Cao Y H. Structure optimization design for underwater glider with blended-wing-body based on Pareto optimal solution[J]. Journal of Unmanned Undersea Systems, 2017, 25(3): 243-249. (in Chinese) |
| [25] |
Deb K, Pratap A, Agarwal S, et al. A fast and elitist multiobjective genetic algorithm:NSGA-Ⅱ[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2): 182-197. DOI:10.1109/4235.996017 |