吳孟龍 常士楠 冷夢(mèng)堯 王 超

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

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

過(guò)冷水滴在空氣流場(chǎng)中的運(yùn)動(dòng)屬于氣粒兩相流動(dòng).按照對(duì)水滴顆粒的不同處理方式，主要有兩種氣粒兩相流動(dòng)研究方法［2］:

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

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

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

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

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

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

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

1.1 控制方程

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

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

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

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

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

1.2 邊界條件

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

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

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

1.3 耗散函數(shù)

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

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

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

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

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

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

2.1 網(wǎng)格劃分

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

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

2.2 計(jì)算條件

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

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

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

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

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

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

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

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

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

選取旋轉(zhuǎn)速度不同，其他條件相同的狀態(tài)1～3的局部水滴撞擊特性計(jì)算結(jié)果進(jìn)行分析.沿旋轉(zhuǎn)帽罩軸向做y=0截面，作局部水滴收集系數(shù)沿弦長(zhǎng)變化曲線如圖4所示.由圖4可以看出，各轉(zhuǎn)速下，旋轉(zhuǎn)帽罩局部水滴收集系數(shù)分布結(jié)果大致相同.隨著轉(zhuǎn)速的增加，水滴撞擊極限逐漸減小.這與文獻(xià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的水滴運(yùn)動(dòng)軌跡，可以看出水滴偏離基準(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 來(lái)流速度對(duì)水滴撞擊特性影響

選取來(lái)流速度不同，其他條件相同的狀態(tài)3～5的水滴撞擊特性計(jì)算結(jié)果進(jìn)行分析.當(dāng)其他飛行條件與氣象條件保持不變，來(lái)流速度減小時(shí)，水滴受到的慣性力減小.慣性力對(duì)水滴運(yùn)動(dòng)軌跡的影響減小，相應(yīng)地黏性阻力與離心力對(duì)水滴的運(yùn)動(dòng)軌跡影響增大.因此，水滴更易偏離原始軌跡，并且撞擊到旋轉(zhuǎn)帽罩表面的趨勢(shì)減弱，從而得到如圖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 水滴平均容積直徑對(duì)水滴撞擊特性影響

選取水滴平均容積直徑不同，其他條件相同的狀態(tài)3、狀態(tài)6與狀態(tài)7的水滴撞擊特性計(jì)算結(jié)果進(jìn)行分析.當(dāng)其他飛行條件與氣象條件保持不變，水滴平均容積直徑減小時(shí)，水滴本身具有的慣性減小.所以，隨著水滴平均容積直徑的減小，慣性力對(duì)于水滴運(yùn)動(dòng)軌跡的影響會(huì)減小，相應(yīng)地水滴受到的黏性阻力與離心力對(duì)水滴運(yùn)動(dòng)軌跡的影響會(huì)增大.因此，得到如圖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

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

為了使旋轉(zhuǎn)帽罩轉(zhuǎn)速對(duì)水滴撞擊特性的影響更加明顯，應(yīng)該盡量減小慣性力對(duì)水滴運(yùn)動(dòng)軌跡的影響，加大黏性阻力與離心力對(duì)水滴運(yùn)動(dòng)軌跡的影響.因此采用來(lái)流速度較小，水滴平均容積直徑較小的狀態(tài)8、狀態(tài)9進(jìn)行比較.由圖8可以看出，隨著旋轉(zhuǎn)帽罩轉(zhuǎn)速的增大，駐點(diǎ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)帽罩的空氣-水滴兩相流流場(chǎng)，繼而得到旋轉(zhuǎn)帽罩的水滴撞擊特性是可行的.

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

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

References)

［1］裘燮綱，韓鳳華.飛機(jī)防冰系統(tǒng)［M］.北京:航空專(zhuān)業(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］周力行.湍流氣粒兩相流動(dòng)和燃燒的理論與數(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ì)算方法研究［J］.航空學(xué)報(bào)，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］常士楠，蘇新明，邱義芬.三維機(jī)翼結(jié)冰模擬［J］.航空學(xué)報(bào)，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ā)動(dòng)機(jī)進(jìn)氣道水滴撞擊特性分析［J］.北京航空航天大學(xué)學(xué)報(bào)，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ā)動(dòng)機(jī)旋轉(zhuǎn)整流帽罩的水滴撞擊特性分析［J］.燃?xì)鉁u輪試驗(yàn)與研究，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ì)算［J］.航空動(dòng)力學(xué)報(bào)，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é)冰縮比試驗(yàn)相似性研究［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)