發(fā)動機進氣帽罩防冰熱載荷的數(shù)值模擬研究

劉華，郭文，楊軍，婁德倉

發(fā)動機進氣帽罩防冰熱載荷的數(shù)值模擬研究

劉華，郭文，楊軍，婁德倉

(中國燃氣渦輪研究院，四川成都610500)

對某發(fā)動機進氣帽罩防冰熱載荷的計算方法進行研究。利用FLUENT軟件對帽罩周圍的流場進行模擬，同時采用FLUENT中的離散相模型計算帽罩周圍的水滴運動軌跡，得到帽罩表面的水滴撞擊效率；利用數(shù)值計算結(jié)果及流動換熱經(jīng)驗關(guān)系式分析帽罩表面的各項熱流，獲得該發(fā)動機帽罩的防冰熱載荷。計算中還比較了不同湍流模型對防冰熱載荷的影響。

帽罩；局部撞擊效率；防冰熱載荷；湍流模型；數(shù)值模擬

1 引言

飛機在結(jié)冰氣象條件下飛行時，許多部件(如機翼、尾翼的前緣，風(fēng)擋玻璃，進氣道前緣及發(fā)動機進氣部件(支撐、整流罩等))都會積冰，另外暴露在飛機外部的測量傳感器(如皮托管、測高度的傳感頭及測溫度的傳感頭等)也會積冰。這些部件的積冰，會影響飛機的安全飛行[1]。

由飛機結(jié)冰問題而引發(fā)的飛行事故屢次發(fā)生，給生命和財產(chǎn)造成巨大損失。國外早在上世紀(jì)30～40年代就已開始重視防冰問題，并取得長足進展。目前，國外已開發(fā)出多種結(jié)冰計算軟件(如美國NASA的LEWICE、加拿大的FENSAP-ICE、法國ONERA的積冰模擬軟件等)[2]，建立了一些大型的冰風(fēng)洞[3]。西方一些航空發(fā)達國家在上世紀(jì)70年代就擁有了較為完備的防冰系統(tǒng)設(shè)計體系，對發(fā)動機防冰能力提出了明確要求：不論發(fā)動機是否設(shè)計了防冰系統(tǒng)，都必須通過特定結(jié)冰環(huán)境下的防冰能力考核驗證。這些國家在結(jié)冰環(huán)境參數(shù)確定、發(fā)動機水?dāng)z取率計算、水撞擊特性計算、防冰表面熱平衡分析、冰聚集量預(yù)估、飛行循環(huán)分析等方面都有成熟的計算軟件支持，可進行詳細、準(zhǔn)確的分析和計算。

目前，我國的防冰系統(tǒng)設(shè)計還處于借鑒國外防冰系統(tǒng)設(shè)計結(jié)構(gòu)的水平，僅能進行較為初步的防冰系統(tǒng)流動與傳熱方案計算。對飛機結(jié)冰的研究，主要是對翼面結(jié)冰過程的數(shù)值模擬，以及水滴撞擊特性、傳質(zhì)傳熱效應(yīng)對結(jié)冰影響的分析等；根據(jù)數(shù)值模擬方法編制程序，計算過冷水的收集系數(shù)、水滴的撞擊極限、冰的種類和幾何外形以及熱防冰系統(tǒng)的能量需求等。近年來對明冰、霜冰的數(shù)值模擬和積冰熱力學(xué)模型的研究越來越多，不過大多數(shù)積冰熱力學(xué)模型以Messinger模型為基礎(chǔ)，基于熱力學(xué)第一定律，對單個控制體內(nèi)的質(zhì)量和能量平衡方程進行分析。

航空發(fā)動機進氣系統(tǒng)結(jié)冰對發(fā)動機的正常工作影響很大，會造成功率/推力損失，引起發(fā)動機損壞，甚至停車[4]，所以必須進行發(fā)動機防冰系統(tǒng)設(shè)計。目前發(fā)動機防冰系統(tǒng)大多采用熱氣防冰，即從壓氣機引出熱空氣，通過防冰引氣管路進入防冰部件內(nèi)部的防冰腔，對防冰部件內(nèi)壁面進行加熱，從而使其外壁表面溫度高于結(jié)冰溫度來防止結(jié)冰。因此，防冰熱載荷計算是發(fā)動機熱氣防冰系統(tǒng)設(shè)計的基礎(chǔ)。本文利用數(shù)值分析軟件并結(jié)合經(jīng)驗關(guān)系式，對某發(fā)動機進氣帽罩進行熱載荷計算。

2 計算模型與計算條件

本文的研究對象是某發(fā)動機錐形進氣帽罩。由于帽罩內(nèi)部結(jié)構(gòu)形狀對帽罩外部熱載荷計算影響不大，故可將其結(jié)構(gòu)簡化，如圖1所示。在帽罩的外部流場計算中，取過帽罩軸線的二維截面結(jié)構(gòu)為計算模型。其計算區(qū)域如圖2所示，從帽罩的駐點向前取3倍帽罩軸向長度為速度進口邊界，向后取5倍軸向長度為壓力出口邊界，向上、向下各取0.4倍軸向長度為壁面邊界，帽罩采用無滑移壁面邊界條件。使用ICEM軟件對計算區(qū)域劃分O型網(wǎng)格。

圖1 帽罩結(jié)構(gòu)示意圖Fig.1 Cowl structure

圖2 計算區(qū)域及網(wǎng)格Fig.2 Calculated region and grids

計算條件參數(shù)[5]為：高度1 500 m，速度73 m/s，溫度257.05 K，壓力84 689 Pa，液態(tài)水含量(LWC) 0.29 g/m3，平均水滴直徑20 μm。

3 帽罩表面水滴撞擊特性計算

3.1 外部流場計算

在利用FLUENT計算流場時，將氣流的流動視為二維定?？蓧赫承粤鲃印?/p>

本文將分別采用Standardk-ε模型、RNGk-ε模型及Realizablek-ε模型進行流場計算，對比分析不同模型下的計算差異。計算中方程用一階迎風(fēng)格式進行離散，壓力速度耦合采用SIMPLE算法。

Standardk-ε模型下計算流場結(jié)果如圖3所示。由圖中可看出，速度在駐點處最小，沿帽罩表面向后呈逐漸增大的趨勢。

圖3 Standardk-ε模型下流場速度云圖Fig.3 Velocity contour of flow-fields with Standardk-εmodel

3.2 水滴軌跡計算

在計算帽罩表面局部撞擊效率之前，首先計算水滴運動軌跡。

在FLUENT中，可采用離散項模型來計算水滴運動軌跡?？諝饬鲌鲎鳛檫B續(xù)相，水滴作為空氣流場中的離散相，兩相之間存在動量、質(zhì)量及能量交換。當(dāng)水滴穿過空氣流場運動時，其軌跡可通過當(dāng)?shù)乜諝庾饔迷谒紊系母鞣N平衡作用力計算獲得。不考慮水滴對空氣流場的影響，水滴運動方程為：

式中：u為空氣速度，up為水滴速度，μ為空氣動力粘度，ρp為水滴密度，dp為水滴直徑，Re為水滴相對空氣運動的雷諾數(shù)，CD為拽力系數(shù)。

在離散項模型中通過求解方程(1)獲得水滴運動軌跡，Standardk-ε模型下水滴軌跡計算結(jié)果見圖4。由圖中可見，一部分水滴撞擊在帽罩表面上，另一部分水滴則發(fā)生繞流現(xiàn)象。這是因為水滴的慣性比空氣大得多，其繞過帽罩運動軌跡的曲率比空氣流線小，從而導(dǎo)致一部分水滴撞擊在帽罩表面上。從圖中還可獲得水滴撞擊范圍，即上、下兩條與帽罩表面相切的水滴軌跡之間的帽罩表面。

圖4 Standardk-ε模型下水滴軌跡圖Fig.4 Trajectories of water droplets with Standardk-εmodel

3.3 帽罩表面水滴局部撞擊效率計算

局部撞擊效率β用于表征撞擊水量在水滴撞擊范圍沿表面的分布情況，定義為微元表面上實際水撞擊效率與最大可能水撞擊效率的比值，即：

式中：dy0為相鄰兩條水滴軌跡的起始縱坐標(biāo)的差值，ds為這兩條水滴軌跡的撞擊在帽罩表面的撞擊點之間的弧長。

由已計算出的水滴軌跡通過上式即可計算得到帽罩表面水滴局部撞擊效率，采用3種湍流模型所獲得的局部水滴撞擊效率的計算結(jié)果如圖5所示。橫坐標(biāo)S/C表示撞擊點與駐點之間的弧長和帽罩軸向長度的比值（駐點往上的弧長為正，往下為負）。由圖5可知，在駐點附近，不同湍流模型計算的水滴局部撞擊效率略有不同，其它位置數(shù)值差別較小。由于在計算設(shè)置時dy0相同，因此β的大小取決于弧長ds，即取決于水滴在空氣流場中的運動情況。所以，湍流模型的選取對計算結(jié)果有影響。

圖5 水滴局部撞擊效率Fig.5 Local impingement efficiency of water droplets

4 防冰熱載荷的計算

4.1 熱流分析

水滴撞擊帽罩的范圍為防冰熱載荷計算區(qū)域。水滴撞擊帽罩表面后，一部分水蒸發(fā)，另一部分水沿帽罩表面向后流動，并在其表面形成水膜。將帽罩分為若干個微元段，取其中一個微元段進行熱平衡分析，建立能量守恒方程和連續(xù)方程。

微元段的連續(xù)方程為：

式中：m˙in為流入該微元段的質(zhì)量流量密度；m˙imp為水滴撞擊在該微元段上的質(zhì)量流量密度，且m˙imp= v∞β·LWC，v∞為水滴來流速度；m˙out為流出該微元段的質(zhì)量流量密度；m˙evap為該微元段上蒸發(fā)的質(zhì)量流量密度。

當(dāng)持續(xù)增加熱載荷時，由能量守恒方程可得出該微元段防冰熱載荷q的計算公式：

式中：r為附面層恢復(fù)系數(shù)，hair為空氣對流換熱系數(shù)，ve為附面層空氣速度，cp,air為空氣比熱。

(2)對流換熱熱流密度qair

式中：Twater為該微元段的水膜溫度，Trec為氣流恢復(fù)溫度。

(3)水膜流動帶入與帶出的熱流密度qwaterin、qwaterout

(1)空氣對附面層的氣動加熱熱流密度qaero

式中：cp,water為水在Twater下的比熱。

(4)撞擊水的動能轉(zhuǎn)化成的熱流密度qdropin

(5)對撞擊水的加熱熱流密度qdroph

式中：T∞為水滴溫度。

(6)水蒸發(fā)熱流密度qevap

式中：i為水的蒸發(fā)潛熱。

4.2 對流換熱系數(shù)

對流換熱系數(shù)是熱氣防冰系統(tǒng)計算中的一個重要參數(shù)，空氣對附面層的氣動加熱熱流密度、對流換熱熱流密度和水蒸發(fā)熱流密度的大小都與其有關(guān)。文中帽罩表面的對流換熱系數(shù)由換熱準(zhǔn)則公式計算得到。

對于帽罩前緣(自駐點的角度θ小于80°)區(qū)域，其表面對流換熱系數(shù)采用文獻[6]中的圓柱表面強制對流換熱經(jīng)驗公式進行計算：

式中：ReD是以帽罩最前端當(dāng)量球體直徑D為特征尺度的來流雷諾數(shù)，λ為空氣的導(dǎo)熱系數(shù)，Pr為空氣的普朗特數(shù)。

對于帽罩后方區(qū)域，其表面對流換熱系數(shù)用文獻[7]中平板表面強制對流換熱經(jīng)驗公式進行計算：

式中：ReS是以計自駐點的壁面長度S為特征尺度的當(dāng)?shù)乩字Z數(shù)。

4.3 防冰熱載荷計算

根據(jù)上述公式，可求解各微元段的各項熱流密度值及防冰熱載荷。由于二維帽罩軸對稱，故只取自駐點向上的防冰區(qū)域進行分析計算。要使帽罩表面水膜不結(jié)冰，所加防冰熱載荷應(yīng)使帽罩表面水膜溫度保持在結(jié)冰點以上。利用上述條件，得出帽罩的防冰熱載荷計算結(jié)果。圖6、圖7分別為水膜溫度275.15 K時對流換熱系數(shù)和Standardk-ε模型下各項熱流密度的分布圖，表1為水膜溫度275.15 K時三種湍流模型計算的帽罩防冰熱載荷，圖8為不同湍流模型下熱載荷密度分布，圖9為帽罩在Stan?dardk-ε模型下不同水膜溫度所需的熱載荷。

由圖6可看出，帽罩表面駐點的換熱系數(shù)比其它區(qū)域大。

由圖7可得出，帽罩的防冰熱載荷密度、對流換熱熱流密度和水蒸發(fā)熱流密度在駐點處較大，這是因為駐點處對流換熱系數(shù)大；防冰熱載荷密度從駐點沿壁面往后逐漸減小，對流換熱熱流密度、水蒸發(fā)熱流密度的分布趨勢一致；撞擊水的加熱熱流密度在駐點處較大，從駐點沿壁面向后逐漸減小，這與水滴局部撞擊效率的分布趨勢一致。空氣對附面層的氣動加熱熱流密度在駐點處較小，沿駐點向后逐漸增加，這是因為沿駐點往后氣流速度逐漸增加。撞擊水的動能轉(zhuǎn)化成的熱流密度很小。水膜流動帶出與帶入的熱流密度的差值沿駐點往后逐漸減小至0。從圖中還可看出，除撞擊水動能轉(zhuǎn)化成的熱流密度對防冰熱載荷密度分布影響不大外，其它熱流密度對防冰熱載荷密度分布影響較大。

表1 不同湍流模型下的防冰熱載荷Table 1 Thermal loads calculated with different turbulence models

圖6 對流換熱系數(shù)分布Fig.6 Distribution of convective heat transfer coefficients

圖7 Standardk-ε模型下各項熱流密度分布Fig.7 Distribution of heat fluxes calculated with Standardk-εmodel

圖8 不同湍流模型下熱載荷密度分布Fig.8 Distribution of total heat flux calculated with different turbulence models

圖9 Standardk-ε模型下的熱載荷Fig.9 Thermal loads calculated with Standardk-εmodel

由圖8可看出，采用不同湍流模型獲得的熱載荷密度分布有所差異。利用計算得到的各微元段的防冰熱載荷密度乘以相應(yīng)微元段面積并求和，就可得到帽罩的總防冰熱載荷，結(jié)果見表1?？梢姡牧髂Ｐ偷倪x取對計算結(jié)果的影響較大。

由圖9可看出，帽罩表面水膜溫度越高，對防冰加熱熱載荷的需求也越大。故可根據(jù)實際需要來確定最小的熱載荷，在保證防冰的前提下，避免造成熱量的過多消耗和浪費。

5 結(jié)論

(1)湍流模型的選取對防冰熱載荷的計算結(jié)果影響較大。

(2)帽罩表面各項熱流密度中，撞擊水的動能轉(zhuǎn)化成的熱流密度所占比例較小。

(3)帽罩表面水膜溫度越高，所需的防冰熱載荷就越大。

[1]裘燮綱，韓鳳華.飛機防冰系統(tǒng)[M].北京：北京航空教材編審組，1985.

[2]陳維建，張大林.飛機機翼結(jié)冰過程的數(shù)值模擬[J].航空動力學(xué)報，2005，20(6)：112—119.

[3]戰(zhàn)培國.結(jié)冰風(fēng)洞研究綜述[J].實驗流體力學(xué)，2007，21 (3)：92—96.

[4]曾海軍.航空發(fā)動機進氣系統(tǒng)結(jié)冰、冰吸入及結(jié)冰保護的合格審定[J].中國民航飛行學(xué)院學(xué)報，2006，17(3)：16—19.

[5]卜雪琴，林貴平.基于CFD的水收集系數(shù)及防冰表面溫度測試[J].北京航空航天大學(xué)學(xué)報，2007，33(10)：55—58.

[6]朱劍鋆，董葳.防冰部件表面流動換熱與溫度計算分析[J].燃氣渦輪試驗與研究，2011，24(1)：15—18.

[7]楊世銘，陶文銓.傳熱學(xué)[M].3版.北京：高等教育出版社，2006.

Numerical Investigation of Anti-Icing Thermal Loads for Engine Inlet Cowl

LIU Hua，GUO Wen，YANG Jun，LOU De-cang
(China Gas Turbine Establishment，Chengdu 610500，China)

Calculating methods of anti-icing thermal loads of an engine inlet cowl are investigated.FLU?ENT software is used to simulate the flow-field of cowl and the trajectories of water droplets are calculated by employing Discrete Phase Model of FLUENT.The local impingement efficiency of water drops is gained with the trajectories.Then the anti-icing thermal loads of cowl is achieved by analyzing the heat fluxes in?volved in the energy exchange on the surface with the results of numerical investigation and empirical rela?tion formulas.Thermal loads of different turbulence models are also compared.

cowl；local impingement efficiency；anti-icing thermal loads；turbulence models；numerical simulation

V244.1+5

A

1672-2620(2012)01-0044-05

2011-03-07；

2011-12-06

劉華(1986-)，女，江西人，碩士研究生，主要從事航空發(fā)動機空氣系統(tǒng)設(shè)計和零件熱分析研究。