大粒径过冷水溢流结冰的翼型气动影响分析

李焱鑫，张 辰，刘 洪，王福新

（上海交通大学 航空航天学院，上海 200240）

0 引 言

2002年12月21日，台湾复兴航空一架ATR－72货机，在飞行中遭遇严重结冰天气，飞机在积冰后升阻特性大幅下降，最终坠毁。事故调查报告［1］表明，失事飞机的飞行环境中包含大粒径过冷水滴（SLD）。SLD的粒径大于50mm，超出了适航条例25部附录C的规定。飞机在遭遇SLD结冰环境时，过冷水滴在机翼表面会继续向后方流动，并在飞机除冰装置后方冻结。当除冰装置开启时，机翼前缘受保护区的积冰被清除，机翼表面在未保护区形成脊状的溢流积冰［2－3］。

鉴于SLD结冰严重的危害性，2005年美国联邦航空局（FAA）在工作报告［4］中，将SLD结冰列为重点研究方向，并在适航咨询通告［5］中特别提出，需要对SLD引起的特殊冰型进行审定。Broeren等［6］通过实验比较了多种积冰外形的气动性能损失，发现在实验条件下脊状冰使翼型的最大升力系数由1.85下降到0.52，失速迎角由18.1°下降至5.6°，比流向冰和角状冰对升阻的破坏都要大。Edward等［7］进一步分析指出，SLD结冰的危害性在于，溢流冰脊使翼型流场提前分离。

脊状冰的分离效应与翼型几何外形密切相关。Lee［8］分析了脊状冰对不同翼型气动性能的影响，NACA 23012m和NLF 0414的最大升力系数（CLmax）分别下降了83.7%和54.8%，这表明脊状冰的影响与翼型形状有关。超临界翼型广泛应用于大型客机的设计中，但前人并未对此进行单独分析。Bragg等［2，9］指出脊状冰生长在翼型的主要升力产生区，使翼型的最大升力系数下降。超临界翼型的主要升力产生区，与溢流结冰位置有较大重合，易形成吸力峰引发气动分离，因此本文着重分析了溢流积冰对超临界翼型的影响。

目前研究结冰影响的方法主要为风洞实验，但风洞实验数据获取不全面，具有一定局限性。数值分析方法具有高效、全面的特点，得到了广泛应用。在翼型结冰后的数值气动分析领域，陈科等［10］应用SST模型对机翼前缘的复杂积冰引起的流动分离进行了分析；Mortensen［11］比较了多种湍流模型，认为 S－A模型用于评价前缘积冰的流动分离具有较高精度。上述两位学者均未对溢流脊状冰进行分析。Dunn等［12］使用NSU2D程序分析了带脊状冰的NACA 23012m的升力曲线，计算结果在大迎角条件下误差较大，不能准确预测出带冰翼型的CLmax和失速迎角αstall。综合分析前人工作，本文发展了一种更适合溢流结冰气动分析的数值方法。

我国面临大型民机设计的关键阶段，为保证安全性，迫切需要分析SLD溢流结冰对机翼气动性能的影响，为我国的民用飞机设计和适航审定提供指导。为此，本文采用数值计算的方法，对溢流冰脊的气动性能进行了分析。本文介绍了验证算例的计算条件，确定了用于溢流冰脊计算的湍流模型，在此基础上验证了计算网格的有效性；将不同冰脊高度的计算结果与实验结果对比，确定了方法的有效性；分析了机翼超临界翼型和平尾翼型的流场特征；分析了SLD溢流结冰条件下，大型客机机翼超临界翼型和平尾翼型各项气动系数的变化。

1 研究方法

1.1 计算模型与网格

本节以Lee［8］的实验数据和实验条件为依据，在CFD软件中对带脊状冰的NACA 23012m的升力曲线进行计算，并与实验结果对比。由于翼型气动性能受脊状冰高度变化的影响较大［13］，本文在确定数值方法时选择了不同的冰脊高度k／c。SLD结冰时具有溢流积冰的特点，脊状冰位置x／c选定为机翼保护区后方x／c＝0.1。冰脊形状为前端面半圆，Lee和Bragg［14］的实验表明此形状的冰脊对翼型的气动破坏作用最明显。计算条件参见表1。

表1 算例计算条件Table 1 Computational conditions of the case

1.1.1 湍流模型

本文基于N－S方程求解带冰翼型流场，脊状冰会引发翼型流场的气动分离，需要选择合适的湍流模型对流场进行分析。Mortensen［11］比较了不同湍流模型在前缘结冰算例计算中的精度，包括S－A模型、重整化 RNG k－ε 模型、可实现性 Realizable k－ε 模型、标准k－ε模型，文献计算结果表明S－A模型在以上几种湍流模型中最优。雷诺应力模型［15］（RSM）并没有采用上述几种模型的涡粘各向同性假设，而是直接求解雷诺应力分量，理论上应更为精确，因此本文采用RSM模型作为计算用湍流模型。

1.1.2 网格验证

考虑到带冰翼型表面较复杂本文生成非结构网格。为验证网格的无关性，生成尺度不同的三套网格，并对网格进行自适应加密验证。

在本文的模拟中，网格区域为20倍弦长的圆，翼型位于区域中心。为捕捉近壁面流场特征，对初始近壁面网格进行了加密处理，在初步计算时，选择了不同的网格加密尺度，见表2。

表2 网格信息Table 2 Grids details

图1是不同网格尺度下的计算升力曲线与实验结果的对比。在图中，三种网格在低角度下均与实验结果符合较好，但在接近失速迎角时，加密后的网格3与实验结果符合较好。

图1 不同网格计算结果对比Fig.1 Comparison of results on different grids

在网格3基础上进行初步计算，得到带冰翼型流场的基本信息，按照速度梯度的变化对原始网格进行自适应加密，得到自适应网格。此时网格区域含185259个非结构网格，共107860个网格节点，如图2。从图3中可见，在使用自适应网格对带冰原始网格进行加密后，自适应网格对失速迎角的计算结果与网格3一致，在不同冰脊高度k／c＝0.0056、0.0083、0.0137时，两套网格在最大升力系数计算上的误差分别为0.02、0.01和0.01。因此，使用网格3对溢流结冰进行计算时，已具有相当精度，进一步进行自适应加密后，计算结果变化很小，网格无关性可得到验证。

图2 计算网格Fig.2 Computational grids

图3 自适应网格计算结果对比Fig.3 Comparison between original grids and adaption grids

1.2 算例计算结果

本文应用结合RSM模型的计算方法，将计算结果与实验进行对比，并分析计算误差。如图4（a），当k／c＝0时，CLmax的计算值比实验值偏大约0.05，失速角的计算值比实验值高约0.6°；如图4（b），当k／c＝0.0056时，CLmax的计算最大误差值发生在－4°，为0.13；如图4（c），当k／c＝0.0083时，CLmax计算最大误差值发生在4°，为0.06；如图4（d），当k／c＝0.0137时，CLmax计算最大误差值发生在－4°，为0.09（图中右侧坐标为升力系数计算绝对误差）。冰脊高度的增长，未使计算误差产生明显增加。

在图4中RSM模型计算所得最大升力系数和失速迎角与实验结果符合较好。表3给出了RSM模型的最大升力系数、失速迎角的计算结果对比。综合分析发现，RSM模型计算的最大升力系数相对误差不超过5.8%，不同迎角下升力系数的绝对误差小于0.2；除k／c＝0.0056外，失速角的预测误差小于1°。可认为基于雷诺应力的流场计算方法，能够应用于带脊状冰的翼型气动性能损失评估。

表3 最大升力系数和失速迎角结果对比Table 3 Comparison between the computation and the experiment on CLmaxandαstall

图4 计算结果与实验结果对比Fig.4 Comparison between the computation and the experiment results

2 大型客机翼型结冰影响分析

本节采用上文确定的研究方法，开展大型客机机翼超临界翼型和平尾翼型的SLD溢流结冰气动分析。由于机翼翼型沿展向的变化较大，本文截取了机翼和平尾在展向不同区域的四个典型翼型截面作为代表性翼型进行带冰气动力分析，分别命名为AF－1、AF－2、AF－3和 AF－T，如图5。

图5 飞机翼型及冰脊形状Fig.5 Aircraft airfoils and ridge ice shape

Bragg指出［2］，雷诺数对带冰翼型最大升力系数的影响非常有限。考虑到结冰多发于起降阶段，来流条件的设置参考民用客机实际的飞行环境，脊状冰位置x／c＝0.1，见表4。

表4 翼型计算条件Table 4 Computational conditions of airfoil

2.1 流场特征分析

2.1.1 机翼超临界翼型流场特征

如图6，以翼型AF－3为例，对比超临界翼型未结冰和带脊状冰两种情况在迎角14°时的流线：当AF－3未结冰时，流动未发生分离；当 AF－3带有k／c＝0.0137的脊状冰时，冰脊导致流动产生了分离，分离区自冰脊处开始生长，在x／c＝0.4处再附着；AF－1／2有类似特征。

图6 AF－3流线图（α＝14°）Fig.6 Streamlines for AF－3（α＝14°）

图7给出了翼型AF－3对应的表面压力分布：在x／c＝0～0.1范围，带冰翼型上表面压力系数（负值）大于干净翼型；在x／c＝0.1处冰脊使上表面压力系数（负值）减小；在冰脊后x／c＝0.1～0.4的范围，上表面压力系数（负值）逐渐增大；在x／c＝0.4～1.0范围内，压力系数变化很小。

图7 AF－3表面压力分布对比（α＝14°，k／c＝0.0137）Fig.7 Surface pressure distributions for AF－3（α＝14°，k／c＝0.0137）

综合分析，冰脊导致机翼超临界翼型表面发生流动分离。从压力分布来看，与干净翼型相比，带冰翼型上表面压力系数较大而下表面较小，上下表面压差小于干净翼型。

2.1.2 平尾翼型流场特征

由于平尾翼型与超临界翼型的几何结构不同，受脊状冰影响有所差异。如图8，以翼型AF－T为例，对比平尾未结冰与带脊状冰两种情况在14°时的流线特征：图8a）中，干净翼型在该条件下完全分离，已发生失速；但当 AF－T带有k／c＝0.0137的脊状冰时，虽然冰脊同样引发了分离，但对整个流场而言，冰脊反而起到了整流的作用，分离区明显减小。

图8 AF－T流线图（α＝14°）Fig.8 Streamlines for AF－T（α＝14°）

图9给出了翼型对应的表面压力分布：在x／c＝0～0.25范围内，带冰翼型上表面的压力系数（负值）小于干净翼型；在x／c＝0.25～1.0范围内，带冰翼型上表面压力系数持续增大，而干净翼型基本不变。

图9 AF－T表面压力分布对比（α＝14°，k／c＝0.0137）Fig.9 Surface pressure distributions for AF－T（α＝14°，k／c＝0.0137）

综合分析，对于平尾翼型，在大迎角下干净翼型已发生分离，但在结冰翼型中由于冰脊的存在，使分离区明显减小。从压力分布看，在x／c＝0～0.25范围内，带冰翼型上下表面压差大于干净翼型，x／c＝0.25～1.0两者压差相若。

2.2 气动参数分析与讨论

2.2.1 升力曲线分析

从结冰翼型的流场特征来看，超临界翼型与平尾翼型受溢流结冰的影响有所不同。为进一步评估结冰的气动影响，对翼型结冰后的升力系数的变化进行分析。

图10给出了大客翼型在不同冰脊高度下的的升力曲线。从变化趋势来看，脊状冰高度k／c＝0.0137时，AF－1／2／3 和 AF－T的CLmax分别下降65.9%、56.1%、57.7% 和 38.5%。AF－1、AF－2、AF－3的失速迎角分别提前8.2°、5.7°和6.6°，而 AF－T的失速迎角延后2.9°。

图10 翼型升力曲线Fig.10 Lift coefficient curve of airfoils

图11给出了大型客机各段翼型不同冰脊高度的CLmax变化图，由该图分析脊状冰高度对CLmax的影响：AF－1在积冰高度k／c由0增长到0.0056时CLmax从1.76 减小到0.97，k／c 达到0.0137 时CLmax减小到0.60，主要的升力破坏发生在结冰初始阶段；翼型 AF－2、AF－3与 AF－1保持着类似的趋势；AF－T在积冰高度k／c由0增长到0.0056时CLmax从1.19减小到0.85，k／c达到0.0137时CLmax减小到0.72，升力降幅较前三类翼型相对平稳。

图12给出了飞机各段翼型不同冰脊高度的失速迎角变化图，由该图分析脊状冰高度对失速迎角的影响：翼型AF－1～AF－3的失速迎角随冰脊高度增加而提前，在k／c达到0.0056时，失速迎角降幅不足2°，随着k／c达到0.0137，失速迎角产生4°～6°不等的降幅，失速角前移主要发生在冰脊形成阶段；翼型AF－T的失速迎角随冰脊高度的增加而延后，在k／c达到0.0137时，失速角比未结冰时延后2.9°，且失速角延后幅度基本与冰脊高度增加呈线性关系。

图11 最大升力系数对比Fig.11 Comparison of maximum lift coefficient

图12 失速迎角对比Fig.12 Comparison of stalling angles

2.2.2 阻力曲线分析

从结冰事故报告［1］分析，在飞机发生结冰后，飞机受到的阻力显著增加，同时原本起稳定作用的纵向操纵参数将会减小甚至反号，因此有必要对结冰后超临界翼型和平尾翼型阻力曲线和力矩曲线的变化进行分析。

图13是翼型阻力曲线图，分别以 AF－3、AF－T为例，分析超临界翼型和平尾翼型结冰后阻力曲线的变化。在机翼超临界翼型和平尾翼型中，溢流积冰使翼型的阻力系数大幅增加，增长幅度随迎角的增加而变大。

2.2.3 力矩曲线分析

图14是翼型力矩曲线图，分别以 AF－3、AF－T为例，分析超临界翼型和平尾翼型结冰后力矩曲线的变化。在AF－3翼型中，结冰后力矩系数下降，随溢流冰高度增加力矩非线性变化趋势增加。在AF－T翼型中，结冰后力矩系数的变化显著小于超临界翼型，结冰后力矩曲线存在线性变化趋势。

图13 翼型阻力曲线Fig.13 Drag coefficient curve of airfoils

图14 翼型力矩曲线Fig.14 Pitching moment curve of airfoils

3 结 论

（1）RSM模型适合于分析溢流结冰的流动分离，可准确预测带冰翼型的最大升力系数与失速迎角。

（2）超临界翼型受到SLD溢流结冰影响时，翼型表面压差减小，流场提前发生分离，最大升力系数和失速迎角显著下降；平尾翼型受溢流结冰影响较小，在大迎角下带冰翼型流场产生的分离区小于干净翼型，失速迎角提前。超临界翼型和平尾翼型结冰后阻力系数显著增加。结冰后超临界翼型的力矩系数显著变化，平尾翼型的力矩变化较小。

（3）本文所做SLD溢流结冰对超临界翼型和平尾翼型气动影响的分析，对大型客机的设计和适航审定具有一定指导意义。

［1］FAA.Occurrence investigation report：in－flight icing encounter and crash into the sea，Transasia airways flight 791［R］.ASC－AOR－05－04－001，2005.

［2］BRAGG M B，BROEREN A P，BLUMENTHAL L A.Icedairfoil aerodynamics［J］.Progress in Aerospace Sciences，2005，41（5）：323－362.doi：10.1016／j.paerosci.2005.07.001

［3］LYNCH F T，KHODADOUST A.Effects of ice accretions on aircraft aerodynamics［J］.Progress in Aerospace Sciences，2001，37（8）：669－767.doi：10.1016／S0376－0421（01）00018－5

［4］Ice Protection Harmonization working group（IPHWG）.Task 2 Working group report on supercooled large droplet rulemaking［R］.2005.

［5］FAA.Advisory Circular：aircraft ice protection［Z］.20－73A.2006.

［6］BROEREN A P，BRAGG M B，ADDY H E，et al.Effect of high－fidelity ice accretion simulations on the performance of a full－scale airfoil model［R］.AIAA 2008－434.doi：10.2514／6.2008－434

［7］EDWARD A W，BROEREN A P，BRAGG M B.Aerodynamics of scaled runback ice accretions［J］.Journal of Aircraft，2008，45（2）：591－603.doi：10.2514／1.32274

［8］LEE S.Effects of supercooled large droplet icing on airfoil aerodynamic［D］.Urbana：University of Illinois，2001.

［9］BROEREN A P，BRAGG M B.Effect of airfoil geometry on performance with simulated intercycle ice accretions［R］.AIAA 2003－0728.doi：10.2514／1.4734

［10］CHEN K，CAO Y H，AN K W，et al.Application of hybird grid to analyzing complex iced airfoil aerodynamic performance［J］.ACTA Aeronautica et Astronautica Sinica，2007，28（s）：87－91.（in Chinese）陈科，曹义华，安克文，等.应用混合网格分析复杂积冰翼型气动性能［J］.航空学报，2007，28（s）：87－91.

［11］MORTENSEN K.CFD simulations of an airfoil with leading edge ice accretion［D］.Copenhagen：Technical University of Denmark，2008.

［12］DUNN T A，LOTH E，BRAGG M B.Computational investigation of simulated large－droplet ice shape on airfoil aerodynamics［J］.Journal of Aircraft，1999，36（5）：836－843.doi：10.2514／2.2517

［13］BROEREN A P，EDWARD A W，BUSCH G T，et al.Aerodynamic simulation of runback ice accretion［J］.Journal of Aircraft，2010，47（3）：924－939.doi：10.2514／1.46475.

［14］LEE S，BRAGG M B.Effects of simulated－spanwise ice shapes on airfoils：experimental investigation［R］.AIAA Paper 99－0092，1999.

［15］LAUNDER B E，SPALDING D B.The numerical computation of turbulent flows［J］.Comput.Methods in Appl.Mech.Eng.，1974，3（2）：269－289.doi：10.1016／0045－7825（74）90029－2.