结冰条件下大型飞机翼面分离流场结构及空气动力学特性研究

魏 扬, 李 杰, 李 明, 李 哲, 张 恒

(1.空军工程大学航空工程学院，西安，710038；2.西北工业大学航空学院，西安，710072；3.中国空气动力研究发展中心空气动力学国家重点实验室，四川绵阳,621000)

结冰作为引发飞行失控(Loss of Control，LOC)环境因素中最重要的因素[1]，长期以来一直受到研究人员的高度关注。飞机结冰会导致飞机的飞行性能下降和操稳特性恶化，严重时会危害飞行安全甚至造成毁灭性灾难。尽管人们对于结冰问题的研究从上世纪就已经开始，但由于结冰现象的复杂性与随机性，迄今为止，国际上对结冰的发生、发展及灾害的诱发、形成、演化等本质规律还缺乏细致深入的认识。基于现有有限认知能力发展起来的结冰防护技术，还不能完全消除结冰对飞行安全的危害，因此国内外由于结冰导致的飞行事故仍时有发生。据美国国家运输安全委员会(National Transportation Safety Board，NTSB)统计，1998～2007年，结冰有关的飞行事故共计565起，并造成229人遇难[2]。

国外很早便系统地开展针对结冰对飞机空气动力学特性的研究工作。总体来看，发展趋势是：从二元翼型积冰[3-10]研究发展到三维机翼积冰研究[11-14]、从单段机翼积冰研究发展到多段机翼(缝翼)积冰研究[15-16]、从前缘积冰研究[17-20]发展到冰脊[21-26]、回流冰[27-30]等复杂冰形的研究，研究的方法不断拓展，飞行试验、结冰风洞试验和基于CFD的数值模拟计算成为主流的研究手段。在此期间取得了不少重要的试验和理论成果，相应的试验设备、试验及测试技术、结冰动力学数值计算的基本理论以及计算软件(如LEWICE[31]、FENSAP-ICE[32]等)也得到发展。国内针对结冰后气动特性的影响研究起步较晚，大多借助于CFD数值模拟[33-37]的手段，而大型的冰风洞等结冰试验设备相较国外建造较晚，相关的报道并不多见。

目前，结冰对飞机气动特性的影响研究常见于简单机翼或翼身组合体模型，结冰后流场数值模拟方法精度较低，流场结构分析不够精细，涉及的冰形和翼型有限，对失速分离流场细节和气动特性变化过程的机理认识不足，同时缺乏全面的大型飞机结冰致灾气动数据和系统的机理分析研究。这在一定程度上限制了对结冰条件下大型飞机气动特性影响规律的本质认识。

为此，本文采用基于RANS的数值模拟方法，对构建的典型大型客机机翼结冰构型进行了数值模拟，分析了该背景飞机结冰复杂空间流场结构及气动特性，为大型飞机结冰条件下飞行动力学研究提供了数据支撑，也为翼面防除冰装置设计提供参考。

1 数值模拟方法

1.1 控制方程及其离散

不考虑体积力及外部热源，直角坐标系下三维非定常可压缩N-S方程组为：

(1)

式中：

式中：Q为守恒变量；F、G、H分别为直角坐标系下x，y，z3个方向上的无粘通量；Fv、Gv、Hv分别对应3个方向粘性通量。采用Jameson有限体积法对控制方程进行空间离散，对无粘项空间离散采用三阶MUSCL格式，对粘性项空间离散采用二阶中心格式。

采用隐式LU-SGS伪时间子迭代格式进行时间推进，该方法在经典隐式LU-SGS格式中引入伪时间导数项，借助伪时间方向的子迭代技术(τ-TS迭代)，使得时间离散精度达到二阶。

1.2 湍流模拟方法

基于雷诺平均(RANS)方法对背景飞机进行基本气动力计算工作。该方法将满足流体动力学方程的瞬时运动分解为平均运动和脉动运动两部分，脉动部分对平均运动的影响通过湍流模型描述；依靠对湍流的理论分析、试验数据或直接数值模拟结果，对雷诺应力做出各种假设，建立附加方程组以描述湍流平均量。RANS方法总体计算所需的网格数量相对较少，对计算资源的需求较易满足，因而在工程领域中得到了广泛应用。由于本文涉及全机的流场计算，考虑到计算资源的限制，网格数量不能太大，因此在本文中采用RANS方法。湍流模型采用Spalart-Allmaras(SA)模型。

2 背景飞机的几何建模及网格划分

2.1 背景飞机结冰构型三维几何数模构建

综合考虑知识产权和研究对象的工程适用性，完成了背景飞机几何模型构建工作。对A320、B737等大型客机几何外形和机翼结冰情况开展了综合调研，形成了具有典型客机特征的背景飞机模型。确定了背景飞机基本几何参数；进行了超临界机翼设计，建立了干净构型的三维数模见图1。根据类似机型结冰风洞的实验数据，形成了典型结冰环境下的机翼重度和中度角状冰模型，见图2。

图1 背景飞机的三维几何模型

图2 重度、中度角状冰的三维数模

2.2 背景飞机结冰构型网格划分

流场计算所用网格是利用先进的网格生成工具ICEM软件生成的点对接多块网格。针对机翼前缘角状冰的不规则几何形状，在保证分离区域网格质量的前提下，确定了能够准确描述冰形几何特征的网格拓扑结构。基于多块面搭接网格思路，根据结冰位置特点对计算区域进行了合理分区，形成了适用于背景飞机机翼结角状冰后流动分析的基本网格拓扑形式。在该网格拓扑基础上对背景飞机带冰构型复杂外形生成了高质量网格，见图3。为了说明网格疏密对数值模拟结果的影响，构建了密、中、疏3套网格进行无关性验证，以气动力量值作为网格无关判定依据。当半模网格总量约为3.0E+7时，数值计算结果保持较高精度，同时能够保证近壁面y+≦1。

图3 背景飞机整机网格、机翼前缘带冰网格及平尾舵偏网格

飞机半模的计算域为一个长方体，其尺寸为：上下边界距离(计算域的高度)为飞机长度的10倍，入口边界到出口边界的距离(计算域的长度)为飞机长度的20倍，展向方向距离(计算域的宽度)为飞机长度的15倍。全模网格可由半模对称得到。

3 背景飞机机翼结冰构型的复杂流场结构分析

将生成的网格输入到CFL3D开源流体力学计算程序中进行流场计算。计算域采用远场自由来流条件，对称面采用对称边界条件，物面采用绝热、无滑移和法向零压力梯度的固壁条件。

图4给出了Ma=0.2、平尾安装角无变化、无舵偏条件下背景飞机干净构型(Clean)、机翼前缘中度结冰构型(IceB)和重度结冰构型(IceA)翼面分离的流动形态。从图可知，背景飞机机翼前缘结冰后，翼面分离始发位置、分离沿展向和弦向的发展变化过程都与干净无冰状态存在本质区别。干净无冰状态在较大迎角下(α=12°)出现分离，且分离始发位置位于内翼，内外翼分离流动发展梯次较为明显，机翼具备良好的失速分离特性。而在重度结冰状态下，即使迎角较小(α=4°)，内外翼分离就几乎同时出现，始发位置均位于机翼前缘，且外翼的展向分离特征极为明显。随着迎角增加，分离流动的发展并不存在明显的展向梯次，而是沿弦向向后迅速推进，并且在较大迎角下，附着流动区域也显示出了较为明显的展向流动趋势，在α=8°时，翼面流动就几乎完全分离，失速分离特性相对干净无冰状态全面恶化。对于中度结冰情况而言，相同迎角下分离流动发展的趋势与重度结冰定性一致，但分离区域相对较小。失速分离特性较重度结冰没有定性差异。

图4 背景飞机翼面分离流动形态对比

图5给出了分别在Ma=0.2、Ma=0.3时，重度结冰条件下翼面分离流动演化的过程。可以看到在低速、小迎角条件下，马赫数改变对翼面分离形态的影响并不显著。但在失速点附近，马赫数的增加明显加剧了翼面流动分离的趋势。这是由于在失速点附近，随着来流马赫数增加，外部流动与分离泡回流区域之间的速度差异更大，流向剪切效应更强，分离泡结构能够承受剪切效应的迎角范围相对更小，宏观表现为流动分离的趋势更强。内翼分离流动沿弦向的扩展速显著增加，且无论在流动分离区域还是附着区域，外翼的展向横流效应均更加明显，但分离形态并未发生定性改变。

图5 马赫数变化对重度结冰条件下翼面分离流动演化过程的影响

图6给出了侧滑角影响下，背景飞机机翼/平尾带冰构型失速分离流动的特征(α表示迎角、β表示侧滑角)。从图6可以看到，在侧滑角的影响下，两侧机翼分离流动形态存在定性差异，加剧了结冰条件下横侧气动力的非线性趋势。大迎角下，迎风面机翼呈现弦向分离特征，背风面机翼呈现展向分离特征。迎风面和背风面机翼流动特征的差异导致滚转力矩非线性。

图6 大迎角机翼/平尾同时带冰构型不同迎角/侧滑角表面极限流线对比情况

图7以速度等值线的形式给出了重度和中度结冰条件下，机翼展向不同站位的空间流场基本结构。由图可知，对于重度结冰条件，在较小迎角下翼面上就形成了典型的分离泡结构，呈现典型的全翼展同时分离特征。随着迎角增加，再附位置不断后移，回流强度不断增加。在失速点附近，分离泡体积迅速膨胀，形成大范围回流区域，表明分离泡的弦向发展过程和再附决定了翼面的基本分离形态。其中内翼和外翼中段位置显示出了较高的回流强度，虽然内翼分离区域范围较外翼小，但回流强度却明显大于外翼。表明机翼/机身、机翼/短舱等部件之间的气动干扰削弱了分离泡的弦向扩展过程，同时抑制了横向流动，降低了当地流动分离的趋势。对于中度结冰情况，分离泡的基本形态与重度结冰一致，较小迎角下的回流区域大小与重度结冰差异并不明显，但在失速点附近，内翼分离流动强度低于重度结冰。

图7 背景飞机机翼重度和中度结冰翼面分离流动空间结构

4 背景飞机机翼结冰后的气动特性分析

根据背景飞机机翼重度结冰和中度结冰条件下纵向基本气动力计算所得数据，对结冰条件下背景飞机气动特性的变化规律进行了分析研究。如图8所示，为背景飞机在Ma=0.2，零舵偏，零安装角条件下机翼重度结冰气动力的计算结果。

图8 背景飞机机翼重度结冰气动力计算结果示例(Ma=0.2，零舵偏，零安装角)

低速条件下的机翼中度和重度结冰计算结果表明：对于背景飞机而言，机翼结冰主要影响飞机失速点附近的气动特性，将导致失速迎角普遍提前2°～4°不等，最大升力系数下降10%～20%左右，升力线斜率降低，纵向静稳定度下降，力矩非线性特征提前出现，之后低头力矩迅速增加。失速后普遍存在升力系数缓慢爬升的现象，这与回流区域的相对稳定发展有关，与国外相关参考文献[38～39]中得到的结论一致。在文献[38]中，NASA针对通用运输机模型(Generic Transport Model，GTM)在Langley实验室的12 inch低速风洞中开展了带冰条件下三维全机构型风洞测力试验。从试验结果来看，在机翼前缘带典型角状冰的前提下，失速点后方也出现了升力缓慢爬升的现象，与本文结论定性一致。其中GTM飞机模型与本文中的背景飞机模型气动布局较为相似。文献[39]中，NASA针对带典型角状冰条件下二维翼面进行了风洞测力试验，得到了类似结论。

计算结果还反映了结冰构型与干净无冰构型之间的气动特征差异明显大于重度结冰与中度结冰，这是由翼面分离流动之间的定性差异所决定的。

针对背景飞机机翼/平尾同时结冰构型，在不同迎角/侧滑角下对全机气动力和力矩特性进行了计算分析。通过与中国空气动力研究与发展中心的风洞试验数据进行对比，验证了数值方法刻画结冰后全机气动特性特别是失速特性的能力。如图9所示，为背景飞机机翼前缘和平尾前缘重度结冰条件下，在12°大侧滑角的风洞试验测力结果与CFD计算结果的对比。从图9可以看到，大侧滑条件下，计算得到的宏观气动力和力矩与试验结果吻合良好，能够反映失速之后全机气动特性的变化趋势。

图9 重度结冰条件下12°侧滑角CFD计算结果与风洞试验值对比情况(Ma=0.1，Re=0.51E+6)

图10给出了背景飞机机翼/平尾重度结冰构型下在0°、6°及12°侧滑角的风洞试验气动测力结果。可以看到，随着侧滑角增加，全机阻力上升，纵向静稳定度下降，但升力特性基本不变。过失速条件下，侧力、偏航力矩和滚转力矩量值均以近似线性形式降低。大侧滑角下，气动力和力矩的变化趋势与中等侧滑角一致，横侧气动力的非线性特征更加显著。

图10 重度结冰构型下不同侧滑角风洞试验测力结果(Ma=0.1，Re=0.51E+6)

5 结论

本文针对典型大型飞机构型，研究了机翼结冰条件下的复杂流场结构及气动特性变化情况。主要的研究内容及结论包括：

1)研究了机翼结冰条件下背景飞机的气动力特性，计算结果表明机翼结冰主要影响背景飞机失速点附近的气动特性，将导致失速迎角普遍提前2°～4°不等，最大升力系数下降10%～20%左右，升力线斜率降低，纵向静稳定度下降，力矩非线性特征也提前出现，全机失速特性全面恶化。

2)针对背景飞机机翼结冰条件下大迎角失速分离复杂空间流场结构及特性进行了研究。计算结果表明：结冰将导致翼面分离形态由内翼始发分离变为全翼展同时分离。分离始发大幅提前、梯次完全消失、流动非定常特征显著是全机失速特性恶化、气动边界缩小的直接原因。分离泡结构的弦向发展和再附过程决定了翼面的基本分离形态。

3)针对背景飞机全机结冰构型下的流场及气动特性相互耦合特征进行了研究。在不同迎角/侧滑角下对全机气动力和力矩特性进行了计算分析。随着侧滑角增加，全机阻力上升，纵向静稳定度下降，但升力特性基本不变。大侧滑角下，横侧气动力的非线性特征更加显著。