基于歐拉方法的2 維翼型冰晶結冰數值計算

譚 燕

（中國民用航空飛行學院航空發動機維修培訓中心，四川廣漢618307）

符號表

0 引言

自20 世紀90 年代以來，全球已經發生超過100起由于冰晶造成的航空發動機功率損失事件[1-4]。航空發動機冰晶結冰問題已經受到各國重點關注，其中美國國家航空航天局（NASA）和加拿大國家研究委員會（NRC）在該方面的研究走在世界最前沿[5-8]。由于試驗方法的局限性，學者開始利用數值手段進行機理研究。冰晶結冰數值模擬基本與過冷液滴結冰模擬一致，大致可以分為流場計算、粒子軌跡計算和冰形生長計算。其中，粒子軌跡計算分為拉格朗日方法[9-10]和歐拉方法[11-12]。Villedieu 等[9]開發了ONERA 2D 工具包，考慮了冰晶球度對軌跡的影響，以及冰晶反彈等因素的影響，并用NASA-NRC 的試驗結果與計算結果進行了對比；Zhang 等[10]建立了冰晶反彈破碎模型和混合相模型，通過用戶定義函數（User Defined Function，UDF） 方式利用 Fluent 軟件模擬了NACA0012 冰晶結冰過程；而Norde 等[11]和Nilamdeen等[12]則采用計算量更少的歐拉方法計算了NACA0012冰晶結冰過程，并均與試驗進行了對比?？傮w而言，冰晶結冰研究成果相對于過冷液滴結冰的少很多[1]。

本文采用歐拉方法對某對稱楔形翼型結冰過程進行數值計算，數學模型中考慮了相變、粒子反彈等問題；采用NASA-NRC 試驗結果驗證了本文方法的可行性，并分析了壓力、CLWC/CTWC等因素對結冰的影響。

1 模型

由于本文多相流計算中微粒濃度PL＜0.1（Particulate Loading，PL=αdρd/αcρc），可以認為載體相對分散相的影響是單向的。因此，本文通過基于Spalart-Allmaras 湍流模型[13]計算出空氣流場，將計算結果作為已知條件用于冰晶控制方程的計算，利用Messinger 模型進行冰形生長計算。

1.1 粒子運動和相變控制方程

由于混合相是冰晶結冰的必要條件[1，5，9]，體積分數α 由固相αi和液相αw組成，故動量守恒方程由原來1 個增加至2 個

雖然基于歐拉方法的冰晶結冰計算步驟與過冷液滴結冰計算步驟相同，但冰晶在運動過程中存在質量、能量、相態的變化，因此，原有的質量和能量控制方程[14]不再適合，必須在守恒方程中考慮傳質和傳熱問題。

當T＜Tm時αw=0，冰晶完全為固相

質量守恒方程（式（3））等號右邊為升華造成的冰晶質量損失，而能量方程（式（4））等號右邊分別為對流和升華造成的能量轉化。

當T=Tm時α=αi+αw，冰晶粒子內部為固態，而外部為液態，混合相粒子通過換熱獲得的熱量Q˙conv用于表面液態蒸發和冰晶融化

因此，通過混合相粒子能量守恒（式（5））可以推導出冰晶融化率[11]

冰晶固態相質量損失主要由融化造成

冰晶液態相損失，除了新增融化質量外（等號右邊第1 項），還有由于表面蒸發造成的質量損失（等號右邊第2 項）

能量方程為

當T＞Tm時αi=0，冰晶粒子完全融化成水，質量αi=0 損失主要由蒸發造成

而對流換熱和蒸發引起的能量變化為

上述質量、動量和能量方程中需要求解未知變量為α、u 和T，但無法確定粒子直徑dp變化。由于粒子數量密度是不變的，因此可以通過粒子數量密度守恒來反映dp變化

上述控制方程可通過流線迎風伽遼金有限元方法進行求解[15]。

1.2 冰形生長模型

本文冰形生長模型采用經典Messinger 模型，該模型以控制體積建立質量和能量守恒方程[11]。進入控制體積的質量主要是冰晶融化部分、過冷液滴前一個控制體流入質量，而離開控制體積的質量主要是液膜蒸發，結冰、下一個控制體積流出質量（如圖1 所示）。因此控制體積的質量方程為

冰層質量守恒方程為

圖1 Messinger 模型

控制體積能量守恒方程為

當進入控制體積的液態水質量小于結冰質量m˙f時，即，說明冰晶和過冷液滴在表面立刻結冰，形成霜冰。此時，在控制體內不存在液膜，液態水不會進入下一個控制體積（），同時也沒有由于蒸發造成的質量（），冰層升華造成的能量損失將取代式（16）的蒸發能量損失于是冰層質量守恒方程和能量方程為

2 算例分析

本文以NRC試驗[5，16]中的對稱楔形翼型作為研究對象，該翼型前部夾角為40°，后部夾角為20°，弦長220.92 mm（如圖2（a）所示）。該試驗設備可產生液滴直徑=40 μm，冰晶直徑=100～200 μm，在本文計算中冰晶尺寸取其中間值，即150 μm。整個流場計算模型（154064 個六面體）如圖2（b）所示。進口距翼型前緣約1.27 m，由于該距離過短，同時受試驗條件的限制，無法完全模擬來流在壓氣機中逐漸增溫增壓過程，純冰晶在該距離短時間內無法形成混合相。為了模擬壓氣機工作環境下所形成的混合相，NRC 在試驗過程中補充了40 μm 液滴。

圖2 對稱楔形翼型實物[16]與網格模型

2.1 對比驗證

為了驗證本文數值方法的準確性，首先以NASA-NRC 風洞試驗[16]中的第139 號試驗工況（見表1）作對比。

表1 NASA-NRC 第139 號試驗工況[16]

本文采用上述數值方法計算后，得到了豐富的流場數據（如圖3 所示），通過對比可見，數值結果與試驗結果基本相同，但二者存在一定偏差，這是由于本文未考慮風洞壁面效應，出口采用壓力遠場所致。

圖3 流場結果

粒子軌跡計算結果如圖4 所示。從圖中可見，盡管計算工況中總溫為12.5 ℃，但濕球溫度為負值，蒸發過程比較強烈，將吸收大量熱量，導致液滴溫度逐漸降低，駐點處液滴溫度由最初的10 ℃降低至-8℃左右，最終形成結冰條件。而在蒸發過程中，撞擊駐點處的液滴在較為強烈的蒸發作用下，其直徑由初始的40 μm 變為39.3 μm。在該工況下暴露402 s 后，計算結果與試驗結果較為一致，但由于采用單步計算，二者尚存在一定差異，如圖5 所示。

圖4 軌跡計算結果

圖5 NASA-NRC 的第139 號試驗結果對比

2.2 壓力對結冰影響分析

為了研究總壓對結冰的影響，進行如下工況的計算（表1 中Model A 工況實際為NASA-NRC 試驗的第543 號 試驗工 況[16]，由于公開文獻僅提供定性結論，沒有太多定量結論，本文沒有進行試驗數據對比），所有模型計算攻角為0°，冰晶尺寸為150 μm，液滴尺寸為40 μm，結冰時間均為60 s。

翼型表面流場中粒子（液滴和冰晶）直徑對比如圖6 所示。從中可見，Model A 中翼型表面的液滴和冰晶的直徑明顯要小于其他2 個模型的。在相同工況下，總壓降低會導致蒸發速率增大和濕球溫度進一步降低，因此，受蒸發速率影響，在低壓環境下（Model A）無論冰晶還是液滴，其質量損失最大。

表2 計算工況I kPa

圖6 翼型表面流場中粒子直徑對比

翼型表面粒子收集速率和冰層厚度結果對比如圖7 所示。從圖中可見，液滴/冰晶最有可能撞擊翼型前緣（z=±2 mm）位置，在前緣位置3 種工況粒子收集情況相差不大（圖7（a））。但在翼型迎風面其他部位，在低壓環境下（Model A），粒子收集速率高于高壓環境（Model B 和Model C）下的，這是因為粒子尺寸越大，撞擊表面后發生反彈或破碎的概率越高，造成表面收集質量速率下降。而從冰層厚度對比結果中（圖7（b））可見壓力造成的結冰情況差異較大，在低壓環境下更容易結冰（Model A 中駐點冰層厚度為4.6 mm，而試驗數據為4.5 mm，數值結果與試驗結果基本一致）。一方面，低壓環境（Model A）表面質量收集速率要高（圖7（a））；另一方面，低壓環境濕球溫度更低（Model A/B/C 的濕球溫度分別為-3、-1 和1℃），更利于冰層生長。

圖7 翼型表面質量收集速率和冰層厚度

綜上所述，壓力越低，濕球溫度越低，蒸發作用越強烈，越容易形成結冰條件，這也解釋了為什么發動機內部結冰通常發生在高空環境。

2.3 CLWC/CTWC比例對結冰影響分析

本文還進行工況Ⅱ（見表3）的計算。濕球溫度Twb=-3℃，攻角為0°，冰晶尺寸DMM=150 μm，液滴尺寸DMV=40 μm，總含 水 量CTW=4 g/m3（CTW=CIW+ CLW），結冰時間均為120 s。

翼型表面結冰對比如圖8 所示，翼型表面液膜厚度如圖9 所示，翼型表面液膜速度如圖10 所示。從圖中可見，Model D 不含液態水（CLW=0 g/m3），冰晶也未形成冰水混合相，干燥翼型表面將無法黏附冰晶形成冰層；Model E 中液態水含量較少，僅僅在迎風面形成液膜（圖9、10），冰層也僅僅出現在迎風面；隨著液滴水含量增加，翼型表面的液膜開始流向背風面，并在背風面形成冰層（Model F、Model G 和Model H）。此外，從駐點處冰層厚度對比可見（圖8），當液滴水占據總含水量一半時（CLW/CTW=0.5）冰層最厚；液態水較少時，液膜厚度較薄，吸附冰晶能力較弱；而冰晶較少時，且液態水較多時，雖然液膜較厚，同樣不利于結冰（圖9）。

表3 計算工況Ⅱ g/m3

圖8 翼型表面結冰

圖9 翼型表面液膜厚度

圖10 翼型表面液膜流動速度

3 結束語

本文采用數值方法研究了對稱楔形翼型冰晶結冰問題。由于NASA-NRC 僅公開極為有限的定量試驗結果，大多為定性試驗結果，因此本文采用部分定量試驗結果對數值方法進行了部分驗證，同時進行了壓力、CLWC/CTWC對結冰影響的機理分析，計算結果與試驗定性結果基本相符。本文方法可為后續發動機結冰機理研究或防冰系統設計提供一定參考。