某层流验证机层流翼段气动改进设计

魏自言，李杰，张恒，杨钊

西北工业大学 航空学院，西安 710072

降低飞机的排放及化石燃料的消耗一直是各国的研究工作重点，美国NASA(National Aeronautics and Space Administration)根据“新一代航空运输系统”计划，将亚声速固定翼民用运输机的发展划分成了3个阶段，其中第2阶段即2020—2025年计划将飞机燃油消耗水平相比波音777或GE90燃油消耗降低40%，而第3阶段计划进一步降低至少70%[1]。对于现代典型民用宽体客机如B-747等，其摩擦阻力约占总阻力的40%，而机翼、尾翼及机身所产生的摩阻占比分别为38%、25%、20%，而对于较小的民用客机如MD-80等，其摩阻占总阻力的比例更大[2]。Schrauf指出若长航程民用客机减少3%的摩阻，则每年每架飞机能节省约15 000美元的燃油[3]。因此，减少飞机的摩阻是降低飞机消耗、排放的重要途径。

由于层流边界层所产生的摩阻相比湍流边界层要小，因此维持或扩大机身、机翼表面的层流区来减少湍流边界层所带来的摩阻是减少飞机摩阻的一个重要、直接的手段，也是各国研究客机减阻的主要研究方向，就是使用自然层流控制(Nature Laminar Flow Control，NLFC)技术，即通过对翼型本身进行优化设计，利用顺压梯度分布延迟转捩，使得其能够维持较大的自然层流区域，称为自然层流翼型设计。相比于依靠吹吸气或鼓包的手段改变局部流动获得层流区的层流流动控制(Laminar Flow Control,LFC)技术，NLFC技术更加根本且可靠，仅管使用NLFC技术需要在飞机设计初期就加以考虑，但可以减少使用LFC技术带来的复杂吹吸气装置及缝道结构[4]。

仅管自然层流翼型在小迎角状态相比传统翼型可以获得更大的升阻比，但是层流翼型在低速大迎角状态下较传统翼型更容易产生流动分离甚至失速。因此，由于层流翼型的失速特性使得其目前在实际工程中没有得到大范围的应用[5]。

对于分离流动地捕捉一直是计算流体力学的重要研究领域。分离流动种类繁多、形态复杂，而目前工程中常用的雷诺平均(Reynolds Averaged Navier-Stokes, RANS)方法由于经验性较强，对复杂分离流动的模拟能力较弱，难以为大迎角失速分离问题提供较为满意的结果。而较为先进的湍流模型如雷诺应力输运模型(Reynolds Stress Model, RSM)其直接预测雷诺应力张量的每个分量，原理上能够更好地预测流动分离过程，但需要额外求解和雷诺应力相关的6个输运方程，引入了大量额外未知量，导致其鲁棒性较差[6]。

近年来越来越多的学者使用大涡模拟(Large Eddy Simulation, LES)、直接数值模拟(Direct Numerical Simulation, DNS)方法对分离流动进行研究，这2种方法在分离流动的捕捉上较RANS方法有较好的普适性且流动细节地捕捉也具有较好的精度，但LES与DNS计算所需的网格规模与雷诺数(Re)相关，其关系约为Nxyz∝Re13/7、Re9/4，因此目前的研究主要集中在低雷诺数问题计算中。而对于在工程上、尤其在航空领域动辄几百上千万雷诺数的实际问题，若要使用这2类方法进行计算，并得到一个较为可靠的计算结果，其计算的时间成本与计算资源需求十分巨大[7]。

为了保证获得较好分离特征的同时节约计算成本，在近壁面区域使用RANS方法，而在远离壁面的区域使用LES方法的RANS/LES混合方法逐渐受到重视。NASA在报告中指出RANS/LES方法将在学术范畴和工程问题中都将获得广泛的应用[8]。其中，DES (Detached Eddy Simulation)类方法自1997年提出至今，经历了推广、研究和改进等多个发展阶段，到目前为止，DES方法已经衍生出了Delay Detached Eddy Simulation(DDES)、Improved-DDES(IDDES)等改进方法，并成功应用于各类研究和实际问题之中。仅管目前对于DES类方法亚格子模型的构造尚存在多种方法，但大部分学者所提出新的亚格子模型构造形式，其目的主要是为了提高DES类方法在复杂流场中流动细节的捕捉精度，并没有对模型进行本质上的修改，因此DES类方法已经属于较为完备的一类方法[9]。

为了对高亚声速层流翼型气动性能进行深入研究及改进设计，本文针对中国航空工业集团第一飞机设计研究院专门特殊设计的某层流验证机进行了IDDES计算及设计修改。该层流验证机中央层流段为某高亚声速层流翼型。仅管该层流验证机在高亚声速状态下具有良好的气动性能，但根据后续全机风洞实验结果表明，中央层流翼段在低速大迎角状态时失速较早，降低了全机的可用升力系数，说明该高亚声速层流翼型在低速大迎角时较外翼段失速较早，因此推迟该高亚声速层流翼型的低速失速迎角十分重要，是提高层流验证机低速气动性能及未来对高亚声速层流飞机设计的关键。本文首先介绍了层流翼型设计的难点及所使用数值方法当前的发展现状。之后,对剪切层自适应 (Shear Layer Adaptive，SLA) IDDES方法的构造进行了介绍。然后,主要对标准翼型验证算例及实验模型进行计算验证。最后通过SLA-IDDES方法对风洞模型的中央层流翼段进行计算分析,并针对层流翼段低速气动性能进行了翼型修改设计。

1 SLA-IDDES方法

1.1 IDDES方法

本文中所用到的SLA-IDDES是由基于k-ωSST(Shear Sress Transport)模型的DDES方法改进得到的。Menter和Kuntz所提出的k-ωSST模型结合了k-ε与k-ω模型，其在边界层和自由剪切层区域使用k-ω模型，而在远离壁面的全湍区域使用k-ε模型，在k-ωSST模型中采用和Wilcoxk-ω一致的k方程和雷诺应力计算方式，而ε方程通过ε=kω转换后与ω方程融合[10]。

k-ωSST 模型守恒形式的k方程与ω方程如下：

(1)

(2)

DDES方法对k-ωSST模型中的k方程进行了修改，引入了特征长度lhyb。

(3)

式中：τij与Sij的计算方式为

DDES方法中的特征长度lhyb定义为

(4)

(5)

式(4)中亚格子尺度Δmax的定义为Δmax=max(Δx,Δy,Δz)，其中Δx、Δy、Δz代表当地网格在3个方向的网格尺度，当网格与壁面距离较近时，DDES方法退化为k-ωSST模型。其中混合函数FSST与文献[10]中混合函数F2的构造一致。

Menter和Kuntz使用k-ωSST中的F1对CDES进行了重构[10]

(6)

当计算网格足够密时，理论上DDES方法可以用做Wall-Modelled LES(WMLES)方法，然而DDES方法在边界层内部产生了对数层不匹配的问题。Shur等提出了IDDES方法，通过引入新的特征长度及亚格子尺度改善了此类问题[12]。

CDESΔ

(7)

Δ=min[max(Cwdw,CwΔmax,Δy),Δmax]

(8)

式中：Δy为壁面法向距离;dw为网格点到壁面的距离;Cw=0.15。混合经验函数定义为

(9)

(10)

(11)

fb=min{2exp(-9α2),1.0}

(12)

(13)

其中：Ωij为旋度张量；hmax为流场中当地网格3个方向最大值，经验系数Cdt1=8,Cdt2=3与经验函数fb保证在距离壁面0.5Δmax

另一个经验函数fe是为了缓解将DDES方法用于WMLES模型产生的LLM问题，其定义为

fe=fe2·max((fe1-1.0),0)

(14)

(15)

fe2=1.0-max(ft,fl)

(16)

(17)

(18)

(19)

式中：κ为Karman常数取值为0.41;常数Cl和Ct根据不同背景模型取值不同，对于两方程SST模型，Cl=5.0和Ct=1.78。

1.2 SLA亚格子尺度

如前所述，原始DES方法中所使用的亚格子尺度[13]较为保守，因此Shur等[14]借鉴了文献[13]中提出的亚格子尺度并进行了改进，对于一个六面体网格，亚格子尺度定义为

(20)

式中：

In=nω×rn

(21)

(22)

(23)

同时为了将SLA模型与IDDES方法结合，依据文献[15]中给出的建议，将亚格子尺度改为

Δ=min[max(Cwdw,CwΔmax,Δy),ΔSLA]

(24)

2 计算验证

为验证SLA-IDDES方法针对翼型小分离流动的捕捉精度，首先选取相对较为简单的三维A-Airfoil翼段进行模拟计算，验证低速大迎角状态下，翼段后缘区域分离位置捕捉精度。后续对层流验证机中央层流翼段进行低速与高速计算，并与实验进行对比，确保方法及实验准确、可靠。本文中所有计算使用课题组发展的基于多块结构网格求解器，无黏项使用五阶WENO格式，黏性项采用二阶中心差分，时间项采用双时间步隐式推进[16]。

2.1 A-Airfoil计算验证

由于A-Airfoil翼段试验数据丰富，在EUROVAL[17]、 ECARP[18]和LES-FOIL[19]等文献中已经做了大量研究，其中包含了分离位置较难捕捉的低速大迎角下翼型后缘小分离流动状态。本节使用粗、中、密3套O型网格对A-airfoil在迎角α=13.3°、Ma=0.13、Re=2.1×106的状态下进行计算验证，计算的来流湍流度为0.03%与风洞试验保持一致。网格节点参数如表1所示，中等数量网格图如图1所示。

图1 中等数量网格

表1 3套网格分布信息

为验证网格无关性，使用SLA-IDDES方法对3套网格均计算稳定且到相同时间步后，取无量纲时间跨度为15ΔT，对流场结果进行平均，其中ΔT=c/U∞。将3套网格计算所得表面压力系数(Cp)与试验进行对比，Cp对比如图2所示。

从图2可以看出，粗网格在前缘负压峰值位置与实验相差较大，其翼型中段与后缘位置Cp的趋势基本相吻合，而中等网格与密网格仅管也在负压峰值位置略有偏差，但趋势及数值上与实验更为接近。而密网格与中等网格计算结果接近，表明中等网格已经满足求解精度的需要。对于前缘负压峰值的偏差，有可能是由于翼型在低速大迎角状态下仍然能保持一定的层流区，且根据文献[20]中所指出的前缘区域存在层流分离泡的现象，由于前缘区域没有明显分离产生，RANS/LES方法在该区域采用RANS求解，而使用全湍计算很难捕捉层流区域的层流分离泡，因此前缘负压峰值与实验结果产生了一定偏差，但偏差相对较小，且本文主要关注翼型的分离流动，对分离流动及翼型改进设计进行研究，因此密网格计算结果与试验结果在前缘区域压力分布存在的偏差尚可接受。

图2 3套网格表面Cp与实验对比

依据试验[17]，A-Airfoil约在弦长82%处分离起始，而在密网格上使用SLA-IDDES模型进行计算得到的分离位置在弦长80.5%处，分离点附近速度型如图3所示，其中C为当地声速。计算所得分离位置与试验相比略微提前但相差不大，表明该方法对翼型后缘小分离流动的捕捉精度较好。

图3 密网格分离点位置速度型

2.2 试验模型计算验证

由于RANS方法对RANS/LES计算的准确性有很重要的影响，尤其部分边界层区域的流动仍是由RANS方法决定，为了保证后续DES计算的正确性，首先使用满足本文所提到的SLA-IDDES方法求解精度的网格，通过RANS方法对实验模型进行计算验证，以确保后续计算分析的准确性。

风洞试验模型如图4所示，低速试验在FL-8风洞进行，试验中采用1∶3.25的缩比模型，试验雷诺数大概在150万量级，马赫数Ma=0.2。该低速试验中采用腹撑式天平进行测力，并通过图4(a)中所示的对称天平测量试验来定量扣除天平装置气动外形对飞机流场和气动力的影响。高速试验在FL-2进行，风洞试验中采用1∶7 的全金属通气缩比模型，试验雷诺数在300万量级。

图4 风洞试验及模型

全机网格如图5所示，整体网格数量约为3 000万。飞行雷诺数为4.93×106，采用k-ωSST模型对缩比前后的模型分别进行计算。计算结果如图6所示，图中高雷诺数代表飞行雷诺数条件下的计算结果，即未缩比模型的计算结果，低雷诺数则代表风洞试验雷诺数条件下的计算结果，即对缩比模型的计算结果。

图5 试验模型网格

对比3条曲线可以看出，计算结果与试验结果在线性段内吻合较好，而非线性段整体趋势较为近似，表明所使用的网格满足计算精度。而高雷诺数的计算与低雷诺数计算结果在线性段内相差不大，在非线性段则有所差异。表明该飞机升力系数在线性段范围内对雷诺数不敏感。

为了实现高低速协调设计，所使用的网格需要进行试验模型高速计算以对比试验结果。该飞机设计巡航马赫数Ma=0.7，计算雷诺数约为Re=3×106，与风洞试验保持一致，采用同样的网格对该状态进行计算并与试验结果进行对比，对比结果如图7所示。

从图6及图7的对比结果可以看出，所使用的网格计算得到的结果在高低速条件下对宏观气动力的捕捉精度较好，与试验结果基本吻合。因此网格基本满足了后续中央层流翼段翼型改进设计的精度要求。

图6 Ma=0.2时计算与试验结果

图7 Ma=0.7计算与试验结果

3 层流翼型计算及改进设计

3.1 原始构型计算分析

通过第2节的验证，表明所使用的SLA-IDDES方法能够捕捉到较难捕捉的小分离流动，同时也验证了所使用的网格能够准确地捕捉宏观气动性能。仅管如此，为对层流翼型进行高低速条件下的气动性能改进设计，需要对层流翼型的流场进行高精度计算，捕捉翼型在较大迎角状态下空间流动的特征，以针对性的对翼型进行改进设计。

为捕捉中央翼段的流动细节，同时优化计算量、加速计算，仅保留部分机身及半展长中央翼段，所截取的模型如图8所示。

图8 计算模型示意

利用2.2节所使用的网格，仅保留关注部分，并将网格整体进行适当加密，使得RANS/LES计算时LES区可以捕捉到足够多的细节用于分析，加密后的网格如图9所示，网格尺度在层流翼段范围内维持x-y方向长度基本一致，z方向则根据IDDES方法的计算需求进行加密，保证网格能够合理激活WMLES，以在边界层附近获得较为精确且细节丰富的复杂流动[21]，尤其是层流翼型容易触发前缘分离，分离后所脱出的尾迹可能对翼段上表面有较大的影响。而后缘之后的尾迹流动对中央翼段影响较小，因此对后缘没有进行过多的加密，以节省计算资源。加密后模型的网格数量约为2 200万。

图9 层流翼段网格

根据第2节的试验模型验证算例结果可以看出，其升力系数从α=8°开始进入非线性段，因此在α=8°之前无论是外翼段还是中央翼段理论上应该不存在明显分离，故对于中央层流翼段的IDDES计算分析可以从α=8°之后进行，因此仅对中央层流翼段在Ma=0.2，8°～14°范围内进行IDDES计算，平均流场采用2.1节中A-Airfoil验证算例流场平均方法得到。

对中央层流翼段在Ma=0.2，迎角8°～14°范围内的表面及空间流动进行分析。图10为中央层流翼段在α=8°、10°、12°及14°下平均后流线及表面Cp云图。可以明显看出，层流翼段在α=8°时，其分离主要在翼段前缘与机身结合的位置，其产生的原因主要是由于机身边界层及翼段边界层在交界处的相互干扰所导致[22]，而非翼型本身所导致的分离。而且层流翼段后缘区域有明显的横流，其主要是由于翼身结合处流动挤压所导致的，与层流翼型本身无关。这个挤压会随着迎角逐渐增大、前缘分离流动的产生而逐渐弱化[23]。

图10 原始层流翼段平均表面Cp及流线

为了更清晰地分析层流翼段前缘区域的流动特征随迎角增大的变化规律，采用Q准则方法对SLA-IDDES方法计算所得到的同一瞬时时刻流场中的涡进行显示。图11为层流翼段在α=8°、10°及12°时Q=1的云图。从图11的Q云图可以看出，中央层流翼段8°迎角时，除翼身结合处产生的马蹄涡在前缘区域诱导产生了流动分离外，整个表面基本维持附着流动的状态，仅后缘处存在分离趋势。在α=10°时，前缘处则发生了明显的分离流动，但分离后流动出现再附，因此除前缘部分，中央层流翼段的流动基本为附着流，且后缘也不存在分离趋势。并且α=12°时，其流动特征与10°类似，分离区域仍维持在前缘区域。

图11 原始模型不同迎角下瞬时空间等值面云图(Q=1)

因此，对于该层流翼段，其低速状态下迎角由8°逐渐增大至14°后，由仅有翼段与机身交界处存在明显小分离，逐渐发展为整个前缘部分均存在明显分离，最终层流翼段上表面有接近一半的区域产生了分离。对比α=10°与12°状态，从瞬时涡量场可以看出，仅管α=10°时翼段上表面前缘区域已经全部分离，但是翼身结合处前缘位置流动分离所扬起的涡较远离翼根处更高，表明翼根处分离更严重，而对于α=12°，其上表面翼身结合处前缘位置流动分离所扬起的涡与靠近中间位置的分离涡扬起高度基本接近，表明此时的分离已经完全由翼型本身的分离特性所主导，而非边界层掺混、干扰所致。

同样从Cp云图可以发现，α=8°～12°的前缘区域有明显的负压峰值区，这个负压峰值区与分离区域基本一致，而在α=14°状态下，前缘区域的负压峰值区较α=12°有所减小，压力分布变为了杂乱无章的分布形态，此时层流翼段上表面的表面及空间流动是有较为复杂的分离与回流。因此，该层流翼型的最大升力系数应处于α=12°附近。

3.2 翼型改进设计及计算分析

从上述分析可以得知，该层流翼段流动分离的主要原因存在于层流翼型的前缘区域，由于层流翼型普遍前缘半径较小，容易触发流动分离，而该层流翼段所使用的层流翼型已经较好地满足了高速飞行，改变最大厚度与弯度等参数会使得高速气动性能随之改变，尤其可能会带来阻力增大或巡航设计点改变，因此为了保证高速气动性能改变较小，同时推迟低速大迎角状态下的前缘分离，修改层流翼型的前缘区域形状是较为合适的途径。

进而对层流翼型的前缘部分进行修改，增加了前缘半径，对翼型0%～4%站位之间的下表面进行略微修改，使得驻点位置前移，低速大迎角状态下前缘区域更不容易产生流动分离，其余部分与原层流翼型保持一致，翼型修改前后对比如图12 所示。

图12 改进前后翼型对比图

同样对中央层流翼段局部模型在α=8°及14°进行SLA-IDDES计算，以分析修改后空间流动的变化，计算条件仍与3.1节保持一致。图13为修改后模型的平均表面流线及Cp云图。

图13 改进后层流翼段平均表面Cp及流线

从α=8°～14°的流线及表面Cp云图可以看出,相比原始模型α=8°时就在翼身结合处产生了前缘局部分离，改进后的模型在α=8°～10°状态，整体仍然维持附着流动，结合区域的前缘部分也没有产生明显分离，三维效应所导致的分离推迟到α=12°时出现。修改后的模型在α=14°时产生了明显的前缘分离，但分离区域相比原始模型在α=14°时较小。

图14为改进后模型同一瞬时时刻Q=1云图。改进后的中央层流翼段在α=12°时产生翼身结合处的前缘小分离，α=13°时发展为接近1/3的前缘区域存在分离，而α=14°时整个前缘区域及靠近中部区域产生流动分离。对比原始构型表面Cp及Q等值面可以看出，改进前后整体流动形态较为类似。翼身结合区前缘部分的小分离流动，从原始构型的α=8°推迟到改进后构型的α=12°。而且改进后的翼型在α=14°尽管其上表面仍然产生了较大范围的流动分离，但相比改进前的模型，其分离区有所较小，表明修改前缘半径对于大分离流动也能起到缓解流动分离的作用。

图14 改进模型不同迎角下瞬时空间等值面云图(Q=1)

为验证改进是否对有效，对层流验证机中央层流翼段的气动性能是否有所提升，将改进前后的中央翼段SLA-IDDES计算结果与2.2节中算例验证的原始构型中央层流翼段RANS结果进行提取对比。

图15为中央层流翼段改进前后CL对比图，可以看出，仅管RANS计算结果较SLA-IDDES计算结果整体趋势较为近似，但RANS方法计算结果偏乐观，其原因可能是RANS方法对小分离捕捉效果较差。从改进前后中央翼段的SLA-IDDES计算结果可以看出原始中央层流翼段在α=12°获得最大升力系数。而对层流翼型的改进设计不仅使得中央翼段大迎角下升力系数有所增加，最大升力系数对应的迎角向右移动了2°，而且CL线性段也相应延长了2°，表明对中央层流翼段前缘半径及下表面修型的改进设计较为有效。

图15 翼型改进前后中央翼段CL对比图

3.3 改进构型高速计算

仅管中央层流翼段的低速性能得到了一定提升，但是对于高速飞行状态时的影响仍然需要探究，仅管理论上翼型前缘及其附近下表面形状的小范围、小幅度调整修形只会稍微改变前缘下表面局部的压力分布特征，对中央层流翼段巡航状态的气动特性几乎无影响。因此仅需对全机在巡航状态下的设计点进行计算对比验证，设计点状态为Ma=0.7，α=1.5°，计算雷诺数仍保持300万量级，与2.2节高速部分计算验证保持一致。

图16为修改前后中央层流翼段对称面处的Cp云图，可以看出修改后的压力分布与原始翼型在设计点处的压力分布整体基本没有变化，仅在前缘区域极小范围内略有变化。图17为修改前后，模型在设计点的表面Cp及流线。可以看出表面流动基本一致，前缘的修改无论对中央层流翼段流动、Cp还是对机身附近的影响都没有明显地变化，则其在巡航设计点的气动特性无明显变化。因此，修改得到的层流翼型也满足了高速巡航的设计要求。

图16 设计点修改前后翼型Cp

图17 设计点压力分布及表面流线

4 结 论

本文通过基于5阶WENO格式的SLA-IDDES方法，对A-Airfoil翼型进行了小分离流动捕捉精度验证，同时对某层流验证机全机进行了RANS计算，利用SLA-IDDES方法对层流验证机的中央层流翼段进行了计算分析，并对中央层流翼段的翼型前缘区域所导致的流动分离进行了分析与翼型改进设计，推迟了前缘分离，增大了失速迎角，同时增加了最大升力系数。并且改进后的翼型在全机巡航马赫数设计点影响不大，基本与原构型保持一致。文章主要结论有以下3点：

1) 使用SLA-IDDES方法对大迎角下的分离流动捕捉较好，能够较为精细地捕捉到流动分离细节。

2) 对于层流翼型，其由于前缘半径较小，在低速大迎角状态时容易触发前缘流动分离，尤其在翼身结合区由于三维效应所引起的分离流动。通过修改前缘半径可以有效地延迟分离出发的迎角。

3) 层流翼型前缘分离对前缘半径及驻点位置参数敏感，增大前缘半径及前移驻点位置可明显抑制低速大迎角下层流翼型翼段的前缘流动分离，并且对高速巡航状态的压力分布影响较小，表面流动也无明显差异。