混合相態(tài)冰晶積冰的數(shù)值研究

郭琪磊，?？〗?，安 博，桑為民，*，周 峰

(1. 西北工業(yè)大學 航空學院，西安 710072；2. 中國民用航空飛行學院 航空工程學院，廣漢 618307；3. 中國商飛上海飛機設計研究院，上海 201210)

0 引 言

冰晶積冰是飛行安全的重要威脅。傳統(tǒng)觀點認為發(fā)動機結冰是由于過冷水滴撞擊到發(fā)動機進氣道前緣、整流罩、支板以及導流葉片表面進而引起外部結冰現(xiàn)象[1]。然而自從上世紀90年代以來，在超過海拔7 000 m且不存在過冷水滴的高空發(fā)生了多起大涵道比渦扇發(fā)動機功率損失事件，研究表明事故原因是由發(fā)動機所吸入冰晶在內涵道部件表面發(fā)生部分融化，進而結冰所導致的[2-3]。

當飛機為躲避低空中發(fā)生降雨的區(qū)域而在其上空巡航時，發(fā)動機會吸入大量冰晶。冰晶在發(fā)動機高溫內流的作用下，會部分融化而形成冰水混合相。冰水混合相粒子會進一步在壓氣機葉片等發(fā)動機內部結構表面形成液態(tài)水膜，后續(xù)吸入的冰晶依附在這些潤濕表面并將表面溫度冷卻至冰點以下。伴隨著冰晶的持續(xù)吸入，最終會在內部結構表面累積成冰。與過冷水滴撞擊所引發(fā)的發(fā)動機外部結冰不同，在發(fā)動機整個低壓壓氣機至高壓壓氣機一級靜子葉片區(qū)域均可能發(fā)生冰晶結冰現(xiàn)象[2-4]。冰晶積冰會導致發(fā)動機喘振、熄火、功率損失，甚至會由于冰脫落而造成內部結構損傷[3]。

國內外學者分別采用地面試驗和數(shù)值模擬方法對冰晶結冰機理進行深入研究。在試驗研究方面，Struck[5]、Currie[6]等認為可以用濕球溫度Twb作為結冰是否穩(wěn)固的判斷依據(jù)，即當Twb＜ 0℃時，潤濕表面可以形成穩(wěn)固的積冰。Currie[7]用黏附效率描述積冰過程，并探討了融化率（Melt Ratio, MR）和撞擊角度對黏附效率的影響。Al-Khalil等[8]分別對霜冰和明冰條件下的冰晶結冰進行了實驗，并發(fā)現(xiàn)冰晶撞擊會對冰晶結冰起到削弱作用。Knezevici等[9]則進一步發(fā)現(xiàn)由于冰晶侵蝕效應，大尺寸冰晶較小尺寸冰晶結冰量有所降低。Hauk等[10-11]研究了球形及不規(guī)則外形冰晶的撞擊特性和影響因素。

然而，模擬航空發(fā)動機冰晶結冰試驗具有成本高、周期長、結果不具備普適性等特點。因此數(shù)值模擬方法也成為了研究冰晶結冰的重要手段。Villedieu等[12]采用Lagrange方法對冰晶形狀、傳熱傳質、相變以及冰晶撞擊過程中的黏附、破碎、反彈、飛濺等現(xiàn)象建立數(shù)學模型。Trontin等[13]在此基礎上，考慮冰晶侵蝕效應影響，對冰晶黏附效率模型加以改進，并根據(jù)試驗數(shù)據(jù)校正數(shù)值模型中經驗參數(shù)；隨后對過冷大液滴、冰晶及混合相等條件下結冰情況進行了數(shù)值仿真[14]。Norde等[15]采用歐拉法計算粒子軌跡，建立了冰晶撞擊和侵蝕模型，并針對混合相傳熱傳質過程的特點，改良了經典Messinger結冰熱力學模型。Currie[16]、Bartkus[17]、Baumert[18]等分別根據(jù)冰風洞試驗數(shù)據(jù)的校準并完善了其冰晶結冰數(shù)值計算軟件。

而國內對冰晶結冰研究仍處于起步階段。Zhang等[19]綜合考慮液滴的飛濺以及冰晶的破碎和反彈等現(xiàn)象，建立了冰晶撞擊模型，并在此基礎上提出了混合相態(tài)結冰熱力學模型。姜飛飛等[20]對冰晶的傳熱傳質過程進行離散化處理，計算分析了冰晶的粒子半徑、冰晶溫度、冰晶速度等參數(shù)變化，獲得了冰晶在低壓壓氣機內涵通道內的運動軌跡及與碰撞特性。卜雪琴等[21]考慮了冰晶黏附效應，分別對霜冰與明冰條件下二維NACA0012翼型表面結冰情況進行了數(shù)值研究，結果表明冰晶黏附效應對混合相下結冰量及冰形均有較大影響。

本文分析了混合相態(tài)下結冰表面的傳熱傳質過程，通過改進的Messinger結冰熱力學模型，建立了混合相態(tài)冰晶積冰熱力學數(shù)值模型。通過與Cox冰風洞NACA0012翼型冰晶積冰實驗對比，驗證了本文數(shù)值模型的正確性，并在此基礎上分析了融化率對冰晶積冰過程的影響機制，以及環(huán)境溫度、馬赫數(shù)等參數(shù)對積冰形態(tài)和積冰生長率的影響。

1 混合相態(tài)下積冰模型

1.1 混合相態(tài)下結冰熱力學模型

1.1.1 質量守恒

混合相態(tài)下結冰熱力學行為需要對經典Messinger模型進行擴展，忽略冰層內部、冰層與壁面之間的熱傳導，基于控制體積法建立如下所示的質量和能量平衡方程[21]。

控制體內液態(tài)水質量傳遞如圖1(a)所示，流入控制體的質量流量有：撞擊到翼面的液滴和冰晶已融化部分的質量流量，流入控制體的溢流水質量流量。而流出控制體的質量流量為：凍結成冰的質量流量（為負時表示冰融化為水進入控制體），蒸發(fā)的質量流量，以及流出控制體的溢流水質量流量。因此可建立控制體內質量平衡方程：

控制體內結冰的質量流量應包含以下三部分：控制體內凍結成冰的質量流量，撞擊并依附在壁面上的冰晶尚未融化的質量流量，及由于升華現(xiàn)象而損失的質量流量，即

圖1 控制體內(a)質量守恒；(b)能量守恒Fig. 1 (a) Mass conservation and (b) Energy conservation in control volume

上述式（1-2）中，

式中，LWC、IWC分別為流場中過冷液滴的液態(tài)水含量與冰晶含量，V∞,d與V∞,ic分別為自由來流的水滴速度和冰晶速度，βd、βic分別為過冷水滴和冰晶的收集系數(shù)，Sw為控制體的底面積，ηic為冰晶融化比，ηic=LWCic/IWC，LWCic為部分融化冰晶中液態(tài)水含量。

式（1）中當前控制體流入的溢流水質量流量為上一個控制體流出的溢流水質量流量，即上一個控制體,out；而當前控制體流出的溢流水質量流量為：

式中，ρw為水的密度，hf為控制體內液膜高度，Uf為控制體內液態(tài)水溢流速度，Zw為沿翼型展向控制體寬度。

最后，質量平衡方程中蒸發(fā)的質量流量m˙ev為 ：

式中，hc為對流換熱系數(shù)，cp,air為空氣比熱容，Pvap,w為壁面處飽和蒸氣壓，Pvap,∞為環(huán)境空氣飽和蒸氣壓，Ps,w為控制體外靜壓，HR 為相對濕度。

1.1.2 能量守恒

控制體內能量傳遞機理如圖1(b)所示?？刂企w內熱量損失項有：因對流損失的熱量conv，蒸發(fā)散熱ev，壁面收集液滴顯熱c,d，壁面收集的冰晶顯熱c,ic（由冰晶未融化部分的固態(tài)顯熱c,ic,i和 已融化部分的液態(tài)顯熱c,ic,w兩部分組成），及流出控制體溢流水的顯熱out。而控制體內熱量增加項有：壁面收集液滴動能ke,d，壁面收集冰晶動能ke,ic（同樣由冰晶未融化部分的固態(tài)動能ke,ic,i和已融化部分的液態(tài)動能ke,ic,w組成），結冰熱f，水膜與冰層間溫差引起的顯熱i，以及流入控制體溢流水的顯熱in。因此可建立控制體內能量平衡方程：

上式中進入控制體的質量項動能分別為：

式中，Uimp,d與Uimp,ic分別為水滴和冰晶粒子撞擊到結冰表面時的撞擊速度。

控制體內存在液態(tài)水的凍結與蒸發(fā)，以及冰的升華三種相變過程，式（8）中相變潛熱項為：

式中，Lf為結冰潛熱，Lev為蒸發(fā)潛熱，Lsub為 升華潛熱。

能量守恒方程(8)中對流換熱項為：

式中，Ts為壁面溫度，Tinf為自由來流溫度。

控制體內顯熱傳遞，主要有水滴和冰晶粒子溫度變化引起的顯熱傳遞和流入/流出控制體的水引起的顯熱傳遞。能量守恒方程中所有溫差引起的顯熱傳遞可以表示為：

式中，Tm為融化溫度，cp,w和cp,ic分別為水和冰的比熱容。

1.2 結冰熱力學模型求解

對于二維翼型，駐點所在控制體流入的溢流水質量流量為零，且此控制體內的液態(tài)水等分為兩部分，分別沿翼型表面向上下游溢流。值得注意的是，前一個控制體的溢流水流出質量流量等于下一個控制體的流入質量流量，即。

首先假定壁面溫度等于融化溫度，即Ts=Tm。根據(jù)能量守恒方程可得凍結成冰的質量流量f，進而根據(jù)f判斷結冰狀態(tài)[15]：

1.3 冰晶黏附模型

通常情況下認為過冷水滴在撞擊結冰表面后會全部發(fā)生黏附參與結冰過程。而冰晶與過冷水滴存在很大不同，冰晶撞擊表面后可能發(fā)生破碎、反彈和黏附等結果。Nilamdeen等[22]基于Euler方法定義了黏附系數(shù)，并指出冰晶撞擊動力學行為受冰晶粒徑尺寸、撞擊速度以及液膜厚度等參數(shù)影響。文獻[22]假定在霜冰區(qū)域冰晶全部反彈，在液膜區(qū)域冰晶全部黏附，即黏附系數(shù)分別為εst=0和εst=1；而在明冰區(qū)域內，冰晶黏附系數(shù)εst則與冰晶撞擊速度和液膜厚度等參數(shù)有關，遵循以下關系式：

式中，hf為液膜高度，hf,max為計算最大液膜厚度，vn為冰晶的法向速度分量。參數(shù)χ用于控制撞擊速度閾值vc，即當所有冰晶粒子法向速度分量大于速度閾值時（vn＞vc），冰晶全部發(fā)生反彈而不發(fā)生黏附，即黏附系數(shù)εst=0。假定hf=hf,max、εst=Ψ，χ可由式（25）推得，具體表達式如下所示。

其中，Ψ為極小正值，vc為定義反彈邊界的撞擊速度閾值，，dp為冰晶直徑。

2 數(shù)值結果與分析

2.1 數(shù)值模型驗證

本文基于FENSAP-ICE結冰數(shù)值計算軟件，湍流模型采用Spalart-Allmaras一方程模型，采用歐拉法計算液滴和冰晶撞擊特性。冰晶粘附模型采用NTI Bouncing Model[22-23]，通過改進的Messinger結冰熱力學模型，采用單步法對積冰增長及結冰表面水膜流動進行求解。

根據(jù)文獻[8]中Cox冰風洞實驗結果進行數(shù)值模型驗證工作。該實驗以NACA0012翼型為研究對象，弦長為0.914 4 m，攻角為0°。液滴當量直徑為20 μm，冰晶當量直徑為150或200 μm，結冰時間均為600 s，詳細結冰條件參見表1。

網格劃分采用結構化網格，翼型附近網格劃分如圖2所示，全局網格數(shù)量為321 600，展向分布5層網格，且近壁處網格劃分滿足y+≈1原則，由文獻[15]可知，本文所采用的結構化網格滿足收斂性要求。計算域高度與文獻[8]中Cox冰風洞試驗段尺寸保持一致。

表1 混合相結冰條件Table 1 Mixed-phase icing conditions

圖2 數(shù)值計算網格Fig. 2 Mesh for numerical simulation

如圖3所示，取本文數(shù)值模型計算的冰形與文獻[8]中實驗結果進行對比。可以看出，無論是工況run19中霜冰結冰條件，還是工況run 10中明冰結冰條件，本文數(shù)值結果與實驗所得冰形均較為一致，進而驗證了本文數(shù)值方法的正確性?；旌舷鄳B(tài)積冰問題中冰晶侵蝕作用主要受冰晶粒徑和撞擊速度以及黏附效率影響響[23]，圖3（b）中在考慮冰晶侵蝕作用后，雖然結冰范圍較實驗結果略大，卻明顯抑制了翼型前緣處冰形生長，可以更好地吻合試驗結果。

圖3 工況(a) run 19和(b) run 10條件下冰形對比Fig. 3 Comparison of ice shape: (a) run19; (b) run 10

2.2 參數(shù)分析

2.2.1 融化率

為研究融化率(MR = LWC/TWC)對冰晶結冰過程的影響，在保證總水含量TWC = 1.4 g/m3前提下，本文分 別選取IWC/LWC = 0.1/1.3、0.4/1.0、0.7/0.7、1.0/0.4、1.3/0.1共五組工況進行對比研究，其余結冰參數(shù)與表1中run 10所示參數(shù)保持一致。

圖4為不同融化率下冰形對比?？芍斠簯B(tài)水含量LWC占主導時，結冰范圍較大，且結冰范圍隨LWC下降而不斷縮小。這是由于此時壁面所收集到的液態(tài)水在氣動力作用下，從前緣駐點沿翼型上下表面向下游溢流，最終全部凍結；而冰晶粒徑尺度相對較大，其運動軌跡不易受到氣動力影響，更加集中地收集于翼型前緣。當冰水含量IWC占主導時，結冰范圍基本保持不變，但結冰厚度顯著降低。

圖4 不同融化率下冰形對比Fig. 4 Comparison of ice shape at different melt ratio

圖5為不同融化率下翼型前緣液膜厚度。當IWC/LWC = 1.3/0.1時，由于液態(tài)水含量過小，液滴撞擊到翼型表面后會在低溫作用下即刻結冰，故此時水膜厚度為0 μm。隨著LWC不斷增大，液膜的厚度和潤濕范圍均隨之增大，可以黏附更多的冰晶在結冰表面換熱積冰。當IWC/LWC = 0.1/1.3時，液膜厚度和潤濕范圍均達到最大值，其中液膜厚度峰值約為10.7 μm，此時潤濕范圍對應圖5中的積冰極限。此外由圖4、圖5可知，隨著LWC不斷降低，雖然冰晶的侵蝕效應隨之減弱[14]，但與此同時結冰表面液膜厚度變小，不足以黏附更多的冰晶，導致總體結冰量逐漸減小。

圖5 不同融化率下液膜厚度Fig. 5 Comparison of film thickness at different melt ratios

圖6為不同融化率下翼型前緣駐點處結冰厚度。為確定翼型前緣駐點處達到最大結冰厚度時的融化率，增加了IWC/LWC = 0.5/0.9、0.6/0.8兩組工況。由圖6可知，在run 10結冰條件下，當IWC/LWC =0.5/0.9時前緣駐點結冰厚度達到最大值，約為9.3 mm。綜上所述，混合相冰晶積冰若達到最大結冰厚度，不僅要有足夠的冰晶含量，以保證冰晶經過相變換熱使結冰表面溫度降至冰點以下，也需要有足夠的液態(tài)水含量，足以黏附冰晶在結冰表面進行換熱。

圖6 不同融化率下前緣駐點處積冰厚度Fig. 6 Ice thickness at the stagnation point at different melt ratio

2.2.2 環(huán)境溫度

為探究環(huán)境溫度對于冰晶結冰過程影響，本文分別選取環(huán)境溫度T為-4℃、-7℃、-10℃。由表2可知，在保證相對濕度（Relative Humidity，RH）一定情況下，不同溫度下翼型表面粒子收集質量差異較小，隨溫度逐漸升高，粒子收集質量略有增加。

環(huán)境溫度的變化也直接影響了濕球溫度Twb變化，三種情況下濕球溫度Twb與環(huán)境溫度近似，同樣為-4℃、-7℃、-10℃。濕球溫度變化會影響粒子的相變傳熱過程，進而影響結冰過程。如圖7所示，翼型表面所收集的液態(tài)水并未在撞擊處完全凍結，而是在駐點附近（約-0.03＜Y＜ 0.03）形成溢流水且壁面溫度約為0℃，從而形成混合態(tài)積冰條件。同時可以發(fā)現(xiàn)，隨著環(huán)境溫度升高，液膜的厚度和潤濕范圍均隨之增大。

表2 不同環(huán)境溫度下濕球溫度及粒子收集質量Table 2 Wet-bulb temperature and particle collection mass at different ambient temperatures

圖7 不同溫度下下壁面溫度和液膜厚度對比Fig. 7 Comparison of wall temperature and film thickness

圖8為不同溫度下翼型前緣駐點積冰生長量對比。可以看出隨著溫度降低，翼型前緣駐點處結冰量和積冰速率均有明顯增加。此外，可知在結冰初期積冰速率均較快，相對應的增長量曲線較為陡峭；而隨著積冰的生長，積冰速度逐漸穩(wěn)定，積冰增長量曲線也趨于平緩。

圖8 不同溫度下駐點處積冰增長量對比Fig. 8 Comparison of tip growth at different temperature

2.2.3 馬赫數(shù)

本文取馬赫數(shù)Ma= 0.16、0.24、0.48，以研究其對混合相態(tài)結冰過程影響，相對應的液滴與冰晶收集系數(shù)和駐點積冰生長率分別如圖9（a）、（b）所示。由圖9（a）可知，隨著馬赫數(shù)增大，翼型表面液滴收集系數(shù)與潤濕極限均隨之增大，其中前緣駐點處液滴收集系數(shù)從0.575增加至0.711；而冰晶收集系數(shù)雖略有增加但變化并不顯著，前緣駐點處冰晶收集系數(shù)僅從0.391增加至0.415。發(fā)生該現(xiàn)象的原因是液滴質量相較于冰晶更小，因而在高氣動力作用下更容易改變其運動軌跡。

圖9 不同馬赫數(shù)下收集率和駐點處積冰增長量對比Fig. 9 Comparison of collection efficiency and tip growth at different Mach number

而從圖9（b）中可得知，當馬赫數(shù)較?。ㄈ鏜a=0.16）時，前緣駐點處積冰生長量隨Ma數(shù)增大而顯著增大；當馬赫數(shù)較大（如Ma= 0.32）時，繼續(xù)增大馬赫數(shù)，前緣駐點積冰生長量并沒有明顯增長。此外，由圖9（b）可知，在經歷初期較快結冰速率后，隨著積冰生長，結冰表面的積冰速度同樣表現(xiàn)出逐漸平緩的趨勢。

3 結 論

本文分析了混合相態(tài)下結冰表面的傳熱傳質過程，通過改進的Messinger結冰熱力學模型，建立了混合相態(tài)冰晶積冰熱力學數(shù)值模型。通過與文獻中Cox冰風洞下NACA0012翼型積冰實驗結果對比，驗證了本文數(shù)值模型的正確性，并在此基礎上分析和研究了融化率對冰晶積冰過程的影響機制，以及環(huán)境溫度、馬赫數(shù)等參數(shù)對積冰形態(tài)和積冰生長率的影響。具體結論如下：

1）混合相冰晶積冰若達到最大結冰厚度，不僅要有足夠的冰晶含量，以保證冰晶經過相變換熱使結冰表面溫度降到冰點以下，也需要有足夠的液態(tài)水含量，足以黏附冰晶在結冰表面進行換熱。因此針對不同結冰環(huán)境，存在最嚴酷結冰融化率；

2）環(huán)境溫度直接影響了粒子濕球溫度變化，進而影響粒子傳熱傳質和相變過程。隨環(huán)境溫度升高，液膜的厚度和潤濕范圍也隨之增大；

3）當環(huán)境溫度降低或馬赫數(shù)增大，翼型前緣駐點處結冰量和積冰速率均有明顯增加；隨著積冰的生長，積冰速度逐漸穩(wěn)定，積冰增長量曲線也趨于平緩。