冰脊对翼型气动特性影响的数值模拟研究

周莉，徐浩军，杨哲，蔡军

（1.空军工程大学工程学院，陕西西安 710038；2.中国人民解放军驻陕飞公司军代表室，陕西汉中 723213）

引言

飞机穿越含有过冷水滴的云层时，机体碰撞过冷水滴后会在其部件表面出现结冰。对于机翼来说，水滴易撞击在翼剖面的前缘，使得机翼前缘结冰现象较为普遍，因此，出现了大量关于机翼前缘结冰过程及其对飞机性能影响的研究文献［1-3］。但当过冷水滴的尺寸超出FAR25附录C所给的2～50 μm结冰包线范围时，在机翼上表面防冰区域后部会形成突起的展向冰脊，造成飞机气动性能恶化，严重时甚至会引发飞行事故。

近几年来，国外对大尺度过冷水滴的碰撞特性以及冰脊的形成进行了一系列试验与数值模拟研究。Tan等［4］通过风洞试验研究了大尺度过冷水滴的碰撞、喷溅和分裂特性，为带冰脊翼型气动特性的分析奠定了基础。Raimund等［5］采用欧拉法对大尺度过冷水滴的碰撞特性进行了数值模拟，所得收集效率的分布与试验结果吻合较好。此外，Lee等［6］和Pan等［7］分别运用试验手段以及雷诺数平均的N-S方程研究了展向冰脊对翼型气动特性的影响，使人们对冰脊危害的认识不断深入。相比之下，国内关于冰脊的研究工作开展得较少，仅文献［8］对冰脊高度和弦向位置对气动特性的影响进行了仿真，但所计算的迎角范围较小，没有充分给出带冰脊翼型的失速特性，也没有考虑冰脊几何形状以及雷诺数对翼型性能的影响。针对冰脊对翼型气动的危害以及前期研究的不足，本文采用数值模拟方法，系统研究了不同冰脊几何形状、弦向位置、冰脊高度以及不同雷诺数对翼型气动性能的影响。

1 数值计算方法

1.1 网格划分

本文计算采用O型结构网格，选取具有代表性的前向1/4圆、后向1/4圆、半圆和前向斜坡来模拟不同形状的冰脊，其中斜坡面的基底长度与高度比为3。考虑到国外针对NACA23012翼型已开展了大量结冰方面的研究工作，有较为丰富的试验数据，本文选用该翼型进行数值模拟。图1为NACA23012翼型在未结冰以及分别加入以上各种冰脊后的翼型网格。

图1 计算网格

1.2 方法验证

本文对未结冰的NACA23012翼型进行数值模拟，并将计算结果与风洞试验结果对比，以此来验证本文数值计算方法的可行性和准确性。

采用无反射压力远场作为远场边界条件，计算时给定无穷远处来流马赫数及其方向。翼型表面为固壁边界，流动无滑移、无穿透。针对SSTk-ω在近壁面区有良好的精度和算法稳定性的特点，选用SSTk-ω湍流模型进行计算［9］。

计算验证状态为:Ma=0.2，Re=1.8 ×106，α=5°。经过流场计算，表1中给出了本文所得升力系数、阻力系数的计算结果。与文献［6］的试验值进行比较，结果表明本文升力系数及阻力系数的计算结果与试验值较为接近。由图2中给出的NACA23012翼型数值计算与风洞试验所得压力分布的对比情况同样可以看出，本文计算结果与试验值吻合较好。因此，运用本文的计算方法模拟翼型外流较为可靠。

表1 计算气动力系数与试验值的比较

图2 NACA23012翼型计算与试验压力系数对比

2 带冰脊翼型的数值模拟

在不同的气象条件以及飞行条件下，不同弦向位置将出现不同形状、不同高度的冰脊，对飞机的气动特性会产生不同程度的影响。通过数值模拟，得到不同条件下的翼型气动系数。计算基准条件为:Ma=0.12，Re=10.5 ×106，冰脊高度为k/c=1.39%，冰脊弦向位置为x/c=0.1。其中，k为冰脊高度，x为冰脊在弦向的坐标，c为翼型弦长。

2.1 冰脊几何形状的影响

在未结冰以及冰脊几何形状分别为前向1/4圆、后向1/4圆、半圆以及前向斜坡四种类型时，通过数值计算得到不同迎角下的升力系数、阻力系数如图3所示。

由图3（a）可以看出，与未结冰情况相比，结冰脊后最大升力系数下降、失速迎角减小，且不同冰脊形状对升力系数的影响程度并不相同。冰脊为前向1/4圆和前向斜坡时对升力系数的影响较为接近、且最为严重；冰脊为半圆时对升力系数的影响最小；而冰脊为后向1/4圆时对升力系数的影响程度介于半圆和前向1/4圆之间。对比具有相同钝迎风面形状的前向1/4圆和前向斜坡冰脊的升力系数可知，虽然下游形状差别较大，但两者的升力系数曲线几乎一致，可见钝迎风面冰脊的下游形状对性能的影响并不是很重要；而通过具有相同流线型迎风面形状的后向1/4圆和半圆冰脊的升力系数比较、以及具有相同下游形状的前向1/4圆和半圆冰脊的升力系数比较可知，流线型迎风面冰脊的下游形状以及冰脊迎风面的形状对性能的影响较大，且冰脊迎风面形状的影响程度更大。

图3（b）所示为不同冰脊形状下的阻力系数变化。可以看出，结冰脊后阻力系数均增大，且与不同形状的冰脊对升力系数的影响趋势一致。

图3 冰脊几何形状对气动系数的影响

2.2 冰脊弦向位置的影响

冰脊选用前向1/4圆，冰脊弦向位置分别为x/c=0.02，0.10和0.20时，不同迎角下的升力系数、阻力系数曲线如图4所示。

由图4（a）的计算结果可以看出，x/c=0.10时的最大升力系数最小，并且最早失速，可见此处对冰脊最为敏感，因此x/c=0.10处被认为是 NACA23012翼型最危险的冰脊位置。图4（b）中，阻力系数随迎角的变化曲线同样也说明在冰脊高度一定的情况下，x/c=0.10处的气动特性恶化最为严重。

图4 冰脊弦向位置对气动系数的影响

图5为α=0°时未结冰翼型上表面的压力云图。可以发现，最危险的冰脊位置与此处流速（压力）大小及逆压梯度大小有关。这是因为在流速越大（压力越小）处结冰会消耗掉更多的边界层能量，同时，当结冰区越靠近最大逆压梯度区域，则冰脊引起的下游气流分离越严重，所产生的分离气泡越大，导致翼型气动性能损失越大。而x/c=0.10恰好同时处于NACA23012翼型的最小压力与最大逆压梯度区域，因此冰脊在此处对翼型气动特性的影响最为明显，即此处为冰脊的危险位置。

图5 未结冰翼型上表面的压力云图

图6给出了迎角α=0°时冰脊分别在x/c=0.02，0.10和0.20处的压力系数分布曲线。可以看出，计算所得Cp分布与文献［6］的试验值基本吻合，但是并没有充分捕捉到冰脊后的分离区。对比图6可得出以下结论:冰脊在x/c=0.02时对Cp分布的影响最小，这是因为x/c=0.02处于顺压梯度范围，所以冰脊所引起的分离泡能够较容易地重新附着到翼型表面，导致分离泡较小；与前面压力云图的分析结果一致，x/c=0.10时冰脊上的Cp变化量最大，即冰脊影响最为严重；冰脊在x/c=0.20时，Cp分布变化程度介于x/c=0.02和x/c=0.10之间，这主要是因为虽然冰脊处于逆压梯度区导致分离泡较大，但是x/c=0.20并不在最低压力范围，导致此处冰脊的影响并不是最严重。

图6 冰脊在不同弦向位置处的压力系数分布

2.3 冰脊高度的影响

冰脊选用前向1/4圆，当冰脊高度分别为k/c=0（未结冰），0.83%和1.39%时，升力系数、阻力系数曲线如图7所示。

图7（a）中，与未结冰情况相比，带有冰脊时的升力系数曲线斜率减小，且最大升力系数及失速迎角明显降低。k/c=1.39%时，最大升力系数由未结冰时的1.69降为0.23，失速迎角由17.7°减小到3.93°。k/c=0.83%时，最大升力系数为0.41，略大于k/c=1.39%时的0.23，失速迎角为5.98°，也高于k/c=1.39%时的3.93°。可见，冰脊越高时最大升力系数及失速迎角越小。

由图7（b）可以看出，与未结冰时的阻力系数相比，k/c=0.83%和1.39%时阻力系数均增大，并且当冰脊越高时，阻力系数的增加幅度越大。

图7 冰脊高度对气动系数的影响

综上分析可知，冰脊使翼型的气动性能下降明显，且冰脊越高，影响越大，即翼型气动特性越差。

图8给出了迎角α=0°、冰脊高度分别为k/c=0，0.83%和1.39%时的压力系数分布曲线。

图8 不同冰脊高度的压力系数分布

可以看出，计算所得Cp分布与文献［6］的试验值基本吻合，与图6相似，并没有很好地捕捉到冰脊后的分离区，但是也充分反映出了不同冰脊高度下的翼型压力系数分布:冰脊的存在使得翼型前缘上表面的压力增加，即前缘吸力峰被拉低，导致升力减小；同时，随着冰脊高度的增加，下表面的压力会减小，从而造成升力的进一步损失。

2.4 雷诺数的影响

冰脊选用前向1/4圆，雷诺数分别为Re=3.5×106，7.5 ×106，10.5 ×106时，计算所得升力系数、阻力系数曲线如图9所示。

图9 雷诺数对气动系数的影响

从图9可以看出，三种不同的雷诺数下升力系数及阻力系数变化并不明显，可见雷诺数对带冰脊翼型的气动特性影响较小。虽然对于未结冰翼型来说，当雷诺数增大时转捩点前移，导致紊流区增大，从而延迟分离，使最大升力系数和失速迎角增大。但是当翼型带有冰脊后，在冰脊处会发生气流分离，此时对失速起决定性作用的已不再是附面层的气流特性，因此改变雷诺数对气动特性影响不大。

3 结论

本文针对冰脊对翼型气动特性的严重影响，运用数值模拟方法系统研究了不同冰脊特性条件下翼型气动性能的变化。计算结果表明:

（1）飞机在结冰脊后，最大升力系数下降、失速迎角减小，且阻力系数增大。

（2）冰脊迎风面的形状对翼型性能的影响较大，并且当迎风面不同时，冰脊下游形状对翼型性能的影响程度也不同:流线型迎风面冰脊的下游形状影响翼型性能，而钝迎风面冰脊的下游形状对性能的影响并不明显。

（3）不同冰脊弦向位置对气动特性的影响不同，其中，x/c=0.10为最危险冰脊位置。

（4）冰脊高度也会影响翼型气动特性，在一定冰脊高度范围内，气动特性的恶化程度随冰脊高度的增加而增大。

（5）雷诺数对带冰脊翼型气动性能的影响并不明显。

［1］孙志国，朱程香，付斌，等.二维翼型结冰数值计算［J］.航空动力学报，2010，25（7）:1485-1490.

［2］李林，王立新，彭小东.结冰对民机飞行性能的影响研究［J］.飞行力学，2004，22（3）:12-16.

［3］钟长生，杜亮，洪冠新.飞机结冰引起的飞行动力学问题探讨［J］.飞行力学，2004，22（3）:64-68.

［4］Tan S C，Papadakis M，Miller D，et al.Experimental study of large droplet splashing and breakup［R］.AIAA-2007-904，2007.

［5］Raimund H，Wagdi G H.Eulerian modeling of in-flight icing due to supercooled large droplets［J］.Journal of Aircraft，2008，45（4）:1290-1296.

［6］Lee S，Bragg M B.Experimental investigation of simulated large-droplet ice shapes on airfoil aerodynamics［J］.Journal of Aircraft，1999，36（5）:844-850.

［7］Pan J P，Loth E.Reynolds-averaged Navier-Stokes simulations of airfoils and wings with ice shapes［J］.Journal of Aircraft，2004，41（4）:879-891.

［8］易贤.飞机结冰的数值计算与结冰试验相似准则研究［D］.绵阳:中国空气动力研究与发展中心，2007.

［9］王福军.计算流体动力学分析:CFD软件原理与应用［M］.北京:清华大学出版社，2004.