卜雪琴，李皓，黃平，林貴平

北京航空航天大學 航空科學與工程學院，北京 100083

從現(xiàn)代飛機飛行開始，飛機結(jié)冰問題就逐漸引起人們的重視。早期的研究主要關(guān)注過冷水滴結(jié)冰，隨著研究地不斷深入，飛機結(jié)冰研究也逐漸由過冷水滴結(jié)冰擴展到冰晶結(jié)冰、冰晶和過冷水滴同時存在的混合相結(jié)冰[1-3]。Mason等[4]總結(jié)發(fā)現(xiàn)渦扇發(fā)動機的動力損失與冰晶密切相關(guān)。冰晶結(jié)冰的機理尚不明確，目前普遍認為冰晶引發(fā)結(jié)冰需要液態(tài)水的存在[5]。盡管純冰晶條件下，冰晶撞擊到機翼、尾翼等冷表面會彈開而無法結(jié)冰，但撞擊到傳感器等熱表面或在發(fā)動機內(nèi)部熱環(huán)境下運動時，冰晶會(部分)融化成液態(tài)水，液態(tài)水使得冰晶黏附凍結(jié)，同時冰晶融化吸熱使得熱表面、熱環(huán)境溫度下降。傳感器結(jié)冰會導致傳感器精度下降。發(fā)動機內(nèi)部核心部件冰晶結(jié)冰嚴重者可能引發(fā)發(fā)動機動力下降、喘振、熄火等危險，同時冰層脫落可能造成發(fā)動機的機械損傷。研究表明，冰晶結(jié)冰時經(jīng)常伴隨著總溫傳感器異常[4]。此外，在同時含有冰晶及過冷水滴的混合相條件下，由于液態(tài)水的存在，結(jié)冰不僅可能發(fā)生在機翼、尾翼等冷表面，造成飛機氣動外形的改變，降低升力、增大阻力、影響飛機穩(wěn)定性，而且可能發(fā)生在傳感器、發(fā)動機內(nèi)部壓氣機等熱表面，使傳感器指示異常，對發(fā)動機構(gòu)成威脅。由于冰晶的存在，冰晶結(jié)冰傳熱傳質(zhì)模型和過冷水滴結(jié)冰模型有顯著區(qū)別。另外，由于冰晶形狀各異，可能是扁平盤狀、柱狀等非球形，其運動以及在熱環(huán)境中運動的融化相變也是一個復雜過程。

冰晶/混合相結(jié)冰的試驗研究需要先進的冰風洞作為支持，存在的難點包括冰晶的制備與運輸、真實結(jié)冰氣象條件的模擬等，具有成本高、周期長的特點[6-8]。而采用數(shù)值模擬能大大降低人力、物力等資源消耗，節(jié)約時間成本，因此也是冰晶/混合相結(jié)冰研究的重要手段。目前國外針對冰晶/混合相結(jié)冰數(shù)值模擬問題已經(jīng)展開廣泛研究，在已有結(jié)冰軟件基礎(chǔ)上進行了改進和補充，使其能夠應用到冰晶結(jié)冰問題中來，已經(jīng)取得了不錯的進展[9-12]。而國內(nèi)研究相對較少，Zhang等[13]使用FLUENT軟件研究了混合相結(jié)冰問題。袁慶浩等[14]對航空發(fā)動機內(nèi)部冰晶結(jié)冰文獻進行了調(diào)研和總結(jié)，綜述了冰晶結(jié)冰和過程水結(jié)冰的區(qū)別、冰晶結(jié)冰試驗和計算方法等。姜飛飛等[15]計算分析了在渦扇發(fā)動機內(nèi)涵道內(nèi)運動過程中冰晶粒子的傳熱傳質(zhì)情況，得到了冰晶粒子的溫度、半徑、液態(tài)水質(zhì)量分數(shù)的變化結(jié)果。

本文暫不考慮冰晶的非球形特性，認為是球形，旨在建立混合相條件下的結(jié)冰熱力學模型和數(shù)值計算方法，并開展二維翼型混合相條件下的結(jié)冰模擬?；贛essinger模型進行擴展，考慮冰晶黏附效應，建立了混合相的結(jié)冰熱力學模型。首先采用FLUENT軟件計算空氣流場，然后基于單向耦合假設(shè)求解歐拉法下的冰晶運動控制方程，獲得冰晶撞擊特性，最后分析結(jié)冰過程中的傳熱傳質(zhì)，通過FLUENT提供的用戶自定義函數(shù)(UDF)編程求解混合相熱力學模型，實現(xiàn)混合相結(jié)冰預測。將數(shù)值計算結(jié)果與文獻試驗數(shù)據(jù)進行對比，分析計算方法的有效性。

1 混合相結(jié)冰熱力學模型

1.1 質(zhì)量守恒

圖1 控制體質(zhì)量守恒Fig.1 Mass balance of control volume

建立質(zhì)量守恒方程：

(1)

(2)

同時，還有如下關(guān)系：

(3)

(4)

(5)

式中，U∞,ic和U∞,d分別為遠場的冰晶速度和水滴速度；β為收集系數(shù)；ηic為冰晶融化比，ηic=LWCic/IWC，LWCic表示部分融化冰晶中液態(tài)水的含量，IWC表示空氣流場中的冰晶含量；A為控制體的底面積；LWCd為空氣流場中過冷水滴的液態(tài)水含量。

(6)

式中：hc為對流換熱系數(shù)；ρ為空氣密度；cp為比熱容；Pr為普朗特數(shù)；hm為對流傳質(zhì)系數(shù)；Sc為施密特數(shù)，物理意義為動量擴散與質(zhì)量擴散之比，即

(7)

式中：v和μ分別表示空氣的運動和動力黏度；D為擴散系數(shù)。

利用氣體狀態(tài)方程可以推導出對流傳質(zhì)的質(zhì)量流量計算式[18]，當表面有溢流水時為蒸發(fā)傳質(zhì)，表面僅為冰時為升華傳質(zhì)：

(8)

式中：Le為路易斯數(shù)，Le=Sc/Pr；MW為相對分子質(zhì)量；Ts為壁面溫度；pv,s,sat為壁面平衡溫度下的飽和蒸汽壓；pT為總壓；TT為總溫；pe為邊界層外邊界的空氣壓力；rh為相對濕度，本文中取值為1。

1.2 能量守恒

建立能量守恒方程：

(9)

能量守恒方程式(9)中，每一項都有特定的計算式，首先有進入控制體各質(zhì)量項的動能：

圖2 控制體內(nèi)能量守恒Fig.2 Energy balance of control volume

(10)

(11)

式中：Uimp,d和Uimp,ic分別為水滴和冰晶粒子撞擊到結(jié)冰表面時候的速度。

控制體內(nèi)涉及水的凍結(jié)、水的蒸發(fā)和冰的升華3種相變過程(水的蒸發(fā)和冰的升華取其一，表面有溢流水則為蒸發(fā)散熱，表面只有冰時則為升華散熱)，因此存在相變潛熱項：

(12)

(13)

(14)

式中：Lf為結(jié)冰潛熱，取值333 kJ/kg；Lev為蒸發(fā)潛熱，取值2 500 kJ/kg；Lsub為升華潛熱，取值2 833 kJ/kg。

對流換熱引起的能量傳遞：

(15)

式中：Tinf為自由流溫度。

(16)

(17)

(18)

(19)

(20)

式中：Tm為融化溫度，即0 ℃；cp,w為水的比熱容，設(shè)為常數(shù)4 174 J/kg·K；cp,ic為冰的比熱容，設(shè)置為常數(shù)2 102 J/kg·K。

2 冰晶黏附模型

通常情況下認為過冷水滴撞擊表面后會全部停留在表面參與結(jié)冰過程，而冰晶存在很大不同，試驗表明冰晶撞擊表面后是否黏附與液態(tài)水的存在有很大關(guān)系[19]。冰晶撞擊表面可能的結(jié)果有：反彈、破碎和黏附[20-22]。Baumert等[23]的冰晶撞擊NACA0012翼型表面的試驗表明反彈后的粒子除駐點附近極少能再次撞擊到壁面，對結(jié)冰過程影響不大，因此主要關(guān)注冰晶撞擊后的黏附效應。

Trontin等[24]根據(jù)冰風洞試驗結(jié)果，包括純冰晶部分融化條件結(jié)冰試驗以及同時含有過冷水滴和冰晶的混合相結(jié)冰試驗，分析液態(tài)水含量對冰晶黏附撞擊表面可能性的影響，建立起符合統(tǒng)計學規(guī)律和冰風洞試驗的經(jīng)驗公式。Trontin等[24]指出根據(jù)加拿大NRC (National Research Council) 的試驗結(jié)果，相同LWC/TWC時(TWC為冰水總含量)，部分融化的純冰晶條件下的黏附系數(shù)要高于混合相條件下的黏附系數(shù)，因此引入了εs,ic和εs,d，分別代表只考慮融化冰晶且無液態(tài)過冷水滴的影響以及考慮有液態(tài)水滴和(部分融化)冰晶共同的影響，黏附系數(shù)ε取兩者中較大者：

ε=max(εs,ic,εs,d)

(21)

式中：

(22)

εs,d=Kd(Yd+ηicYic)

(23)

其中：ηic為冰晶融化比，根據(jù)Currie等[25]的純冰晶環(huán)境下的試驗，這里取0<ηic<0.2；Yd為撞擊到表面的過冷水滴質(zhì)量占撞擊到表面的總質(zhì)量的比值；Yic是撞擊到表面的冰晶質(zhì)量占撞擊到表面的總質(zhì)量的比值；Yd和Yic兩者之和為1；Kic和Kd是常數(shù)。

(24)

(25)

考慮到表面溫度越低時，冰晶黏附越少，特別是溫度較低的干態(tài)冰表面時，冰晶不會黏附。于是，黏附模型中的常數(shù)Kd和表面溫度有關(guān)：

(26)

可以發(fā)現(xiàn)，式(23)中的Yd+ηicYic代表了撞擊到表面的所有液態(tài)水(包括過冷水滴和融化的冰晶)占總質(zhì)量的比值。因此當混合相中冰晶不發(fā)生融化時，ηic=0，則ε=max(εs,ic,εs,d)=εs,d=KdYd。

本文計算的混合相結(jié)冰算例中環(huán)境溫度較低，認為冰晶在運動過程中不發(fā)生融化，LWCic=0，式(21)的黏附系數(shù)ε則簡化為ε=KdYd。

考慮冰晶黏附效應后，撞擊冰晶中實際黏附在結(jié)冰表面參與結(jié)冰過程的冰晶質(zhì)量流量為

(27)

3 模型求解方法

3.1 整體計算方法

結(jié)冰熱力學模型用C語言編寫，并通過UDF掛靠到FLUENT中完成編譯和執(zhí)行。首先計算空氣流場，在空氣流場的基礎(chǔ)上計算水滴和冰晶的流場，獲得水滴和冰晶的收集系數(shù)，隨后把計算得到的剪切力、對流換熱系數(shù)、壓強、收集系數(shù)等作為輸入?yún)?shù)計算結(jié)冰熱力學模型，獲得最終的結(jié)冰量。計算流程如圖3所示。

圖3 結(jié)冰計算流程Fig.3 Icing calculation flowchart

3.2 結(jié)冰熱力學模型求解方法

對于二維翼型，認為駐點處流入控制體的質(zhì)量流量為0，駐點所在控制體內(nèi)的水平均分為兩部分分別向上下表面溢流；流入流出控制體的水之間存在約束關(guān)系，即前一個控制體流出水的質(zhì)量流量等于后一個控制體流入水的質(zhì)量流量：

(28)

(29)

混合相結(jié)冰熱力學模型求解的具體方法如下：將結(jié)冰表面從駐點處分為上表面和下表面兩部分，首先計算駐點所在的控制體，然后分別針對上表面和下表面按順序逐個求解駐點下游的控制體，直到表面水溢流極限處，完成全部控制體的求解。對每一個控制體來說，結(jié)冰過程熱力學模型的求解通過迭代過程來實現(xiàn)。定義流入流出控制體能量的凈值Q*=進入控制體的熱量-流出控制體的熱量，即能量守恒方程左右相減的差值，通過迭代使得Q*為0即滿足控制體能量守恒，完成能量守恒方程的求解。

(30)

式中：t是結(jié)冰時間；ρice是冰的密度；A是控制體的面積(二維情況下為長度)。

4 計算結(jié)果及對比驗證

取Al-Khalil在Cox冰風洞做的混合相結(jié)冰試驗作為驗證算例[19]。該試驗采用了NACA0012翼型，弦長0.914 4 m，水滴平均直徑為20 μm，冰晶平均直徑分別為150 μm和200 μm。結(jié)冰試驗條件如表1所示。其中算例1和算例2環(huán)境溫度較高，算例2比算例1增加了冰晶含量，其他參數(shù)完全相同；算例3和算例4環(huán)境溫度較低，兩者的液態(tài)水含量和冰晶含量不同，算例3的液態(tài)水含量少冰晶含量多，而算例4的液態(tài)水含量多冰晶含量少。

表1 結(jié)冰計算條件Table 1 Calculation conditions of icing

4.1 空氣流場計算

空氣流場的計算劃分結(jié)構(gòu)網(wǎng)格、翼型周圍網(wǎng)格如圖4所示。網(wǎng)格數(shù)量為24 477個，計算域半徑10 m，保證遠場邊界條件的成立。選取了底層網(wǎng)格高度0.1 mm、0.05 mm、0.005 mm和0.001 mm進行網(wǎng)格無關(guān)性計算，最終確定底層網(wǎng)格為0.005 mm的網(wǎng)格能達到要求，此時壁面附近的y+≈1。空氣為理想氣體，采用轉(zhuǎn)捩剪切應力輸運(SST)湍流模型，假設(shè)表面等效沙粒高度為1 mm，使用SIMPLE算法，為了提高計算精度，除壓力采用標準格式外，其他參數(shù)全部采用二階迎風格式。

圖4 計算網(wǎng)格Fig.4 Computation mesh

由于計算參數(shù)相近，4個算例的計算結(jié)果十分相近，算例1中翼型周圍空氣速度和溫度分布如圖5所示。由于攻角為0°，速度和溫度分布均為上下對稱，駐點處速度為0 m/s，溫度最高，較遠場處升高2 ℃，但仍低于0 ℃，這也說明冰晶在運動過程中不可能會發(fā)生融化，前面假設(shè)冰晶運動過程中不融化是合理的。

4.2 冰晶/水滴撞擊特性計算

混合相中既有過冷水滴又有冰晶，本文有如下假設(shè)：① 水滴或冰晶粒子在空氣中的含量小(小于0.7 g/m3)，各自以及混合相的容積分數(shù)(α=Vpriticle/Vmixture)在10-6或10-7量級，可認為水滴和冰晶互不影響，空氣和水滴以及空氣和冰晶之間作用均是單向耦合，忽略粒子對空氣流場的影響；② 環(huán)境溫度較低，忽略水滴、冰晶運動過程中的變形、蒸發(fā)、升華、冰晶融化、粒子碰撞融合等現(xiàn)象；③ 遠場處粒子的初始速度和周圍空氣相同。

圖5 空氣速度和溫度分布云圖Fig.5 Contours of air velocity and temperature distribution

水滴粒子的撞擊特性計算采用歐拉法，基于利用課題組在FLUENT基礎(chǔ)上二次開發(fā)的用戶自定義標量的方法，詳見文獻[26]，在此僅作簡單描述。在空氣流場計算收斂的基礎(chǔ)上，自定義了包括水滴容積分數(shù)、水滴x、y方向速度的3個標量；阻力系數(shù)采用球形阻力系數(shù)，對水滴運動的質(zhì)量、動量方程進行描述，利用FLUENT求解器進行求解，采用二階迎風格式；為了不讓水滴在撞擊壁面處堆積使得計算發(fā)散，壁面采用自定義的排水邊界。局部水滴收集系數(shù)計算式為

(31)

式中：αd和α∞,d分別是壁面當?shù)睾妥杂蓙砹髦兴稳莘e分數(shù)；ud是壁面當?shù)厮嗡俣仁噶?；n是壁面單位法向量。

為了盡可能接近真實情況，水滴撞擊特性的計算采用7種不同直徑的水滴計算結(jié)果取加權(quán)平均：

(32)

式中：pd,i為第i種直徑水滴所占的液態(tài)水含量比重；βd,i是相應水滴直徑的局部水收集系數(shù)。水滴直徑d分布如表2所示。

表2 水滴直徑分布Table 2 Distribution of droplet diameters

算例1的計算結(jié)果如圖6所示，圖中s表示弧長，c表示弦長。其他算例計算結(jié)果相近。水滴直徑顯著影響撞擊極限和局部收集系數(shù)大小。

圖6 水滴局部收集系數(shù)Fig.6 Local collection efficiency of droplets

冰晶的撞擊特性計算中，由于本文對比文獻中的試驗溫度較低，故不考慮冰晶運動過程中的融化。由于文獻試驗中采用的冰晶粒子接近球形，故本文針對冰晶的撞擊特性計算方法和水滴撞擊特性計算的方法類似，冰晶運動中的阻力系數(shù)采用了Clift和Gauvin的球形粒子阻力系數(shù)模型[27]。冰晶粒子的局部撞擊系數(shù)計算式為

(33)

算例2和算例3中冰晶的局部收集系數(shù)結(jié)果如圖7所示，冰晶的直徑越大，冰晶撞擊翼型的范圍越廣。當量直徑為200 μm和150 μm的冰晶收集系數(shù)峰值相差不大。算例3和算例4的飛行條件以及冰晶當量直徑一樣，其冰晶收集系數(shù)結(jié)果也一致。

圖7 冰晶局部收集系數(shù)Fig.7 Local collection efficiency of ice crystals

4.3 黏附效應及結(jié)冰冰形計算

算例1不含冰晶，因此無冰晶黏附。算例2～算例4冰晶黏附系數(shù)和局部撞擊系數(shù)的計算結(jié)果如圖8所示。從黏附系數(shù)結(jié)果可以看出，僅算例2有冰晶黏附，算例3和算例4無冰晶黏附，說明算例3和算例4中冰晶撞擊壁面后全部反彈，不參與結(jié)冰過程。

無冰晶的算例1表面溫度和溢流水結(jié)果如圖9(a)所示。算例2有冰晶黏附是因為此狀態(tài)下環(huán)境溫度相對較高，從圖9(b)可看出算例2表面收集的水并沒有在撞擊處全部凍結(jié)，駐點附近(-0.03

圖8 冰晶局部黏附系數(shù)Fig.8 Local adhesion efficiency of ice crystals

存在冰晶黏附，黏附系數(shù)小于0.2。算例3和算例4無冰晶黏附的原因是這兩個狀態(tài)下環(huán)境溫度相對較低，從它們的表面溫度和溢流水結(jié)果即圖9(c) 和圖9(d)中看到，撞擊水在撞擊處能夠全部凍結(jié)而無溢流水，使得表面形成干態(tài)冷表面，此時黏附系數(shù)模型中Kd=0，沒有形成冰晶黏附條件。換句話說，算例3和算例4中的冰晶量并沒有對圖10是本文冰形預測結(jié)果與文獻[18]試驗結(jié)果的對比。為了體現(xiàn)黏附模型的作用，本文預測結(jié)果展示了是否考慮黏附模型兩種情況。不考慮冰晶黏附效應時，認為撞擊到機翼表面的冰晶全部黏附并參與結(jié)冰過程，忽略了冰晶可能出現(xiàn)的反彈飛濺。從圖中可看到對于不含冰晶的算例1，計算結(jié)果與試驗結(jié)果相差不大，其余算例中無黏附模型的冰形計算結(jié)果均遠大于試驗結(jié)果?？紤]黏附效應后，計算結(jié)果得到極大改善。這說明冰晶的黏附效應是影響混合相結(jié)冰的重要因素，數(shù)值模擬時必須考慮。

算例1和算例2因為表面存在溢流水，在溢流區(qū)表面溫度為273.15 K(圖9)，結(jié)冰為明冰。從圖10(b)可以看到混合明冰結(jié)冰氣象條件下所結(jié)冰偏向楔形狀而不是槽狀，這和Baumert等[23]混合相結(jié)冰冰形一般是楔形冰的結(jié)論一致。另外，算例1計算結(jié)果較試驗結(jié)果偏小，而算例2偏大。分析原因如下，算例1與算例2的結(jié)冰條件不同之處在于算例2中增加了冰晶含量，從圖10(a) 和圖10(b)中可看到算例2的試驗結(jié)冰量反而較算例1的略有減小，試驗現(xiàn)象說明了混合相條件下除了冰晶黏附效應還存在冰晶侵蝕效應，這說明本文液態(tài)水結(jié)冰量預測值偏小，加上本文計算未考慮冰晶侵蝕效應，因此算例1冰形預測值偏小，而算例2預測值偏大。冰晶侵蝕效應機理還不明確，目前尚未有合適的模型。

圖9 表面溫度和溢流水結(jié)果Fig.9 Results of surface temperature and runback water

圖10 冰形計算結(jié)果對比驗證Fig.10 Comparison and validation of computational results of ice shape

算例3和算例4表面無溢流水，表面溫度低于273.15 K(圖9)，結(jié)冰為霜冰。從圖10(c)和圖10(d)可看出霜冰條件下的仿真結(jié)果與試驗結(jié)果吻合較好。

5 結(jié) 論

本文針對混合相結(jié)冰問題，建立了混合相結(jié)冰熱力學模型和冰晶黏附模型，利用FLUENT軟件及其用戶自定義函數(shù)開展了混合相結(jié)冰模擬，和試驗結(jié)果對比說明了方法的有效性。主要結(jié)論如下：

1) 基于Messinger模型，增加冰晶有關(guān)質(zhì)量和能量相，建立了適用于混合相結(jié)冰的結(jié)冰熱力學模型。

2) 分別計算了NACA0012翼型在霜冰條件和明冰條件下的結(jié)冰情況，仿真結(jié)果與試驗結(jié)果具有較好的一致性，說明了本文計算模型和方法對混合相結(jié)冰的適用性。混合明冰結(jié)冰氣象條件下所結(jié)冰偏向楔形狀而不是槽狀。

3) 黏附效應是影響混合相結(jié)冰的關(guān)鍵因素，添加黏附效率經(jīng)驗公式能有效改善仿真計算的準確性。本文的黏附模型是經(jīng)驗公式，還有待從機理上進一步開展黏附模型的試驗研究，以獲得更精確的黏附模型。

4) 非球形冰晶在運動過程中的傳熱傳質(zhì)模型有待于增加到現(xiàn)有歐拉模型中，從而擴展對熱環(huán)境下的非球形冰晶運動及黏附結(jié)冰的計算。

5) 冰晶的侵蝕效應是區(qū)別于過冷水滴結(jié)冰的重要方面，其理論模型有待深入研究。