基于歐拉方法的2 維翼型冰晶結(jié)冰數(shù)值計算

譚 燕

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

符號表

0 引言

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

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

1 模型

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

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

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

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

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

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

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

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

冰晶固態(tài)相質(zhì)量損失主要由融化造成

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

能量方程為

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

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

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

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

1.2 冰形生長模型

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

冰層質(zhì)量守恒方程為

圖1 Messinger 模型

控制體積能量守恒方程為

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

2 算例分析

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

圖2 對稱楔形翼型實物[16]與網(wǎng)格模型

2.1 對比驗證

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

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

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

圖3 流場結(jié)果

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

圖4 軌跡計算結(jié)果

圖5 NASA-NRC 的第139 號試驗結(jié)果對比

2.2 壓力對結(jié)冰影響分析

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

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

表2 計算工況I kPa

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

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

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

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

2.3 CLWC/CTWC比例對結(jié)冰影響分析

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

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

表3 計算工況Ⅱ g/m3

圖8 翼型表面結(jié)冰

圖9 翼型表面液膜厚度

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

3 結(jié)束語

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