王 超 常士楠 楊 波 黎明中

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

在飛機(jī)的防／除冰過程中形成冰脊的現(xiàn)象一直是危害飛行安全的重大隱患之一.1994年10月31日，美國一架ATR-72商用飛機(jī)在雨中失事［1－2］，機(jī)載人員全部遇難；2002年12月，復(fù)興航空公司一架貨機(jī)在澎湖地區(qū)失速墜海；2006年6月，我國一架軍用運(yùn)輸機(jī)因結(jié)冰在安徽失事，機(jī)上40人全部遇難，以上三起事故均為冰脊所引發(fā).形成冰脊的原因主要有兩種：一是由于大尺度過冷水滴撞擊結(jié)冰所致；另一種是對于常規(guī)尺寸的水滴而言，雖然機(jī)翼采取了熱防護(hù)措施，但由于機(jī)翼壁面提供的熱量不足以使結(jié)冰融化后的液態(tài)水完全蒸發(fā)，此時液態(tài)水會在氣流的吹拂下發(fā)生溢流，在機(jī)翼未設(shè)熱防護(hù)的區(qū)域形成冰脊.后一種情況往往更具有隱蔽性，不容易被人們發(fā)現(xiàn)，從而造成重大損失.

自20世紀(jì)90年代以來，研究人員針對冰脊問題開展了大量的研究［3－4］，但研究內(nèi)容多關(guān)注于冰脊對機(jī)翼的氣動方面的影響，而對冰脊的形成特點(diǎn)和影響因素則缺少深入而系統(tǒng)的研究.經(jīng)過10多年的研究，人們還只能夠通過冰風(fēng)洞試驗來得到對應(yīng)氣象條件下冰脊的形狀、位置及大小等信息，精確描述冰脊形成的數(shù)理模型仍未建立，至今還沒有實(shí)現(xiàn)冰脊形成過程的數(shù)值計算［5］.本文針對均勻防護(hù)熱流條件下冰脊的生長特點(diǎn)和影響因素展開研究，通過對帶有冰脊的機(jī)翼的氣動特性進(jìn)行計算分析，研究不同條件下冰脊對機(jī)翼的氣動性能的影響.

1 流場與局部水收集系數(shù)計算

1.1 流場計算

對于二維非定常Navier-Stokes方程，針對某一控制體Ω寫成如下積分形式［6］：

式中，V為微元控制體；U為質(zhì)量、動量和能量各守恒變量；S為翼型的壁面邊界；F為無粘通量；Fv為粘性通量；n為控制體的外法線向量；Re為翼型的特征雷諾數(shù).對于該控制方程的求解借助流體動力力學(xué)軟件Fluent完成.圍繞翼型的結(jié)構(gòu)網(wǎng)格如圖1所示，計算域網(wǎng)格節(jié)點(diǎn)總數(shù)為107 430個，網(wǎng)格單元數(shù)是106 444個，對機(jī)翼的附面層網(wǎng)格進(jìn)行了加密處理.

圖1 NACA0012翼型網(wǎng)格圖

1.2 對流換熱系數(shù)計算

實(shí)際的結(jié)冰表面并不是平整光滑的，而是存在大小不等的粗糙度.對于結(jié)冰表面粗糙度的計算，本文采用等效沙粒粗糙度模型［7］.等效沙粒粗糙度高度ks的計算公式為

式中，QLWC為液態(tài)水含量（g／m3）；T∞為環(huán)境溫度（K）；u∞為來流速度（m／s）；c為翼型的弦長（m）.

結(jié)冰翼型壁面的流動非常復(fù)雜，流動存在一個從層流到湍流變化的復(fù)雜過程，轉(zhuǎn)捩點(diǎn)位置的判定對流動類型的區(qū)分很重要.為此，引入粗糙度雷諾數(shù)的概念［8］：

式中，Rek為粗糙度雷諾數(shù)；νa為空氣的運(yùn)動粘度（m2／s）；uk為 對 應(yīng) 粗 糙 度ks處 的 空 氣 速度（m／s）.

當(dāng)Rek≤600認(rèn)為流動為層流，層流區(qū)域中的對流換熱系數(shù)為

式中，λ為空氣的導(dǎo)熱系數(shù)（W·m－1·k－1）；ue為邊界層外邊界速度（m／s）.

當(dāng)Rek＞600認(rèn)為流動為湍流，流動區(qū)域中的對流換熱系數(shù)為

式中，Cpa為空氣的定壓比熱容（J·kg－1·K－1）；Cf為壁面摩擦系數(shù)，其數(shù)值可以通過經(jīng)驗公式得到［9］：

式中s表示距離前緣駐點(diǎn)的弧長（m）.

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

局部水收集系數(shù)是表征機(jī)翼表面收集水滴率的參數(shù)，對于水滴撞擊結(jié)冰計算具有重要的意義.本文采用拉格朗日軌跡追蹤法求解水滴的運(yùn)動軌跡，進(jìn)而確定機(jī)翼壁面上局部收集系數(shù)的分布.水滴運(yùn)動的軌跡方程為

式中，mw為水滴質(zhì)量（kg）；d為水滴直徑（m）；ρa(bǔ)和ρw分別為空氣和水滴的密度（kg／m3）；g為重力加速度（m／s2）；u為當(dāng)?shù)氐臍饬魉俣龋╩／s）；uw為水滴的速度（m／s）；Rew＝｜u－uw｜d／νa為水滴相對于氣流運(yùn)動的雷諾數(shù)；Cd為阻力系數(shù)，表達(dá)式為［10］

本文采用四階Runge-Kutta法對式（7）進(jìn)行求解，得到各個水滴的軌跡后就可以計算翼型壁面各控制體的局部水收集系數(shù).

2 控制體內(nèi)熱質(zhì)平衡分析與算法驗證

2.1 質(zhì)量平衡分析

微元控制體內(nèi)質(zhì)量的流入與流出如圖2所示，m（i）cap為 微元體 收集到 的水量；m（i－1）in為從上一控制體流入的水量；m（i）e為蒸發(fā)或升華到空間的水量；m（i）f為凍結(jié)成冰的水量，當(dāng)防冰熱載荷不足時，該項不為0，需要加以考慮；m（i）out為流向下一控制體的水量；上述質(zhì)量單位為kg／（m2·s）.

圖2 微元體內(nèi)質(zhì)量平衡示意圖

2.2 能量平衡分析

微元控制體內(nèi)熱量的流入與流出如圖3所示，q（i）wv為水滴撞擊到結(jié)冰表面時的動能轉(zhuǎn)化成的比熱流；q（i）cov和q（i）v分別為由于對流散失到空間的比熱流和附面層摩擦產(chǎn)生的比熱流；q（i－1）in和q（i）out分別為流入和流出控制體的液態(tài)水所攜帶的比熱流；q（i）w為加熱收集水滴的比熱流；q（i）f和q（i）e分別為液態(tài)水結(jié)冰時釋放和蒸發(fā)所需的比熱流.以上熱量單位均為W／m2，各項熱流的計算表達(dá)式可參閱文獻(xiàn)［10］，此處不再詳述.容易求得微元體內(nèi)的熱量平衡關(guān)系式為

當(dāng)采用熱氣進(jìn)行防冰時，由于防冰腔結(jié)構(gòu)以及防冰系統(tǒng)內(nèi)部傳熱的影響，翼型表面不同位置所獲得的供熱熱流是不同的，但是對于電熱防冰系統(tǒng)，加熱熱流可控，可以實(shí)現(xiàn)對不同防冰位置處供給相同的加熱熱流.因此，假定為翼型壁面的熱防護(hù)區(qū)域內(nèi)的每個微元控制體提供一恒定且相同的熱流密度q（i）s，如圖3所示.對于熱防護(hù)區(qū)域之外的微元控制體，則認(rèn)為q（i）s＝0.此條件類似于在防護(hù)區(qū)域內(nèi)布置連續(xù)分布的電加熱層，且各處的加熱熱流密度近似相等的情況.根據(jù)以上分析，微元體內(nèi)的熱量平衡關(guān)系式改寫為

圖3 微元體內(nèi)熱流平衡示意圖

2.3 算法驗證

由于冰脊的形成和生長是在機(jī)翼的整個防冰過程進(jìn)行的，因此，為更準(zhǔn)確描述冰脊的特征，開發(fā)了多步結(jié)冰程序.為便于與文獻(xiàn)［11］中的冰風(fēng)洞試驗結(jié)果進(jìn)行對比，進(jìn)而驗證本文算法的正確性，選取NACA0012翼型為計算對象，分別模擬了環(huán)境溫度為－6.1，－10，－19.4℃的3種典型冰形，具體計算條件列于表1中.

表1 結(jié)冰計算條件

在所選取計算條件下的結(jié)冰計算結(jié)果與文獻(xiàn)中的試驗結(jié)果的對比如圖4所示.從圖4a和圖4b可以看出，本文的計算結(jié)果在結(jié)冰厚度及冰角的生長方向和大小等方面都與試驗結(jié)果取得了較好的吻合，對于霜冰冰形而言，如圖4c所示，除在結(jié)冰下表面的局部位置計算結(jié)果略小于試驗冰形外，其余大部分計算結(jié)果與試驗結(jié)果幾乎重合，因此，利用本文發(fā)展的計算程序可以進(jìn)行復(fù)雜條件下的熱質(zhì)耦合計算.

圖4 3種典型的結(jié)冰冰形比較

3 冰脊的生長模擬

翼型表面熱防護(hù)范圍由直徑為40μm的過冷水滴的撞擊范圍確定，環(huán)境溫度分別為－13.35℃和－19.4℃，其他計算條件同表1.圖5以橫坐標(biāo)為距駐點(diǎn)弧長、縱坐標(biāo)為結(jié)冰厚度的形式，展開放大給出了兩種結(jié)冰溫度條件下冰脊隨著結(jié)冰時間的累積情況.當(dāng)環(huán)境溫度為－13.35℃時，如圖5a所示，翼型上、下表面的冰脊主要是在貼近熱防護(hù)極限的位置處開始生長；隨著時間增加，冰脊在寬度和高度方向上均有增長，且在高度方向上的增長較為明顯；就冰脊寬度方向上的增長而言，翼型下表面的冰脊向著駐點(diǎn)方向增加，而機(jī)翼上表面冰脊增長不明顯；當(dāng)環(huán)境溫度為－19.4℃時，如圖5b所示，在機(jī)翼駐點(diǎn)處出現(xiàn)了顯著的結(jié)冰，直到結(jié)冰時間為240s，在機(jī)翼下表面熱防護(hù)極限外開始有冰脊出現(xiàn)，機(jī)翼上表面的冰脊形成不明顯.綜合兩種情況下冰脊的生長特點(diǎn)可看出，當(dāng)結(jié)冰時間大于240s時，冰脊在高度方向的增長速度顯著減小.冰脊的生長，開始是由于溢流水量隨著時間不斷累積，而等到冰脊增長到一定高度后，氣流中的水滴則可直接撞擊到冰脊表面，如圖6所示，從而造成冰脊在高度上的不斷增加.

圖5 冰脊隨結(jié)冰時間的累積

圖6 冰脊導(dǎo)致水滴軌跡變化的示意圖

4 冰脊生長的影響因素分析

影響飛機(jī)結(jié)冰的因素，如環(huán)境溫度、氣流速度等也必然會對冰脊產(chǎn)生影響，探求不同的結(jié)冰條件下冰脊的特點(diǎn)，對避免冰脊的產(chǎn)生、保障飛機(jī)性能和飛行安全具有重要的意義.本節(jié)除所要研究的結(jié)冰條件外，其他計算條件同表1.

1）環(huán)境溫度對冰脊的影響.

環(huán)境溫度對機(jī)翼表面的水形態(tài)有直接的影響，進(jìn)而對冰脊產(chǎn)生影響.圖7給出了加熱功率恒定時冰脊隨環(huán)境溫度的變化情況.

圖7 不同環(huán)境溫度對冰脊的影響（q（i）s＝20kW／m2）

從圖7可以看出，當(dāng)環(huán)境溫度從－4.4℃降低到－7.8℃時，在熱防護(hù)區(qū)域外形成的冰脊寬度有所減小，但冰脊高度顯著增加；當(dāng)環(huán)境溫度從－7.8℃降低到－13.35℃時，機(jī)翼上表面熱防護(hù)區(qū)域外的冰脊高度和寬度減小，但是在駐點(diǎn)處開始有結(jié)冰形成；當(dāng)溫度降低到霜冰條件時，如圖7b所示，相對于駐點(diǎn)處的結(jié)冰高度而言，機(jī)翼熱防護(hù)區(qū)域外冰脊的高度相對較小.環(huán)境溫度較低時，壁面與氣流之間的溫差增大，從而導(dǎo)致?lián)Q熱量增加，亦即液態(tài)水的凍結(jié)速率變大，液態(tài)水的溢流速度和溢流量降低，當(dāng)供給加熱熱流恒定時，在機(jī)翼駐點(diǎn)處就出現(xiàn)了顯著的結(jié)冰.

2）飛行速度對冰脊的影響.

圖8給出了不同飛行速度條件下冰脊的特征，從圖中可以看出，當(dāng)飛行速度在一定范圍內(nèi)增加時（如Ma＝0.141～0.304），翼型上下表面的冰脊體積會增加，但是當(dāng)速度繼續(xù)增加到某一數(shù)值時（如Ma＝0.416），翼型下表面的冰脊消失，只在上表面有少量存在，但不足以對機(jī)翼的氣動性能造成較大的危害.這是因為飛行速度的增加加速了結(jié)冰表面液態(tài)水蒸發(fā)，同時氣動加熱量也增加了，并最終使得冰脊體積減小甚至消失.

圖8 Ma對冰脊的影響（t＝－13.35℃，q（i）s＝20kW／m2）

3）液態(tài)水含量和水滴直徑對冰脊的影響.

圖9和圖10分別給出了不同液態(tài)水含量和水滴直徑條件下冰脊的變化情況.液態(tài)水含量和水滴直徑對冰脊的影響規(guī)律相同，即隨著兩者的增大，翼型上、下壁面處冰脊的體積也迅速增加.液態(tài)水含量和水滴直徑增加時，水滴的撞擊范圍和收集系數(shù)增加，同時，結(jié)冰表面的溢流水量也增加，從而導(dǎo)致冰脊的迅速增長.

圖9 QLWC對冰脊的影響（t＝－13.35℃，q（i）s＝20kW／m2）

4）加熱功率對冰脊的影響.

圖10 dMV對冰脊的影響（t＝－13.35℃，q（i）s＝20kW／m2）

從圖11中不同加熱功率對冰脊的影響可以看出，隨著加熱功率的增加，駐點(diǎn)處的結(jié)冰范圍和厚度迅速減小，直至消失；當(dāng)加熱功率從5kW／m2增加到10kW／m2時，機(jī)翼下表面的冰脊體積急劇增加，上表面未出現(xiàn)冰脊；當(dāng)加熱功率從10kW／m2增加到20kW／m2時，機(jī)翼下表面的冰脊體積急劇減小，而機(jī)翼上表面開始有冰脊出現(xiàn)；當(dāng)加熱功率繼續(xù)增加到30kW／m2時，機(jī)翼下表面的冰脊消失.冰脊隨著加熱功率增加出現(xiàn)上述變化的原因源于結(jié)冰表面液態(tài)水的變化.結(jié)冰表面的液態(tài)水量主要取決于結(jié)冰融化的水量、蒸發(fā)量和收集水量，由于翼型收集的水量恒定，結(jié)冰融化的水量和蒸發(fā)量都隨加熱功率的增加而增加.當(dāng)加熱功率從5kW／m2增加到10kW／m2時，單位時間內(nèi)結(jié)冰融化的水量大于液態(tài)水的蒸發(fā)量，因此，在氣流的吹拂作用下，熱防護(hù)區(qū)域外的溢流量增加，導(dǎo)致冰脊體積變大；當(dāng)加熱功率增加到20kW／m2時，單位時間內(nèi)積冰融化的水量小于液態(tài)水的蒸發(fā)量，溢流水量減少，從而導(dǎo)致冰脊體積減??；隨著加熱功率的繼續(xù)增加，機(jī)翼表面的液態(tài)水量不斷降低，并最終使冰脊消失.

圖11 加熱功率對冰脊的影響（t＝－13.35℃，q（i）s＝20kW／m2）

5 冰脊對機(jī)翼氣動特性的影響分析

通過以上對冰脊的計算分析可知，冰脊形成的位置和大小受諸多因素的影響，表現(xiàn)出較大的不確定性.基于上述研究，選取4種典型溫度及4種不同的防護(hù)熱流條件下形成的冰脊，計算冰脊對機(jī)翼氣動特性的影響.圖12給出了加熱功率為20kW／m2時，不同環(huán)境溫度條件下冰脊對機(jī)翼升力系數(shù)的影響，對比看出，霜冰條件下形成的冰脊對機(jī)翼的氣動性能具有較大破壞性，結(jié)合圖7b也可以發(fā)現(xiàn)，在機(jī)翼的氣動駐點(diǎn)區(qū)內(nèi)有形狀很不規(guī)則結(jié)冰形成，嚴(yán)重破壞了機(jī)翼的氣動外形.

圖12 冰脊對機(jī)翼氣動特性的影響（q（i）s＝20kW／m2）

結(jié)合圖11和圖13可以看出，溫度為－13.35℃時，隨著加熱功率的增大，冰脊尺度變小，其對機(jī)翼氣動性能的破壞性也顯著降低.

圖13 冰脊對機(jī)翼氣動特性的影響（t＝－13.35℃）

6 結(jié) 論

通過對不同條件下的機(jī)翼表面冰脊的生長過程和影響因素以及帶有冰脊的機(jī)翼的氣動特性的計算分析，得出以下主要結(jié)論：

1）對于明冰條件，冰脊主要在熱防護(hù)區(qū)域的邊緣位置處形成和發(fā)展，并隨著結(jié)冰時間的累積在寬度和高度方向上迅速增加；而霜冰條件下更容易在熱防護(hù)區(qū)域內(nèi)容易出現(xiàn)結(jié)冰.

2）環(huán)境溫度、飛行速度、液態(tài)水含量、水滴直徑和加熱功率對冰脊的形成和發(fā)展有不同程度的影響.當(dāng)飛機(jī)進(jìn)入明冰或混合冰結(jié)冰范圍的云層時，可以增加機(jī)翼的防冰加熱功率或是加大飛行速度，進(jìn)而防止冰脊的形成和發(fā)展.

3）霜冰條件下，機(jī)翼駐點(diǎn)處的結(jié)冰對機(jī)翼的氣動性能破壞性較大，當(dāng)飛機(jī)進(jìn)入該類結(jié)冰條件的云層時，應(yīng)當(dāng)盡量增大防／除冰系統(tǒng)的加熱功率，避免在加熱區(qū)域內(nèi)形成結(jié)冰.

（References）

［1］Pereira C M.Status of NTSB aircraft icing certification-related safety recommendations issued as a result of the 1994 ATR-72accident at Roselawn，IN ［R］.AIAA 1997-0410，1997

［2］Bond T H，Miller D R，Potapczuk M G.Overview of SLD engineering tools development［R］.AIAA 2003-386，2003

［3］Bragg M B.Aircraft aerodynamic effects due to large droplet Ice Accretions［R］.AIAA 96-0932，1996

［4］李延.表面粗糙度與飛機(jī)結(jié)冰的相關(guān)性研究［D］.北京：北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，2011

Li Yan.Study on relativity between icing surface roughness and aircraft icing［D］.Beijing：School of Aeronautic Science and Engineering，Beijing University of Aeronautics and Astronautics，2011（in Chinese）

［5］Alegre N.Full-scale runback ice：accretion and aerodynamic study［D］.Cranfield University，2010

［6］蔣勝矩，李鳳蔚.基于N-S方程的翼型結(jié)冰數(shù)值模擬［J］.西北工業(yè)大學(xué)學(xué)報，2004，22（5）：559－560

Jiang Shengju，Li Fengwei.A better numerical simulation of ice accretion on airfoil［J］.Journal of Northwestern Polytechnical University，2004，22（5）：559－560（in Chinese）

［7］Shin J，Berkowitz B，Chen H，et al.Prediction of ice shapes and their effect on airfoil performance［R］.AIAA 91-0264，1991

［8］von Donehoff A，Horton E.A low-speed experimental investigation of the effect of a sandpaper type of roughness on boundary-layer transition［R］.NACA Report No.1349，1956

［9］Beaugendre H，Morency F，Habashi W G.Roughness implementation in fensap-ice：model calibration and influence on ice shapes［J］.Journal of Aircraft，2003，40（6）：1212－1215

［10］Gent R W，Dart N P，Cansdale J T.Aircraft icing［J］.Philosophical Transactions：Mathematical，Physical and Engineering Sciences，2000，358（1776）：2873－2911

［11］Shin J，Bond T H.Results of an icing test on a NACA 0012 airfoil in the NASA Lewis icing research tunnel［R］.AIAA-92-0647，1992