某层流验证机中央翼段失速分离特性优化及分析

唐松祥，李杰，张恒，牛笑天

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

减小飞行阻力一直是民用飞机设计中重要的一环，在目前绿色航空发展的大背景下，减阻这一目标更是被给予了更高的关注。目前，减阻具有较多的技术手段，其中层流机翼设计和控制技术是一种行之有效的方法。层流区的控制可有效提高机翼的升阻特性、降低燃油消耗率及有害气体排放。有研究指出当雷诺数保持相同时，层流边界层的摩擦阻力可达湍流边界层的1/10[1]。

作为结合自然层流技术和层流控制技术的混合层流控制技术，美国国家航空航天局(National Aeronautics and Space Administration，NASA)在早期已将其应用于波音757中[2-3]。近年来有关层流控制方面也得到了广泛的研究[4-10],而层流控制需要增加的额外装置在一定程度上阻碍了其应用。

自然层流技术则无需额外复杂的层流控制设备，其通过翼型控制压力梯度等流动参数达到增大层流区的效果。因而对层流翼型的设计成为自然层流技术的关键。国内外学者针对层流翼型的优化设计及边界层转捩开展了不少工作[11-16]。Han等[11]通过一种代理模型优化了跨声速条件下的自然层流翼型，并证实了该算法模型在跨声速这一特定条件下翼型优化的合理性。Xu等[12]通过试验及数值模拟探索了超临界层流翼型的边界层特性，并为翼型设计提供了一定思路。马晓永等[13]利用了NURBS(Non-Uniform Rational B-Spline)方法和CST(Classfunction/Shapefunction Transformation)方法对三维层流翼型进行了优化对比，结果表明NURBS方法在层流翼型的设计上具有更好的实用性。Cella等[14]提出了层流翼型的设计思路，并结合一跨声速层流翼型进行了相关设计及试验工作，说明了该设计思路的准确性。张彦军等[15]利用转捩预测技术与计算流体力学的手段设计超临界翼型，并利用试验总结了转捩特性与来流参数的相关规律，为超临界层流翼型设计提供了一定的参考。许联铭等[16]设计了前掠层流机翼，其设计方法使机翼在跨声速条件下获得了50%～60%的层流区长度。

上述自然层流相关文献的工作大多是基于计算流体力学开展的，说明模拟计算在层流翼型设计上具有足够的准确性。对于自然层流机翼的设计，合理的数值方法至关重要。常用的湍流模型如Spalart-Allmaras (S-A)模型、k-ω剪切应力运输(Shear Stress Transport，SST)模型等基于全湍的湍流模型均无法模拟层流转捩这一过程。目前工程上常用Menter等[17-18]提出的γ-Reθ转捩模型预测转捩。对于复杂分离流动的机理，由于传统的雷诺平均(Reynolds-Averaged Navier-Stokes，RANS)框架下流场在近壁面区域具有较大的湍流黏性，使分离流场预测的准确性受到了一定限制，不少学者甚至开始尝试将该模型与分离涡模拟(Detached Eddy Simulation，DES)方法结合[19-22]，以此对具有层流转捩的分离流场进行预测。但从文献[19-20]的研究结果看，基于雷诺平均RANS框架的γ-Reθ转捩模型在层流转捩的预测上具有良好的准确性。

本文研究对象为中国航空工业集团公司第一飞机设计研究院设计的某特殊构型层流机翼技术验证机，该无人机采用创新布局形式，将试验验证段与飞行平台较好地结合起来，整体布局方案具备性价比高、飞行试验效率高的特点，旨在验证不同飞行高度和马赫数时层流在真实大气飞行条件下的流动状态和减阻效果，是一项具有开创性意义的工作。针对该验证机低速条件下中央翼段的失速分离特性进行优化修形，旨在提高低速飞行状态下的失速特性。首先基于CFL3D求解器利用RANS方法对Aerospatiale-A翼型进行数值方法及网格拓扑在低速条件下转捩及分离现象准确性的验证，然后结合试验验证高速条件下层流验证机构型的计算准确性，最后针对层流验证机中央翼段翼型进行优化修形，以期提高其在低速、水平面高度飞行状态下的气动特性，并保证其在高速巡航状态下依旧具有较好的气动特性。

1 计算方法

湍流模型采用基于k-ωSST模型的四方程γ-Reθ转捩模型，该模型利用参数γ对标准两方程k-ωSST模型中湍动能方程的生成项与破坏项进行修正模拟转捩过程。标准形式的两方程k-ωSST模型湍动能方程为[23]

(1)

式中：ρ为密度;k为湍动能;t为时间;uj为各方向速度;xj为空间坐标;τij为切向应力张量;Sij为应变率张量;β*为模型常数，β*=0.09;ω为比湍流耗散率;μ为动力黏度;σk为模型变量;μt为涡黏系数。

四方程转捩模型湍动能方程具体形式为[24]

min(max(γeff,0.1),1.0)ρβ*ωk+

(2)

式中：γeff为间歇因子。

利用Aerospatiale-A翼型分析湍流模型和网格拓扑的准确性。

采用的计算程序CFL3D以格心格式的有限体积法离散控制方程。对于空间离散，无黏项采用Roe平均通量差分分裂格式(FDS)，单元界面上差值模板为三阶MUSCL(Monotonic Upstream-Centered Scheme for Conversation Laws)格式，采用min-mod限制器防止在间断处出现数值震荡。黏性项采用中心差分格式。时间推进方式采用近似因子分解(AF)隐式时间推进算法。

2 Aerospatiale-A翼型验证算例

2.1 计算模型

利用Aerospatiale-A翼型分析标准两方程k-ωSST模型和四方程γ-Reθ转捩模型对层流翼型气动参数的计算准确性。

选择计算迎角α为13.3°，在该迎角下上翼面同时具有转捩与后缘分离现象[25]。计算雷诺数Re=2.1×106，计算马赫数Ma=0.15，湍流强度Tu=0.2%，湍流黏性比为10，模型展向长度为1倍弦长。

表1为计算网格参数，计算网格示意图参见图1。翼型上翼面网格节点数约为下翼面的1.6倍。

表1 计算网格参数

图1 Aerospatiale-A翼型绕流计算网格

2.2 计算结果与讨论

对于Aerospatiale-A翼型，目前可用F1和F2试验结果[25]进行参照，F1与F2数据在不同风洞下得到。F2试验数据包含了流场中的部分速度测量数据，相较于F1试验数据更为全面，因此与计算结果作对比的试验数据均为F2试验数据。

图2(a)和图2(b)分别为采用表1中网格2的标准k-ωSST和不同网格下γ-Reθ转捩模型在非定常求解的翼面沿展向平均后的压力系数(Pressure Coefficient，Cp)和摩擦阻力系数(Friction Coefficient，Cf)。可知采用转捩模型时，计算网格2与计算网格3的Cp与Cf预测结果吻合良好，计算网格1在预测Cp时前缘峰值更大；而在Cf分布中，其对转捩位置的预测更加靠前，说明计算网格1无法准确预测翼面参数，计算网格2 已达到了良好的网格收敛性，后续将针对计算网格2计算结果进行分析。在Cp与Cf的预测中，3种γ-Reθ转捩模型均较标准k-ωSST模型获得了与试验结果匹配更好的结果。对于图2(a)所示的Cp，标准k-ωSST模型在前缘附近预测的Cp有较大程度的失真，在试验值Cp最大位置误差约8.7%，而γ-Reθ转捩模型在该处Cp误差约2.9%，且在整个翼段上标准k-ωSST模型预测的Cp均较γ-Reθ转捩模型的低，因此γ-Reθ转捩模型在Cp预测上具有更好的准确性。对于图2(b)所示的Cf，由于标准k-ωSST模型不具备捕捉转捩现象的能力，其在整个翼段上的Cf预测均有较大误差，而γ-Reθ转捩模型能准确地捕捉转捩现象，从而在整个翼段的Cf预测上具有更好的准确性。

图2 不同计算模型压力与表面摩阻系数分布

图3为采用计算网格2时不同弦向站位下标准k-ωSST和γ-Reθ转捩模型在非定常求解的翼面沿展向平均后的速度型分布，其中横坐标为当地速度Uw与来流速度Uinf的比值，纵坐标为翼面法向高度yw与翼型弦长c的比值。可看出标准k-ωSST模型预测的附面层均偏厚，相比较而言γ-Reθ转捩模型预测的速度型更能匹配F2试验结果。从图3(d)可看出在x/c=0.825时翼面开始出现分离，而γ-Reθ转捩模型的预测结果显示该模型在翼面分离的预测上具有一定的延迟。而从图3(e) 可知在充分发展的分离区中γ-Reθ转捩模型预测的速度型同样具有良好的准确性。

图3 上翼面不同站位的速度型分布

对于同时具有转捩和流场分离现象的Aerospatiale-A翼型二维翼段低速状态分析，在使用的网格拓扑结构下四方程γ-Reθ转捩模型在压力系数、摩擦阻力系数和边界层速度型分布的预测上均能较为合理地吻合F2试验结果，为γ-Reθ转捩模型在低速条件下翼型优化提供了有力的依据。

3 某传统翼型验证算例

3.1 风洞试验及计算模型

研究对象为中国航空工业集团公司第一飞机设计研究院设计的某特殊构型层流验证机。利用某传统翼型验证高速条件下层流验证机中央翼段数值模拟的准确性，该层流验证机在高速巡航时飞行马赫数为0.7，升力系数为0.19，迎角约1°～2°。试验采用的风洞为风雷FL-60风洞。采用1∶7 的缩比模型，试验压力为101.651 kPa，马赫数为0.7，雷诺数约8.0×106，来流湍流强度为0.4%。图4为层流转捩验证风洞试验照片。

图4 层流转捩验证风洞试验

数值模拟计算条件与试验条件相同，翼型与翼段数模如图5所示，中央翼段的计算网格拓扑与图1所示Aerospatiale-A翼型网格拓扑相似，网格节点数相同。

图5 中央翼段转捩验证模型

3.2 试验与计算结果对比

图6所示为迎角1°时风洞试验结果与计算Cf云图的对比。从图6(a)的试验热敏成像图可看出中央翼段的转捩区呈现出锯齿状特征，其中较亮的区域为层流区。而试验结果与图6(b)计算结果所呈现的光滑转捩区存在一定差异，该锯齿状形态考虑为用于试验的翼段表面涂层存在一定的如颗粒等使表面粗糙的物质或是来流污染物等因素提前诱导了转捩，因而实际转捩位置取层流区能达到的较大且延展向较平整的弦向距离。

图6 迎角1°时试验与计算转捩位置对比

表2为不同迎角下试验和计算转捩位置与弦长c之比的对比数据，其中计算网格2和计算网格3分别为中央翼段与表1所示的计算网格2与计算网格3节点分布相同的计算网格。需要说明的是，由于计算网格1无法准确模拟低速条件下的Cp与Cf分布，故不再对其进行验证。从对比数据可看出在高速条件下两套网格预测的转捩位置吻合良好，计算所得的转捩位置在-2°与-1°迎角下与计算结果存在较大差异，而随着迎角增加，计算误差逐渐减小，而迎角达2°时误差仅有2%，考虑后续对翼型优化时高速条件下主要保证其巡航状态具有较好的气动特性，而计算预测的巡航状态转捩位置具有较好的合理性，因而计算模型可用于后续层流翼型的优化。

表2 试验与计算转捩位置对比

图7为迎角1°时翼段上翼面在对称面附近的试验Pressure-Sensitive Paint (PSP)压力系数云图，其中横坐标为距前缘弦向距离与弦长之比，纵向范围是以中央翼段展向50%站位为中心，两侧距离为20%弦向长度的区域。从图7可知，翼段对称面附近区域的Cp分布较为均匀。图8给出了翼段对称面处的试验与计算Cp曲线。可看出计算网格2与计算网格3的计算结果吻合良好，说明计算网格2具有良好的计算准确性。但从1°迎角的计算与试验结果对比可知，除前缘部分外计算结果均高于试验结果。为分析该原因，给出了0°迎角的计算结果，0°迎角的计算Cp与迎角1°的试验Cp吻合较好，说明试验迎角存在一个系统误差，而前缘Cp存在差异的原因可能为实际加工的试验翼段前缘弯度不匹配或测压装置使实际模型与计算模型存在一定差异。>

图7 迎角1°时试验PSP压力系数云图

图8 试验与计算压力系数对比

4 中央翼段层流翼型优化修形及分析

4.1 计算模型及计算条件

图9为层流验证机的数模，其中中央翼段翼型为针对验证机在高速巡航状态飞行时设计的翼型。由于主要分析中央翼段的气动特性，并从翼型的角度对其气动特性进行优化修形，故需保证中央翼段近物面的计算网格与图1所示的Aerospatiale-A翼型计算网格相同。

图9 层流验证机的数模

针对低速与高速飞行状态对层流验证机的绕流流场进行数值计算，分析中央翼段在低速及高速条件下的气动特性。低速计算状态为马赫数Ma=0.20、海拔高度H=0、雷诺数Re=4.93×106；高速计算状态为Ma=0.7、H=8 000 m、Re=7.86×106、湍流强度Tu=0.4%，湍流黏性比为10。

4.2 计算结果对比及分析

由于计算中采用的数模沿对称轴对称，图10给出了采用原始层流翼型时层流验证机不同迎角下半模上翼面Cp云图及流线图。可看出在8°迎角时中央翼段前缘在与机身连接处已出现了小尺度分离；随迎角增加至10°，前缘分离已由翼身连接处向机翼中段延伸；当迎角增加到12°后，分离逐渐向后缘发展；而在迎角达14°时分离区几乎已占据半个上翼面。

图10 低速条件下原始翼型中央翼段上翼面压力系数及流线分布

图11为高速巡航状态下中央翼段为原始层流翼型的上翼面Cp、流线与Cf分布情况。可知在Ma=0.7的巡航状态(升力系数保持在0.19)下，上翼面的流线保持了较好的附着状态，且从上翼面Cf分布可知，中央翼段保持了较好的层流特性。

图11 高速条件下原始翼型中央翼段上翼面压力系数、流线与摩擦阻力系数分布

针对低速情况下上翼面存在的分离问题，通过对翼型进行适当优化修形减缓低速大迎角情况下的分离发展。对于翼型外形的调整主要集中在翼型前缘下表面位置，通过减小前缘下表面位置的曲率变化改善前缘驻点附近气流的速度特性以延缓分离。翼型上表面不做修形以维持其高速层流特性不变。

图12所示为修形前后翼型对比，可看出改进前的翼型前缘具备尖劈特征、钝度很小、极易触发分离，因此只需对前缘外形进行微量调整，消除尖劈特征，即可获得较为显著的失速特性改善效果。改进后的翼型下表面曲率过渡更加平缓，上表面曲率仅因曲率过渡而存在略微的调整，而除前缘附近的曲率变化外，改进翼型与原始翼型的差异非常小。

图12 修形前后翼型对比

图13所示为改进中央翼段翼型后低速条件下不同迎角上翼面的压力系数Cp分布及流线，可看出当迎角为8°时，相比于图10相同计算状态下中央翼段采用原始翼型的结果，中央翼段前缘在与机身连接处出现的小尺度分离已得到了较好的抑制。迎角增大到10°时，该处分离亦没有出现，说明减小前缘下表面的曲率变化对低速前缘分离的初始发展有良好的抑制作用。在迎角增大到12°后，中央翼段前缘与机身连接处出现了较小的分离，从流线形态上来看，该分离强于原始翼型在8°迎角时该处的分离，说明改进翼型可将初始分离推迟2°～4°迎角。当迎角增加到14°时分离已向翼段中央发展，相较于原始翼型在14°迎角时的分离形态，改进翼型后的翼段整体分离程度较弱，且翼段中央部分分离程度低于翼段与机身结合处，说明分离区还未由结合处向翼段中央完全发展，较原始翼型的翼段分离看，翼型改进已较大程度地改善了低速大迎角情况下的前缘分离特性，对气动特性有较大提升。

图13 低速条件下改进翼型后中央翼段上翼面压力系数及流线分布

图14为高速条件计算得到的巡航状态下原始翼型上翼面Cp、流线与Cf分布。与图11相同条件下中央翼段采用原始翼型数据的计算结果进行对比，可见改进翼型后中央翼段上表面的流线几乎与原始模型一致，说明改进翼型并未使翼段上翼面的流态在高速巡航状态下恶化，且从Cp云图与Cf云图对比可看出改进翼型后，中央翼段上表面的Cp和Cf分布几乎与原始计算模型的Cp和Cf相同，说明改进翼型保持了原始模型在高速巡航状态下的气动特性。

图14 高速条件下改进翼型中央翼段上翼面压力系数、流线与摩擦阻力系数分布

图15为高速巡航状态下中央翼段Cp在50%展向站位处的曲线对比。从对比曲线可知，中央翼段翼型改进前后，除前缘附近的Cp分布有较小差异外，整体Cp分布几乎一致，说明改进翼型较好地保持了中央翼段的气动特性。

图15 高速条件下中央翼段剖面压力系数分布

图16为低速状态下中央翼段翼型改进前后全机升力曲线与升阻极曲线对比,图中CL和CD分别为升力系数和阻力系数。可知改进翼型后低速计算状态下全机的升力整体高于采用原始翼型的全机升力，虽升力曲线的非线性段均出现在10°迎角后，但随迎角增加，改进翼型对全机升力的贡献逐渐明显。对于采用原始翼型的全机升力曲线，其拐点出现在迎角为16°时，而采用改进翼型后拐点被推迟到了迎角为18°时，且最大升力系数提高了近0.15，进一步说明中央翼段采用改进翼型后可较好地提升全机的气动性能。从升阻极曲线可知改进翼型在低升力时气动特性略低差于原始翼型，但在升力系数高于1.2后其气动特性较原始翼型有明显的提升。

图16 低速条件下全机气动特性曲线对比

5 结 论

研究的特殊构型层流验证机机身通过中央层流验证翼段连接，该中央翼段对全机气动特性有较大影响。通过标准Aerospatiale-A翼型验证算例对计算模型及基准网格分布的准确性进行评估，并利用风洞试验对高速巡航状态的计算精度进行考量，在保证计算精确性的前提下利用数值计算对中央翼段的翼型进行优化修形，在兼顾高速巡航状态气动特性的前提下合理提升低速状态的气动特性。

1) 在对Aerospatiale-A翼型的验证计算中，对比了标准两方程SST模型和四方程γ-Reθ转捩模型在压力分布、摩擦阻力系数与上翼面不同站位速度型分布上的预测情况。标准SST模型预测的上翼面压力分布整体低于试验值，且由于其无法预测转捩过程，其上翼面摩擦阻力系数的模拟结果也与试验值相差较大，而四方程γ-Reθ转捩模型的计算结果与试验结果吻合良好。在速度型的预测上，标准SST模型预测的附面层厚度在不同站位上均较厚，整体而言四方程γ-Reθ转捩模型对速度型的预测更符合实际。

2) 在高速巡航状态的验证中，利用试验对转捩位置和翼面压力分布进行了对比验证。从转捩位置看，虽在迎角较小时计算与试验结果存在一定差异，但在巡航角度附近转捩位置的预测与试验结果匹配较好。从压力分布预测看，计算结果具有与试验结果匹配较好的趋势。总体而言，在高速巡航状态下的计算结果较为合理。

3) 在低速条件下将改进翼型应用于中央翼段后，前缘分离从迎角为8°延缓到迎角为12°，虽在迎角为14°时翼面也出现了较为严重的分离，但相比于原始翼型在迎角为14°时的流线，其整体分离程度依旧较小。结合低速状态的升力和升阻极曲线，改进翼型明显整体提升了低速状态下的全机气动特性。对于高速巡航状态，改进翼型仅在前缘附近与原始翼型的压力分布有略微差异，整体而言其并未对物面参数产生影响，说明新翼型满足优化目标。