基于Myers模型的三維結冰數值仿真

雷夢龍，常士楠，楊波

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

近年來，大量的飛行事故導致飛機的結冰研究成為一個熱點。根據最近的飛行數據，9%的飛行事故是由于飛機在嚴酷天氣下結冰所致，所有氣象因素帶來的飛行事故中，12%是飛機結冰導致的，且92%是在飛行中發生結冰[1]。

飛機上結的冰，就其外形可以分為明冰、霜冰、混合冰3種。霜冰與機翼外形較為接近，對機翼的氣動性能影響較小。明冰是在氣溫較高、液態水含量較大時，撞擊到氣動面上的水滴沒有立即全部結冰，而是產生水膜沿機翼表面流動，形成溢流冰。嚴重情況下可能形成冰角，破壞氣動性能[2]。目前對機翼結冰問題的研究方式主要分為理論、實驗和數值仿真，數值仿真作為一種經濟快速的模擬方法，一直是結冰領域最常用的方式之一。

機翼表面結冰問題的數值仿真過程中，由于復雜的結冰機理及界面處氣-液-固耦合傳熱傳質特性，準確預測明冰表面的水膜流動一直是一個難點[3]。早期對飛機結冰問題的研究中，通常使用Messinger結冰模型[4]求解溢流水和結冰量，如現在仍在使用的LEWICE[5]和ONERA[6]等。Messinger模型引入了凍結系數，在一個控制體內，撞擊到機翼表面的水一部分會凍結成冰，未凍結的水沿氣流方向全部流入位于下游的控制體。Messinger模型相對于實際結冰過程過于簡陋：① 沒有考慮到空氣壓力梯度和剪切力對水膜流動的影響，它假設某個控制體中沒有凍結的水全部流入下游的控制體中，與實際情況不符，尤其是在有回流流場的冰角位置；② 在水膜結冰過程中，沒有考慮空氣、水膜、冰層之間的導熱；③ Messinger模型在提出之初是用于計算二維結冰問題的，雖然被擴展到三維結冰情況，但在三維情況下，Messinger模擬更為困難，每個控制體的上游和下游都不止一個控制體，水流流入和流出的方向和流量難以準確確定。國內應用Messinger模型計算三維結冰的文獻已經有很多，如文獻[7]計算了MS-317后掠翼的結冰，文獻[8]計算了發動機進氣道的結冰，文獻[9]計算了轉子的結冰。

為了克服Messinger模型的缺點，水膜的流動必須考慮在內，現有的商業軟件中，FENSAP-ICE就考慮了水膜流動[10]。此外，Myers提出了適用于任意三維表面的水膜流動模型[11]，并基于此流動模型，提出了Myers結冰模型[12]。Myers模型在計算明冰的結冰增長時，考慮了水膜、冰層之間的導熱，在計算水膜流動時，考慮到了空氣剪切力、空氣壓力梯度和曲面特性[13]。目前對Myers模型的研究還較少，其中，比較深入的有文獻[14-15]。

對機翼結冰的數值仿真通常是單步法，在一次計算空氣和水滴流場的基礎上，認為機翼結冰后的冰形對空氣和水滴流場沒有影響，不再根據新的冰形進行流場計算。單步法需要的計算量更少，且算法更簡單，因此大部分結冰計算都使用單步法計算結冰[16]。Myers模型計算結冰[17]時，需要考慮此處的結冰厚度、水膜厚度等信息，但在冰形改變后重新計算結冰時，機翼上的結冰厚度和水膜厚度會初始為0，因此Myers模型在使用多步法計算時需要復雜的算法。

為了充分考慮水膜流動對結冰的影響，提高結冰計算的準確性，本文采用了Myers結冰模型，并在其基礎上加以改善，采用單步法和多步法計算結冰冰形。為了驗證結冰模型的正確性，本文將NACA0012平直翼和GLC-305后掠翼的計算結果與商業軟件以及實驗結果對比，并提供了三維發動機進氣道的計算結果。本文分析了仿真結果的合理性，并討論了Myers結冰計算模型的優勢和不足之處。

1 Myers水膜流動模型

水膜在曲面上的流動受到空氣壓力梯度、重力分量、空氣剪切力等影響，為了簡化水膜流動模型，需要進行以下假設：

1) 水的黏度、密度在不發生相變時，隨溫度變化很小，在計算水膜流動時將其視為常數。

2) 忽略機翼表面粗糙度對水膜流動的影響。

圖1為任意三維曲面上的貼體坐標系示意圖，貼體坐標系用(s1，s2，η)表示，其中s1、s2與曲面相切，η是曲面在該點的法線，不同控制體的s1、s2和η的方向不同。

曲面可以用第1類和第2類基本量描述，對于不同類型的曲面，曲面基本量的取值不同，實際計算中，可以假設曲面為平板來簡化計算[18]，如ICECREMO軟件就是以平板假設簡化水膜流動模型。Myers等以潤滑理論簡化水膜在平板上流動的Navier-Stokes控制方程，得到[19]

(1)

(2)

(3)

式中：u為s1方向的速度；v為s2方向的速度；μw為水的黏度；p為空氣壓力；ρw為水的密度；ε為水膜高度和長度的比值；Re為雷諾數；gs1、gs2、gη分別為重力在s1、s2和η方向的分量。根據潤滑理論，O(ε,ε2Re)的數值大小可以忽略。

水膜在三維表面上s1和s2方向單位長度的體積流量為水的流動速度在高度上的積分，即

(4)

(5)

式中：h為水膜厚度；b為冰層厚度；As1、As2為空氣剪切力在s1、s2方向的分量。水膜流動的質量守恒方程為

(6)

(7)

式中：LWC為液態水含量；V∞為空氣遠場的速度；β為水滴收集系數；me為水膜的蒸發質量;ρi為冰的密度。

2 Myers結冰增長模型

飛機表面的結冰厚度和水膜厚度相比于機翼的整體尺寸通常很小，因此，Myers模型在分析機翼表面、冰層、水膜之間的導熱時，將結冰過程視為準靜態過程，忽略s1、s2方向的導熱，只考慮與機翼表面垂直的η方向導熱，并且忽略機翼與其表面冰層或水膜之間的接觸熱阻，認為機翼表面溫度與冰層或水膜底部溫度相同。

圖 2為機翼表面的冰層結構，從下到上依次為冰層、水膜、空氣。冰層的底部溫度為機翼表面溫度Ts，頂部溫度為水的凝固點Tm，水膜底部溫度為Tm，頂部溫度為Th，空氣溫度為Ta。通常一個控制體內的熱流分為向空氣的對流換熱熱流Qc、蒸發熱流Qe、升華熱流Qs、水滴動能Qk、結冰釋放的潛熱Ql、水滴顯熱Qd、氣動加熱Qa、水膜流動帶入的熱流Qin、水膜流出帶走的熱流Qout、冰層內的導熱熱流Qcond。

Qc=hcv(Tm-Ta)

(8)

Qe=mele

(9)

Qs=msls

(10)

(11)

Ql=micelf

(12)

Qd=mimpcw(Td-Tm)

(13)

(14)

Qin=mincw(Tin-Tm)

(15)

Qout=moutcw(Tout-Tm)

(16)

(17)

式中：hcv為對流換熱系數；me為蒸發質量；ms為升華質量；mimp為撞擊的水滴質量；mice為結冰質量；vd為水滴速度；le、ls、lf分別為蒸發、升華、結冰潛熱；cw為水的比熱；Td為水滴溫度；r為恢復系數；ue為邊界層速度；ca為空氣的比熱；Tin為流入溫度；Tout為流出溫度；ki為冰的導熱系數;min和mout為流入和流出質量。對于霜冰和明冰，熱流項的個數以及其中一些參數的取值略有不同。

2.1 霜冰結冰模型

結冰開始階段，過冷水滴撞擊到機翼表面時，由于機翼表面溫度低于水的凝固點，結冰類型為霜冰。結霜冰時，機翼表面的熱流平衡方程為

Qk+Qa+Ql-Qc-Qd-Qs-Qcond=0

(18)

在每一個時間步長對式(18)迭代求解可得出霜冰結冰情況下機翼表面的溫度Ts。

霜冰的結冰量為撞擊到機翼表面的水量，此結冰增長微分方程為

(19)

將結冰過程視為準靜態過程，忽略s1、s2的導熱和冰層本身的儲熱能力，根據控制體內的熱流平衡，可得霜冰表面的導熱微分方程為

(20)

式(20)的右側是空氣向機翼表面方向傳熱的熱流密度，將其寫為E0d-E1d·Ts的形式，即

E0d=LWC·V∞·β·[cw(Td-Tm)+ciTm]+

(21)

E1d=LWC·V∞·β·ci+hcv

(22)

解式(20)，可得出冰層內的溫度分布為

(23)

霜冰頂部溫度為Tm時對應的冰厚為臨界厚度bf，超過bf，冰的表面會出現水膜，霜冰轉化為明冰。令式(23)中T=Tm、η=b，可得臨界冰厚bf為

(24)

若bf>0且結冰厚度b≤bf，冰層的表面溫度小于水的凝固點Tm，此時結冰類型為霜冰；當結冰厚度升至b>bf時，冰層的表面溫度等于Tm，結冰類型為明冰；若bf<0，則無論結冰厚度為多少，冰層表面溫度都小于Tm，結冰類型為霜冰。

2.2 明冰結冰模型

明冰結冰時，水膜除向空氣導熱外，還向冰層導熱，機翼表面每個控制體內的熱流與霜冰情況有所不同，其熱流平衡方程為

Qk+Qa+Ql+Qin-Qc-Qd-Qe-

Qcond-Qout=0

(25)

在每一個時間步長內解此方程，可得出明冰結冰情況下的機翼表面溫度Ts。

Myers模型認為，只要冰厚b≤bf，結冰類型就是霜冰。但是實際情況下，若水的流量過大，結冰釋放的潛熱不足以使水膜全部結冰時，即使冰厚b≤bf，也可能產生水膜，繼續向后流動。因此，本文在計算結冰時，將明冰分為兩種類型。

當臨界冰厚bf>0且冰厚b

根據控制體內的熱流平衡，水膜與空氣接觸面的導熱微分方程為

(26)

式(26)右側是空氣向水膜方向傳熱的熱流密度，將式(26)的右側寫為E0f-E1f·Tm的形式，可以得出水膜內η方向的溫度分布，即

mincwTin-moutcwTout

(27)

E1f=LWC·V∞·β·cw+hcv+mincw-moutcw

(28)

導熱微分方程式(26)的解為水膜內的溫度分布，即

(29)

冰層的溫度分布為

(30)

忽略s1、s2方向的導熱，只考慮η方向的導熱，水膜的結冰速率微分方程為

訊問筆錄通常是由監察機關調查人員制作，內容通常是一問一答的記錄雙方的語言，筆錄制作本身可能帶有傾向性，訊問筆錄不能完全真實地反映訊問過程中發生的情況；此外，被調查人的傷痕也并不意味著他一定受到了刑訊逼供，造成傷痕的可能性很多，比如打架斗毆，甚至不排除被調查人以自傷的方式來誣陷調查人員，況且非法取證手段多樣，并不一定會在身體上留下痕跡，若是采取精神摧殘或者惡劣訊問環境的方式，很難證明當時的真實情況。錄音錄像不僅可以反映訊問時雙方的問答，還可以記錄當時所發生的一切情況，可以起到證明是否有非法取證行為的作用。

(31)

本文對Myers模型的改進主要是對結冰類型判斷標準的修改。Myers模型認為，水滴撞擊到機翼表面發生結冰，若結冰厚度b≤bf，結冰類型為霜冰，結冰厚度增長至b>bf后，結冰類型轉變為明冰。此推論是假設結冰過程為準靜態得出的，但是實際結冰情況下，即使結冰厚度b≤bf，若上游流入和水滴撞擊的流量Q足夠大，也有可能不會全部凍結，從而使結冰類型由霜冰轉化為明冰。因此，本文對結冰類型判斷依據進行修改，結霜冰時，判斷流量Q是否超過最大結冰質量，若超過最大結冰質量，則水不會全部凍結，冰的表面會產生水膜，結冰類型轉變為第2類明冰。令式(31)中的水膜厚度h=0 m，此時的結冰量即是當前位置的最大結冰量。

為了比較本文對Myers模型的改進對結冰冰形的影響，本文提供了改進的Myers模型與原Myers模型在兩種不同條件下的計算結果對比，兩個計算結果都使用單步法。

表1給出了NACA0012翼型的環境條件，圖3是Case 1的結冰冰形計算結果，圖 4是Case 2的結冰冰形計算結果。從圖3和圖4可以明顯看出，本文改進的Myers模型與原Myers模型計算結果的區別主要體現在結冰的上下極限附近，這兩個位置的結冰量主要來自上游流入的水膜。由于其結冰厚度沒有超過臨界結冰厚度，原Myers模型將這兩個位置處的冰作為霜冰計算，積聚在這里的水膜會全部凍結，形成冰角，造成此處的結冰量過大；而本文改進的Myers模型會通過計算這兩個位置處的水膜流量大小，將其分類為完全凍結的霜冰和部分凍結的明冰區別處理，因此結冰冰形更符合實際情況。由于采用單步法，結冰冰形會與多步法計算結果產生較大差距，多步法的計算結果將在4.1節給出。

表1 NACA0012翼型結冰條件Table 1 Icing condition of NACA0012 airfoil

圖 5和圖 6分別是Case 1和Case 2的水膜厚度計算結果對比。以機翼駐點為分界線，機翼上部和下部各有3個水膜厚度峰值。在駐點處，水滴大量撞擊，水膜厚度逐漸增加，形成第1個厚度峰值，在空氣壓力和剪切力作用下，水膜流速加快，厚度減小；水膜流量繼續累積，厚度再次增加，形成第2個厚度峰值，水滴撞擊量開始小于凍結量，水膜厚度再次減小；最后，空氣壓力梯度和剪切力減小，使水膜流速減緩，厚度增加，形成第3個厚度峰值。Case 1相比于Case 2，由于氣溫較高，水膜覆蓋的區域更大。原Myers模型對明冰的判斷標準過于簡單，把上下結冰極限處的冰歸類為霜冰，使水膜全部凍結為冰，因此水膜覆蓋范圍比改進的Myers模型小，結冰厚度更大。

3 數值計算方法

數值計算的第1步是確定s1、s2和η的方向，本文中，η的方向為網格的法向，指向機翼外部；s1的方向為網格表面空氣流動方向在網格面上的投影，在此方向上壓力梯度和空氣剪切力很大，因此，計算水膜流動時，s1的方向為水膜流動的主要方向；s2的方向為s1和η方向的垂向。

(32)

式中：t為時間；A為網格的面積；n為網格邊的法向，指向網格外，方程右邊中括號內的各項依次為流入流出項、凍結項、撞擊項。

在結冰初始階段，結冰厚度為0 m，此時的結冰類型為霜冰，霜冰結冰的數值計算式為

(33)

無論是第1類明冰還是第2類明冰，冰層頂部溫度都等于水的凝固點Tm。明冰結冰量的數值計算式為

(34)

在單步法的計算中，需要在從計算開始到計算結束的時間內，以一定的時間步長Δt為間隔，計算結冰量和水膜厚度，最終得出結冰的冰形；在多步法的計算中，需要每隔一定的時間間隔更新冰形及其流場，并重復上述計算步驟。

4 計算條件及結果

本文使用的流場數據均由流體計算軟件Fluent得到。流場數據需要包含機翼表面的坐標信息、空氣壓力、空氣剪切力、水滴撞擊系數、對流換熱系數等數據。然后，用本文改進的Myers模型結冰計算程序讀取數據并求解。

4.1 NACA0012平直翼

本文使用的是弦長為0.533 4 m的NACA0012平直翼，水滴直徑為20 μm，速度分別為60.1、102.8 m/s，空氣液態水含量分別為1.3、0.55 g/m3，空氣溫度分別為266.45、262.04 K，結冰時間分別為480、420 s，如表 1所示。NACA0012翼型對比使用的商業軟件計算結果與實驗結果來自文獻[20]。

兩組計算條件均為明冰結冰類型，不同的是Case 1的氣溫高，液態水含量高，結冰時間長，相比于Case 2會產生更多水膜在氣流作用下流動，形成明顯的溢流冰。在結冰計算中，結冰造成的機翼形狀變化會對流場產生不可忽視的影響，因此在整個過程中，本文使用多步法根據冰形變化重新計算流場。

圖 7為Case 1的仿真結果與實驗結果及ONERA軟件仿真結果的對比。在結冰的冰角處，本文的結冰極限與實驗結果及ONERA軟件仿真結果非常接近。上結冰極限的下游處沒有結冰是因為此處的空氣回流使得水膜在壓力梯度的作用下不會向下游位置流動。因此，上冰角的位置可能是由于流場計算軟件對空氣流場計算結果的差異所致。在下結冰極限位置處，本文的計算結果與實驗結果十分吻合。在結冰的上半部分，本文的結冰量略小于實驗結果，原因可能是Myers模型在計算水膜流量時，把機翼表面作為光滑表面。整體上本文計算的結冰量與實驗結果相近。

圖 8為Case 2的仿真結果與實驗結果及LEWICE軟件計算結果的對比。相比于實驗結果和LEWICE軟件計算結果，本文計算的結冰冰形較為圓潤，原因可能與Case 1相同，是由于沒有考慮粗糙度對水膜流動的影響所致。本文計算的上冰角的位置以及下結冰極限與實驗結果吻合很好。

本文提供了單步法和多步法(5步)的計算結果。多步法的上下結冰極限處的結冰厚度較低，原因可能是因為隨著結冰的進行，冰形對流場的影響逐漸突顯出來，在上下結冰極限處，受到結冰冰形的影響，水滴收集量相比于沒有結冰的機翼要小，因此結冰厚度也減小。在水滴的主要撞擊位置，由于空氣-水滴流場受到結冰冰形的影響，多步法的結冰厚度稍高于單步法的結冰厚度。單步法相對于多步法的優勢不明顯，原因可能是計算網格在每一步的網格重構時會產生誤差，隨著結冰過程的進行，這些誤差累積起來，最終對冰形產生較大影響。

4.2 GLC-305后掠翼

本文的目的是對三維機翼的結冰進行仿真。為此，選取三維GLC-305后掠翼。表 2為GLC-305后掠翼的計算條件，大氣壓力為101 300 Pa，水滴直徑為20 μm。圖 9為GLC-305后掠翼的幾何外形，RLE為可變機翼前緣，MAC是平均氣動弦長，截面 A、截面 B與機翼前緣垂直，截面 C為機翼翼根,圖中1 in=25.4 mm。圖 10為GLC-305后掠翼的表面網格。本文對比使用的實驗結果和LEWICE計算結果來源于文獻[21]。

表2 GLC-305后掠翼的結冰條件Table 2 Icing condition of GLC-305 swept wing

GLC-305后掠翼的結冰計算條件為明冰條件，從計算結果中可明顯觀察到明冰的特征。圖 11為翼根處的結冰厚度云圖，圖 12為翼梢處的結冰厚度云圖。翼梢處的結冰厚度比翼根處更大，原因是翼梢處的水滴收集系數和對流換熱系數比翼根處大，因此結冰情況更嚴重。

本文將計算結果與實驗結果以及LEWICE軟件計算結果對比。圖 13(a)為截面A處的結冰冰形對比，圖 13(b)為截面B處的結冰冰形對比，圖 13(c)為截面C處的結冰冰形對比。在截面A處，本文計算的冰角不明顯，但結冰量與實驗結果大致相同。在截面B處，本文計算的上冰角稍小于實驗結果，但是對于上下結冰極限以及結冰量總體趨勢，本文與實驗結果符合得很好。在截面C處，本文的計算結果與實驗結果十分接近。在截面A和截面C中，LEWICE過度預測了冰角的大小，并錯誤預測了冰角的位置。

4.3 發動機進氣道

發動機進氣道結冰也是威脅飛機飛行安全的重要問題。發動機進氣道在進口部分通常做成翼型形式，其內表面結冰范圍和結冰強度通常比外表面大得多，原因是內表面氣動特性惡化，速度場分布不均勻以及氣流局部分離。

本文計算了發動機進氣道在明冰結冰情況下的結冰冰形，結冰條件如表 3所示，空氣壓力為101 325 Pa，水滴直徑為20 μm。發動機進氣道的外形和表面網格如圖 14所示。

本文選取的結冰條件較為嚴重，速度為150 m/s，溫度為263 K，液態水含量為2 g/m3，屬于明冰結冰條件。圖 15為發動機進氣道表面的水滴收集系數(β)云圖，圖 16為發動機進氣道的結冰厚度云圖，圖 17為發動機進氣道頂部截面的結冰外形，圖 18為底部截面的結冰外形。可以看出，本文中發動機進氣道處的冰形有明顯的溢流冰特征。結合圖17和圖18可以發現，發動機進氣道的結冰主要集中在進氣道的內側，而在外側結冰量就少了很多。原因是水滴主要撞擊在內側，在外側撞擊量少很多，造成進氣道內部的水滴收集量相對外部較大。

表3 發動機進氣道的結冰條件Table 3 Icing condition of engine inlet

5 結 論

本文以Myers模型為基礎，修改了其結冰類型判斷標準，得出了一套考慮三維水膜流動，以及水膜、冰層和空氣間傳熱傳質的結冰計算算法。本文將NACA0012平直翼和GLC-305后掠翼的計算結果與文獻中的實驗及仿真結果進行對比，得出以下結論。

1) Myers模型充分考慮了水膜在空氣壓力、空氣剪切力作用下沿著機翼表面的流動。而且，Myers模型在計算結冰速率時，考慮了已產生的結冰厚度和水膜厚度對結冰速率造成的影響。因此，本文的計算結果相比于Messinger模型計算結果，結冰外形與實驗結果符合良好。

2) 本文計算明冰結冰問題時，對于環境溫度很低、液態水含量小的計算條件，計算得出的冰形與實驗結果吻合很好。但在計算空氣溫度較高、液態水含量較大的結冰工況時，對機翼的上冰角的計算結果與實驗有一定的差距。

3) 本文計算了發動機進氣道的結冰冰形，得出了合理的計算結果。本文使用的發動機進氣道外形復雜，三維特征明顯，表明本文的三維計算程序有良好的通用性。

Myers模型雖然考慮了水膜、冰層內部復雜的傳熱和流動過程，但是為了簡化計算，Myers模型對水膜流動和結冰過程作出了大量的簡化假設，這些假設在簡化計算的同時會使結冰計算結果產生誤差。在未來的工作中，必須將這些被忽略的內容加入結冰計算模型中。