航空发动机冰晶结冰研究进展

苏龙伟，申世才，田晓平，吴 悠

（中国飞行试验研究院，西安 710089）

0 引言

航空发动机是飞行装置的核心部件，其安全问题对飞行十分重要。在20 世纪90 年代中期，几架商用飞机的喷气发动机在高空经历了频繁的功率损失，引起了航空科研人员的普遍重视。2013 年，GE 公司发生了9 次由吸入冰晶导致的发动机推力损失事件[1]，如俄罗斯1 架波音747-8HVF 货机在海拔12000 m 时误入含有冰晶的空域，发动机吸入冰晶造成核心机结冰，导致2号发动机喘振停车，1号发动机推力减小。

Mason 等[2]通过对20 世纪90 年代以来在海拔6 km 以上的高空发生的超过100 起不明原因的发动机功率损失事件进行分析，选取了46 个典型事件，重点对1架通勤运输飞机发动机故障和1架大型运输机推力损失事件进行分析，最终确定冰晶在压气机内结冰是导致该类故障的主要原因。

2014 年美国联邦航空局FAA[3]、欧洲航空安全局EASA[4]先后将航空发动机冰晶结冰纳入适航要求（FAR33 部附录D）。中国大型民用客机C919 进行自然结冰试飞取证时也曾遭遇冰晶结冰气象，给飞行安全造成隐患。在中国某民用大涵道比发动机在后续适航取证时，已明确要求进行冰晶结冰验证，因此有必要充分认识和研究冰晶结冰的内在机理，为中国后续冰晶结冰试航取证提供支持。

本文整理了近年来冰晶结冰的相关文献，对比了冰晶结冰与传统过冷水滴结冰的区别，重点总结了国内外在地面试验和数值模拟方面取得的重要成果，并对未来研究趋势进行展望，以期为中国冰晶结冰风洞研究和数值模拟研究提供思路。

1 冰晶结冰与过冷水滴结冰的区别

早在20 世纪四五十年代，国外就已经针对过冷水滴结冰开展了相关研究，制定了飞机自然环境结冰飞行标准，并逐步应用到飞机的设计当中[5]，而冰晶结冰则是在2006 年才开始受到航空界关注[1，6]，研究相对匮乏。冰晶结冰和过冷水滴结冰在结冰位置、结冰条件和结冰机理方面均存在差异，因此有必要将冰晶结冰与过冷水滴结冰进行区分。

过冷水滴结冰是低温云层中的过冷水滴撞击到低于0 ℃的飞机表面后凝结产生结冰，主要发生在机翼表面、发动机旋转帽罩、进气道前缘、和空速管等位置[7-9]。冰晶结冰则主要发生在发动机低压压气机全段、高压压气机前几级静子叶片及空速管等探测部件表面，发生冰晶结冰的潜在位置如图1所示。

图1 发生冰晶结冰的潜在位置

过冷水滴结冰一般发生在海拔7 km 以下的对流云层下部，飞行速度不高于600 km/h。而冰晶结冰主要发生在亚太沿海地区海拔7 km 以上的高空[1]，这部分区域经常发生强对流和风暴天气[10]，发动机推力损失事件发生区域如图2 所示。受夏季洋流季风的影响，高空冰晶区域逐渐向外扩散，当飞机经过雷暴区时，为躲避降雨区，往往选择在降雨区上部飞行，从而导致发动机吸入大量冰晶颗粒，发生冰晶结冰。

图2 发动机推力损失事件发生区域（红黑点）[10]

冰晶进入压气机内部高温环境时，会部分融化形成冰水混合物，湿润发动机叶片，后续进入的冰晶黏附在叶片上形成积冰，并降低叶片温度；随着冰晶不断吸热融化，叶片温度降低区域向发动机内部延伸，后续进入的冰晶持续黏附造成结冰区域不断扩大。发动机产生冰晶结冰，轻则导致流道截面积变小，内部流场不稳定，发动机推力降低，功率下降；重则随着积冰的脱落，打伤发动机叶片，造成叶片损伤，甚至损坏整台发动机[8，10-12]。过冷水滴结冰和冰晶结冰的主要区别见表1。

表1 过冷水滴结冰和冰晶结冰的主要区别

上述关于冰晶结冰的介绍仅停留在定性层面，并不能准确描述其内在机理，且未分析总温、总压、冰晶直径、形状和液态水含量等因素对结冰的影响，因此有必要深入的开展冰晶结冰机理研究。

2 冰晶结晶机理

过往研究表明[1-2，6]，飞机驶入含有高浓度冰晶区域时，可能会发生严重的发动机结冰故障，开展冰晶结冰试验研究和数值仿真研究，可以从源头上解决航空发动机冰晶结冰故障。

2.1 国内外试验研究进展

2.1.1 静态部件冰晶结冰试验

加拿大国家研究委员会（National Research Council，NRC）从2007 年以来就针对冰晶结冰开展了相关研究[13]，在已有的冰粒生成系统的基础上开发了全新的冰晶结冰试验台架，模拟高空中的冰晶结冰过程。该系统[14-15]主要包括冰粒生成系统、压缩机试验管道系统和出口管道，NRC 可视化试验台冰粒试验系统如图3 所示。冰粒生成系统由冰进料机，研磨机和冰粒注入设备组成，冰进料机通过调节水流调整获取的冰水含量，研磨机通过改变系统参数调整颗粒大小，喷射系统将冰粒输送到试验台，该系统可以控制冰粒的速度和方向，喷嘴处安装有液滴雾化器，雾化后的水滴直径为20~40 μm。

图3 NRC可视化试验台冰粒试验系统

压缩机试验管道系统由120 cm×37 cm×28 cm 的入射管道和过渡测试管道组成，过渡测试管道如图4所示。管道内安装探针测试来流压力和温度，过渡测试管道长约70 cm，垂直偏移量为20 cm，翼形测试件分别固定在上下两侧的铝板上，通过调整铝板温度间接模拟旁路气流温度对压气机内叶片结冰的影响。

图4 过渡测试管道[15]

针对小颗粒冰晶结冰，NASA Glenn 中心推进系统实验室（Propulsions Systems Lab，PSL）对其全尺寸高空发动机试验装置进行改造[16]，PSL 实验室结冰试验设备如图5 所示。通过喷淋装置喷出液滴水，随后构建外部冷却环境，液态水滴与周围冷空气对流传热和蒸发冷却，部分小水滴凝结成为小冰晶，大水滴仍旧保持液态，利用冰晶和液态水的混合相在试验件表面进行结冰试验，该设备可产生的冰粒直径为15 μm左右。

图5 PSL实验室结冰试验设备

2006年Mason[1]首次确认高空冰晶会进入发动机内部并产生结冰，导致发动机推力损失。2007 年，MacLeod[13]利用NRC 结冰试验设施，将加热铝板置于冰晶条件下，在短时间内产生快速的冰晶堆积，进一步证明了冰晶结冰的存在。2009 年，波音公司联合NRC实验室与NASA合作[15，17]，在前述试验的基础上，使用S 型几何管道模拟高压压气机和低压压气机之间的叶片范围，在流场中掺入冰粒撞击几何测试表面，模拟发动机内部冰晶结冰过程，证明存在最佳的结冰条件，结冰效果与液态水含量（Liquid Water Content，LWC）占总水含量（Total Water Content，TWC）的比值MR有关，结冰效果与液态水含量占比的关系如图6 所示。当LWC占比较大时，缺乏足够的冰晶冷却叶片表面造成结冰减少，当LWC占比较小时，没有足够的液态水湿润叶片黏附冰晶，同样导致结冰效果降低。

图6 结冰效果与液态水含量占比的关系

确定存在冰晶结冰现象后，探究影响冰晶结冰的因素成为科研人员关注的重点。研究表明，LWC/TWC[18]，湿球温度（Wet-bulb Temperature，Twb）[14，21]，压力[14，21]，粒径[18，26]，冰晶形状[21-23]等都是影响冰晶结冰的关键因素。

Currie等[18]用黏附效率表征冰晶在发动机内的结冰量，研究MR对冰晶结冰效果的影响，试验表明，当MR=10%~20%时，黏附效率达到最大的0.4~0.5；并证明黏附效率与冰晶入射角度相关，当冰晶倾斜入射时，黏附效率与TWC、马赫数和粒径密切相关，当冰晶垂直入射时，上述因素相互独立；Currie 基于此试验数据，忽略冰层与壁面之间的热传导、相变和溢流等过程，开发了冰晶结冰的数值模拟方法，并预测了球形表面的结冰量。

研究表明[13]，湿球温度决定结冰是否牢固。2010年，NASA 与NRC 在发动机高空试验研究设备（Research Altitude Test Facility，RATFac）上开展试验测试[13]，在总压为45～93 kPa 下进行了一系列翼型结冰试验，试验发现，在压力为45 kPa，即Twb<0 ℃时，产生明显的明冰；在最高压力为93 kPa，即Twb>0 ℃，形成易脱落的疏松冰。分析表明，气压降低，会增加测试件表面液态水的蒸发冷却，导致湿球温度降低，结冰更加坚固，因此低压情况下更易形成稳定的明冰。但冰晶含量较小时，即使Twb<0 ℃，由于没有足量的冰晶，黏附效率也较差，而当Twb>0 ℃时，即使含有足量的冰晶，由于湿球温度较高，产生的结冰也较为稀松[25]。湿球温度对冰晶融化过程也存在影响[26]，湿球温度增加会加速冰晶的融化过程。

2016 年，研究人员[16，19]在PSL 试验设备上开展全液滴和冰晶的结冰试验，并发现结冰时湿球温度几乎不变或仅有少量增长。

Bartkus 等[20]根据能量平衡和湿球温度大小提出了2 种结冰原理：冻结式结冰和融化式结冰。冻结式结冰发生在Twb<0 ℃时，结冰主要由冰水混合相中的过冷水滴冻结产生，冰层较为坚固；融化式结冰发生在Twb>0 ℃时，结冰主要由冰水混合相中未融化的冰晶颗粒堆积产生，冰层较为稀松。

2018 年，NRC 联合NASA 在RATFac 试验设备上利用NACA 0012 翼型进一步开展混合相态结冰机理试验[21]，探究湿球温度、马赫数、总水含量和粒径对结冰的影响。当Twb>3 ℃时，结冰速率最大，但由于冰层中含有较多的液态水，结冰脱落明显，随湿球温度降低，结冰脱落量减小，冰形更加稳固，Twb降至0 ℃时结冰量最大；试验发现，马赫数增加导致结冰速率减小，这可能是侵蚀加剧的原因。

在湿球温度相同时，增加马赫数会降低冰晶融化率，见表2、3，这是因为马赫数增大，冰晶融化时间减小导致；试验验证了Tsao等[27]提出的最小结冰阈值假说，探究了TWC对结冰的影响，当TWC<3 g/m3时，进入管道的冰晶都被融化，几乎不会发生结冰，随TWC增加，结冰量近似线性增加，当TWC=10 g/m3时，结冰量达到最大。

表2 Ma=0.25时MR（T0=15 ℃，P0=5 Pa，small PSD）

表3 Ma=0.40时MR（T0=15 ℃，P0=5 Pa，small PSD）

粒径对冰晶结冰效果有重要影响。Knezevic等[26]利用压气机内管道排放槽开展冰晶结冰模拟，发现在湿球温度、冰晶流量相同时，冰晶粒径越小，结冰效果越明显；Currie 等[18]研究了大小不同的冰晶在半球形表面的结冰效果，证明较小直径的冰晶可以产生更多的结冰，不同尺寸结冰效果如图7 所示；Struk等[21]分别在马赫数为0.25 和0.40 时进行不同尺寸的冰晶结冰试验，结果表明小尺寸冰晶可以产生较好的尖锐箭头形结冰。主要有2 种因素解释该现象，一方面，在质量流量相同时，颗粒直径越大，对积冰的撞击和侵蚀效果越明显，导致积冰量越少；另一方面，颗粒直径越小，越容易被结冰表面的水膜黏附，减弱了了冰晶粒子的侵蚀作用，而冰晶直径增大时，水膜黏附作用减弱，侵蚀加剧，结冰量减少。

图7 不同尺寸结冰效果[18]

上述试验证明冰晶粒子会对结冰表面产生侵蚀作用，但冰晶侵蚀常与多种过程耦合，如冰晶粒子撞击后的反弹、飞溅等，这使研究变得复杂。为单独研究冰晶对结冰表面的侵蚀过程，NASA 开展了一系列概念验证测试[28-30]，将冰晶侵蚀过程与其他物理过程解耦[20]，首先利用含有过冷水滴或冰水混合相的气流在翼形表面产生结冰，其次设置特定的喷嘴参数和流动参数，将表面附有冰层的翼形暴露在含有冰晶的流场中，单独研究冰晶对结冰表面的侵蚀作用。为定量评估侵蚀过程，定义侵蚀系数E为被侵蚀冰的质量通量与撞击的冰晶的质量通量的比值

式中：ṫ为冰层厚度增加速率；ρacc为冰层密度；β为壁面冰晶收集系数；U∞为来流速度；TWC为总水含量。

试验研究了空气流速、粒径直径、TWC和总温对侵蚀的影响，结果表明：空气流速越大，侵蚀系数越大；增大粒径直径，会加剧冰晶侵蚀效果；同时总温降低可能会造成喷头冻结，使得撞击到测试件表面的TWC变小，导致侵蚀系数偏大。

2.1.2 动态部件冰晶结冰试验

前述内容均是静态的测试件表面冰晶结冰试验，无法反映发动机运转过程中关键部位的影响，如冰晶撞击转子叶片后导致的破碎、冰晶沿程的快速热力学变化、冰粒经过压气机转子叶片的离心、弹跳等物理过程。为此，NRC 在美国联邦航空管理局（Federal Aviation Administration，FAA）的资助下研发了冰晶环境轴向压气机（Ice-Crystal Environment Modular Axial Compressor Rig，ICE-MACR），用于模拟更真实的发动机高空冰晶结冰环境[33-35]，ICE-MACR 设备如图8所示。该系统不仅配备2 个轴向压缩机来模拟涡扇的风扇级（stage 1）和助推级（stage 2），还配备了更先进的仪器测量系统测量和观察发动机内部结冰情况。首先，冰粒在入口处随气流吹入，在压缩级与旋转叶片发生碰撞，离心，在扩展段（extender duct）吸热融化，产生冰水混合相，其次，仪器测试段（instrument segment）对环境参数进行详细测量，最后利用试验段（test article）中周向排列的15个NACA 0012翼型叶片进行结冰测试。该压气机结冰系统轴向长度与全尺寸发动机相等，径向方向通过全尺寸发动机缩放满足风洞测量需求。

图8 ICE-MACR设备[34]

2021 年，Chow 等[32]利 用ICE-MACR 在IWT 风 洞中开展冰晶结冰试验，在试件表面加装热通量计，研究冰晶与结冰表面的热量传递和温度变化过程；文献[34]研究了单段压气机的冰晶沿程变化，得到了冰粒经过单段压缩级后的粒径尺寸分布；文献[35]研究了湿球温度对液态水含量的影响，通过改变湿球温度观测粒子的沿程融化率，线性拟合得到湿球温度与颗粒融化率的关系。

相比单段压缩级，2 段压缩级不仅对冰晶产生额外的加压和加热，还会产生更大的轴向速度和离心速度。文献[36]分别研究了冰晶在单段和2段压气机转子叶片上的碰撞与破碎过程，探究了叶片旋转产生的离心力对颗粒分布和流场的影响，证明随转子转速增加，颗粒尺寸渐进线逐渐向最小尺寸靠拢，得到冰粒破碎后的尺寸分布，并提供了压缩级后径向方向的冰粒尺寸和TWC大小分布；文献[37]分别研究了采用单段和双段压缩级后结冰的效果，量化结冰严重程度，并对比了结冰效果，单段和双段压缩级结冰如图9所示。

图9 单段和双段压缩级结冰[37]

试验表明，冰晶融化比（MR）是影响结冰效果的关键，存在1 段MR范围可以产生最佳的冰晶结冰。2022 年，Mason 等[38]利用ICE-MACR 研究TWC对结冰的影响，探讨MR与TWC之间的关系，通过观察总水含量在TWC阈值附近变化时积冰的生长和脱落过程，模拟真实高空环境下的冰晶结冰过程。结果表明，产生最佳冰晶结冰的最小MR值不随TWC变化，但该过程中结冰速率与TWC密切相关。

需要说明的是，目前针对冰晶结冰所开展的试验研究并不充分，如上游吸积过程对下游结冰的影响，决定结冰的各因素之间的相互影响研究[31]，触发结冰的关键因素研究[37]，冰晶结冰相似准则研究，需要进一步发展测量技术，如喷雾冷凝装置的防冻技术[31]、液态水含量和总水含量的测试方法[39]、结冰生长速率测量技术[40]等，并更新风洞试验设备，进一步深入开展冰晶结冰机理研究。

2.2 数值模拟

飞机结冰研究中，数值模拟占据着重要地位。相比试验手段，数值模拟计算成本低，节约人力和财力，同时规避了试验可能遇到的风险。冰晶结冰本质上是传质传热的过程，涉及冰晶粒子融化相变、碰撞、弹开、飞溅、黏附和冻结等多种物理现象。根据冰晶吸入发动机的过程，可将冰晶结冰数值模拟分为3 个阶段，分别是空气流场计算、冰晶运动轨迹和撞击特性计算、结冰计算。

空气流场计算主要分为面元法和求解Euler方程或者Navier-Stokes 方程的方法，第1 种方法主要用于求解简单物体周围的流场分布，随着计算机技术的发展，目前主要采用CFD 方法求解空气流场。以往的过冷水滴结冰计算认为空气流场是水滴-空气两相流，水滴尺寸小、含量低，因此忽略过冷水滴对空气的作用，但冰晶颗粒直径在200 μm左右，含量甚至可以达到9 g/m3，因此是否可以忽略冰晶对空气流场的作用是最先需要解决的问题。Rios 等[41-42]对冰晶运动过程中的动力学特性进行研究，将颗粒流动分为稀疏离散相流动和稠密离散相流动，认为冰晶会影响气流分布，因此应当采用稀疏离散相双向偶和的方法。而Nilamdeen 等[43]则认为尽管冰晶粒径大，含量高，但冰粒对空气流场的作用依然较小，并针对空气、水滴和冰晶的三相流动开展了单向耦合计算研究。综上所述，冰晶与流场是否需要双向耦合计算仍需进一步研究，但为简化起见，目前主流研究仍旧采用单向耦合的方法。

冰晶运动轨迹和撞击特性计算主要有2 种方法，分别是Lagrange 法和Euler 法。拉格朗日法将粒子视为离散相，建立每个粒子的控制方程，并求解其运动轨迹，欧拉法则将粒子看作连续相，引入粒子容积分数的概念，求解粒子的连续方程和动量方程。运动轨迹和撞击特性计算需要考虑多种因素，如冰粒运动过程中的融化相变和传热，冰粒的大小、形状，冰粒撞击到壁面后的弹开、破碎、飞溅以及冰粒对结冰表面的侵蚀等过程。受限于试验数据的缺乏，目前开展的冰晶运动轨迹计算主要针对某一些特定因素建立相关模型，对其他因素则进行简化。

Villdedieu 等[22]用颗粒球形度定义非球形冰晶粒子，并模拟冰粒运动和相变的过程，根据冰粒温度（Ts）与融化温度（Tm）的相对大小判断冰晶的融化程度，将粒子沿轨迹的相变分为3个阶段：

（1）当Ts

（2）当Ts=Tm时，颗粒逐渐融化，并在外表面形成水膜，整个冰核和水膜温度不变，颗粒吸收的热量与冰晶融化潜热和水膜蒸发潜热平衡。

（3）当Ts>Tm时，颗粒为球形液滴，液滴从外界吸收热量且表面蒸发，直至达到平衡。

对流与蒸发过程引起冰晶运动过程中的质量和能量传递，Villdedieu 定义了粒子表面的对流和水膜蒸发模型，其中对流是颗粒与其周围空气之间温差驱动的结果，与努塞尔数相关，液膜蒸发是液膜表面和周围气流的蒸汽浓度差驱动的结果，与舍伍德数相关，其表达式为

式中：为对流换热量为蒸发质量；Nu为努塞尔数；Sh为舍伍德数。

Trontin 等[47]利用冰晶融化比、液态水含量和总水含量表征冰晶撞击壁面后的黏附和弹开过程，定义黏附系数为

式中：εs,c为纯冰晶气象时的黏附系数；εs,d为混合相气象时的黏附系数。

式中：ηm为冰晶融化比；φd为撞击到壁面的过冷水质量占总水质量的比值；φic为撞击壁面的冰晶质量占总水质量的比值；φd与φic的和为1；Kd与Kc均为常数。

该式表示，在纯冰晶条件下，若冰晶粒子全部融化，则ηm为1，黏附系数达到最大值1，若冰晶融化率为0，则黏附系数为最小值0。

该模型通过考虑液态水含量占总水含量的比值间接计算液膜厚度对黏附系数的影响，计算较为方便，但并未考虑粒子速度、尺寸和形状的影响，因此在模拟中的实用性有限，仍旧具有较大的改进空间。

Norde 等[23，44-45]基于欧拉法对Villdedieu 定义的冰粒形状、融化和相变模型进行了计算验证，建立MooseMBIce方法对3维机翼进行了积冰模拟，采用单向耦合计算，阻力系数采用文献[46]中定义的Ganser模型，即阻力系数CD取决于颗粒球度和颗粒横向球度，并简化了Trontin 的碰撞模型[22，47]，认为冰粒只会在表面反弹或黏附，用冰晶法相动能和表面能量之比模拟冰晶的弹开和破碎过程，定义冲击数L以区分黏附和弹开，计算冰晶撞击反弹对结冰的影响。

式中：eσ为表面能量；un为冰晶撞击速度。

2019 年，Norde[48]将1 个发动机的3 级几何结构（定子-转子-定子）进行建模，分别采用拉格朗日法和欧拉法模拟冰晶在压气机内部的运动和融化过程，结果表明，当假定冰晶完全沉积在叶片表面时，2 种方法计算的撞击质量流量基本一致，仅在融化率方面略微不同，并研究了颗粒尺寸对冰晶融化过程的影响；但当考虑冰晶的反弹和破碎等因素时，欧拉法在处理凹面时由于轨迹发散导致结果不同，拉格朗日法计算时结果难以收敛。

最初的结冰计算均是基于Messinger 结冰热力学模型[7，49-50]，通过建立结冰表面的质量和热量平衡方程计算结冰量，该模型适用于计算过冷水滴结冰，但未考虑冰层及水膜内的温度分布梯度和传热等非稳态特性，在后续研究中，学者基于试验结果提出了基于Messigner 热力学模型的改进模型，在模型中增加了有关冰晶的质量和能量传递项方程[48]。

控制体内质量守恒方程如图10所示，即

图10 控制体内质量守恒方程

式中：为流场中冰晶融化进入控制体的部分；ṁc,d为流场中液态水滴进入控制体的部分为上一控制体溢流到该控制体的部分为控制体中蒸发带走的水蒸气为控制体中液态水结成冰的部分；为本控制体溢流到下一控制体的部分。

控制体内结冰的质量守恒方程为

式中：为流场中未融化冰晶撞击到水膜上黏附的部分；为冰层中升华的部分为总结冰量。

控制体内能量守恒方程如图11所示，即

图11 控制体内能量守恒方程

式中：Qke,d为流场中液滴碰撞带来的动能；Qke,ic,w为流场中冰晶融化的液滴碰撞带来的动能；Qke,ic,i为流场中未融化的冰晶碰撞带来的动能；Qin为上一个控制体溢流水带来的显热；Qf为控制体内液态水结冰的潜热；Qi为冰层与水膜的热量传递显热；Qconv为控制体对流换热带走的能量；Qev为控制体内液滴蒸发的潜热；Qc,d为控制体内液态水的显热；Qc,ic,w为控制体内冰晶融化的液态水的显热；Qc,ic,i为控制体内未融化冰晶的显热；Qout为流向下一控制体的显热。

对2 维结冰壁面而言，驻点所在控制体流入的溢流水质量为0，且此控制体内的液态水分别向上表面和下表面溢流。对其他控制体，控制体内溢流进的水的质量等于前一个控制体溢流出的水的质量[59]。

假定壁面温度等于融化温度，即Tm=Ts，根据能量守恒可以得到结成冰的质量ṁf，进而判断结冰状态。

（2）当+<0 时，进入控制体内的液态水都未结冰，表面为湿润状态，此时= 0,Qi= 0，质量方程为

随研究不断深入，研究人员基于试验结果开发了结冰数值模拟软件，如美国的LEWICE、GlennICE，加拿大的FENSAP-ICE，并基于结冰软件开展冰晶结冰研究。Bourgault 等[51]利用FENSAP-ICE 进行基于欧拉法的结冰计算；Nilamdeen等[43，52]考虑了冰晶反弹的动力学效应，利用FENSAP-ICE 软件进行空气、水滴和冰晶的三相流结冰模拟，并与Cox 风洞试验结果进行对比，该模型成功预测了实际工况下的结冰趋势和方向，但由于未考虑冰层密度变化和冰晶的飞溅过程，预测结冰量大于实际结冰量；文献[53]基于大量试验的结果，改进了过冷大水滴运动撞击模型，对比了LEWICE 3D 和GleenICE 结冰软件的计算结果；Wright等[54]利用GlennICE 软件开展混合相结冰模拟，考虑粒子球形度的影响，利用冰晶与空气的换热和温度变化模拟冰晶的相变过程，利用回弹系数表征冰晶黏附与反弹量的大小，利用侵蚀经验公式模拟侵蚀的效果；Bidwell 等[55]使用LEWICE 3D 的过冷大水滴模型模拟了粒子的破碎和飞溅过程，并与试验结果进行对比；2021 年，Nilamdeen 等[56]利用ANSYS 流体仿真模块对ICE-MACR 冰风洞结冰设施进行建模，在FENSAP-ICE 软件中开展空气、水和冰晶三相流结冰模拟，并与试验结果进行对比，证明该方法可以准确地预测流场流动、冰晶运动和结冰过程，为进一步更新结冰软件提供指导。

中国冰晶结冰研究相比国外起步较晚，由于缺乏冰晶结冰风洞和试验数据的不足，研究仍处在初步阶段。张丽芬等[12]建立了混合相结冰模型，考虑了冰晶粒子的破碎、反弹作用及表面水膜的飞溅过程，并将仿真结果与风洞试验进行比较。姜飞飞等[57]研究了冰晶在压气机内的传质和传热情况，得到了冰晶在不同条件下的沿程变化和运动轨迹；黄平等[58]基于欧拉法建立了冰晶粒子的运动相变模型，模拟了环境温度、粒径大小、初始球形度和气流相对湿度对冰晶粒子融化的影响；卜雪琴等[59]对Messinger 结冰热力学模型进行了改进，研究了机翼表面的传质传热过程，模拟了2 维机翼表面的混合相结冰过程；马乙楗等[60]基于拉格朗日法建立了冰晶运动-传热传质耦合的数值计算方法，更新了冰晶撞击模型和黏附系数计算方法，研究了冰晶形状对融化过程和运动相变特性的影响，研究了来流总温与液态水含量对黏附特性的影响，得到了冰晶运动相变与黏附特性的影响规律；郭向东等[61]发展了基于欧拉法的两相流传质传热耦合计算方法，研究了典型大型结冰风洞内的冰晶运动和相变过程。谭燕[62]基于欧拉法对某楔形翼型进行冰晶结冰数值模拟，分析了压力与液态水含量对结冰情况的影响。郭琪磊等[62]利用FENSAP-ICE 结冰软件，研究了混合相态条件下的结冰过程，分析了环境温度、马赫数对收集系数、冰型的影响，明确了融化率对结冰的作用机制。

3 研究展望

综上所述，冰晶结冰研究刚刚起步，由于试验数据的缺乏，冰晶结冰的数值模拟仍在初级阶段，需进一步开展冰晶结冰机理试验，研究冰晶进入发动机后结冰过程，并开发更加准确的数值模拟方法。针对冰晶/混合相结冰，还有以下方面需要重点关注：

（1）冰晶碰撞模型。目前数值建模中几乎都对冰晶与表面撞击过程进行简化处理，认为冰晶在撞击表面弹开或者黏附，并未详细考虑冰晶碰撞到叶片后的破碎、飞溅等现象，因此对冰晶撞击过程的详细描述，是今后研究的关键方向。

（2）完善叶片表面的结冰热力学模型。目前叶片表面结冰热力学模型均是基于Messinger 模型的改进模型，并未细致考虑结冰过程中冰层与结构的热量传递和冰层及液膜内的温度梯度，忽略了冰晶撞击、破碎、弹开过程中的质量和能量变化，因此并不能反映真实的结冰过程，需进一步完善。

（3）结冰探测技术。目前结冰风洞测量手段单一，且存在精度不准和稳定性不足的情况，需进一步开发更精准的结冰环境测量技术和结冰探测技术。

（4）进一步开展风洞试验研究。数值模型需要大量的试验数据进行模型修正，受限于风洞试验能力不足，很多建立的模型难以通过试验数据进行模型修正和验证，因此需要进一步开展冰晶结冰机理试验研究。

4 结束语

与传统过冷水滴结冰不同，冰晶结冰主要对航空发动机产生危害，深入研究其结冰机理，发展压气机防除冰技术，对飞行安全具有重要意义。本文详细阐述了冰晶结冰的危害、与传统过冷滴结冰的区别，国内外试验研究和数值仿真建模所取得的重要进展。中国冰晶结冰研究刚刚起步，仅建立了初步的冰风洞试验室[64]，但大涵道比民用发动机的科研试飞和试航取证均需在充分认识冰晶结冰机理的基础上开展，深入开展航空发动机的冰晶结冰研究迫在眉睫。