基于歐拉方法的機(jī)翼表面水滴撞擊特性數(shù)值模擬

劉 鑫，邢玉明，杜 佶

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

飛機(jī)結(jié)冰是指飛機(jī)在飛行時(shí)其外表面上水分積聚凍結(jié)成冰的現(xiàn)象。它影響飛機(jī)的氣動(dòng)外形、飛行品質(zhì)及飛行安全，給飛機(jī)飛行帶來(lái)了極大的危害。如：飛機(jī)風(fēng)擋結(jié)冰會(huì)影響駕駛員視野；測(cè)溫、測(cè)壓感頭結(jié)冰會(huì)導(dǎo)致儀表指示失真；機(jī)（尾）翼結(jié)冰將影響氣動(dòng)外形，增大飛機(jī)阻力，減小升力，影響全機(jī)操縱性、穩(wěn)定性；螺旋槳、直升機(jī)旋翼結(jié)冰會(huì)降低升力效率；軸流式壓氣機(jī)進(jìn)氣部件結(jié)冰能使發(fā)動(dòng)機(jī)熄火等。據(jù)統(tǒng)計(jì)，大約有9%的飛行事故是由結(jié)冰造成的。長(zhǎng)期以來(lái)，結(jié)冰機(jī)理和飛機(jī)的防／除冰系統(tǒng)一直是國(guó)內(nèi)外飛機(jī)系統(tǒng)設(shè)計(jì)的重要研究課題［1］。

結(jié)冰源于云層中的過(guò)冷水滴，即低于0℃的液態(tài)水滴，過(guò)冷水滴是極不穩(wěn)定的，受到擾動(dòng)后就可能轉(zhuǎn)化為冰。當(dāng)過(guò)冷水滴撞擊到?jīng)]有采取防冰措施的飛行器表面后，它可能立即結(jié)成明冰或向后溢流再結(jié)成明冰。研究中，可以通過(guò)冰風(fēng)洞試驗(yàn)、飛行試驗(yàn)和數(shù)值計(jì)算等途徑獲取結(jié)冰量或結(jié)冰區(qū)域。但由于試驗(yàn)花費(fèi)昂貴、周期長(zhǎng)且受限于設(shè)備條件，很難模擬結(jié)冰包線中一些特殊狀態(tài)。因此美法英等國(guó)都在積極研發(fā)飛機(jī)結(jié)冰的數(shù)值計(jì)算方法和軟件。

進(jìn)行結(jié)冰計(jì)算的第一步需要計(jì)算水滴在模型表面上的撞擊特性和收集系數(shù)。撞擊特性主要通過(guò)兩種方法來(lái)計(jì)算：一種是Lagrangian方法；另一種是Eulerian方法。多數(shù)結(jié)冰計(jì)算軟件采用Lagrangian方法，在所得流場(chǎng)計(jì)算結(jié)果的基礎(chǔ)上計(jì)算過(guò)冷水滴的運(yùn)動(dòng)軌跡，用來(lái)確定水滴在模型表面上的撞擊點(diǎn)，從而進(jìn)一步計(jì)算水滴收集系數(shù)。Eulerian法由于在計(jì)算主流場(chǎng)的同時(shí)，也計(jì)算得到每一個(gè)網(wǎng)格內(nèi)過(guò)冷水滴的質(zhì)量系數(shù)（或體積系數(shù)），因此不需要計(jì)算過(guò)冷水滴在流場(chǎng)中的軌跡。Lagrangian方法對(duì)二維或幾何形狀簡(jiǎn)單的表面應(yīng)用起來(lái)比較簡(jiǎn)便，但對(duì)多段翼型或發(fā)動(dòng)機(jī)進(jìn)口部件等三維復(fù)雜外形，由于粒子釋放位置較難確定，很難準(zhǔn)確確定撞擊極限，且計(jì)算時(shí)用時(shí)較長(zhǎng)。而Eulerian方法則不涉及這些問(wèn)題［2－4］。

本文應(yīng)用Eulerian方法計(jì)算水滴收集系數(shù)，采用Fluent商業(yè)軟件計(jì)算空氣流場(chǎng)，即將水滴視為擬流體通過(guò)Eulerian方法實(shí)現(xiàn)水滴流場(chǎng)的求解。

1 數(shù)學(xué)模型

1.1 基本假設(shè)

在Eulerian方法中引入容積分?jǐn)?shù)α的概念。它表示單位體積的混合流體中某一相所占的體積，第q相的容積分?jǐn)?shù)可以表示為αq，q相的體積可以定義為

在結(jié)冰條件下，由于水滴的容積分?jǐn)?shù)較小，一般在10－6量級(jí)，可認(rèn)為空氣和水滴之間是單向耦合的，即空氣流場(chǎng)影響水滴，而水滴對(duì)空氣流場(chǎng)的影響可以忽略［5－6］。這使得空氣相和水滴相可以分離求解：先得到空氣流場(chǎng)，再使用空氣流場(chǎng)結(jié)果計(jì)算水滴相。

對(duì)Eulerian法中空氣和水滴的控制方程可以進(jìn)行以下假設(shè)：

（a）過(guò)冷水滴在云層中均勻分布，且以球形存在，不凝聚、不分解、不變形、不破裂；（b）過(guò)冷水滴在云層中，不與空氣進(jìn)行熱、質(zhì)交換，即過(guò)冷水滴不蒸發(fā)、物性參數(shù)不變化；（c）忽略作用在水滴上的空氣阻力對(duì)空氣動(dòng)力的影響；（d）作用在水滴上的作用力僅有阻力和重力；（e）作用在水滴上的阻力是穩(wěn)態(tài)的；（f）空氣的紊流脈動(dòng)不影響水滴的運(yùn)動(dòng)。

1.2 控制方程

過(guò)冷水滴穩(wěn)態(tài)控制方程如下［5］：

其中，u是水滴的速度矢量，ua是空氣的速度矢量，F(xiàn)是作用于水滴上的除阻力之外的力。K為空氣－水滴交換系數(shù)，由下式給出：

式中，μa為空氣的動(dòng)力黏度，dd是水滴直徑，f為阻力函數(shù)，對(duì)于球形物體的阻力函數(shù)可以用下式進(jìn)行計(jì)算：

式中，CD為阻力系數(shù)，Re為相對(duì)雷諾數(shù)，其具體計(jì)算公式如下：

式中，ρa(bǔ)為空氣密度。

1.3 邊界條件

本文借助Fluent軟件計(jì)算空氣流場(chǎng)，其中進(jìn)口采用速度進(jìn)口邊界條件，出口采用壓力出口邊界條件，無(wú)滑移壁面邊界條件，使用k－e湍流模型計(jì)算空氣流場(chǎng)。將水滴設(shè)為擬流體，通過(guò)UDS模塊求解水滴相控制方程。

水滴相的入口采用速度入口邊界條件，速度大小與空氣速度相等，出口采用壓力出口邊界條件。對(duì)于水滴相的壁面邊界條件需要特殊處理［7］，固體壁面對(duì)于空氣來(lái)說(shuō)只是普通的固壁，但是對(duì)于水滴相來(lái)說(shuō)卻相當(dāng)于一個(gè)單向出口，這個(gè)出口只容許水滴流出計(jì)算域而不容許流入計(jì)算域。這是由于水滴在撞擊到壁面時(shí)會(huì)積聚在壁面上而不會(huì)從壁面上繞流經(jīng)過(guò)，從而再一次進(jìn)入計(jì)算域中，故而水滴在撞擊到壁面時(shí)只能流出計(jì)算域而不能再次進(jìn)入。水滴在壁面附近的運(yùn)動(dòng)情況如圖1所示，當(dāng)水滴撞擊到壁面上時(shí)，即水滴速度ud與壁面微元單位法向向量點(diǎn)積小于等于零時(shí)，水滴的體積分?jǐn)?shù)與水滴速度保持不變；在非撞擊區(qū)域，壁面微元單位法向量與水滴速度ud點(diǎn)積大于零時(shí)，水滴體積分?jǐn)?shù)取零，水滴速度取下游網(wǎng)格單元速度值。

圖1 水滴相壁面邊界條件

1.4 局部水滴收集系數(shù)計(jì)算

局部水滴收集系數(shù)β為當(dāng)?shù)乇砻嫖⒃獙?shí)際水滴收集量與最大可能收集量的比值，其計(jì)算公式為：

式中，Wβ為實(shí)際水滴收集量，Wβmax為最大可能水滴收集量，ud為當(dāng)?shù)厮嗡俣仁噶浚琻為部件表面的法向方向矢量，αx為當(dāng)?shù)厮误w積分?jǐn)?shù)，V0為水滴初始速度，ρ為水滴密度，α0為初始水滴體積分?jǐn)?shù)，可由遠(yuǎn)場(chǎng)液態(tài)水含量計(jì)算得出。

2 計(jì)算模型及計(jì)算條件

2.1 計(jì)算模型

本文以NACA0012機(jī)翼作為計(jì)算模型討論不同氣象條件及飛行速度對(duì)水滴撞擊特性的影響。計(jì)算模型中，翼型弦長(zhǎng)1 000mm，在網(wǎng)格劃分過(guò)程中，在圓柱周圍劃分10倍遠(yuǎn)場(chǎng)網(wǎng)格［8］，遠(yuǎn)場(chǎng)長(zhǎng)度方向適當(dāng)延伸。計(jì)算網(wǎng)格如圖2－3所示。

圖2 流場(chǎng)計(jì)算網(wǎng)格

圖3 NACA0012計(jì)算網(wǎng)格

2.2 計(jì)算條件

本文關(guān)注不同氣象條件及飛行速度對(duì)飛行器表面水滴撞擊特性的影響，為便于對(duì)比分析，并從分析中獲得一般規(guī)律，選取計(jì)算條件如表1所示。

表1 飛行速度及氣象條件表

3 計(jì)算結(jié)果及分析

利用上述計(jì)算方法對(duì)表1中所列計(jì)算條件下NACA0012機(jī)翼的水滴撞擊特性進(jìn)行了數(shù)值模擬計(jì)算。水滴在機(jī)翼空氣流場(chǎng)中運(yùn)動(dòng)，其運(yùn)動(dòng)軌跡受到空氣黏性阻力的影響。由于水滴具有慣性，因此其保持原有運(yùn)動(dòng)形式的趨勢(shì)撞擊到機(jī)翼表面。水滴的慣性及受到的阻力與水滴的平均容積直徑、來(lái)流速度等有關(guān)。

選取狀態(tài)1，取距翼根50mm截面處的計(jì)算結(jié)果與文獻(xiàn)［9］中結(jié)果作對(duì)比，水滴收集系數(shù)分布曲線如圖4所示。由圖4可知計(jì)算結(jié)果與文獻(xiàn)結(jié)果走勢(shì)符合，因此證明本文計(jì)算方法是合理的。

3.1 飛行速度對(duì)水滴撞擊特性的影響

選取飛行速度不同而其他條件相同的狀態(tài)4，6，7對(duì)水滴撞擊特性進(jìn)行分析，計(jì)算結(jié)果如圖5所示。由計(jì)算結(jié)果可知，當(dāng)其他飛行參數(shù)和氣象條件不變時(shí)，飛行速度減小，水滴受到的慣性力減小，慣性對(duì)水滴的影響減小，使得阻力對(duì)水滴的運(yùn)動(dòng)軌跡影響增大。因此，水滴撞擊到機(jī)翼表面的趨勢(shì)有所減弱，從而可知機(jī)翼水滴局部水收集系數(shù)逐漸減小。

圖4 狀態(tài)1條件下機(jī)翼水滴收集系數(shù)分布曲線

圖5 不同飛行速度下機(jī)翼水滴收集系數(shù)分布曲線

3.2 水滴平均容積直徑對(duì)水滴撞擊特性的影響

選取水滴平均容積直徑不同而其他條件相同的狀態(tài)1，2，3對(duì)水滴撞擊特性進(jìn)行分析。隨著水滴平均容積直徑的增大，水滴本身的慣性也相應(yīng)增大，慣性力對(duì)水滴的影響增加，黏性阻力對(duì)水滴軌跡的影響減小，直徑大的液滴由于慣性作用，在機(jī)翼壁面附近不能像空氣流場(chǎng)一樣很輕松地繞流，而是直接撞擊在壁面上，故而，有圖6中所示的計(jì)算結(jié)果，即水滴收集系數(shù)隨著水滴平均容積直徑的增大而增大，且撞擊范圍增大。值得注意的是，在水滴平均容積直徑超過(guò)一定數(shù)值時(shí)（一般認(rèn)為是50μm），水滴將會(huì)產(chǎn)生沿程參數(shù)變化，并在壁面附近產(chǎn)生飛濺、反彈等物理現(xiàn)象，使得撞擊特性發(fā)生變化［10］。

圖6 不同平均容積直徑下機(jī)翼水滴收集系數(shù)分布曲線

3.3 不同液態(tài)水含量對(duì)水滴撞擊特性的影響

選取液態(tài)水含量不同而其他狀態(tài)相同的狀態(tài)3，4，5對(duì)水滴撞擊特性進(jìn)行分析，由于在控制方程中，液態(tài)水含量的值可以換算出相應(yīng)的體積分?jǐn)?shù)α值，當(dāng)水滴流場(chǎng)處于穩(wěn)態(tài)時(shí)，即時(shí)間導(dǎo)數(shù)項(xiàng)為零，方程中體積分?jǐn)?shù)項(xiàng)可以消去，故而液態(tài)水含量對(duì)水滴撞擊特性影響不大，計(jì)算結(jié)果（見圖7）也證明了理論分析的結(jié)果，即當(dāng)液態(tài)水含量變化時(shí)，機(jī)翼附近水滴收集系數(shù)基本不變。

圖7 不同液態(tài)水含量下機(jī)翼水滴收集系數(shù)分布曲線

4 結(jié)論

本文基于歐拉方法計(jì)算了不同飛行速度和氣象條件下NACA0012機(jī)翼水滴撞擊特性，使用Fluent商用軟件計(jì)算空氣流場(chǎng)，并通過(guò)編制UDS程序求解水滴相控制方程，對(duì)水滴相撞擊邊界條件進(jìn)行了特殊處理使之更好地反映實(shí)際物理過(guò)程。從計(jì)算結(jié)果中可以得出：隨著飛行速度的增加，局部水滴收集系數(shù)明顯上升，水滴撞擊范圍加大；隨著水滴平均容積直徑的增加，水滴收集系數(shù)上升明顯，水滴撞擊范圍加大；而液態(tài)水含量對(duì)水滴撞擊特性影響不大。

由于水滴平均容積直徑對(duì)水滴收集系數(shù)影響很大，且在水滴直徑達(dá)到一定數(shù)值時(shí)（一般為50μm以上），水滴參數(shù)沿程會(huì)發(fā)生改變，并會(huì)在壁面上產(chǎn)生飛濺、反彈等現(xiàn)象。故而，在過(guò)冷大水滴條件下的撞擊特性還有待深入研究。

［1］ 裘燮綱.韓鳳華.飛機(jī)防冰系統(tǒng)［M］.北京：航空專業(yè)教材編審組出版，1985：44－58.

［2］ Mark G Potapczuk.LEWICE／E：An Euler Based Ice Accretion Code［R］.AIAA 92－0037，1992.

［3］ Hedde T，Guffond D.Improvement of the ONERA 3D icing code comparison with 3Dexperimental shapes［R］.AIAA 93－0169，1993.

［4］ Morency F，Heloise Beaugendre，Wagdi G Habashi.FENSAP－ICE：a study of ice shapes on droplet impingment［R］.AIAA 2003－1223，2003.

［5］ Suikno Wirogo，Shashidhar Sriramblatla.An eulerian method to calculate the collection efficiency on two and three dimensional bodies［R］.AIAA 2003－1073，2003.

［6］ Chirag Bhargava，Eric Loth，Mark Potapczuk.Aerodynamic simulations of the NASA glenn icing research tunnel［R］.AIAA 2003－566，2003.

［7］ 楊勝華，林貴平，申曉斌.三維復(fù)雜表面水滴撞擊特性計(jì)算［J］.航空動(dòng)力學(xué)報(bào)，2010，25（2）：284－290.

［8］ ANSYS Inc.ANSYS fluent12.0User’s manual［M］.New Hampshire：ANSYS Inc，2009.

［9］ Gray A Ruff，Brian M Berkowitz.User’s manual for the NASA lewis ice accretion prediction code（LEWICE）［R］.NASA，1990：1－27.

［10］ Brahimi M T，Tran P，Chocron D，et al.Effect of supercooled large droplets on ice accretion characteristics［R］.AIAA 97－0306，1997.