吳孟龍 常士楠 冷夢堯 王 超

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

飛機在含有大量過冷水滴的云層中飛行時，其迎風(fēng)面會發(fā)生結(jié)冰現(xiàn)象.發(fā)動機進口旋轉(zhuǎn)帽罩處于飛機的迎風(fēng)面，也會發(fā)生結(jié)冰現(xiàn)象.旋轉(zhuǎn)帽罩結(jié)冰會降低發(fā)動機入口氣流品質(zhì)，使發(fā)動機性能降低.并且，旋轉(zhuǎn)帽罩表面結(jié)冰脫落，被吸入發(fā)動機內(nèi)部，會導(dǎo)致發(fā)動機損毀，造成重大飛行事故［1］.因此十分有必要對旋轉(zhuǎn)帽罩進行結(jié)冰研究.飛機結(jié)冰數(shù)值模擬研究一般分為空氣流場求解、水滴軌跡及撞擊特性求解、結(jié)冰表面熱平衡分析以及結(jié)冰冰形計算4個部分.本文主要關(guān)注于旋轉(zhuǎn)帽罩水滴軌跡及撞擊特性分析.

過冷水滴在空氣流場中的運動屬于氣粒兩相流動.按照對水滴顆粒的不同處理方式，主要有兩種氣粒兩相流動研究方法［2］:

1)把氣體當(dāng)作連續(xù)介質(zhì)，而將顆粒視為離散體系，即拉格朗日法;

2)把氣體與顆粒都看成共同存在且相互滲透的連續(xù)介質(zhì)，把顆粒作為擬流體，即歐拉法.

目前，國內(nèi)外對于靜止部件的水滴撞擊特性的求解既有采用拉格朗日法［3］，也有采用歐拉法［4-6］.而對于旋轉(zhuǎn)部件的水滴撞擊特性，國內(nèi)外研究較少，并且多數(shù)采用拉格朗日法［7-8］.

本文在歐拉法求解靜止部件水滴撞擊特性的基礎(chǔ)上，提出一種旋轉(zhuǎn)帽罩水滴撞擊特性數(shù)值模擬方法，實現(xiàn)了歐拉法求解旋轉(zhuǎn)帽罩表面的水滴撞擊特性.通過對不同轉(zhuǎn)速旋轉(zhuǎn)帽罩水滴撞擊特性進行數(shù)值模擬研究，發(fā)現(xiàn)旋轉(zhuǎn)帽罩轉(zhuǎn)速越大，其表面水滴撞擊極限越小.但是，在典型飛行條件與氣象條件下，離心力對水滴運動軌跡的影響遠(yuǎn)小于慣性力對其的影響，旋轉(zhuǎn)帽罩轉(zhuǎn)速對水滴撞擊特性影響較小.

1 水滴相數(shù)學(xué)模型

基于歐拉法在靜止部件水滴撞擊特性的成功應(yīng)用［9-12］，本文提出了一種旋轉(zhuǎn)帽罩水滴撞擊特性求解方法.本文借助于Fluent軟件，采用S-A湍流模型實現(xiàn)旋轉(zhuǎn)帽罩空氣流場的求解.而將水滴視為擬流體，采用歐拉法，借助于Fluent UDS模塊，實現(xiàn)旋轉(zhuǎn)帽罩水滴流場的求解.

結(jié)冰條件下，由于水滴的體積分?jǐn)?shù)一般在10-6量級，可認(rèn)為空氣和水滴之間是單向耦合的，即空氣流場影響水滴的運動，而水滴的運動對空氣流場沒有影響［9，13］.這使得空氣相和水滴相可以分離求解:先得到空氣流場，再計算水滴流場.并且對水滴進行一些合理的假設(shè):①水滴是球形的，無變形或破裂;②水滴之間無碰撞、聚合，撞擊到壁面后無飛濺;③水滴和周圍空氣無質(zhì)量和熱量的傳遞;④只考慮作用在水滴上的空氣阻力.

1.1 控制方程

基于以上假設(shè)，穩(wěn)態(tài)水滴相控制方程為

式(1)右邊項為耗散函數(shù)，主要作用為消除數(shù)值計算過程中產(chǎn)生的奇異值，使計算可以順利進行并得到合理的數(shù)值，在1.3節(jié)中將會詳細(xì)闡述.式中，α是水滴體積分?jǐn)?shù);ρd和ud分別為水滴的密度和速度;ua為空氣速度;K為水滴慣性系數(shù)，其表達(dá)式為

其中μa為空氣的動力黏性系數(shù).基于假設(shè)①認(rèn)為水滴為球形，因此可應(yīng)用球形水滴的阻力系數(shù)函數(shù)f，其表達(dá)式為

其中，水滴的阻力系數(shù)CD為

相對雷諾數(shù)Red計算公式為

1.2 邊界條件

旋轉(zhuǎn)帽罩空氣-水滴流場數(shù)值模擬中，空氣相與水滴相的邊界條件處理方式不同.其中，空氣相入口采用速度入口邊界條件，出口采用出流邊界條件，壁面采用動壁面中的旋轉(zhuǎn)、無滑移壁面邊界條件.

水滴相的入口采用速度入口邊界條件，速度大小與空氣速度保持一致;出口則采用出流邊界條件.水滴相的壁面邊界條件需要進行特殊處理.這是由于水滴撞到壁面后，會積聚到壁面上，從而從流場中離開.因此，在水滴撞擊區(qū)域，如圖1所示情況a和b，壁面微元單位法向量n與水滴速度ud點積n·ud≤0時，水滴的體積分?jǐn)?shù)與水滴速度保持不變;在非撞擊區(qū)域，如圖1所示情況c，n·ud＞0時，水滴的體積分?jǐn)?shù)取零，水滴速度取下游網(wǎng)格單元的速度值［9-10，12］.

圖1 水滴相壁面邊界條件Fig.1 Wall boundary condition of droplet flowfield

1.3 耗散函數(shù)

由于水滴相控制方程屬于雙曲型方程［10］，數(shù)值迭代過程會出現(xiàn)局部數(shù)值奇異，因此常采用一些耗散函數(shù)，將局部奇異值耗散，使計算能夠得到收斂.本文采用的耗散函數(shù)［9］為

其中，αp為當(dāng)前網(wǎng)格的水滴體積分?jǐn)?shù);αNi為周圍網(wǎng)格的水滴體積分?jǐn)?shù);a為經(jīng)驗常數(shù)，a取值過大，會使得計算結(jié)果誤差增大，a取值過小，則會使計算結(jié)果發(fā)散.本文根據(jù)文獻［9］中a的取值，并通過與Fensap進行算例計算結(jié)果比較，取a=0.001.本文通過編寫UDF函數(shù)，改變材料項中的耗散值，從而將耗散函數(shù)添加到水滴相控制方程.

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

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

式中，n為壁面單元的單位法向量;V∞為自由來流的水滴速度;α為壁面單元的水滴體積分?jǐn)?shù);α0為入口處單元的水滴體積分?jǐn)?shù).

2 計算模型與計算條件

2.1 網(wǎng)格劃分

由于橢圓型帽罩較易結(jié)冰［14］，因此本文采用橢圓型旋轉(zhuǎn)帽罩作為計算數(shù)模，其中長軸半徑為0.09 m，短軸半徑為0.06 m.在網(wǎng)格劃分過程中，對旋轉(zhuǎn)帽罩尾部進行適當(dāng)?shù)难由?在旋轉(zhuǎn)帽罩前緣、尾部以及周圍劃分10倍遠(yuǎn)場網(wǎng)格.所劃網(wǎng)格如圖2所示.

圖2 旋轉(zhuǎn)帽罩計算域網(wǎng)格Fig.2 Computational mesh of the spinner

2.2 計算條件

本文關(guān)注重點在于旋轉(zhuǎn)速度對于旋轉(zhuǎn)帽罩水滴撞擊特性的影響，因此，取不同旋轉(zhuǎn)速度的飛行及氣象條件作為計算條件，分析研究旋轉(zhuǎn)速度對旋轉(zhuǎn)帽罩水滴撞擊特性的影響.參照典型民用航空發(fā)動機旋轉(zhuǎn)帽罩轉(zhuǎn)速，選取旋轉(zhuǎn)帽罩轉(zhuǎn)速與其他飛行條件如表1所示.

表1 飛行及氣象條件Table 1 Flight and weather conditions

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

利用上述計算方法對表1中所列計算條件下的旋轉(zhuǎn)帽罩水滴撞擊特性進行了數(shù)值模擬計算.

水滴在旋轉(zhuǎn)帽罩空氣流場中運動，其運動軌跡受到空氣的黏性阻力與離心力的影響.由于水滴具有慣性，因此會保持原有運動形式，撞擊到旋轉(zhuǎn)帽罩表面.而水滴受到的黏性阻力與離心力則會使水滴偏離原始軌跡，阻礙水滴撞擊旋轉(zhuǎn)帽罩.水滴的慣性、受到的黏性阻力與離心力跟水滴的平均容積直徑、來流速度以及旋轉(zhuǎn)帽罩旋轉(zhuǎn)速度有關(guān).由于旋轉(zhuǎn)帽罩表面的水滴撞擊特性取決于水滴的運動軌跡.因此，本文對不同旋轉(zhuǎn)速度、來流速度以及水滴平均容積直徑的計算結(jié)果進行了分析.

本文首先取狀態(tài)3計算結(jié)果與成熟的商業(yè)結(jié)冰計算軟件Fensap進行對比.沿旋轉(zhuǎn)帽罩軸向做y=0截面，作局部水滴收集系數(shù)沿弦長變化曲線如圖3所示.

圖3 狀態(tài)3，y=0截面局部水滴收集系數(shù)分布曲線Fig.3 Case 3，cross-section y=0，local water droplets collect coefficient distribution curves

由圖3可以看出，本文計算結(jié)果與Fensap計算結(jié)果十分吻合，因此可以證明本文計算方法是合理的.由于本文采用的耗散函數(shù)與Fensap所采用的耗散函數(shù)不同，導(dǎo)致帽罩后部存在一定差異.

3.1 旋轉(zhuǎn)帽罩轉(zhuǎn)速對水滴撞擊特性影響

選取旋轉(zhuǎn)速度不同，其他條件相同的狀態(tài)1～3的局部水滴撞擊特性計算結(jié)果進行分析.沿旋轉(zhuǎn)帽罩軸向做y=0截面，作局部水滴收集系數(shù)沿弦長變化曲線如圖4所示.由圖4可以看出，各轉(zhuǎn)速下，旋轉(zhuǎn)帽罩局部水滴收集系數(shù)分布結(jié)果大致相同.隨著轉(zhuǎn)速的增加，水滴撞擊極限逐漸減小.這與文獻［7，15］采用拉格朗日法求解旋轉(zhuǎn)帽罩水滴撞擊特性所得結(jié)論相同.

圖4 狀態(tài)1～3，y=0截面局部水滴收集系數(shù)分布曲線Fig.4 Case 1 ～3，cross-section y=0，local water droplets collect coefficient distribution curves

圖5所示為狀態(tài)3的水滴運動軌跡，可以看出水滴偏離基準(zhǔn)線，隨氣流發(fā)生旋轉(zhuǎn)，但其偏轉(zhuǎn)角度較小，即水滴受到旋轉(zhuǎn)帽罩的旋轉(zhuǎn)影響較小.

圖5 狀態(tài)3，速度78 m/s，轉(zhuǎn)速5000 r/min水滴軌跡圖Fig.5 Case 3，velocity=78 m/s，rotate speed=5000 r/min，the scheme of droplet trajectory

3.2 來流速度對水滴撞擊特性影響

選取來流速度不同，其他條件相同的狀態(tài)3～5的水滴撞擊特性計算結(jié)果進行分析.當(dāng)其他飛行條件與氣象條件保持不變，來流速度減小時，水滴受到的慣性力減小.慣性力對水滴運動軌跡的影響減小，相應(yīng)地黏性阻力與離心力對水滴的運動軌跡影響增大.因此，水滴更易偏離原始軌跡，并且撞擊到旋轉(zhuǎn)帽罩表面的趨勢減弱，從而得到如圖6所示旋轉(zhuǎn)帽罩表面局部水滴收集系數(shù)明顯減小，水滴撞擊極限范圍也明顯減小的結(jié)果.

圖6 狀態(tài)3～5，y=0截面局部水滴收集系數(shù)分布曲線Fig.6 Case 3 ～5，cross-section y=0，local water droplets collect coefficient distribution curves

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

選取水滴平均容積直徑不同，其他條件相同的狀態(tài)3、狀態(tài)6與狀態(tài)7的水滴撞擊特性計算結(jié)果進行分析.當(dāng)其他飛行條件與氣象條件保持不變，水滴平均容積直徑減小時，水滴本身具有的慣性減小.所以，隨著水滴平均容積直徑的減小，慣性力對于水滴運動軌跡的影響會減小，相應(yīng)地水滴受到的黏性阻力與離心力對水滴運動軌跡的影響會增大.因此，得到如圖7所示，隨著水滴平均容積直徑的減小，旋轉(zhuǎn)帽罩局部水滴收集系數(shù)明顯減小，水滴撞擊極限范圍也明顯減小的結(jié)果.

圖7 狀態(tài)3，6，7，y=0截面局部水滴收集系數(shù)分布曲線Fig.7 Case 3，6，7，cross-section y=0，local water droplets collect coefficient distribution curves

由以上計算結(jié)果與分析可知，當(dāng)來流速度較小，水滴平均容積直徑較小時，水滴運動軌跡受黏性阻力與離心力的影響較大.

為了使旋轉(zhuǎn)帽罩轉(zhuǎn)速對水滴撞擊特性的影響更加明顯，應(yīng)該盡量減小慣性力對水滴運動軌跡的影響，加大黏性阻力與離心力對水滴運動軌跡的影響.因此采用來流速度較小，水滴平均容積直徑較小的狀態(tài)8、狀態(tài)9進行比較.由圖8可以看出，隨著旋轉(zhuǎn)帽罩轉(zhuǎn)速的增大，駐點處局部水滴收集系數(shù)略有減小，并且水滴撞擊極限也有所減小.

圖8 狀態(tài)8和9，y=0截面局部水滴收集系數(shù)分布曲線Fig.8 Case 8 and 9，cross-section y=0，local water droplets collect coefficient distribution curves

4 結(jié)論

1)歐拉法求解旋轉(zhuǎn)帽罩的空氣-水滴兩相流流場，繼而得到旋轉(zhuǎn)帽罩的水滴撞擊特性是可行的.

2)在飛行條件與氣象條件保持不變的情況下，隨著旋轉(zhuǎn)帽罩轉(zhuǎn)速增加，水滴撞擊極限稍有減小.

3)在典型飛行條件與氣象條件下，離心力對水滴運動軌跡的影響遠(yuǎn)小于慣性力對其的影響.旋轉(zhuǎn)帽罩轉(zhuǎn)速對水滴的運動軌跡影響較小，水滴撞擊特性變化不大.

References)

［1］裘燮綱，韓鳳華.飛機防冰系統(tǒng)［M］.北京:航空專業(yè)教材編審組，1985:50-53 Qiu Xiegang，Han Fenghua.Aircraft anti-icing system［M］.Beijing:Compilation and Examination Group of Aero Specialized Teaching Materials，1985:50-53(in Chinese)

［2］周力行.湍流氣粒兩相流動和燃燒的理論與數(shù)值模擬［M］.北京:科學(xué)出版社，1994:293-298 Zhou Lixing.Theory and numerical modeling of turbulent gasparticle flows and combustion［M］.Beijing:Science Press，1994:293-298(in Chinese)

［3］楊倩，常士楠，袁修干.水滴撞擊特性的數(shù)值計算方法研究［J］.航空學(xué)報，2002，23(2):173-176 Yang Qian，Chang Shinan，Yuan Xiugan.Study on numerical method for determining the droplet trajectories［J］.Acta Aeronuatica et Astronautica Sinca，2002，23(2):173-176(in Chinese)

［4］常士楠，蘇新明，邱義芬.三維機翼結(jié)冰模擬［J］.航空學(xué)報，2011，32(2):212-222 Chang Shinan，Su Xinming，Qiu Yifen.Ice accretion simulation on three dimensional wings［J］.Acta Aeronautica et Astronautica Sinca，2011，32(2):212-222(in Chinese)

［5］申曉斌，林貴平，楊勝華.三維發(fā)動機進氣道水滴撞擊特性分析［J］.北京航空航天大學(xué)學(xué)報，2011，37(1):1-5 Shen Xiaobin，Lin Guiping，Yang Shenghua.Analysis on three dimensional water droplets impingement characteristics of engine inlet［J］.Journal of Beijing University of Aeronautics and Astronautics，2011，37(1):1-5(in Chinese)

［6］ Bourgault Y，Boutanios Z，Habashi W G.Three-dimensional eulerian approach to droplet impingement simulation using FENSAP-ICE，part 1:model，algorithm，and validation［J］.Journal of Aircraft，2000，37(1):95-103

［7］趙秋月，董威，朱劍鋆.發(fā)動機旋轉(zhuǎn)整流帽罩的水滴撞擊特性分析［J］.燃?xì)鉁u輪試驗與研究，2011，24(4):32-35 Zhao Qiuyue，Dong Wei，Zhu Jianjun.Droplets impinging characteristic analysis of the rotating fairing of aero-engine［J］.Gas Turbine Experiment and Research，2011，24(4):32-35(in Chinese)

［8］ Das K，Hamed A，Basu D.Ice shape prediction for turbofan rotating blades［R］.AIAA 2006-0209，2006

［9］楊勝華，林貴平，申曉斌.三維復(fù)雜表面水滴撞擊特性計算［J］.航空動力學(xué)報，2010，25(2):284-290 Yang Shenghua，Lin Guiping，Shen Xiaobin.Water droplet impingement prediction for three-dimensional complex surfaces［J］.Journal of Aerospace Power，2010，25(2):284-290(in Chinese)

［10］ Slater S A，Young J B.The calculation of inertial particle transport in dilute gas-particle flows［J］.International Journal of Multiphase Flow，2001，27(1):61-87

［11］ Wirogo S，Srirambhatla S.An Eulerian method to calculate the collection efficiency on two and three dimensional bodies［R］.AIAA 2003-1073，2003

［12］ Tong X，Luke E A.Eulerian simulation of icing collection efficiency using a singularity diffusion model［R］.AIAA 2005-1246，2005

［13］ Cao Y，Ma C，Zhang Q，et al.Numerical simulation of ice accretions on an aircraft wing［J］.Aerospace Science and Technology，2012，23(1):296-304

［14］ Linke D A.Systems of commercial turbofan engines［M］.Berlin:Springer，2008:179-183

［15］李西園.旋轉(zhuǎn)部件結(jié)冰縮比試驗相似性研究［D］.北京:北京航空航天大學(xué)，2010 Li Xiyuan.Similarity of icing scaling test for aircraft rotating components［D］.Beijing:Beijing University of Aeronautics and Astronautics，2010(in Chinese)