機(jī)翼混合相冰晶結(jié)冰現(xiàn)象的數(shù)值研究

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

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

中圖分類號(hào)：V211.41"文獻(xiàn)標(biāo)志碼：A

DOI：10.7652/xjtuxb202410015"文章編號(hào)：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

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

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

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

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

侵蝕效應(yīng)是冰晶結(jié)冰區(qū)別于過冷水滴結(jié)冰的獨(dú)有現(xiàn)象，即來流中的冰晶對(duì)已形成的冰層具有侵蝕作用。Struk等［15-16］以及Currie等［17］在研究冰晶結(jié)冰現(xiàn)象時(shí)觀察到了侵蝕效應(yīng)，并發(fā)現(xiàn)侵蝕程度與冰晶的融化率和粒徑密切相關(guān)。在此基礎(chǔ)上，Trontin等［13］基于實(shí)驗(yàn)數(shù)據(jù)提出了冰晶侵蝕半經(jīng)驗(yàn)?zāi)Ｐ停撃Ｐ透鶕?jù)冰晶含量、粒徑等參數(shù)來預(yù)測(cè)冰晶侵蝕效果。

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

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

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

1.1"冰晶運(yùn)動(dòng)計(jì)算方法

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

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

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

0.4，Redgt;1300（1）

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

Rej=ρa(bǔ)ua－ujdμa（2）

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

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

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）

當(dāng)0.01lt;ReD≤1.5時(shí)

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

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

ω=lgReD（7）

當(dāng)1.5lt;ReD≤100時(shí)

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

當(dāng)100lt;ReD≤300時(shí)

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

當(dāng)ReDgt;300時(shí)

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為冰晶半短軸長(zhǎng)度，m；b為冰晶半長(zhǎng)軸長(zhǎng)度，m；deq為冰晶體積等效直徑，m。

1.2"冰晶碰撞黏附模型

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

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

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

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

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

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

將無(wú)量綱參數(shù)用兩個(gè)臨界值將冰晶碰撞行為分為完全彈性碰撞、非彈性碰撞和碰撞破碎3種情況，并定義恢復(fù)系數(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表示撞擊時(shí)液滴的質(zhì)量分?jǐn)?shù)；φi表示撞擊時(shí)冰晶的質(zhì)量分?jǐn)?shù)。ηi、φi、φd、Ki、Kd表達(dá)式為

η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表示冰晶開始結(jié)冰的壁面溫度；Ki和Kd通過風(fēng)洞結(jié)冰實(shí)驗(yàn)擬合得到。

二次碰撞特性為

ξnn=1，ψ≤ψc1

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

ξnt=0.41－ψc2ψ（28）

dmax=ψ211c2ψ（29）

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

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

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

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

uc=2d（31）

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

γ=1u2cln1θ+（32）

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

1.3"結(jié)冰模型和侵蝕模型

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

冰晶侵蝕將導(dǎo)致結(jié)冰冰形產(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表示翼型弦長(zhǎng)，m；Tf表示冰的相變溫度，K；R=8.314J·mol-1·K-1表示理想氣體常數(shù)。B、u0、βIWC0、ηi，0、D0、l0、A1～A5均是由實(shí)驗(yàn)獲得的模型經(jīng)驗(yàn)常數(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)凍結(jié)的質(zhì)量流率，kg·s－1；m·out表示溢流水流出控制體的質(zhì)量流率，kg·s－1；m·imp，i表示撞擊冰晶的質(zhì)量流率，kg·s－1。

假設(shè)冰晶侵蝕造成的質(zhì)量損失與結(jié)冰速率和溢流水的積累速率成正比，得到控制體內(nèi)結(jié)冰的質(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)凍結(jié)的質(zhì)量流率修正值，kg·s－1；m·*out表示溢流水流出控制體的質(zhì)量流率修正值，kg·s－1。

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

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

2.1"網(wǎng)格無(wú)關(guān)性驗(yàn)證

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

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

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

2.2"模型驗(yàn)證

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

3"計(jì)算結(jié)果與討論

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

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

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

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

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

3.2"黏附模型對(duì)冰晶結(jié)冰的影響

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

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

圖7是不同黏附模型結(jié)冰厚度分布?？梢钥闯觯寒?dāng)無(wú)冰晶黏附模型時(shí)，翼型結(jié)冰厚度遠(yuǎn)大于實(shí)驗(yàn)結(jié)果；對(duì)于高IWC工況，NRC黏附模型的駐點(diǎn)結(jié)冰厚度更接近實(shí)驗(yàn)結(jié)果，誤差為7.2%；對(duì)于低IWC工況，NTI模型與NRC模型的駐點(diǎn)結(jié)冰厚度接近，誤差分別為11.6%和6.8%。

3.3"侵蝕模型對(duì)冰晶結(jié)冰的影響

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

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

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

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

4"結(jié)"論

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

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

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

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

參考文獻(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ā)動(dòng)機(jī)中冰晶結(jié)冰的研究進(jìn)展 ［J］. 實(shí)驗(yàn)流體力學(xué)， 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］黃平， 卜雪琴， 劉一鳴， 等. 混合相/冰晶條件下的結(jié)冰研究綜述 ［J］. 航空學(xué)報(bào)， 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ā)動(dòng)機(jī)內(nèi)部冰晶結(jié)冰研究綜述 ［J］. 推進(jìn)技術(shù)， 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］馬乙楗， 柴得林， 王強(qiáng)， 等. 翼面結(jié)冰過程中的冰晶運(yùn)動(dòng)相變與黏附特性 ［J］. 航空學(xué)報(bào)， 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ā)動(dòng)機(jī)內(nèi)相變換熱特性 ［J］. 航空動(dòng)力學(xué)報(bào)， 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］卜雪琴， 李皓， 黃平， 等. 二維機(jī)翼混合相結(jié)冰數(shù)值模擬 ［J］. 航空學(xué)報(bào)， 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］. 空氣動(dòng)力學(xué)學(xué)報(bào)， 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維翼型冰晶結(jié)冰數(shù)值計(jì)算 ［J］. 航空發(fā)動(dòng)機(jī)， 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.

（編輯"陶晴）