摘要：為探究冰晶與過冷水滴共存時的混合相結冰特性，基于FENSAP-ICE結冰計算平臺，以NACA0012翼型為研究對象，采用阻力系數(shù)對非球形冰晶的運動過程進行修正，考慮了冰晶黏附效應和侵蝕效應的影響，建立了歐拉框架下冰晶運動、碰撞、黏附、結冰和侵蝕全物理過程的混合相結冰數(shù)值計算方法。分析了冰晶質(zhì)量濃度（ice water content，IWC）、液態(tài)水質(zhì)量濃度（liquid water content，LWC）以及冰晶黏附和侵蝕效應對結冰形態(tài)和結冰厚度的影響規(guī)律，對比了NTI和NRC 兩種冰晶黏附模型的適用條件。研究結果表明：當總水含量（total water content，TWC）一定時，隨著冰水含量比IWC與LWC之比的增加，結冰覆蓋面積逐漸減小，最大結冰厚度先增大后減小，冰形由冠狀逐漸變?yōu)榧饨菭?，當IWC與LWC分別為0.4g·m－3和1.0g·m－3時，駐點處的結冰厚度達到最大值；對于高IWC工況，NRC黏附模型的駐點結冰厚度更接近實驗結果，誤差為7.2%；對于低IWC工況，NTI模型與NRC模型的駐點結冰厚度接近；侵蝕作用主要影響翼型駐點附近的區(qū)域，IWC或LWC的增大均會加劇侵蝕效應。研究可為飛行器高效防除冰系統(tǒng)的設計提供理論依據(jù)。

關鍵詞：混合相；冰晶結冰；黏附；侵蝕；數(shù)值模擬

中圖分類號：V211.41"文獻標志碼：A

DOI：10.7652/xjtuxb202410015"文章編號：0253-987X（2024）10-0168-10

Numerical Study on Mixed-Phase Ice Crystal Icing on Airfoils

ZHONG Fuhao1， LIU Senyun2， LIU Xiufang1，3， DAI Xinbo2，

CHEN Jiajun1， MIAO Qingshuo1， YI Xian2， HOU Yu1，3

（1. School of Energy and Power Engineering， Xi’an Jiaotong University， Xi’an 710049， China;

2. Key Laboratory of Icing and Anti/De-icing of China Aerodynamics Research and Development Center

（CARDC）， Mianyang， Sichuan 621000， China; 3. MOE Key Laboratory of Cryogenic Technology and

Equipment， Xi’an Jiaotong University， Xi’an 710049， China）

Abstract：To explore the characteristics of mixed-phase icing when ice crystals and supercooled water droplets coexist， a numerical study is conducted utilizing the FENSAP-ICE simulation platform. The NACA0012 airfoil is selected as the subject of the study. The motion of non-spherical ice crystals is adjusted by modifying the drag coefficient， and both ice crystal adhesion and erosion effects are taken into consideration. Within the Euler framework， a comprehensive numerical calculation method is developed to simulate the complete physical process of ice crystal motion， collision， adhesion， accretion， and erosion. The influence of ice water content （IWC）， liquid water content （LWC）， as well as ice crystal adhesion and erosion effects on the shape and thickness of icing is analyzed. The applicability conditions of two ice crystal adhesion models， namely NTI and NRC， are compared. The research findings suggest that， when maintaining a constant total water content （TWC）， an increase in the ratio of IWC to LWC leads to a gradual reduction in the ice coverage area. Additionally， the maximum ice thickness initially increases and then decreases， accompanied by a transition in ice shape from crown to angular. The maximum ice thickness at the stagnation point is achieved when IWC and LWC are 0.4g·m－3 and 1.0g·m－3， respectively. Under high IWC conditions， the NRC adhesion model demonstrates a stagnation point ice thickness that closely aligns with experimental results， with an error of 7.2%. Conversely， under low IWC conditions， both the NTI and NRC models yield similar stagnation point ice thicknesses. The erosion effect primarily affects the region near the airfoil’s stagnation point， and an increase in IWC or LWC intensifies this erosion effect. This study provides a theoretical foundation for the design of an efficient anti-icing system for aircrafts.

Keywords：mixed phase; ice crystal icing; adhesion; erosion; numerical simulation

當飛機穿越含有過冷水滴或冰晶的云層時，可能面臨嚴重的結冰現(xiàn)象。早期的研究認為飛行器結冰起因于過冷水滴結冰，過冷水滴結冰是低空云層中的過冷水滴撞擊至飛行器迎風表面凝結產(chǎn)生的結冰現(xiàn)象［1］。然而，自20世紀90年代以來，海拔7km以上，在過冷水滴難以存在的高空，依然發(fā)生了多起不明原因的發(fā)動機功率損失事件，這引起了航空科研人員的廣泛關注。對多起發(fā)動機推力損失事件的調(diào)查結果表明：發(fā)動機吸入冰晶所產(chǎn)生的黏附結冰是導致該類故障的主要原因［2］。后續(xù)深入研究發(fā)現(xiàn)，純冰晶以固相的形式撞擊冷表面會發(fā)生反彈，不會發(fā)生黏附結冰現(xiàn)象，液態(tài)水的存在是冰晶黏附結冰的必要條件［3］。由于液態(tài)水潤濕了構件表面，來流中的冰晶在撞擊表面過程中才可能發(fā)生黏附結冰現(xiàn)象。冰水混合相的形成可分為以下兩種情況：一是對流云層中存在既有冰晶又有水滴的混合相冰晶結冰氣象條件，會引起機翼、尾翼等冷表面結冰；二是當冰晶被吸入發(fā)動機內(nèi)流道后，在熱氣流作用下部分融化，形成冰水混合相，可能導致發(fā)動機核心部件結冰。機翼結冰會導致飛機氣動性能惡化，發(fā)動機核心部件結冰則會引發(fā)喘振、失速、燃燒室熄火等嚴重后果，二者均會對飛行安全構成重大威脅［4-5］。

混合相冰晶結冰的實驗研究需要在結冰風洞中進行，由于結冰氣象條件的地面模擬和冰晶制備難度較大，制約了混合相冰晶結冰地面模擬實驗平臺的建設［6］。數(shù)值模擬是目前開展混合相冰晶結冰研究的重要手段［7-10］。根據(jù)混合相冰晶結冰發(fā)生的物理過程，其數(shù)值研究可分為冰晶運動過程數(shù)值模擬、冰晶黏附結冰特性數(shù)值模擬和冰晶侵蝕效應數(shù)值模擬3個方面。

冰晶運動計算是黏附結冰特性數(shù)值模擬的基礎。冰晶在氣流中的運動屬于離散兩相流動的研究范疇，通常采用拉格朗日方法或歐拉方法進行模擬。拉格朗日方法通過追蹤單個冰晶粒子的運動，能夠精確獲取冰晶在任意時刻的空間位置和運動狀態(tài)。Villedieu等［10］采用拉格朗日方法對翼型外部流動中的冰晶運動過程進行了數(shù)值模擬。與之相對，歐拉法將冰晶視為連續(xù)介質(zhì)，通過求解相應的輸運方程來獲得冰晶的體積分數(shù)和動量分布，此方法在處理冰晶與復雜幾何體相互作用時具有較低的計算成本。Norde等［11］和Nilamdeen等［12］利用歐拉法建立了冰晶在氣流中的運動模型，分別針對翼型外部流動和發(fā)動機內(nèi)部流動進行了數(shù)值模擬。

冰晶黏附結冰與液態(tài)水含量、冰晶粒徑、撞擊速度等多個參數(shù)有關，黏附機理異常復雜［6］。目前，研究中主要采用兩種類型的冰晶黏附模型。Trontin等［13］基于冰晶粒子的平均融化率、水滴和冰晶的質(zhì)量分數(shù)建立了冰晶黏附半經(jīng)驗模型，并通過加拿大National Research Council（NRC）的實驗數(shù)據(jù)對模型中的參數(shù)進行了擬合，從而發(fā)展了NRC冰晶黏附模型。Nilamdeen等［14］基于歐拉法建立了Newmerical Technologies International（NTI）冰晶黏附模型，該模型與冰晶粒徑、撞擊速度以及液膜厚度等參數(shù)密切相關。

侵蝕效應是冰晶結冰區(qū)別于過冷水滴結冰的獨有現(xiàn)象，即來流中的冰晶對已形成的冰層具有侵蝕作用。Struk等［15-16］以及Currie等［17］在研究冰晶結冰現(xiàn)象時觀察到了侵蝕效應，并發(fā)現(xiàn)侵蝕程度與冰晶的融化率和粒徑密切相關。在此基礎上，Trontin等［13］基于實驗數(shù)據(jù)提出了冰晶侵蝕半經(jīng)驗模型，該模型根據(jù)冰晶含量、粒徑等參數(shù)來預測冰晶侵蝕效果。

近年來，隨著我國航空事業(yè)的飛速發(fā)展，混合相冰晶結冰現(xiàn)象受到越來越多的國內(nèi)學者關注。西北工業(yè)大學的張麗芬等［18］、中國空氣動力研究與發(fā)展中心的馬乙楗等［19］基于拉格朗日法確定了冰晶和液滴的運動軌跡，建立了機翼混合相結冰模型。上海交通大學的姜飛飛等［20］、西安交通大學的陳佳軍等［21-22］基于拉格朗日方法描述了冰晶在二維壓氣機內(nèi)的運動、碰撞及換熱過程，并分析了冰晶液態(tài)水含量、溫度、速度等參數(shù)的沿程變化規(guī)律。北京航空航天大學的常士楠團隊［23］、卜雪琴等［24］采用歐拉法分別針對航發(fā)進氣系統(tǒng)和翼型建立了冰晶混合相結冰計算模型，探究了氣流條件和冰晶初始參數(shù)對冰晶運動換熱及黏附結冰過程的影響。西北工業(yè)大學的郭琪磊等［25］利用FENSAP-ICE結冰計算軟件研究了馬赫數(shù)和氣流溫度等參數(shù)對混合相冰晶結冰特性的影響規(guī)律。

綜上可知，在冰晶與過冷水滴同時存在的混合相氣象條件下，結冰過程屬于空氣-冰晶-水滴三相混合流動的熱動力學過程，其機理異常復雜。液態(tài)水含量以及黏附和侵蝕效應對結冰特性的影響機理尚不明確，亟需對混合相結冰現(xiàn)象開展更加深入和全面的研究。本文針對過混合相結冰過程，基于FENSAP-ICE結冰計算平臺，考慮了非球形冰晶粒子的運動特性，通過改進的Messinger模型實現(xiàn)了冰晶黏附結冰的計算，并考慮了冰晶侵蝕效應對冰形的影響，建立了歐拉框架下涉及冰晶運動、黏附、結冰和冰晶侵蝕的全物理過程混合相冰晶結冰計算方法。以NACA0012翼型為研究對象，通過與風洞實驗結果比較，驗證了該計算方法的準確性，并在此基礎上分析了混合相中冰水含量比對結冰過程的影響，比較了兩種黏附模型及其適用條件，并研究了侵蝕效應對冰層生長的影響。

1"混合相冰晶結冰數(shù)值方法

1.1"冰晶運動計算方法

本文采用歐拉雙流體模型計算水滴和冰晶的運動換熱和撞擊特性，遵循連續(xù)性方程、動量方程和能量方程［9，26］。

球形與非球形粒子在氣流中運動過程的求解主要區(qū)別在阻力系數(shù)的計算。阻力系數(shù)和雷諾數(shù)分別用Cj和Rej表示，j=d表示液滴，j=D表示冰晶?？紤]到液滴的粒徑較小，因此假設液滴在表面張力的作用下保持為球體。液滴阻力系數(shù)Cd的計算公式［12］如下

Cd=24Red（1+0.15Re0.687d），Red≤1300

0.4，Redgt;1300（1）

液滴或冰晶的雷諾數(shù)Rej表達式為

Rej=ρaua－ujdμa（2）

式中：ρa表示空氣密度，kg·m－3；ua表示氣流速度矢量，m·s－1；uj表示液滴或冰晶速度矢量，m·s－1，j=d 表示液滴，j=D表示冰晶；μa表示空氣動力黏度，N·s·m－2；d表示液滴或冰晶直徑，m，對于非球形冰晶粒子表示當量直徑，即與冰晶粒子具有相同體積的球體直徑。

冰晶為非球形，在冰晶結冰計算過程中，通過引入冰晶長徑比E（顆粒內(nèi)部的最長徑和與它相垂直的最長徑之比）來計算非球形冰晶阻力系數(shù)CD［12］。當冰晶的雷諾數(shù)ReD≤0.01時

CD=8m3ReD1+mReD48+m21440ReDlnReD2（3）

m=12（1－E2）1.5E（1－E2）0.5+

（1－2E2）tan－1［（1－E2）0.5/E］－1（4）

當0.01lt;ReD≤1.5時

CD=AReD（1+10X）（5）

X=－0.883+0.906ω－0.025ω2（6）

ω=lgReD（7）

當1.5lt;ReD≤100時

CD=AReD（1+0.138Re0.792D）（8）

當100lt;ReD≤300時

CD=AReD（1+0.00781Re1.393D）（9）

當ReDgt;300時

CD=1.17（10）

冰晶阻力系數(shù)式（5）、（8）、（9）中的A定義為

A=32fa［φ－（φ2－1）tan－1（1/φ）］（11）

f=（b2－a2）0.5（12）

φ=bdeq（13）

deq=2aE13（14）

式中：a為冰晶半短軸長度，m；b為冰晶半長軸長度，m；deq為冰晶體積等效直徑，m。

1.2"冰晶碰撞黏附模型

冰晶粒子撞擊到翼型表面時，會發(fā)生破碎、黏附、反彈、二次撞擊等行為。為實現(xiàn)冰晶結冰的模擬，需要使用合適的模型判斷冰晶撞擊黏附過程。目前主要有兩種冰晶碰撞黏附模型，NRC冰晶碰撞黏附模型和NTI冰晶黏附模型。

NRC冰晶碰撞黏附模型［10，15］為

ψ=πρdd3u2n/12πeσd2=ρddu2n12eσ（15）

式中：ψ是無量綱參數(shù)，表示冰晶法向動能與表面能的比值；ρd表示冰晶密度，kg·m－3；un表示撞擊速度的法向分量，m·s－1；eσ表示冰晶表面能，J·m－2。表面能與溫度有關，表達式為

eσ=eσ0expQsRT－QsRT0（16）

式中：eσ0表示溫度為253K時的表面能，J·m－2；Qs表示冰晶裂紋產(chǎn)生相關的活化能，J·mol－1·K－1；R為氣體常數(shù)，J·mol－1·K－1；T表示冰晶溫度，K。

將無量綱參數(shù)用兩個臨界值將冰晶碰撞行為分為完全彈性碰撞、非彈性碰撞和碰撞破碎3種情況，并定義恢復系數(shù)ξ

ξ=1，ψ≤ψc1

ξlt;1，ψc1lt;ψlt;ψc2

ξ1，ψ≥ψc2（17）

ψc1=0.5 ； "ψc2=90（18）

黏附率模型為

ε=max［εi，εd］（19）

εi=（Ki－2）η3i+（3－2Ki）η2i+Kiηi（20）

εd=Kd（φd+ηiφi）（21）

式中：ε表示冰晶黏附率；ηi表示混合相中冰晶部分的融化率；φd表示撞擊時液滴的質(zhì)量分數(shù)；φi表示撞擊時冰晶的質(zhì)量分數(shù)。ηi、φi、φd、Ki、Kd表達式為

ηi=βWβIWC（22）

φi=m·imp，im·imp，i+m·imp，d（23）

φd=m·imp，dm·imp，i+m·imp，d（24）

Ki=2.5（25）

Kd=0.6，T≥Tm

0，Tlt;Tm（26）

式中：βIWC表示氣流中冰晶質(zhì)量濃度（ice water content，IWC），kg·m－3；βW表示氣流中冰晶融化部分的液態(tài)水質(zhì)量濃度，kg·m－3；m·imp，i表示撞擊冰晶的質(zhì)量流率，kg·s－1；m·imp，d表示撞擊液滴的質(zhì)量流率，kg·s－1；Tm表示冰晶開始結冰的壁面溫度；Ki和Kd通過風洞結冰實驗擬合得到。

二次碰撞特性為

ξnn=1，ψ≤ψc1

ψc1ψ13，ψgt;ψc1（27）

ξnt=0.41－ψc2ψ（28）

dmax=ψ211c2ψ（29）

式中：ξnn表示碰撞粒子的法向動量轉(zhuǎn)化為二次粒子法向動量的分數(shù)；ξnt表示碰撞粒子的法向動量轉(zhuǎn)化為二次粒子切向動量的分數(shù)；dmax表示二次碰撞粒子碎片最大直徑，m。

NTI冰晶黏附模型［13］認為：冰晶撞擊于液膜超過一定厚度的區(qū)域會全部黏附，撞擊于霜冰區(qū)域會完全反彈；在明冰區(qū)域內(nèi)，當冰晶撞擊法向速度分量大于臨界速度時，冰晶全部反彈，當小于臨界速度時，冰晶部分黏附。

明冰區(qū)域的冰晶黏附率模型為

ε=hfmax（hf）1exp（γu2n）（30）

uc=2d（31）

式中：hf表示液膜厚度，m；uc表示臨界速度，m·s－1；γ是控制參數(shù)，當ungt;uc時，使黏附率為0，γ計算式為

γ=1u2cln1θ+（32）

其中，θ+是使分母不為0、對數(shù)函數(shù)有意義的極小正值。

1.3"結冰模型和侵蝕模型

冰晶混合相結冰模型需要對經(jīng)典Messinger模型進行擴展，忽略冰層內(nèi)部、冰層與壁面之間的熱傳導，本文基于控制體積法［24］建立了質(zhì)量和能量平衡方程。

冰晶侵蝕將導致結冰冰形產(chǎn)生角狀或羽毛狀的粗糙突起，通過侵蝕效率給出冰晶侵蝕模型［13］

Φer=But，iu02βIWCβIWC0－A11－ηi1－ηi，0－A2DiD0A3·

exp－A4QsR1Ts－1Tm·exp－A5QsR1Ti－1Tm（1+（l0κ）2）（33）

式中：Φer表示冰晶侵蝕效率；ut，i表示冰晶撞擊切向速度，m·s－1；Di表示冰晶的平均質(zhì)量直徑，m；Ti表示冰晶溫度，K；κ表示冰層的局部曲率；c表示翼型弦長，m；Tf表示冰的相變溫度，K；R=8.314J·mol-1·K-1表示理想氣體常數(shù)。B、u0、βIWC0、ηi，0、D0、l0、A1～A5均是由實驗獲得的模型經(jīng)驗常數(shù)，具體取值見表1。

通過冰晶侵蝕效率得到冰晶侵蝕質(zhì)量流率

m·er=min［m·i+m·out，m·imp，imin（1，Φer）］（34）

式中：m·er表示冰晶侵蝕的質(zhì)量流率，kg·s－1；m·i表示控制體內(nèi)凍結的質(zhì)量流率，kg·s－1；m·out表示溢流水流出控制體的質(zhì)量流率，kg·s－1；m·imp，i表示撞擊冰晶的質(zhì)量流率，kg·s－1。

假設冰晶侵蝕造成的質(zhì)量損失與結冰速率和溢流水的積累速率成正比，得到控制體內(nèi)結冰的質(zhì)量流率與溢流水流出控制體的質(zhì)量流率的修正公式

m·*i=m·i－m·im·i+m·outm·er=m·i1－m·erm·i+m·out（35）

m·*out=m·out－m·outm·i+m·outm·er=m·out1－m·erm·i+m·out（36）

式中：m·*i表示控制體內(nèi)凍結的質(zhì)量流率修正值，kg·s－1；m·*out表示溢流水流出控制體的質(zhì)量流率修正值，kg·s－1。

2"數(shù)值模型驗證

本文采用FENSAP-ICE結冰數(shù)值計算軟件，湍流模型選擇Spalart-Allmaras單方程模型，考慮了非球形冰晶粒子的運動差異，基于歐拉雙流體模型計算了液滴和冰晶運動過程。本文還對比了NRC冰晶黏附模型和NIT冰晶黏附模型，通過改進的Messinger結冰熱力學模型，采用單步法對積冰增長及結冰表面水膜流動進行求解。

2.1"網(wǎng)格無關性驗證

由于目前只有美國COX結冰風洞公開了關于翼型混合相冰晶結冰的實驗結果［27］，為與實驗進行對比，本文選擇弦長為0.9144m、翼展為1.2192m的NACA0012翼型作為研究對象，對翼型與周圍流場使用ICEM劃分C型結構化網(wǎng)格。由于本文使用Spalart-Allmaras湍流模型，為滿足條件在翼型表面添加附面層，加密翼型的表面網(wǎng)格，滿足y+≈1原則，劃分網(wǎng)格如圖1所示。

用于網(wǎng)格無關性驗證和模型驗證的模擬工況參數(shù)詳見表2。表中：模擬工況001用于網(wǎng)格無關性驗證；模擬工況002～005與美國COX結冰風洞［27］公開的4種實驗工況相同，用于數(shù)值模型實驗驗證。

圖2（a）展示了001模擬工況下，使用不同網(wǎng)格數(shù)計算的駐點結冰厚度。圖2（b）展示了網(wǎng)格數(shù)為157萬時，不同時間的結冰冰形?？梢钥闯?，當網(wǎng)格數(shù)大于157萬時，駐點結冰厚度不再隨網(wǎng)格數(shù)發(fā)生變化，冰形也趨于穩(wěn)定，表現(xiàn)出由于水膜溢流形成的瘤狀凸起。因此，在滿足計算精度的前提下，選取157萬網(wǎng)格計算最為合適。

2.2"模型驗證

取本文數(shù)值模型計算的冰形與文獻［27］中實驗結果、Norde［28］的模擬結果進行對比，結果如圖3所示。可以看出，無論是工況002中過冷液滴結冰條件還是工況003中混合相冰晶結冰條件，本文的模擬結果與實驗結果在冰形形態(tài)上相近，且相較于Norde的模擬更為準確。本文模擬結果與實驗結果的結冰厚度、結冰范圍接近，可以證明本文數(shù)值方法的準確性。

3"計算結果與討論

3.1"液態(tài)水含量與冰晶含量之比對冰晶結冰特性的影響

βLWC和βIWC是影響冰晶結冰冰形的重要因素，其主要影響在于冰晶的黏附作用、侵蝕作用以及水膜溢流。為研究βIWC/βLWC對冰晶結冰的影響，定義總水含量（total water content，TWC，用符號βTWC表示）為LWC與IWC之和，選擇βTWC=1.4g·m－3的模擬工況021～027進行計算，模擬工況如表3所示。

圖4為不同βIWC/βLWC下NACA0012翼型的混合相結冰冰形?？梢钥闯觯寒擳WC一定時，隨著IWC的增大，積結冰覆蓋面積逐漸減??；由于冰晶的黏附、侵蝕以及水膜的流動，冰層的形狀逐漸從冠狀變向尖角狀；當IWC較低時，水膜流動起主要作用，當IWC較高時，冰晶的黏附和侵蝕起主要作用。當IWC為0、LWC為1.4g·m－3時，為水滴結冰，此時壁面所收集的液滴會匯聚形成液膜，液膜在氣流剪切作用下向后溢流并逐漸凍結形成冠狀冰形。隨著LWC減小，液膜溢流范圍逐漸減少，所以結冰區(qū)域也隨之縮減。液膜溢流量的減少使得冰晶的黏附和侵蝕作用變得更加顯著，由此導致的冠狀冰特征逐漸減弱，而冰形開始展現(xiàn)出與翼型表面幾何輪廓相一致的特征。當IWC為1.4g·m－3、LWC為0時，不發(fā)生冰晶結冰現(xiàn)象，這一結果強調(diào)了液態(tài)水在冰晶結冰過程中的必要性。

圖5為βIWC/βLWC對冰晶結冰厚度的影響?？梢钥闯?，在TWC保持恒定的情況下，NACA0012型翼型的駐點結冰厚度和最大結冰厚度均隨IWC的增加呈現(xiàn)出先增后減的趨勢。在較高IWC條件下，翼型駐點與最大結冰厚度處的結冰厚度差逐漸減小，趨向于0，表明最大結冰區(qū)域逐步移至駐點附近。

此現(xiàn)象可以通過水膜溢流與收集系數(shù)解釋：翼型駐點位置的冰晶收集系數(shù)最高，當LWC較高時，翼型表面收集的液滴并不能迅速凍結，水膜在氣流剪切作用下向后溢流導致駐點結冰厚度小于最大結冰厚度；當IWC較高時，只有很少的水膜向后溢流，此時最大結冰處位于駐點。

3.2"黏附模型對冰晶結冰的影響

為探究黏附模型對冰晶結冰過程的影響，本文對表2中的004和005工況進行模擬計算。在本研究中，無黏附模型是假設所有冰晶在碰撞到翼型表面后均不發(fā)生反彈或破碎，完全黏附于表面。

圖6展示了采用不同黏附模型計算得到的結冰冰形與實驗觀測結果的比較。在無黏附模型時，由于未考慮到冰晶的反彈和破碎，模擬得到的結冰冰形與實驗數(shù)據(jù)存在顯著偏差。這一發(fā)現(xiàn)表明，在結冰模擬過程中必須考慮冰晶的反彈和破碎效應，因此引入精確的黏附模型以評估冰晶碰撞和黏附過程是至關重要的。在高IWC條件下，NTI黏附模型能更準確的預測冰形，其模擬結果中因水膜溢流引起的冠狀凸起特征與實驗結果相符。相比之下，NRC模型的預測結果與實驗數(shù)據(jù)存在較大偏差，模擬出的冰形呈現(xiàn)出不真實的雙角狀凸起。然而，在低IWC工況，NTI與NRC兩種黏附模型均能較好地模擬出冰形，其結果均與翼型前緣的實際結冰輪廓相吻合。

圖7是不同黏附模型結冰厚度分布?？梢钥闯觯寒敓o冰晶黏附模型時，翼型結冰厚度遠大于實驗結果；對于高IWC工況，NRC黏附模型的駐點結冰厚度更接近實驗結果，誤差為7.2%；對于低IWC工況，NTI模型與NRC模型的駐點結冰厚度接近，誤差分別為11.6%和6.8%。

3.3"侵蝕模型對冰晶結冰的影響

冰晶侵蝕對翼型上結冰形態(tài)和厚度具有顯著影響。在所采用的冰晶侵蝕模型中，冰晶粒徑和速度均與侵蝕程度呈正相關。為直觀體現(xiàn)侵蝕作用對冰晶結冰的影響，對表2中002～005的4種工況分別采用冰晶侵蝕模型進行模擬。

鑒于工況002中的IWC為0，在此條件下不會發(fā)生冰晶侵蝕，因此以下分析主要集中在其他3組工況下。圖8展示了侵蝕模型對結冰形態(tài)的具體影響。

通過對比圖8（a）和圖8（b）可知，在IWC相同的條件下，當LWC較高時，侵蝕作用對結冰形態(tài)的影響更為顯著。如圖8（b）所示，當LWC為0.3g·m－3時，由于LWC處于較低水平，結冰量相對較小，有無侵蝕模型對結冰模擬的影響不甚明顯。相反，在圖8（a）中，LWC較高，侵蝕作用對結冰形態(tài)的影響更為顯著。同時，不同模型的計算結果在翼型駐點附近顯示出明顯的差異，這表明冰晶侵蝕作用的主要影響區(qū)域集中在翼型駐點附近。

進一步地，對比圖8（a）和圖8（c）可以看出：在LWC相同的條件下，較高的IWC同樣會加劇侵蝕效應對結冰形態(tài)的影響，即LWC和IWC的增大均會加劇侵蝕效應對結冰形態(tài)的影響。這是因為冰晶侵蝕造成的損失與冰晶含量、結冰速率和溢流水的積累速率成正比。當LWC較高時，結冰速率較大，導致侵蝕現(xiàn)象更顯著；當IWC較高時，冰晶沖擊會進一步增加質(zhì)量損失，從而導致結冰形態(tài)發(fā)生顯著變化。

4"結"論

本文結合混合相態(tài)冰晶結冰的熱力學原理，基于FENSAP-ICE平臺對混合相態(tài)冰晶現(xiàn)象進行數(shù)值計算，通過與美國COX結冰風洞的NACA0012翼型結冰實驗結果對比，驗證了本文數(shù)值模型的正確性，并在此基礎上研究了混合相中液態(tài)水含量對冰晶結冰特性的影響規(guī)律，分析了兩種黏附模型的差異，討論了侵蝕效應對結冰的影響，具體結論如下。

（1）當TWC一定時，隨著IWC的增加，結冰覆蓋面積逐漸減小，且冰層的形狀逐漸從冠狀變向尖角狀。在IWC較小的情況下，水膜流動是影響結冰形態(tài)的主要因素；在IWC較大的情況下，冰晶的黏附和侵蝕是影響結冰形態(tài)的主要因素。當TWC相同時，翼型的駐點結冰厚度與最大結冰厚度隨IWC的增大均呈現(xiàn)出先增后減的趨勢。當IWC較大時，翼型駐點處的結冰厚度與最大結冰厚度之間的差異逐漸減小，趨近于0，表明最大結冰區(qū)域逐步移至駐點附近。

（2）考慮冰晶黏附效應時，翼型結冰形狀和駐點結冰厚度與實驗結果的符合性較好。對于高IWC工況，NTI模型能更準確地模擬冰形，其預測的由于水膜溢流形成的冠狀凸起特征與實驗觀測相似。NRC模型在駐點結冰厚度的計算上與實驗數(shù)據(jù)更為接近，其誤差僅為7.2%。在低IWC工況下，NTI模型與NRC模型在計算冰形與駐點結冰厚度方面均顯示出較高的準確性。兩種模型的計算結果均與翼型前緣輪廓接近，且駐點結冰厚度計算誤差分別控制在11.6%和6.8%以內(nèi)。

（3） LWC或IWC增大，均會加劇侵蝕效應對結冰形態(tài)的影響。侵蝕作用的影響主要集中于翼型駐點附近，這是因為該區(qū)域的收集系數(shù)較高、結冰速率較快，所以侵蝕現(xiàn)象更為明顯。

參考文獻：

［1］LYNCH F T， KHODADOUST A. Effects of ice accretions on aircraft aerodynamics ［J］. Progress in Aerospace Sciences， 2001， 37（8）： 669-767.

［2］MASON J G， STRAPP J W， CHOW P. The ice particle threat to engines in flight ［C］//44th AIAA Aerospace Sciences Meeting and Exhibit. Reston， VA， USA： AIAA， 2006： AIAA 2006-206.

［3］沈浩， 韓冰冰， 張麗芬. 航空發(fā)動機中冰晶結冰的研究進展 ［J］. 實驗流體力學， 2020， 34（6）： 1-7.

SHEN Hao， HAN Bingbing， ZHANG Lifen. Research progress of the ice crystal icing in aero-engine ［J］. Journal of Experiments in Fluid Mechanics， 2020， 34（6）： 1-7.

［4］TETTEH E， LOTH E， NEUTEBOOM M O， et al. In-flight gas turbine engine icing： review ［J］. AIAA Journal， 2022， 60（10）： 5610-5632.

［5］黃平， 卜雪琴， 劉一鳴， 等. 混合相/冰晶條件下的結冰研究綜述 ［J］. 航空學報， 2022， 43（5）： 112-130.

HUANG Ping， BU Xueqin， LIU Yiming， et al. Mixed phase/glaciated ice accretion： review ［J］. Acta Aeronautica et Astronautica Sinica， 2022， 43（5）： 112-130.

［6］袁慶浩， 樊江， 白廣忱. 航空發(fā)動機內(nèi)部冰晶結冰研究綜述 ［J］. 推進技術， 2018， 39（12）： 2641-2650.

YUAN Qinghao， FAN Jiang， BAI Guangchen. Review of ice crystal icing in aero-engines ［J］. Journal of Propulsion Technology， 2018， 39（12）： 2641-2650.

［7］VILLEDIEU P， TRONTIN P， GUFFOND D， et al. SLD Lagrangian modeling and capability assessment in the frame of ONERA 3D icing suite ［C］//4th AIAA Atmospheric and Space Environments Conference. Reston， VA， USA： AIAA， 2012： AIAA 2012-3132.

［8］WRIGHT W B. Further refinement of the LEWICE SLD model ［C］//44th AIAA Aerospace Sciences Meeting and Exhibit. Reston， VA， USA： AIAA， 2006： AIAA 2006-464.

［9］HONSEK R， HABASHI W G. FENSAP-ICE： Eulerian modeling of droplet impingement in the SLD regime of aircraft icing ［C］//44th AIAA Aerospace Sciences Meeting and Exhibit. Reston， VA， USA： AIAA， 2006： AIAA 2006-465.

［10］VILLEDIEU P， TRONTIN P， CHAUVIN R.Glaciated and mixed phase ice accretion modeling using ONERA 2D icing suite ［C］//6th AIAA Atmospheric and Space Environments Conference. Reston， VA， USA： AIAA， 2014： AIAA 2014-2199.

［11］NORDE E， SENONER J M，VAN DER WEIDE E T A， et al. Eulerian and Lagrangian ice-crystal trajectory simulations in a generic turbofan compressor ［J］. Journal of Propulsion and Power， 2019， 35（1）： 26-40.

［12］NILAMDEEN S， HABASHI W G， AUB M S， et al. FENSAP-ICE： modeling of water droplets and ice crystals ［C］//1st AIAA Atmospheric and Space Environments Conference. Reston， VA， USA： AIAA， 2009： AIAA 2009-4128.

［13］TRONTIN P， BLANCHARD G， VILLEDIEU P.A comprehensive numerical model for mixed-phase and glaciated icing conditions ［C］//8th AIAA Atmospheric and Space Environments Conference. Reston， VA， USA： AIAA， 2016： AIAA 2016-3742.

［14］NILAMDEEN S， HABASHI W G.Multiphase approach toward simulating ice crystal ingestion in jet engines ［J］. Journal of Propulsion and Power， 2011， 27（5）： 959-969.

［15］CURRIE T C， STRUK P M， TSAO J C， et al. Fundamental study of mixed-phase icing with application to ice crystal accretion in aircraft jet engines ［C］//4th AIAA Atmospheric and Space Environments Conference. Reston， VA， USA： AIAA， 2012： AIAA 2012-3035.

［16］STRUK P M，RATVASKY T P， BENCIC T J， et al. An initial study of the fundamentals of ice crystal icing physics in the NASA propulsion systems laboratory ［C］//9th AIAA Atmospheric and Space Environments Conference. Reston， VA， USA： AIAA， 2017： AIAA 2017-4242.

［17］CURRIE T C， FULEKI D， MAHALLATI A.Experimental studies of mixed-phase sticking efficiency for ice crystal accretion in jet engines ［C］//6th AIAA Atmospheric and Space Environments Conference. Reston， VA， USA： AIAA， 2014： AIAA 2014-3049.

［18］ZHANG Lifen， LIU Zhenxia， ZHANG Meihua.Numerical simulation of ice accretion under mixed-phase conditions ［J］. Proceedings of the Institution of Mechanical Engineers： Part G"Journal of Aerospace Engineering， 2016， 230（13）： 2473-2483.

［19］馬乙楗， 柴得林， 王強， 等. 翼面結冰過程中的冰晶運動相變與黏附特性 ［J］. 航空學報， 2023， 44（1）： 35-46.

MA Yijian， CHAI Delin， WANG Qiang， et al. Phase change and adhesion characteristics of ice crystal movements in wing icing ［J］. Acta Aeronautica et Astronautica Sinica， 2023， 44（1）： 35-46.

［20］姜飛飛， 董威， 鄭梅， 等. 冰晶在渦扇發(fā)動機內(nèi)相變換熱特性 ［J］. 航空動力學報， 2019， 34（3）： 567-575.

JIANG Feifei， DONG Wei， ZHENG Mei， et al. Phase change heat transfer characteristic of ice crystal ingested into turbofan engine ［J］. Journal of Aerospace Power， 2019， 34（3）： 567-575.

［21］CHEN Jiajun， LIU Xiufang， ZHONG Fuhao， et al. Melting and accumulation characteristics of ice crystals in aero-engines ［J］. Transactions of Nanjing University of Aeronautics amp; Astronautics， 2023， 40（6）： 653-662.

［22］ZHONG Fuhao， WEI Zhen， CHEN Jiajun， et al. Dynamic behavior and heat transfer characteristics of non-spherical ice crystals in high-temperature air flow ［J］. Transactions of Nanjing University of Aeronautics amp; Astronautics， 2023， 40（2）： 205-215.

［23］QI Haifeng， CHANG Shinan， YANG Yinglin. Numerical study of mixed phase ice accumulation in aero-engine inlet system ［J］. Applied Thermal Engineering， 2023， 231： 120909.

［24］卜雪琴， 李皓， 黃平， 等. 二維機翼混合相結冰數(shù)值模擬 ［J］. 航空學報， 2020， 41（12）： 195-205.

BU Xueqin， LI Hao， HUANG Ping， et al. Numerical simulation of mixed phase icing on two-dimensional airfoil ［J］. Acta Aeronautica et Astronautica Sinica， 2020， 41（12）： 195-205.

［25］郭琪磊， ?？〗?， 安博， 等. 混合相態(tài)冰晶積冰的數(shù)值研究 ［J］. 空氣動力學學報， 2021， 39（2）： 168-175.

GUO Qilei， NIU Junjie， AN Bo， et al. Numerical simulation of ice crystal icing under mixed-phase conditions ［J］. Acta Aerodynamica Sinica， 2021， 39（2）： 168-175.

［26］譚燕. 基于歐拉方法的2維翼型冰晶結冰數(shù)值計算 ［J］. 航空發(fā)動機， 2020， 46（4）： 30-35.

TAN Yan. Numerical calculation of 2D airfoil ice crystal icing based on Euler method ［J］. Aeroengine， 2020， 46（4）： 30-35.

［27］AL-KHALIL K， IRANI E， MILLER D. Mixed phase icing simulation and testing at the cox icing wind tunnel ［C］//41st Aerospace Sciences Meeting and Exhibit. Reston， VA， USA： AIAA， 2003： AIAA 2003-903.

［28］NORDE E.Eulerian method for ice crystal icing in turbofan engines ［D］. Enschede： University of Twente， 2017： 99-112.

（編輯"陶晴）