基于歐拉方法的機翼表面水滴撞擊特性數值模擬

劉 鑫，邢玉明，杜 佶

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

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

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

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

本文應用Eulerian方法計算水滴收集系數，采用Fluent商業軟件計算空氣流場，即將水滴視為擬流體通過Eulerian方法實現水滴流場的求解。

1 數學模型

1.1 基本假設

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

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

對Eulerian法中空氣和水滴的控制方程可以進行以下假設：

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

1.2 控制方程

過冷水滴穩態控制方程如下［5］：

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

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

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

式中，ρa為空氣密度。

1.3 邊界條件

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

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

圖1 水滴相壁面邊界條件

1.4 局部水滴收集系數計算

局部水滴收集系數β為當地表面微元實際水滴收集量與最大可能收集量的比值，其計算公式為：

式中，Wβ為實際水滴收集量，Wβmax為最大可能水滴收集量，ud為當地水滴速度矢量，n為部件表面的法向方向矢量，αx為當地水滴體積分數，V0為水滴初始速度，ρ為水滴密度，α0為初始水滴體積分數，可由遠場液態水含量計算得出。

2 計算模型及計算條件

2.1 計算模型

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

圖2 流場計算網格

圖3 NACA0012計算網格

2.2 計算條件

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

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

3 計算結果及分析

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

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

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

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

圖4 狀態1條件下機翼水滴收集系數分布曲線

圖5 不同飛行速度下機翼水滴收集系數分布曲線

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

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

圖6 不同平均容積直徑下機翼水滴收集系數分布曲線

3.3 不同液態水含量對水滴撞擊特性的影響

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

圖7 不同液態水含量下機翼水滴收集系數分布曲線

4 結論

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

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

［1］ 裘燮綱.韓鳳華.飛機防冰系統［M］.北京：航空專業教材編審組出版，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］ 楊勝華，林貴平，申曉斌.三維復雜表面水滴撞擊特性計算［J］.航空動力學報，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.