航空發動機冰晶結冰研究進展

蘇龍偉，申世才，田曉平，吳 悠

（中國飛行試驗研究院，西安 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為被侵蝕冰的質量通量與撞擊的冰晶的質量通量的比值

式中：t?為冰層厚度增加速率；ρ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 控制體內質量守恒方程

式中：為流場中冰晶融化進入控制體的部分；m?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，根據能量守恒可以得到結成冰的質量m?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]，但大涵道比民用發動機的科研試飛和試航取證均需在充分認識冰晶結冰機理的基礎上開展，深入開展航空發動機的冰晶結冰研究迫在眉睫。