基于Myers模型的三維結(jié)冰數(shù)值仿真

雷夢龍，常士楠，楊波

北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100083

近年來，大量的飛行事故導(dǎo)致飛機的結(jié)冰研究成為一個熱點。根據(jù)最近的飛行數(shù)據(jù)，9%的飛行事故是由于飛機在嚴(yán)酷天氣下結(jié)冰所致，所有氣象因素帶來的飛行事故中，12%是飛機結(jié)冰導(dǎo)致的，且92%是在飛行中發(fā)生結(jié)冰[1]。

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

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

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

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

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

1 Myers水膜流動模型

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

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

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

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

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

(1)

(2)

(3)

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

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

(4)

(5)

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

(6)

(7)

式中：LWC為液態(tài)水含量；V∞為空氣遠(yuǎn)場的速度；β為水滴收集系數(shù)；me為水膜的蒸發(fā)質(zhì)量;ρi為冰的密度。

2 Myers結(jié)冰增長模型

飛機表面的結(jié)冰厚度和水膜厚度相比于機翼的整體尺寸通常很小，因此，Myers模型在分析機翼表面、冰層、水膜之間的導(dǎo)熱時，將結(jié)冰過程視為準(zhǔn)靜態(tài)過程，忽略s1、s2方向的導(dǎo)熱，只考慮與機翼表面垂直的η方向?qū)幔⑶液雎詸C翼與其表面冰層或水膜之間的接觸熱阻，認(rèn)為機翼表面溫度與冰層或水膜底部溫度相同。

圖 2為機翼表面的冰層結(jié)構(gòu)，從下到上依次為冰層、水膜、空氣。冰層的底部溫度為機翼表面溫度Ts，頂部溫度為水的凝固點Tm，水膜底部溫度為Tm，頂部溫度為Th，空氣溫度為Ta。通常一個控制體內(nèi)的熱流分為向空氣的對流換熱熱流Qc、蒸發(fā)熱流Qe、升華熱流Qs、水滴動能Qk、結(jié)冰釋放的潛熱Ql、水滴顯熱Qd、氣動加熱Qa、水膜流動帶入的熱流Qin、水膜流出帶走的熱流Qout、冰層內(nèi)的導(dǎo)熱熱流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為對流換熱系數(shù)；me為蒸發(fā)質(zhì)量；ms為升華質(zhì)量；mimp為撞擊的水滴質(zhì)量；mice為結(jié)冰質(zhì)量；vd為水滴速度；le、ls、lf分別為蒸發(fā)、升華、結(jié)冰潛熱；cw為水的比熱；Td為水滴溫度；r為恢復(fù)系數(shù)；ue為邊界層速度；ca為空氣的比熱；Tin為流入溫度；Tout為流出溫度；ki為冰的導(dǎo)熱系數(shù);min和mout為流入和流出質(zhì)量。對于霜冰和明冰，熱流項的個數(shù)以及其中一些參數(shù)的取值略有不同。

2.1 霜冰結(jié)冰模型

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

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

(18)

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

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

(19)

將結(jié)冰過程視為準(zhǔn)靜態(tài)過程，忽略s1、s2的導(dǎo)熱和冰層本身的儲熱能力，根據(jù)控制體內(nèi)的熱流平衡，可得霜冰表面的導(dǎo)熱微分方程為

(20)

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

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

(21)

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

(22)

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

(23)

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

(24)

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

2.2 明冰結(jié)冰模型

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

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

Qcond-Qout=0

(25)

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

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

當(dāng)臨界冰厚bf>0且冰厚b

根據(jù)控制體內(nèi)的熱流平衡，水膜與空氣接觸面的導(dǎo)熱微分方程為

(26)

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

mincwTin-moutcwTout

(27)

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

(28)

導(dǎo)熱微分方程式(26)的解為水膜內(nèi)的溫度分布，即

(29)

冰層的溫度分布為

(30)

忽略s1、s2方向的導(dǎo)熱，只考慮η方向的導(dǎo)熱，水膜的結(jié)冰速率微分方程為

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

(31)

本文對Myers模型的改進(jìn)主要是對結(jié)冰類型判斷標(biāo)準(zhǔn)的修改。Myers模型認(rèn)為，水滴撞擊到機翼表面發(fā)生結(jié)冰，若結(jié)冰厚度b≤bf，結(jié)冰類型為霜冰，結(jié)冰厚度增長至b>bf后，結(jié)冰類型轉(zhuǎn)變?yōu)槊鞅?。此推論是假設(shè)結(jié)冰過程為準(zhǔn)靜態(tài)得出的，但是實際結(jié)冰情況下，即使結(jié)冰厚度b≤bf，若上游流入和水滴撞擊的流量Q足夠大，也有可能不會全部凍結(jié)，從而使結(jié)冰類型由霜冰轉(zhuǎn)化為明冰。因此，本文對結(jié)冰類型判斷依據(jù)進(jìn)行修改，結(jié)霜冰時，判斷流量Q是否超過最大結(jié)冰質(zhì)量，若超過最大結(jié)冰質(zhì)量，則水不會全部凍結(jié)，冰的表面會產(chǎn)生水膜，結(jié)冰類型轉(zhuǎn)變?yōu)榈?類明冰。令式(31)中的水膜厚度h=0 m，此時的結(jié)冰量即是當(dāng)前位置的最大結(jié)冰量。

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

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

表1 NACA0012翼型結(jié)冰條件Table 1 Icing condition of NACA0012 airfoil

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

3 數(shù)值計算方法

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

(32)

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

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

(33)

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

(34)

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

4 計算條件及結(jié)果

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

4.1 NACA0012平直翼

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

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

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

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

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

4.2 GLC-305后掠翼

本文的目的是對三維機翼的結(jié)冰進(jìn)行仿真。為此，選取三維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后掠翼的表面網(wǎng)格。本文對比使用的實驗結(jié)果和LEWICE計算結(jié)果來源于文獻(xiàn)[21]。

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

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

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

4.3 發(fā)動機進(jìn)氣道

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

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

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

表3 發(fā)動機進(jìn)氣道的結(jié)冰條件Table 3 Icing condition of engine inlet

5 結(jié) 論

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

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

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

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

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