劉 鑫，邢玉明，杜 佶

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

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

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

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

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

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

1.1 基本假設(shè)

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

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

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

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

1.2 控制方程

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

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

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

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

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

1.3 邊界條件

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

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

圖1 水滴相壁面邊界條件

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

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

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

2 計算模型及計算條件

2.1 計算模型

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

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

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

2.2 計算條件

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

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

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

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

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

3.1 飛行速度對水滴撞擊特性的影響

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

［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］.航空動力學(xué)報，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.