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

譚 燕

（中國(guó)民用航空飛行學(xué)院航空發(fā)動(dòng)機(jī)維修培訓(xùn)中心，四川廣漢618307）

符號(hào)表

0 引言

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

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

1 模型

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

1.1 粒子運(yùn)動(dòng)和相變控制方程

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

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

當(dāng)T＜Tm時(shí)αw=0，冰晶完全為固相

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

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

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

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

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

能量方程為

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

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

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

上述控制方程可通過(guò)流線迎風(fēng)伽遼金有限元方法進(jìn)行求解[15]。

1.2 冰形生長(zhǎng)模型

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

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

圖1 Messinger 模型

控制體積能量守恒方程為

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

2 算例分析

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

圖2 對(duì)稱(chēng)楔形翼型實(shí)物[16]與網(wǎng)格模型

2.1 對(duì)比驗(yàn)證

為了驗(yàn)證本文數(shù)值方法的準(zhǔn)確性，首先以NASA-NRC 風(fēng)洞試驗(yàn)[16]中的第139 號(hào)試驗(yàn)工況（見(jiàn)表1）作對(duì)比。

表1 NASA-NRC 第139 號(hào)試驗(yàn)工況[16]

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

圖3 流場(chǎng)結(jié)果

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

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

圖5 NASA-NRC 的第139 號(hào)試驗(yàn)結(jié)果對(duì)比

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

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

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

表2 計(jì)算工況I kPa

圖6 翼型表面流場(chǎng)中粒子直徑對(duì)比

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

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

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

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

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

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

表3 計(jì)算工況Ⅱ g/m3

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

圖9 翼型表面液膜厚度

圖10 翼型表面液膜流動(dòng)速度

3 結(jié)束語(yǔ)

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