1 引言

1990年以来，在海拔6706m以上的高空已发生超过100起飞机发动机功率损失事件^{[1, 2]}。2006年，一架卡塔尔航空公司的空客A330航班在6437m高空经历雷暴天气，双发同时停车^[3]；2013年，GEnx发动机的运营商报告了9起12km高空发动机推力损失事件^[4]。经调查发现，造成以上故障发生的原因是发动机结冰，但并不是传统上的过冷水结冰。传统上发动机结冰是过冷水造成的外部结冰，极少发生在6km以上的高空，因为高空只含有极少的过冷水滴。研究者从故障发生时飞行员观察到的现象、飞行路线、发生的时间、地点、海拔等方面入手，推测发动机高空功率损失频发的原因。最终由Manson等^[1]，Grzych等^[5]确定上述发动机在高空的功率损失是冰晶在发动机内部结冰导致的。图 1为冰晶在高压压气机静子叶片上结成的冰层。冰晶内部结冰与过冷水滴外部结冰有很大的不同^[6~8]。为方便区别，本文将过冷水导致的结冰称为过冷水结冰，冰晶内部结冰称为冰晶结冰。

冰晶结冰可导致压气机气流通道截面积减小、内部流场紊乱，进而引起发动机功率下降、掉转速（Engine rollback）等现象。积冰脱落可造成低高压压气机叶片断裂、燃烧室熄火等一系列严重后果^{[1, 4, 6, 8, 9]}，图 2为积冰脱落打坏的叶片，日益引起发动机行业各方的重视。而且冰晶结冰的故障在事前很难监测。引起发动机内部结冰的冰晶直径约为40μm，尺寸小于普通冰晶，肉眼无法辨识，气象雷达更无法探测^[4]。

认识到问题的严重性，2006年FAA第一次发出警告，要求飞行员提高意识，避免在强风暴区域飞行，并制定相应的应急处理措施。同时要求发动机供应商设计发动机时必须考虑冰晶结冰问题。随后，FAA在2014年11月公布的适航条例中补充了冰晶结冰气象设计标准（详见FAR33 Appendix D^{[10, 11]}）。目前三大航空发动机公司也积极采取相关措施，如GE通过升级发动机全权数字控制系统（FADEC），探测吸入发动机内部的冰晶，根据探测结果打开放气活门，排出冰晶^{[3, 4]}。但由于还不清楚冰晶结冰的物理机理，这些举措并没有从根本上解决冰晶结冰的问题。

航空发动机冰晶结冰问题已受到国内外重点关注，目前国外研究主要集中在冰晶结冰地面模拟试验以及计算模拟方面；国内仅有少量学者进行计算模拟方面的研究，暂无试验方面的研究。由于冰晶结冰涉及传热、气液固相变、冰晶破碎以及侵蚀等一系列复杂问题，进展相对缓慢，多数研究成果只停留在定性层面。鉴于目前研究成果较分散，并且尚无冰晶结冰方面的综述性文章发表，本文在调研国内外文献资料的基础上，对冰晶结冰的研究进展进行了系统性的总结，并讨论了该领域有待进一步解决的方向性问题。

2 与过冷水结冰的区别

过冷水结冰的研究相对成熟，而冰晶结冰则刚刚起步。为了更好地理解冰晶结冰的现象和机理，有必要将两者进行对比分析。因过冷水滴撞击引起的航空发动机结冰多发生在进气道前缘^[10]、整流罩^{[12, 13]}、支板^{[13, 14]}以及导流叶片^[13]等部位；而根据冰晶结冰事件的调查发现，从整个低压压气机一直到高压压气机一级静子叶片位置均有结冰现象发生。图 3为标有不同结冰位置的典型涡扇发动机结构简图。

由于温度的限制，过冷水滴基本分布在7km以下的大气环境中。因此，过冷水结冰高度范围通常在7km以下。此外，由于过冷水结冰需要低于冰点的环境以释放潜热^[6]，而当飞行速度高于600km/h时，会由于气动加热导致机体表面温度高于0℃，因此过冷水结冰时的飞行速度一般低于600km/h。虽然目前飞机的巡航高度一般均已高于过冷水结冰高度，但在起飞、爬升、空中盘旋和下滑阶段，仍可能会发生比较严重的结冰现象，严重影响飞行安全。目前已研制出多种防除冰方式，如机械、液体、气热和电热等。

研究表明^{[1, 2, 4, 5, 15]}，冰晶结冰多发生在春夏季太平洋西海岸及东南亚地区海拔7km以上的高空，这些地区由于处在热带和亚热带，春夏季多发生强对流雷暴天气。春夏季空气温度的突然升高，导致热空气上升，在对流层遇冷形成降雨，在较强对流天气下，潮湿的空气进一步上升至对流层顶部，上升过程中，温度进一步降低，从而形成微小的冰粒即冰晶。因此故障发生地区，多形成低空有较大降雨，高空含有大面积冰晶的垂直空间分布，如图 4^[16]所示。在风的作用下冰晶区域向外分散，NASA调查发现，冰晶从主风暴区可飘移到56km以外。当飞机经过雷暴区时，由于冰晶尺寸过小且分散，气象雷达无法辨识，只能捕捉到低空发生降雨的区域，飞机为躲避降雨区域，在降雨区域上空进行巡航（如图 4），从而吸入大量的冰晶，冰晶进入到发动机内部冰点以上的环境，部分融化产生冰水混合相，湿润内部结构表面，随后吸入的冰晶粘附在这些湿润的表面并进行冷却，随着冰晶的大量吸入，从而逐渐累积成冰^{[1, 6]}。

综上，将航空发动机过冷水结冰与冰晶结冰的各项差异总结如表 1所示。由此可以看出，冰晶结冰与过冷水结冰存在显著差异。

3 冰晶结冰机理研究

对冰晶结冰机理的研究，是从根本上解决其引发航空发动机故障的关键。一般认为，总温低于10℃并有可见水汽时发生结冰^[9]。在压气机前段有一段总温介于0~10℃的区域^[9]。冰晶进入该区域会部分融化，产生冰水混合相，润湿压气机叶片。随后吸入的冰晶不断粘附在湿润的叶片表面，使得叶片表面温度逐渐降低，形成积冰；同时，由于冰晶的融化吸热也使得该段区域环境温度降低，引起与其相邻区域温度降低，导致低温区域向后延伸，并随着冰晶吸入量的逐渐增大，向压气机内部的延伸越深入，从而结冰的区域也越深入^{[1, 6, 9]}。

这种对冰晶结冰机理的描述仅停留在定性层面，无法精准把握结冰过程的影响因素及其参数确定，如冰水含量、总温、总压等，因此必须开展大量的定量研究。目前主要的研究方式集中在地面试验模拟和仿真计算模拟。

3.1 地面试验模拟研究

由于航空发动机冰晶结冰是近年来才意识到的问题，物理现象复杂，人们对于其结冰机理的认识方面有很大的不足，为此，国外众多学者试图通过试验进行逐步探索。然而，由于发动机内部环境的恶劣，采用传统的探测手段几乎不可能，因此对结冰机理的理解带来巨大挑战；另外，高海拔地区的冰晶尺寸在μm量级，设计生成冰晶的试验装置也是一大难点。这些苛刻的要求都对试验模拟提出了巨大的挑战。

自2009年至今，美国国家航空航天局（NASA）与加拿大国家研究委员会（NRC）在结冰机理的探索上做了一系列重要工作，为冰晶结冰问题的试验研究做出了巨大贡献。这些研究主要围绕NRC高空试验研究设备（Research Altitude Test Facility，RATFac）可视化试验台和NASA推进系统实验室（Propulsions Systems Lab，PSL）全尺寸发动机冰晶结冰试验装置展开。

3.1.1 局部流动模拟地面试验设备

RATFac可视化试验台，属于局部流动模拟，该装置通过研磨装置产生平均直径为45~200μm的冰粒。试验时，直接将冰粒喷射到温暖的环境中，以模拟吸入发动机内部时的环境。其中直径较小的冰粒首先完全融化，而直径较大的冰粒在到达结冰位置前，仍部分保持固态^[17]。试验装置如图 5所示，详见文献[6, 18]。

3.1.2 全尺寸流动模拟地面试验设备

为研究较小尺寸（< 40μm）冰粒的结冰机理，并试图在发动机外部产生关键的冰水混合相，2014年PSL进行了冰晶吸入试验模拟^[19]。与RATFac不同，这是在全尺寸涡扇发动机运转状态下的流动模拟。通过液态水喷嘴装置喷出液滴，随后通过对流传热及蒸发冷却，使得在进入发动机之前，小液滴凝结成冰晶，大液滴保持液态或部分液态，以此产生冰水混合相。该装置可产生尺寸分布在15μm左右的冰粒^[20]。装置简图如图 6~7所示。

3.1.3 冰晶结冰的必要条件

如何在冰点以上的环境下模拟出冰晶结冰的现象，进而探究冰晶结冰的必要条件是试验模拟的首要问题。2010年Mason等^[6]首次在RATFac可视化试验台上模拟出了在高于0℃的环境下冰晶结冰与脱落的现象，如图 8所示为结冰案例。

随后，进一步探究出冰晶结冰过程中，冰水混合相的产生是冰晶结冰的必要条件^[6]。此外，混合相中若没有足够的冰水含量（Ice water content，IWC），则不能使发动机内部结构表面温度降低到冰点以下；反之，若没有足够的液态水含量（Liquid water content，LWC），则不能粘附足够的冰晶来换热。因此结冰的程度，取决于发动机内部冰晶融化量的大小，并存在最佳结冰阶段^[22]。若用LWC占总水含量（Total water content，TWC）的比例，来表征冰晶融化量的大小，则结冰量与融化量的关系如图 9所示，其中TWC=LWC+IWC。

3.1.4 结冰稳固的判断依据

结冰过程中，由于受不同结冰条件的影响，所生成的冰不但在外观形貌和内部结构上有所不同，而且在密度、导热系数、比热容、冰相的内部应力特性及冰在物面的黏附特征等均存在较大差异。结冰是否稳固是影响结冰特性及防/除冰能量需求准确预测的重要指标之一，因此对结冰稳固影响因素的研究是近年来研究者另一个关注的课题。

Struk等^{[22, 23]}为探究在相对稳定的等熵流动下的结冰情况，在Manson试验装置的基础上，改为对称楔形叶型。在来流马赫数为0.2~0.3，总温大于0℃，总压45~93kPa内进行了一系列试验，发现在45kPa时，叶片前缘形成稳固的明冰（Glaze ice），如图 10所示；而在93kPa时，形成易脱落的疏松冰，如图 11。经分析，发现在45kPa的试验中湿球温度（Wet-bulb temperature，T_wb）T_wb < 0℃；而在93kPa的试验中T_wb > 0℃。由此Struk等^[22]认为可以用湿球温度作为结冰是否稳固的判断依据。随后，经Struk等^{[22, 23]}，Currie等^[23]，Knezevici等^[24]的研究，最终确定了只有在T_wb < 0℃的情况下，0℃以上的表面才会形成稳固的结冰区域。

湿球温度是指在绝热条件下，几何表面上的水蒸发所需的潜热完全来自于周围湿空气温度降低所放出的显热，当系统中空气达饱和状态且系统达到热平衡时系统的温度^{[22, 23]}。可见蒸发速率直接影响湿球温度的高低，而蒸发速率又随着压力的降低（海拔升高）而加快^[22]，因此在其他条件相同的情况下，湿球温度随压力的降低而降低，即低压更易诱使T_wb达到0℃以下，这就解释了上述45kPa相比93kPa形成稳固结冰的原因。此外，研究表明总温和湿度对T_wb的影响较小^[23]。

由Manson等^{[1, 2, 5, 6, 16]}提出的冰晶结冰过程可知，吸入的冰晶或冰水混合相越多，结冰量越大，因此无论T_wb在何种条件下，结冰量与IWC，LWC，马赫数等基本均呈正相关，只是结冰是否稳固的区别^{[22, 23]}。另外，考虑到撞击侵蚀的影响，大尺寸冰晶相比小尺寸冰晶结冰量降低^[24]。事实表明^[23]，T_wb < 0℃时，随IWC的增大，由于侵蚀等原因，粗糙稳固冰质向光滑稳固冰质转变^[25~27]；而T_wb > 0℃时，均形成易脱落的疏松冰。值得指出的是，若LWC/TWC较大，由于没有足够的冰晶与壁面换热，使温度降到冰点以下，即使T_wb < 0℃，也有可能产生疏松的结冰区域。

2016年，Struk等^{[20, 28]}在PSL试验设备上，研究了全液滴状态到全冰晶状态下的一系列不同情况，进一步发现结冰时的湿球温度均基本保持不变或只有轻微增长。除以上方面的研究外，由于发动机掉转速的突然性严重威胁飞行安全，文献[29]深入研究了TWC、冰粒尺寸、风扇转速、热源以及海拔对发动机掉转速的影响；文献[30]发现，当冰水流量与空气流量之比小于某临界值时，不会发生掉转速现象。另外，由于侵蚀对冰形及结冰区域有重要影响，文献[31, 32]研究了冰晶的撞击特性及其影响因素。

截止目前，上述研究是冰晶结冰试验方面的重大成果。以上试验成果均在不同程度上揭示了冰晶结冰的物理机制，但由于结冰过程涉及因素众多以及气液固相变机理的复杂，相关研究还需进一步深入。

3.2 计算模拟研究

由3.1节可知，用于模拟航空发动机冰晶结冰的试验设备异常复杂，需要消耗大量人力、物力、财力，试验结果也只是某些特定条件下的情况，没有普适性。而影响冰晶结冰的因素众多，进行各种条件不同状态下的试验模拟是不现实的，因此对计算模拟的研究很有必要。

计算模拟离不开试验数据的校核与修正，由于试验数据的匮乏，只能对比NSAS和NRC所做的试验数据，因此冰晶结冰仿真计算方面相对研究较少。虽然目前国外已有比较完善的结冰软件，如美国的LEWICE^[33]，GlennICE^[34]；加拿大的FEASAP；法国的ONERA^[35]等。但这些软件主要是针对过冷水结冰而设计开发的，因此国外对冰晶结冰的计算研究主要是针对这些软件的改进，使其能够模拟冰晶结冰的情况。国内尚缺乏专门的结冰软件，部分学者利用商业软件Fluent结合用户自定义函数（UDF）进行冰晶结冰模拟^[8]。

冰晶结冰涉及传质、传热、气液固相变、接触、碰撞等一系列问题，物理模型复杂，影响因素众多，图 12为影响冰晶结冰的物理机制。

根据冰晶吸入发动机中经历的过程，可将冰晶结冰过程的模拟大致简化为四个部分：空气流场计算、冰晶运动轨迹和撞击特性计算、结冰计算、冰形计算。

3.2.1 空气流场计算

冰晶在流场中的运动和撞击特性很大程度上决定于空气绕流流场的分布。文献[37]验证了模拟冰晶吸入时空气流与冰晶之间的耦合可忽略，因此计算空气流场时可忽略冰晶的影响。目前空气流场的计算方法主要有两种：（1）面元法，主要用于求解势流模型。（2）利用空间网格求解Euler或Navier-stokes方程^{[10, 38, 39]}。前者主要用于二维计算，求解简单速度快，但无法处理复杂几何形状的流场计算问题；后者由于计算机和CFD的发展，求解雷诺平均的Navier-stokes方程已不再困难^{[10, 38, 39]}。采用CFD方法计算空气流场后，得到空气速度、压力及温度等参数的分布。

3.2.2 冰晶运动轨迹和撞击特性计算

在空气流场计算结果的基础上，针对冰晶建立相应的控制方程，进行运动轨迹和撞击特性计算，获得结构表面的冰晶收集率与水收集率。主要计算方法有两种：Lagrange法和Euler法。Lagrange法根据牛顿第二定律对单个冰粒进行受力分析，建立其运动方程，以此来跟踪冰粒的运动轨迹；Euler法则是将冰晶视作连续相，引入容积分数，建立其连续方程和动量方程^{[10, 38, 39]}。

由于冰晶运动过程中气液固相的转换机理以及冰晶的撞击行为是理解冰晶结冰机制的关键，因此国外学者多用Lagrange法进行计算。冰晶运动轨迹和撞击特性计算需要考虑不同形状和孔隙率的影响、运动过程中冰晶与空气之间的传质传热、相态的变化、以及撞击壁面或水膜时冰晶的粘附、破碎、弹射、飞溅、侵蚀等现象的发生^[35]。因此，学者在这方面的研究方向多致力于建立能描述以上现象的数学模型。

2011年，Lou等^[40]在冰晶为球体以及传热系数与半径成反比的假设下，对冰晶吸入过程中提出了三种包含相变的模型，即裸冰模型（图 13），水膜覆盖模型（图 14）以及疏松孔质模型。建立了前两种直径随时间变化的数学模型，得到进入燃烧室之前冰晶未完全融化的临界直径。由于第三种模型需要考虑空隙率等因素的影响，没有建立相应数学模型。

文献[8]详细论述了冰水混合相的运动轨迹与考虑球形冰晶的破碎、弹射、粘附、飞溅在内的撞击模型；文献[35]系统地论述了冰晶运动轨迹、撞击、结冰计算所需的数学模型，在Mason，Wright，Lou等的基础上通过定义颗粒球度（Particle sphericity）来模拟非球形状态下的几何模型，相应的换热模型和相变模型；通过定义沉积概率来考虑冰晶粘附的现象；通过定义冰晶法向动能与表面能之比来模拟冰晶弹射和破碎的现象。文献[33]使用LEWICE3D软件中的过冷大水滴（Supercooled large droplet，SLD）模型模拟冰晶撞击过程中的弹射和破碎，并与NRC试验结果对比；文献[37]在GlennICE软件中添加了考虑冰粒吸入过程中与空气之间的换热及温度的变化，以此模拟吸入过程中的相变，并考虑了非球状冰粒的影响，利用回弹系数表征冰晶粘附与弹开量的多少，利用低速下的侵蚀经验公式模拟侵蚀效果，但由于试验数据的限制，结果没有进行验证。

虽然运动轨迹和撞击特性计算方面做了很多研究，但由于试验数据的缺乏，目前均是建立在众多简化假设的基础上，与真实情况还有较大差距，并且无法同时模拟撞击过程中冰晶的粘附、破碎、弹射、飞溅、侵蚀等现象。

3.2.3 结冰计算

基于空气流场及撞击表面的冰晶和水的收集率，根据守恒定理求解撞击表面的结冰热力学模型，得到表面温度、结冰量和溢流分布。目前绝大多数软件均是基于Messinger结冰热力学模型，该计算模型最早用来分析过冷水绝热结冰表面的平衡温度，忽略冰层与结构表面之间的传热，该模型综合考虑了结冰过程涉及的主要能量项，通过结冰表面的热量收支平衡来建立结冰热过程的描述，但未考虑冰层及水膜内部的温度梯度分布及传热的非稳态特征^{[10, 38, 39, 41]}。研究表明，Messinger模型可较好模拟过冷水导致的毛冰（Rime ice）的生长，而对于明冰、混合冰，由于溢流水等影响，模拟结果误差较大^[39]。为此，不少学者在NASA和NRC试验的基础上提出了以Messinger结冰热力学模型为基础的冰晶结冰模型改进方法。

由RATFac的一系列试验可知，冰晶导致的结冰类型按湿球温度大致可分为两大类，即T_wb < 0℃时附着稳固的明冰和T_wb > 0℃时易脱落的疏松冰。2014年，Tsao等^[42]在前人过冷水结冰模型的基础上，分别建立了上述两种情况下的冰晶结冰热力学模型来定量描述发动机内部未加热的表面上冰晶结冰的过程。

Nilamdeen等^[43]扩展了用于模拟过冷水结冰的Shallow-water模型，模拟了冰水混合相在未加热不旋转叶型上的结冰；文献[8, 35]在经典Messinger结冰热力学模型的基础上，考虑了冰晶的影响，并通过相关文献试验结果进行对比与修正；文献[44]建立了考虑冰晶在内的疏松孔质的积冰模型，并利用OpenFOAM开源平台进行了模拟。

同运动轨迹与撞击计算类似，虽然结冰计算方面进行了大量研究，但是由于试验数据的限制，目前基本未考虑撞击侵蚀对结冰的影响。

3.2.4 冰形计算

根据结冰计算得到各控制体各时间步的结冰量，再选用合适的冰密度模型，如Macklin经验模型，Jones经验模型等^[10]，即可得到各控制体各时间步的结冰厚度，在原几何表面不断叠加即可得到冰形轮廓^{[10, 38, 39]}。其中每叠加一次，都要重新划分网格，以考虑冰形对流场的影响，再进行下一个时间步的计算，直到设定的求解时间为止^{[10, 38]}。值得指出的是，目前多数模拟为简化计算，忽略了冰形对流场的影响^[38]。

3.2.5 结冰后发动机性能计算

除了对结冰过程的模拟外，文献[45, 46]模拟了发动机掉转速的现象，并提出通过结冰后监测发动机性能的改变来判断发动机内部是否结冰。2015年，Kundu^[47]通过一维瞬态模型模拟气体连续相，拉格朗日模型模拟冰粒离散相，详细论述了内部结冰对含有十级转子的压气机性能的影响。发现随着流量的增加，喘振裕度下降，冰粒尺寸越小，压气机喘振裕度下降得越多。

4 研究展望

冰晶内部结冰近年来受到国内外的广泛关注。随着研究的不断深入，在理论方面初步提出了冰晶结冰的机理；试验方面证实了冰晶结冰现象的发生并找到了一些关键性的表征参数；数值计算方面在一定的假设前提下，模拟出了冰晶结冰的过程，但由于试验数据的限制，很大程度上制约了计算模拟的发展。我国在冰晶结冰方面几乎空白，基本仍停留在过冷小液滴的结/防冰计算，因而应当加强冰晶结冰方面的研究。纵观国内外已有成果，冰晶结冰在试验与计算方面尚有众多基础性问题值得深入细致地研究：

（1）冰晶结冰试验相似准则研究

由于试验条件的限制，针对设计的冰晶结冰气象条件使用全尺寸模型进行试验往往是不现实的。因此需要进行相应的相似准则研究，如几何相似、流场相似、冰晶运动轨迹与撞击特性相似、结冰表面热力学相似等等，从而为缩比试验提供依据。

（2）冰晶运动与撞击模型研究

冰晶运动过程中，伴随着与空气间的传质传热，目前的计算模拟基本忽略了运动过程中物态的变化；撞击到壁面时，可发生弹射、破碎、粘附、飞溅、侵蚀等现象。目前的撞击模型由于试验数据限制，仅考虑其中的部分项，且多数尚未考虑后续吸入的冰晶对已结冰层侵蚀的影响。因此对冰晶运动与撞击过程中现象的完整描述，仍是需要深入研究的关键方向。

（3）冰晶结冰热力学模型研究

结冰热力学模型直接决定了冰形模拟的准确性，冰晶结冰的过程其实是包含气液固相变的移动边界过程，另外，由于气动剪力、重力、表面张力的共同作用，使得结冰表面水膜的流动特性成为影响冰层生长速率的重要因素。虽然不少研究者在试验的基础上提出了基于Messinger结冰热力学模型的改进方法，但多数研究仅限于对热量项在“量”上的修正^[7]，对于冰晶融化冷却与结冰过程中冰晶/水膜与结构表面之间的换热、冰层与水膜内部的温度梯度、水膜的流动特性尚未考虑，因而不能有效表征冰晶结冰过程的传热相变机制。

（4）冰晶结冰微观与宏观多尺度研究

目前的冰晶结冰模拟主要停留在冰形宏观形貌的预测。然而冰晶结冰是一个分子尺度、冰晶尺度、冰形尺度等多尺度问题。随着结冰精细化预测趋势的发展，需要研究微观结构特征与宏观形貌特征间的定量关系，从而能够与宏观预测方法建立有机联系，建立更精细的预测方法。

（5）冰晶结冰防/除冰方法研究

利用中压压气机放气，将脱落的积冰排出机体外是目前GE公司采用的除冰方法，但在巡航时，由于性能和耗油经济性等因素，放气活门处于关闭状态，此时积冰脱落将有很大概率打伤高压压气机叶片，导致燃烧室熄火等。因此在理解冰晶结冰机理的基础上，发展新的防除冰方法，例如在结冰位置增加加热装置、采用超疏水材料涂层等，对于飞行安全至关重要。

5 结束语

随着航空技术的发展及未来飞行安全性能指标的不断提高，深入掌握冰晶结冰过程的热力学、冰晶撞击动力学及空气动力学等综合作用下的复杂物理机制，突破现有认知能力的局限性，深刻把握冰晶结冰的科学本质，对建立结冰条件下相应的科学防护手段和操纵应对策略，具有重要的科学意义。