民用飞机机头冰脱落特性数值模拟

冯丽娟，谭广琨，张 强，易 贤

(1.上海飞机设计研究院总体气动部，上海 201210; 2.中国空气动力研究与发展中心空气动力学国家重点实验室，四川绵阳 621000)

民用飞机机头冰脱落特性数值模拟

冯丽娟1，*，谭广琨1，张 强1，易 贤2

(1.上海飞机设计研究院总体气动部，上海 201210; 2.中国空气动力研究与发展中心空气动力学国家重点实验室，四川绵阳 621000)

建立了带动力条件下飞机表面冰块脱落的数值计算方法，包括全机带动力的空气流场计算方法、冰块运动中气动力和力矩的确定方法以及冰块运动的六自由度方程数值求解方法。采用以上方法对某民用支线飞机巡航和进场状态下的机头冰脱落特性进行了模拟，对比了迎角对冰脱落的影响，分析了冰脱落的速度和轨迹关系，得到了冰块被吸入发动机的概率，为发动机吞冰设计、飞机冰脱落的适航符合性验证以及自然结冰试飞提供了较好的支持。

民用飞机;机头结冰;冰脱落;计算流体力学;数值模拟

0 引言

飞机结冰广泛存在于飞行过程中，是造成危害飞行安全的主要隐患之一［1-3］。飞机结冰主要发生在机翼和尾翼前缘、发动机唇口、风挡和机头等部位［4］。除了机头，以上部位通常会采取防除冰措施［5］。未采取防护措施的机头结冰后，如果发生冰脱落现象，会对飞行安全产生影响。例如，脱落冰块打到后机身或尾翼，可能损坏飞机部件，尤其当脱落冰块的质量大于发动机吞冰能力限制时，被发动机吸入后会打坏发动机叶片，造成发动机严重推力损失甚至停机等后果。如何表明在各种飞行条件下，机头上的结冰脱落后对飞机发动机以及机体的安全不造成影响，满足相关适航条款要求，对于该民机的适航取证是一个难点。如果在风洞中进行冰块投放试验，一方面试验费用较高，无法覆盖大量的飞行条件，同时，现阶段还不能做到完全真实模拟，尤其是无法在试验中评估发动机动力对冰脱落特性的影响。因此，有必要发展冰脱落研究的数值方法。

国外对于冰脱落研究的公开报道不多，主要是用数值手段［6-9］或实验［10］分析冰块的运动轨迹，研究通常采用简化的方块冰形，大部分研究未考虑飞机带动力的情况。国内近年来才有少量的研究。文献［11］研究了冰形外部压力和冰层与飞机蒙皮之间剪切应力对冰脱落的影响，探讨了冰脱落发生的原因和条件;文献［12］采用四自由度模型计算分析了不同尺寸、形状冰粒的脱落轨迹，未考虑三维流场和冰块运动的六自由度;文献［13］计算模拟了全机流场中冰块的六自由度运动，但由于没有考虑到发动机的动力影响，对近发动机区域脱落冰块的运动轨迹的分析与真实情况下脱落冰块的轨迹会存在较大差异。

本文建立了一种考虑飞机发动力影响的、基于六自由度运动方程的飞机冰脱落数值模拟方法，并对采用尾吊式发动机布局的某民用支线飞机的机头冰脱落特性进行了计算模拟，得到了冰块运动轨迹以及冰块被吸入发动机的概率。

1 冰脱落计算思路［13］

脱落冰块的运动是一个非定常过程，同时与飞机流场存在相互干扰作用，考虑到冰块尺寸远小于飞机，因此忽略冰块对飞机流场的干扰效果，只考虑飞机流场对冰块的作用，计算思路为:

(1)选取典型外形的冰块;

(2)计算冰块在各种姿态下的气动特性，建立冰块的气动数据库;

(3)采用CFD方法计算飞机三维流场，获得飞机流场分布;

(4)在上述工作基础上，求解刚体运动六自由度方程，获得脱落冰块的运动轨迹;

(5)根据轨迹计算结果，分析冰块被吸入发动机的概率。

2 流场模拟方法

通过求解附带k-ω SST湍流模型的RANS方程获得飞机绕流流场以及冰块随迎角和侧滑角变化的气动力数据。无量纲参数采用:密度、速度、压力、能量、温度 γ ( γ－1)、长度 L、时间L/和粘性系数。三维直角坐标系下守恒型RANS方程可以表达为:

采用LU-SGS格式的隐式时间推进法求解控制方程，对流项离散采用Jameson中心格式，粘性项采用中心离散［14］。

飞机流场计算考虑了带动力特性的复杂构型模拟，这样能够更为准确地模拟飞机实际的流场特征，使得发动机前方的冰脱落轨迹更加符合实际情况，提高了冰脱落轨迹计算的准确性。

典型的涡轮风扇发动机计算模型如图1所示，对发动机进排气性能利用特定的计算边界条件加以控制，其中风扇入口为进气口，风扇出口(外涵道)和主喷口(内涵道)为排气口。将进气口作为流场出口，给定质量流量，按准一维等熵流理论计算其他变量;排气口作为流场入口，给定流动方向、总温和总压。

图1 涡轮风扇发动机计算模型示意图Fig.1 Turbo fan engine simplified model

3 冰块运动轨迹计算方法

冰块运动轨迹通过数值求解刚体运动的六自由度方程获得［15］。在六自由度方程求解中，以飞机的体轴系为参考来描述脱落物的质心运动速度、质心位移，冰块脱落前的体轴系与飞机的体轴系之间存在俯仰角、偏航角和滚转角三个初始欧拉角。采用四阶龙格-库塔方法求解，具体步骤为:

(1)从脱落瞬间开始，在轨迹模拟的任一时刻，利用冰块的质心坐标从飞机流场数据中读取所在位置的气流速度，再根据冰块姿态求得相对来流迎角和侧滑角，进而从冰块气动数据库中插值得到气动力和力矩;

(2)以该时刻的气动力和力矩作为输入条件，求解刚体运动的六自由度方程，得到冰块在下一时刻的位置和姿态;

(3)根据冰块新的位置和姿态，进行下一时间步的计算，如此反复迭代，直到冰块运动至特定区域。

4 民用飞机机头冰脱落计算分析

4.1 脱落冰块外形的选取

飞机机头结冰外形通常呈圆盘状，如图2所示。

图2 飞机机头结冰外形Fig.2 Ice shape on nose

选取尺寸为1/3圆盘的扇形冰块表征机头圆盘冰破裂后形成的脱落冰块，如图3所示。

4.2 冰块气动数据库的建立

图3 扇形冰块Fig.3 Fan like ice shape

冰块气动数据库在冰块轨迹计算过程中提供气动力和力矩信息，包含三个笛卡尔坐标轴方向的气动力系数和力矩系数。图4为脱落冰块的计算网格，远场边界为10倍冰块长度。

图4 冰块空间网格及远场边界Fig.4 Mesh of ice and far field

图5为冰块气动力随迎角和侧滑角变化的数据库，其中横坐标为迎角，变化区间为(-180°，180°)，不同的曲线族按侧滑角分类，变化区间为(-90°，90°)。

4.3 计算条件

以某尾吊式布局的民用支线飞机作为研究对象，选用的扇形脱落冰块半径为0.12 m，厚度0.04 m，冰块的质量为0.454 kg。

机头脱落点位置:以机头驻点为圆心、半径为0.25 m的圆上选取均匀分布的6个点作为机头冰块的初始脱落位置，如图6所示。

图5 冰块气动数据库Fig.5 Ice slice aerodynamic database

图6 冰块在机头初始脱落位置Fig.6 Initial position of ice slice on nose

冰块初始姿态表示为不同的欧拉角(俯仰角、偏航角、滚转角)的组合。

飞机的计算网格采用结构网格，物面第一层为平均弦长的10-5，远场边界为10倍飞机长度，网格布局注意机身的保形，同时在流场变化较快的地方进行局部加密。

计算了巡航和进场两种状态，具体计算状态见表1。

表1 飞机飞行状态Table 1 Flight condition

4.4 巡航状态冰脱落特性计算分析

图7为巡航状态不同迎角下脱落冰块轨迹的前侧视图，曲线颜色代表冰块质心的移动速度，由于马赫数较大，脱落冰块受气动力影响较重力更明显，冰块群整体位于发动机上侧，冰块运动到发动机唇口截面时的最大速度达到38 m/s。图8为飞机左侧半模冰块轨迹与发动机唇口轴向截面的交点分布图，由图可知，随着迎角的增大，轨迹交点整体上移，其中4.2°迎角时吸入发动机的点数最多，概率为4.32%，随着迎角增大，吸入发动机的冰块数明显减少，迎角6.6°时，所计算的脱落轨迹中，只有一条进入发动机唇口，概率为0.62%，而迎角8.4°则未发现吸入发动机的冰块。

图7 巡航构型冰脱落轨迹Fig.7 Ice shedding trajectories of cruise case

4.5 进场状态冰脱落特性计算分析

图9为进场状态不同迎角下脱落冰块轨迹的前侧视图，与巡航状态相比，由于马赫数明显减小，气动力的作用减弱，因此脱落轨迹向下偏移较多，冰块运动到发动机唇口截面时的最大速度为30 m/s。图10为飞机左侧半模冰块轨迹与发动机唇口轴向截面的交点分布图，在4.3°迎角时轨迹交点均位于发动机下部，没有冰块打入发动机，随着迎角增大，冰块轨迹上扬，但吸入发动机的冰块仍然很少，迎角为7.3°和10.3°时，分别只有一条轨迹打入发动机，概率为0.62%。

图8 巡航构型冰脱落轨迹与唇口轴向截面的交点Fig.8 Footprint of shedding trajectories on engine plane of cruise case

图9 进场构型冰脱落轨迹Fig.9 Ice shedding trajectories of approach case

图10 巡航构型冰脱落轨迹与唇口轴向截面的交点Fig.10 Footprint of shedding trajectories on engine plane of approach case

5 结论

本文建立了冰脱落计算分析方法，并针对某型民用支线飞机带动力条件下机头冰脱落问题进行了计算分析，得到如下结论:

(1)冰块运动轨迹受重力，来流速度和迎角的影响，正迎角下气动力使轨迹上扬，重力则使轨迹下偏，随着来流速度以及迎角的增大，气动力作用影响加大，冰块轨迹升高;

(2)通过对该民机机头冰脱落的轨迹被吸入发动机的概率分析来看，带动力条件下，对于本文所计算的状态，该尾吊式发动机吸入机头脱落冰的概率较低。

(3)本文的数值方法可以为民用飞机冰脱落适航符合性验证提供支持，也可作为评估自然结冰试飞中冰脱落对发动机安全影响的手段。

［1］ Cebeci T，Kafyeke F.Aircraft icing［J］.Annu.Rev.Fluid Mech，2003，35:11-21.

［2］ Kind R J，Potapczuk M G.Experimental and computational simulation of in-flight icing phenomena［J］.Progress in Aerospace Science，1998，34:275-345.

［3］ Frank T L，Abdollah K.Effects of ice accretions on aircraft aer-odynamics［J］.Progress in Aerospace Sciences，2001，37(8):669-767.

［4］ Yi Xian，Zhu Guolin，et al.Numerically simulation of ice accretion on airfoil and research of icing wind tunnel test scale rule［J］.Acta Aerodynamica Sinica，2002，20(4):428-433.(in Chinese)

易贤，朱国林，等.飞机积冰的数值计算与积冰试验相似准则研究［J］.空气动力学学报，2002，20(4):428-433.

［5］ Qiu Xiegang，Han Fenghua.Aircraft anti-icing system［M］.Beijing:Aviation Professional Teaching Material Compliance Group，1985.(in Chinese)

裘燮刚，韩凤华.飞机防冰系统［M］.北京:航空专业教材编审组，1985.

［6］ Luis C de Castro Santos，Ramon Papa，Marcello do Areal S.Ferrari.A simulation model for ice impact risk evaluation［R］.AIAA 2003-0030.

［7］ Wright W B，Keith T G，DeWitt K J.Numerical simulation of icing，deicing and shedding［R］.AIAA-91-0665，1991.

［8］ Michael Papadakis，Hsiung-Wei Yeong，Ian G Suares，et al.Experimental and computational investigation of ice shedding from aircraft surfaces［R］.AIAA 2006-1010，2006.

［9］ Hsiung-Wei Yeong，Michael Papadakis，Koji Shimoi.Ice trajectories and Monte Carlo analyses for a business jet［R］.AIAA 2009-3973.

［10］ Michael Papadakis，Hsiung-Wei Yeong，Koji Shimoi，et al.Ice shedding experiments with simulated ice shapes［R］.AIAA 2009-3972.

［11］Xiao Chunhua，Gui Yewei，Yi Xian，et al.Pressure distribution and stress calculation of glaze ice on iced airfoil［J］.Journal of Aerospace Power，2009，24(7):1457-1463.(in Chinese)

肖春华，桂业伟，易贤，等.结冰翼型表面明冰的压力分布和应力计算［J］.航空动力学报，2009，24(7):1457-1463.

［12］ Tan Meijing，Sun Zhiguo，Zhu Chunling，et al.Four degree of Freedom simulation of ice shedding for aircraft icing［J］.Aeronautical Computing Technique，2012，42(2):100-104.(in Chinese)

檀妹静，孙志国，朱春玲.四自由度飞机结冰脱落模拟研究［J］.航空计算技术，2012，42(2):100-104.

［13］Wang Jiantao，Yi Xian，Xiao Zhongyun，et al.Numerical simulation of ice shedding from ARJ21-700［J］.Acta Aerodynamica Sinica，2013，31(4):430-436.(in Chinese)

王建涛，易贤，肖中云，等.ARJ21-700飞机冰脱落数值模拟［J］.空气动力学学报，2013，31(4):430-436.

［14］ Mou Bin，Xiao Zhongyun，Liu Gang，et al.Low-speed preconditioning based on Roe scheme［J］.Acta Aerodynamica Sinica，2013，31(4):430-436.(in Chinese)

牟斌，肖中云，刘刚，等.基于Roe格式的低速预处理方法研究［J］.空气动力学学报，2010，28(2):125-131.

［15］Xu Mingyou，Ding Songbin.Flight dynamic［M］.Beijing:Science Press，2003.(in Chinese)

徐明友，丁松滨.飞行动力学［M］.北京:科学出版社，2003.

Numerical simulation of ice shedding from civil airplane nose

Feng Lijuan1，*，Tan Guangkun1，Zhang Qiang1，Yi Xian2
(1.General and Aerodynamics Department，Shanghai Aircraft Design Research Institute，Shanghai 201210，China; 2.State Key Laboratory of Aerodynamics，China Aerodynamics Research and Development Center，Mianyang 621000，China)

The shedding ice blocks from a civil airplane nose seriously threat the safety of the civil airplane，and as result it is very important to thoroughly investigate the shedding characteristics.A numerical simulation method of ice shedding was proposed to simulate the trajectories of the shedding ice blocks from airplane nose by establishing the force and moment database of the shedding ice blocks using CFD(Computational Fluid Dynamics)tool，and solving the six degree-of-freedom equations numerically，the effect of engine power was also considered.Employing this method，the trajectories of the shedding ice from a civil airplane nose with aft-mounted engines were simulated，the effect of angle of attack of airplane on the trajectories of shedding ice were studied，the relationship among the velocity，the trajectories of the shedding ice and the probability of shedding ice blocks ingested into the engines was analyzed.The research work is of help for the design of engine ingestion ability，the compliance work of nose ice shedding and the nature icing flight test of the civil aircraft.

civil airplane;airplane nose icing;ice shedding;computational fluid dynamics;numerical simulation

V244.1

A

10.7638/kqdlxxb-2014.0005

0258-1825(2015)04-0536-06

2014-02-17;

2014-05-14

国家自然科学基金(11172314)

冯丽娟*(1976-)，女，硕士，高工，主要从事民机结冰适航取证方面的研究工作.E-mail:fenglijuan@comac.cc

冯丽娟，谭广琨，张强，等.民用飞机机头冰脱落特性数值模拟［J］.空气动力学学报，2015，33(4):536-541.

10.7638/kqdlxxb-2014.0005 Feng L J，Tan G K，Zhang Q，et al.Numerical simulation of ice shedding from civil airplane nose［J］.Acta Aerodynamica Sinica，2015，33(4):536-541.