0 引言

由于现代民机的摩阻占总阻力的40%~55%, 增大机翼表面的层流区域以期降低摩阻成为设计空气动力学近些年的研究热点^[1-3].目前, 自然层流机翼已经在高空长航时无人机^[4]、中小型公务机上得到了广泛的应用^[5-6], 但针对跨音速的大中型民机尚处于理论与试验阶段, 尤其是后掠效应与实现层流对机翼表面压力分布的要求相互矛盾, 更限制了自然层流技术在中大型民机上的实际应用^[7].尽管如此, 自然层流技术依然是提高未来民机性能的发展方向, 而对于大多数在跨声速状态巡航的民机来说, 对跨音速自然层流翼型的深入研究, 必然为真正制造出高性能的层流机翼提供一定的理论依据及技术支持.

维持层流的顺压区长度与激波位置是跨音速自然层流翼型设计的关键, 顺压有利于维持层流, 但过长的顺压区为翼型后缘的压力恢复带来较强的激波, 因此在气动设计中, 需要兼顾这两个关键因素.本文旨在通过压力分布的反设计来协调激波的位置、强度以及顺压区, 并通过考察翼型在不同雷诺数下及不同马赫数下的气动性能, 为跨音速层流翼型的设计提供新方法和新思路.

反设计方法的优点是:在设计人员丰富气动设计经验的基础上, 针对设计要求提出合理压力分布或速度分布, 然后通过求解空气动力学反问题来确定对应的几何外形.近些年发展的特征正交分解(proper orthogonal descomposition, POD)方法是一种高效的数据处理方法, 已经广泛地应用到了诸多领域^[8-9], 它通过模态分析的思想提取到数据的主要特征, 用在翼型反设计的优势是在给定目标压力分布和对已知基准翼型形状修改而形成的翼型设计空间的条件下, 通过基于POD的数据重构方法(Gappy POD)就可以反设计得到与给定的目标压力分布相对应的最优翼型形状.

本文采用了一种基于差量采样解的GappyPOD翼型反设计方法, 将该方法成功应用在跨音速层流翼型的反设计中, 针对两种经典的压力分布形态, 进行了精确高效的反设计, 评估了压力分布形态对转捩发生抑制能力的影响.

1 基于差量采样解的Gappy POD方法

文献[10]详细介绍了POD方法, Gappy POD方法是POD方法用于已知系统中的填补某一缺失数据的元素时的变形, 通过在系统中获取合适的足够数量采样解(快照样本), 在采样解的基础上提取此系统的基模态, 采用Gappy POD方法^[11]可以对该系统中任一缺失数据的元素进行精确填补.

首先定义向量矩阵G及其对应的“标签”向量矩阵n.“标签”向量矩阵n用来识别向量矩阵G中的元素是否缺失.若不缺失n_i^k=1, 若缺失n_i^k=0, 其中n_i^k对应于向量矩阵G中第k个不完全数据向量中的第i个元素的“标签”.假设{φⁱ}_i=1^m是通过POD方法求得的已知采样集{uⁱ}_i=1^p的基向量^[12].g表示该系统中待填补的向量, 由已经得到的m个POD基来填补近似向量$ \mathit{\boldsymbol{\tilde g}}$.那么

$ \mathit{\boldsymbol{\tilde g}} = \mathop \sum \limits_{i = 1}^m {\alpha _i}{\mathit{\boldsymbol{\varphi }}^i}, $

(1)

并通过求解Ma=f线性方程组来确定α_i, 其中:M_ij=(φⁱ, φⁱ)_n; f_i=(g, φⁱ)_n.近似向量$\mathit{\boldsymbol{\tilde g}} $在求出α_i后由公式(1)求得.

本文研究了一种基于差量采样解的Gappy POD方法, 对原始采样解的结构重新定义, 重新定义采样解集合为Δu=[(u¹-u), (u²-u), …, (u^p-u)], 其中:u表示采样解的平均值; uⁱ|_i=1^p表示第i个采样解.通过上述方法求解Δu的基向量及系数, 那么最终修正向量可表示为$\mathit{\boldsymbol{\tilde g}} = \mathit{\boldsymbol{\bar u}}\mathop \sum \limits_{i = 1}^m {\alpha _i}{\mathit{\boldsymbol{\varphi }}^i} $.

2 翼型反设计方法

基于差量采样解Gappy POD方法构建翼型反设计方法分4步:1)根据反设计需要定义合理的设计空间, 利用计算流体力学(computational fluid dynamics, CFD)方法获得设计空间中样本的解, 以此构造快照矩阵; 2)对生成的快照矩阵按文献[10]中的方法求解基向量Φ, 并计算每个基向量的能量; 3)根据所求解问题的需要对低能量模态进行适当截断, 形成对应设计空间的降阶子空间; 4)利用降阶子空间对目标缺失数据进行数据重构.

2.1 CFD数值模拟技术

采用有限体积法求解雷诺平均Navier-Stokes方程, 其积分形式为$\frac{\partial }{{\partial t}}\iiint_V {Q{\text{d}}V} + \iint_S f\cdot{n{\text{d}}S} = 0 $, 在预测转捩时, 采用文献[12-13]提出的基于SST k-ω湍流模型中的完全基于流场当地变量的$ \gamma - \overline {{{\operatorname{Re} }_{\theta t}}} $转捩模型, $ \gamma - \overline {{{\operatorname{Re} }_{\theta t}}} $两方程转捩模型通过γ控制湍动能生成来实现当地转捩, 无量纲守恒形式的γ输运方程为

$ \begin{array}{l} \frac{{\partial \rho \gamma }}{{\partial t}} + \frac{{\partial \rho \gamma {u_i}}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_i}}}[(\mu + \frac{{{\mu _t}}}{{{\sigma _\gamma }}})\frac{{\partial \gamma }}{{\partial {x_i}}}] + {c_{a1}}\rho S{F_{{\rm{length}}}}\\ {F_{{\rm{onset}}}}^{{\rm{0}}{\rm{.5}}}{\gamma ^{0.5}}(1 - {c_{e1}}\gamma ) - {c_{a2}}\rho \mathit{\Omega }{F_{{\rm{turb}}}}\gamma ({c_{e2}}\gamma - 1). \end{array} $

引进变量$\overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} $, 边界层以外的$\overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} $通过经验关联函数Re_θt计算, 边界层之内的$\overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} $通过输运方程从边界层以外联立计算,

$ \begin{array}{l} \frac{{\partial \rho \overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} }}{{\partial t}} + \frac{{\partial \rho \overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} {u_i}}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_i}}}[{\sigma _{\theta t}}(\mu + {\mu _t})\frac{{\partial \overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} }}{{\partial {x_i}}}]\\ + {c_{\theta t}}\frac{\rho }{{{t_{{\rm{scale}}}}}}(\overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} - \overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} )(1 - {F_{\theta t}}). \end{array} $

2.2 快照获取方法

本文的快照矩阵是由翼型数据和压力分布构成的, 每一个翼型数据和对应压力分布组成快照矩阵的列向量, 通过在基准翼型的基础上加入扰动量形成多个快照, 依次构成翼型反设计的设计空间.

对于基准翼型的扰动是通过ΔCST翼型参数化方法实现的, 对于在采用类函数变化(class function/shape function transformation, CST)参数化^[14]翼型时, 定义$ \zeta \left( \psi \right)' = {\psi ^{0.5}}\left( {1 - \psi } \right)\mathop \sum \limits_{i = 0}^n ({A_i} + \Delta {A_i}){K_i}{\psi ^i}{\left( {1 - \psi } \right)^{n - i}} + \psi {\zeta _T}$, 即在初始设计变量基础上引入一个扰动量ΔA_i, 将上式展开得

$ \zeta \left( \psi \right)' = \zeta \left( \psi \right) + {\psi ^{0.5}}\left( {1 - \psi } \right)\mathop \sum \limits_{i = 0}^n \Delta {A_i}{K_i}{\psi ^i}{\left( {1 - \psi } \right)^{n - i}}, $

(2)

引入定义

$ \Delta \zeta \left( \psi \right) = {\psi ^{0.5}}\left( {1 - \psi } \right)\mathop \sum \limits_{i = 0}^n \Delta {A_i}{K_i}{\psi ^i}{\left( {1 - \psi } \right)^{n - i}}, $

(3)

因此, 由式(2)和(3)得到ζ(ψ)'=ζ(ψ)+Δζ(ψ).

2.3 反设计方法

本文的反设计过程为:1)对快照向量进行POD分解, 得到基向量, 用基向量逼近目标压力分布得到反设计翼型; 2)以反设计得到的翼型为初始翼型, 再进行扰动得到扰动翼型, 加上反设计得到翼型, 进行CFD分析获得快照向量, 将这些快照向量加入原始快照向量库, 如不满足条件, 重复第1、2步; 3)满足收敛条件, 结束.

3 跨音速自然层流翼型反设计

本节主要通过两种不同的压力分布形态反设计得到气动外形, 进而在不同雷诺数及马赫数下分析压力分布形态对转捩发生的抑制能力, 同时研究其与阻力发散特性的协调性.

ΔCST参数化上下表面各15个设计变量, 共30个设计变量, 这30个设计变量依次取指定的变量上下限, 共产生60个新翼型, 加上基准翼型共61个.采用文中构建的$\gamma - \overline {{{{\mathop{\rm Re}\nolimits} }_{\theta t}}} $转捩模型进行边界层转捩预测和气动特性分析评估, 计算网格为8万.

对生成的61个翼型形状在如下状态进行CFD计算, 获得采样解:Ma_∞=0.74;C_L=0.53;Re=6.0×10⁶; 飞行环境湍流度Tu=0.2%;湍流黏性比为10.

图 1为目标压力分布反设计结果, 可以看出来, 本文的反设计方法对于不同的压力分布都可以得到令人满意的结果, 记此次反设计得到的翼型为AirFoil_LSC2.图 2为初始基准翼型与两个反设计翼型的比较.选取一个设计状态为Ma_∞=0.74, C_L=0.575, Re=2.0×10⁷, 阻力发散特性较好的传统超临界翼型与上述跨音速自然层流翼型进行对比, 标号分别为:超临界翼型(Airfoil-SC)、Airfoil-LSC1、Airfoil-LSC2, 如图 3所示.

数值模拟3个翼型在雷诺数分别为600万、800万、900万的气动性能.图 4~8给出了各个翼型之间压力分布和阻力发散特性以及层流区随马赫数、雷诺数大小的变化对比.可以看出在雷诺数低于800万时, 随着马赫数的增大, 激波位置后移, 同时上表面转捩位置后移.Airfoil-SC翼型从压力分布上看, 翼型上表面没有明显的顺压梯度, 但依然可以依靠上表面的弱激波维持40%左右的层流区域, 并且层流区域是随着激波位置不断向翼型后缘移动.对于层流理念设计的Airfoil-LSC1和Airfoil-LSC2来说, 翼型上表面层流区域的变化基本与激波位置变化相一致; 从下表面的转捩位置与压力分布的比较来看, 翼型下表面的层流区域的维持是需要压力分布顺压区域的, 顺压区域越长, 层流区越长.从不同雷诺数下的阻力发散特性比较可以看出, 雷诺数在800万以下, 层流理念设计的跨音速层流翼型在层流区大小、阻力系数大小方面有着无可替代的优势, 但阻力发散特性要低于传统超临界翼型, 这是因为自然层流翼型由于顺压梯度的存在, 随着马赫数的增大, 激波不断增强后移, 为激波后翼型压力恢复带来过大的逆压梯度, 因此带来较大的激波阻力(波阻), 过大的波阻削弱了层流摩阻降低的效益, 总阻力的优化在高马赫数时不是特别明显.

从800万雷诺数及900万雷诺数翼型气动特性的比较可以看出, Airfoil_LSC1、Airfoil_LSC2翼型上表面压力分布在20%~40%弦长区域出现了微小的波动(逆压梯度), 这个微小的波动导致了转捩的发生; 下表面的层流区域的变化基本与顺压区吻合, 层流理念设计的翼型, 可以维持层流区域到60%弦长的位置, 微小的压力分布波动并未影响层流区域的维持.从800万雷诺数之下的翼型压力分布也可以观察到这个小波动, 但在800万雷诺数之下, 这个微小的逆压梯度并未造成转捩的发生, 转捩的发生与激波有很大的关系, 但超过800万雷诺数后, 翼型上表面的转捩对顺压梯度要求变得更加苛刻, 需要对压力分布进行更细致的设计.从800万及900万雷诺数的阻力发散特性可以看出, 在低马赫数时, 层流设计的翼型依然能带来气动效率的提高, 但在马赫数大于0.75时, 由于翼型前部顺压梯度带来过大压力恢复梯度, 陡增了波阻, 降低了翼型的阻力发散特性.

综上, 跨音速自然层流设计的翼型典型的问题是在中高雷诺数下阻力发散特性明显下降, 对顺压区域要求更加严苛; 与传统的超临界翼型对比表明, 在中低等雷诺数下顺压设计的层流翼型阻力特性明显, 阻力发散特性不是很差, 转捩区域与激波位置有很大关系, 微小的逆压梯度并未造成转捩的发生; 在中高等雷诺数下, 抑制转捩的发生需要更严苛的顺压梯度, 阻力发散特性变得更差, 在这种情况下, 可以考虑采用严苛的弱顺压梯度、短顺压长度和及早的弱激波恢复等设计理念进行设计.

4 结论

本文基于差量采样解GappyPOD的翼型反设计方法, 对跨声速自然层流翼型进行了反设计, 通过分析设计结果可得到以下结论:

1) 层流理念设计的跨音速翼型在800万雷诺数之下, 层流区域的维持与激波位置有很大的关系, 对顺压区域的要求不太明显, 小的逆压梯度并未造成翼型上下表面转捩的发生; 马赫数小于0.75时, 气动效率明显优于传统的超临界翼型; 马赫数大于等于0.75时, 由于波阻的陡增, 降低了层流带来的气动效率的增益, 因此若考虑阻力发散特性, 需要控制翼型上表面的层流区域, 牺牲设计点的气动效率, 以带来相对较好的阻力发散特性.

2) 层流理念设计的跨音速翼型在雷诺数大于800万时, 层流区域的维持对顺压梯度要求变得严苛, 微小的逆压梯度会造成翼型上表面转捩的发生, 需要对逆压梯度进行更细致的设计, 阻力特性随着马赫数的增大变得更差, 在这种情况下, 需要考虑采用严苛的弱顺压梯度、短顺压长度和及早的弱激波恢复等设计理念进行设计.