機翼防冰過程中冰脊問題的數值分析

王 超 常士楠 楊 波 黎明中

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

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

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

1 流場與局部水收集系數計算

1.1 流場計算

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

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

圖1 NACA0012翼型網格圖

1.2 對流換熱系數計算

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

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

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

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

當Rek≤600認為流動為層流，層流區域中的對流換熱系數為

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

當Rek＞600認為流動為湍流，流動區域中的對流換熱系數為

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

式中s表示距離前緣駐點的弧長（m）.

1.3 局部水收集系數計算

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

式中，mw為水滴質量（kg）；d為水滴直徑（m）；ρa和ρw分別為空氣和水滴的密度（kg／m3）；g為重力加速度（m／s2）；u為當地的氣流速度（m／s）；uw為水滴的速度（m／s）；Rew＝｜u－uw｜d／νa為水滴相對于氣流運動的雷諾數；Cd為阻力系數，表達式為［10］

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

2 控制體內熱質平衡分析與算法驗證

2.1 質量平衡分析

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

圖2 微元體內質量平衡示意圖

2.2 能量平衡分析

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

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

圖3 微元體內熱流平衡示意圖

2.3 算法驗證

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

表1 結冰計算條件

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

圖4 3種典型的結冰冰形比較

3 冰脊的生長模擬

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

圖5 冰脊隨結冰時間的累積

圖6 冰脊導致水滴軌跡變化的示意圖

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

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

1）環境溫度對冰脊的影響.

環境溫度對機翼表面的水形態有直接的影響，進而對冰脊產生影響.圖7給出了加熱功率恒定時冰脊隨環境溫度的變化情況.

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

從圖7可以看出，當環境溫度從－4.4℃降低到－7.8℃時，在熱防護區域外形成的冰脊寬度有所減小，但冰脊高度顯著增加；當環境溫度從－7.8℃降低到－13.35℃時，機翼上表面熱防護區域外的冰脊高度和寬度減小，但是在駐點處開始有結冰形成；當溫度降低到霜冰條件時，如圖7b所示，相對于駐點處的結冰高度而言，機翼熱防護區域外冰脊的高度相對較小.環境溫度較低時，壁面與氣流之間的溫差增大，從而導致換熱量增加，亦即液態水的凍結速率變大，液態水的溢流速度和溢流量降低，當供給加熱熱流恒定時，在機翼駐點處就出現了顯著的結冰.

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

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

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

3）液態水含量和水滴直徑對冰脊的影響.

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

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

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

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

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

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

5 冰脊對機翼氣動特性的影響分析

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

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

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

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

6 結 論

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

1）對于明冰條件，冰脊主要在熱防護區域的邊緣位置處形成和發展，并隨著結冰時間的累積在寬度和高度方向上迅速增加；而霜冰條件下更容易在熱防護區域內容易出現結冰.

2）環境溫度、飛行速度、液態水含量、水滴直徑和加熱功率對冰脊的形成和發展有不同程度的影響.當飛機進入明冰或混合冰結冰范圍的云層時，可以增加機翼的防冰加熱功率或是加大飛行速度，進而防止冰脊的形成和發展.

3）霜冰條件下，機翼駐點處的結冰對機翼的氣動性能破壞性較大，當飛機進入該類結冰條件的云層時，應當盡量增大防／除冰系統的加熱功率，避免在加熱區域內形成結冰.

（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］李延.表面粗糙度與飛機結冰的相關性研究［D］.北京：北京航空航天大學航空科學與工程學院，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方程的翼型結冰數值模擬［J］.西北工業大學學報，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