李焱鑫，張 辰，劉 洪，王福新

（上海交通大學(xué) 航空航天學(xué)院，上海 200240）

0 引 言

2002年12月21日，臺灣復(fù)興航空一架ATR－72貨機(jī)，在飛行中遭遇嚴(yán)重結(jié)冰天氣，飛機(jī)在積冰后升阻特性大幅下降，最終墜毀。事故調(diào)查報告［1］表明，失事飛機(jī)的飛行環(huán)境中包含大粒徑過冷水滴（SLD）。SLD的粒徑大于50mm，超出了適航條例25部附錄C的規(guī)定。飛機(jī)在遭遇SLD結(jié)冰環(huán)境時，過冷水滴在機(jī)翼表面會繼續(xù)向后方流動，并在飛機(jī)除冰裝置后方凍結(jié)。當(dāng)除冰裝置開啟時，機(jī)翼前緣受保護(hù)區(qū)的積冰被清除，機(jī)翼表面在未保護(hù)區(qū)形成脊?fàn)畹囊缌鞣e冰［2－3］。

鑒于SLD結(jié)冰嚴(yán)重的危害性，2005年美國聯(lián)邦航空局（FAA）在工作報告［4］中，將SLD結(jié)冰列為重點研究方向，并在適航咨詢通告［5］中特別提出，需要對SLD引起的特殊冰型進(jìn)行審定。Broeren等［6］通過實驗比較了多種積冰外形的氣動性能損失，發(fā)現(xiàn)在實驗條件下脊?fàn)畋挂硇偷淖畲笊ο禂?shù)由1.85下降到0.52，失速迎角由18.1°下降至5.6°，比流向冰和角狀冰對升阻的破壞都要大。Edward等［7］進(jìn)一步分析指出，SLD結(jié)冰的危害性在于，溢流冰脊使翼型流場提前分離。

脊?fàn)畋姆蛛x效應(yīng)與翼型幾何外形密切相關(guān)。Lee［8］分析了脊?fàn)畋鶎Σ煌硇蜌鈩有阅艿挠绊?，NACA 23012m和NLF 0414的最大升力系數(shù)（CLmax）分別下降了83.7%和54.8%，這表明脊?fàn)畋挠绊懪c翼型形狀有關(guān)。超臨界翼型廣泛應(yīng)用于大型客機(jī)的設(shè)計中，但前人并未對此進(jìn)行單獨分析。Bragg等［2，9］指出脊?fàn)畋L在翼型的主要升力產(chǎn)生區(qū)，使翼型的最大升力系數(shù)下降。超臨界翼型的主要升力產(chǎn)生區(qū)，與溢流結(jié)冰位置有較大重合，易形成吸力峰引發(fā)氣動分離，因此本文著重分析了溢流積冰對超臨界翼型的影響。

目前研究結(jié)冰影響的方法主要為風(fēng)洞實驗，但風(fēng)洞實驗數(shù)據(jù)獲取不全面，具有一定局限性。數(shù)值分析方法具有高效、全面的特點，得到了廣泛應(yīng)用。在翼型結(jié)冰后的數(shù)值氣動分析領(lǐng)域，陳科等［10］應(yīng)用SST模型對機(jī)翼前緣的復(fù)雜積冰引起的流動分離進(jìn)行了分析；Mortensen［11］比較了多種湍流模型，認(rèn)為 S－A模型用于評價前緣積冰的流動分離具有較高精度。上述兩位學(xué)者均未對溢流脊?fàn)畋M(jìn)行分析。Dunn等［12］使用NSU2D程序分析了帶脊?fàn)畋腘ACA 23012m的升力曲線，計算結(jié)果在大迎角條件下誤差較大，不能準(zhǔn)確預(yù)測出帶冰翼型的CLmax和失速迎角αstall。綜合分析前人工作，本文發(fā)展了一種更適合溢流結(jié)冰氣動分析的數(shù)值方法。

我國面臨大型民機(jī)設(shè)計的關(guān)鍵階段，為保證安全性，迫切需要分析SLD溢流結(jié)冰對機(jī)翼氣動性能的影響，為我國的民用飛機(jī)設(shè)計和適航審定提供指導(dǎo)。為此，本文采用數(shù)值計算的方法，對溢流冰脊的氣動性能進(jìn)行了分析。本文介紹了驗證算例的計算條件，確定了用于溢流冰脊計算的湍流模型，在此基礎(chǔ)上驗證了計算網(wǎng)格的有效性；將不同冰脊高度的計算結(jié)果與實驗結(jié)果對比，確定了方法的有效性；分析了機(jī)翼超臨界翼型和平尾翼型的流場特征；分析了SLD溢流結(jié)冰條件下，大型客機(jī)機(jī)翼超臨界翼型和平尾翼型各項氣動系數(shù)的變化。

1 研究方法

1.1 計算模型與網(wǎng)格

本節(jié)以Lee［8］的實驗數(shù)據(jù)和實驗條件為依據(jù)，在CFD軟件中對帶脊?fàn)畋腘ACA 23012m的升力曲線進(jìn)行計算，并與實驗結(jié)果對比。由于翼型氣動性能受脊?fàn)畋叨茸兓挠绊戄^大［13］，本文在確定數(shù)值方法時選擇了不同的冰脊高度k／c。SLD結(jié)冰時具有溢流積冰的特點，脊?fàn)畋恢脁／c選定為機(jī)翼保護(hù)區(qū)后方x／c＝0.1。冰脊形狀為前端面半圓，Lee和Bragg［14］的實驗表明此形狀的冰脊對翼型的氣動破壞作用最明顯。計算條件參見表1。

表1 算例計算條件Table 1 Computational conditions of the case

1.1.1 湍流模型

本文基于N－S方程求解帶冰翼型流場，脊?fàn)畋鶗l(fā)翼型流場的氣動分離，需要選擇合適的湍流模型對流場進(jìn)行分析。Mortensen［11］比較了不同湍流模型在前緣結(jié)冰算例計算中的精度，包括S－A模型、重整化 RNG k－ε 模型、可實現(xiàn)性 Realizable k－ε 模型、標(biāo)準(zhǔn)k－ε模型，文獻(xiàn)計算結(jié)果表明S－A模型在以上幾種湍流模型中最優(yōu)。雷諾應(yīng)力模型［15］（RSM）并沒有采用上述幾種模型的渦粘各向同性假設(shè)，而是直接求解雷諾應(yīng)力分量，理論上應(yīng)更為精確，因此本文采用RSM模型作為計算用湍流模型。

1.1.2 網(wǎng)格驗證

考慮到帶冰翼型表面較復(fù)雜本文生成非結(jié)構(gòu)網(wǎng)格。為驗證網(wǎng)格的無關(guān)性，生成尺度不同的三套網(wǎng)格，并對網(wǎng)格進(jìn)行自適應(yīng)加密驗證。

在本文的模擬中，網(wǎng)格區(qū)域為20倍弦長的圓，翼型位于區(qū)域中心。為捕捉近壁面流場特征，對初始近壁面網(wǎng)格進(jìn)行了加密處理，在初步計算時，選擇了不同的網(wǎng)格加密尺度，見表2。

表2 網(wǎng)格信息Table 2 Grids details

圖1是不同網(wǎng)格尺度下的計算升力曲線與實驗結(jié)果的對比。在圖中，三種網(wǎng)格在低角度下均與實驗結(jié)果符合較好，但在接近失速迎角時，加密后的網(wǎng)格3與實驗結(jié)果符合較好。

圖1 不同網(wǎng)格計算結(jié)果對比Fig.1 Comparison of results on different grids

在網(wǎng)格3基礎(chǔ)上進(jìn)行初步計算，得到帶冰翼型流場的基本信息，按照速度梯度的變化對原始網(wǎng)格進(jìn)行自適應(yīng)加密，得到自適應(yīng)網(wǎng)格。此時網(wǎng)格區(qū)域含185259個非結(jié)構(gòu)網(wǎng)格，共107860個網(wǎng)格節(jié)點，如圖2。從圖3中可見，在使用自適應(yīng)網(wǎng)格對帶冰原始網(wǎng)格進(jìn)行加密后，自適應(yīng)網(wǎng)格對失速迎角的計算結(jié)果與網(wǎng)格3一致，在不同冰脊高度k／c＝0.0056、0.0083、0.0137時，兩套網(wǎng)格在最大升力系數(shù)計算上的誤差分別為0.02、0.01和0.01。因此，使用網(wǎng)格3對溢流結(jié)冰進(jìn)行計算時，已具有相當(dāng)精度，進(jìn)一步進(jìn)行自適應(yīng)加密后，計算結(jié)果變化很小，網(wǎng)格無關(guān)性可得到驗證。

圖2 計算網(wǎng)格Fig.2 Computational grids

圖3 自適應(yīng)網(wǎng)格計算結(jié)果對比Fig.3 Comparison between original grids and adaption grids

1.2 算例計算結(jié)果

本文應(yīng)用結(jié)合RSM模型的計算方法，將計算結(jié)果與實驗進(jìn)行對比，并分析計算誤差。如圖4（a），當(dāng)k／c＝0時，CLmax的計算值比實驗值偏大約0.05，失速角的計算值比實驗值高約0.6°；如圖4（b），當(dāng)k／c＝0.0056時，CLmax的計算最大誤差值發(fā)生在－4°，為0.13；如圖4（c），當(dāng)k／c＝0.0083時，CLmax計算最大誤差值發(fā)生在4°，為0.06；如圖4（d），當(dāng)k／c＝0.0137時，CLmax計算最大誤差值發(fā)生在－4°，為0.09（圖中右側(cè)坐標(biāo)為升力系數(shù)計算絕對誤差）。冰脊高度的增長，未使計算誤差產(chǎn)生明顯增加。

在圖4中RSM模型計算所得最大升力系數(shù)和失速迎角與實驗結(jié)果符合較好。表3給出了RSM模型的最大升力系數(shù)、失速迎角的計算結(jié)果對比。綜合分析發(fā)現(xiàn)，RSM模型計算的最大升力系數(shù)相對誤差不超過5.8%，不同迎角下升力系數(shù)的絕對誤差小于0.2；除k／c＝0.0056外，失速角的預(yù)測誤差小于1°?？烧J(rèn)為基于雷諾應(yīng)力的流場計算方法，能夠應(yīng)用于帶脊?fàn)畋囊硇蜌鈩有阅軗p失評估。

表3 最大升力系數(shù)和失速迎角結(jié)果對比Table 3 Comparison between the computation and the experiment on CLmaxandαstall

圖4 計算結(jié)果與實驗結(jié)果對比Fig.4 Comparison between the computation and the experiment results

2 大型客機(jī)翼型結(jié)冰影響分析

本節(jié)采用上文確定的研究方法，開展大型客機(jī)機(jī)翼超臨界翼型和平尾翼型的SLD溢流結(jié)冰氣動分析。由于機(jī)翼翼型沿展向的變化較大，本文截取了機(jī)翼和平尾在展向不同區(qū)域的四個典型翼型截面作為代表性翼型進(jìn)行帶冰氣動力分析，分別命名為AF－1、AF－2、AF－3和 AF－T，如圖5。

圖5 飛機(jī)翼型及冰脊形狀Fig.5 Aircraft airfoils and ridge ice shape

Bragg指出［2］，雷諾數(shù)對帶冰翼型最大升力系數(shù)的影響非常有限?？紤]到結(jié)冰多發(fā)于起降階段，來流條件的設(shè)置參考民用客機(jī)實際的飛行環(huán)境，脊?fàn)畋恢脁／c＝0.1，見表4。

表4 翼型計算條件Table 4 Computational conditions of airfoil

2.1 流場特征分析

2.1.1 機(jī)翼超臨界翼型流場特征

如圖6，以翼型AF－3為例，對比超臨界翼型未結(jié)冰和帶脊?fàn)畋鶅煞N情況在迎角14°時的流線：當(dāng)AF－3未結(jié)冰時，流動未發(fā)生分離；當(dāng) AF－3帶有k／c＝0.0137的脊?fàn)畋鶗r，冰脊導(dǎo)致流動產(chǎn)生了分離，分離區(qū)自冰脊處開始生長，在x／c＝0.4處再附著；AF－1／2有類似特征。

圖6 AF－3流線圖（α＝14°）Fig.6 Streamlines for AF－3（α＝14°）

圖7給出了翼型AF－3對應(yīng)的表面壓力分布：在x／c＝0～0.1范圍，帶冰翼型上表面壓力系數(shù)（負(fù)值）大于干凈翼型；在x／c＝0.1處冰脊使上表面壓力系數(shù)（負(fù)值）減??；在冰脊后x／c＝0.1～0.4的范圍，上表面壓力系數(shù)（負(fù)值）逐漸增大；在x／c＝0.4～1.0范圍內(nèi)，壓力系數(shù)變化很小。

圖7 AF－3表面壓力分布對比（α＝14°，k／c＝0.0137）Fig.7 Surface pressure distributions for AF－3（α＝14°，k／c＝0.0137）

綜合分析，冰脊導(dǎo)致機(jī)翼超臨界翼型表面發(fā)生流動分離。從壓力分布來看，與干凈翼型相比，帶冰翼型上表面壓力系數(shù)較大而下表面較小，上下表面壓差小于干凈翼型。

2.1.2 平尾翼型流場特征

由于平尾翼型與超臨界翼型的幾何結(jié)構(gòu)不同，受脊?fàn)畋绊懹兴町?。如圖8，以翼型AF－T為例，對比平尾未結(jié)冰與帶脊?fàn)畋鶅煞N情況在14°時的流線特征：圖8a）中，干凈翼型在該條件下完全分離，已發(fā)生失速；但當(dāng) AF－T帶有k／c＝0.0137的脊?fàn)畋鶗r，雖然冰脊同樣引發(fā)了分離，但對整個流場而言，冰脊反而起到了整流的作用，分離區(qū)明顯減小。

圖8 AF－T流線圖（α＝14°）Fig.8 Streamlines for AF－T（α＝14°）

圖9給出了翼型對應(yīng)的表面壓力分布：在x／c＝0～0.25范圍內(nèi)，帶冰翼型上表面的壓力系數(shù)（負(fù)值）小于干凈翼型；在x／c＝0.25～1.0范圍內(nèi)，帶冰翼型上表面壓力系數(shù)持續(xù)增大，而干凈翼型基本不變。

圖9 AF－T表面壓力分布對比（α＝14°，k／c＝0.0137）Fig.9 Surface pressure distributions for AF－T（α＝14°，k／c＝0.0137）

綜合分析，對于平尾翼型，在大迎角下干凈翼型已發(fā)生分離，但在結(jié)冰翼型中由于冰脊的存在，使分離區(qū)明顯減小。從壓力分布看，在x／c＝0～0.25范圍內(nèi)，帶冰翼型上下表面壓差大于干凈翼型，x／c＝0.25～1.0兩者壓差相若。

2.2 氣動參數(shù)分析與討論

2.2.1 升力曲線分析

從結(jié)冰翼型的流場特征來看，超臨界翼型與平尾翼型受溢流結(jié)冰的影響有所不同。為進(jìn)一步評估結(jié)冰的氣動影響，對翼型結(jié)冰后的升力系數(shù)的變化進(jìn)行分析。

圖10給出了大客翼型在不同冰脊高度下的的升力曲線。從變化趨勢來看，脊?fàn)畋叨萲／c＝0.0137時，AF－1／2／3 和 AF－T的CLmax分別下降65.9%、56.1%、57.7% 和 38.5%。AF－1、AF－2、AF－3的失速迎角分別提前8.2°、5.7°和6.6°，而 AF－T的失速迎角延后2.9°。

圖10 翼型升力曲線Fig.10 Lift coefficient curve of airfoils

圖11給出了大型客機(jī)各段翼型不同冰脊高度的CLmax變化圖，由該圖分析脊?fàn)畋叨葘Lmax的影響：AF－1在積冰高度k／c由0增長到0.0056時CLmax從1.76 減小到0.97，k／c 達(dá)到0.0137 時CLmax減小到0.60，主要的升力破壞發(fā)生在結(jié)冰初始階段；翼型 AF－2、AF－3與 AF－1保持著類似的趨勢；AF－T在積冰高度k／c由0增長到0.0056時CLmax從1.19減小到0.85，k／c達(dá)到0.0137時CLmax減小到0.72，升力降幅較前三類翼型相對平穩(wěn)。

圖12給出了飛機(jī)各段翼型不同冰脊高度的失速迎角變化圖，由該圖分析脊?fàn)畋叨葘κ儆堑挠绊懀阂硇虯F－1～AF－3的失速迎角隨冰脊高度增加而提前，在k／c達(dá)到0.0056時，失速迎角降幅不足2°，隨著k／c達(dá)到0.0137，失速迎角產(chǎn)生4°～6°不等的降幅，失速角前移主要發(fā)生在冰脊形成階段；翼型AF－T的失速迎角隨冰脊高度的增加而延后，在k／c達(dá)到0.0137時，失速角比未結(jié)冰時延后2.9°，且失速角延后幅度基本與冰脊高度增加呈線性關(guān)系。

圖11 最大升力系數(shù)對比Fig.11 Comparison of maximum lift coefficient

圖12 失速迎角對比Fig.12 Comparison of stalling angles

2.2.2 阻力曲線分析

從結(jié)冰事故報告［1］分析，在飛機(jī)發(fā)生結(jié)冰后，飛機(jī)受到的阻力顯著增加，同時原本起穩(wěn)定作用的縱向操縱參數(shù)將會減小甚至反號，因此有必要對結(jié)冰后超臨界翼型和平尾翼型阻力曲線和力矩曲線的變化進(jìn)行分析。

圖13是翼型阻力曲線圖，分別以 AF－3、AF－T為例，分析超臨界翼型和平尾翼型結(jié)冰后阻力曲線的變化。在機(jī)翼超臨界翼型和平尾翼型中，溢流積冰使翼型的阻力系數(shù)大幅增加，增長幅度隨迎角的增加而變大。

2.2.3 力矩曲線分析

圖14是翼型力矩曲線圖，分別以 AF－3、AF－T為例，分析超臨界翼型和平尾翼型結(jié)冰后力矩曲線的變化。在AF－3翼型中，結(jié)冰后力矩系數(shù)下降，隨溢流冰高度增加力矩非線性變化趨勢增加。在AF－T翼型中，結(jié)冰后力矩系數(shù)的變化顯著小于超臨界翼型，結(jié)冰后力矩曲線存在線性變化趨勢。

圖13 翼型阻力曲線Fig.13 Drag coefficient curve of airfoils

圖14 翼型力矩曲線Fig.14 Pitching moment curve of airfoils

3 結(jié) 論

（1）RSM模型適合于分析溢流結(jié)冰的流動分離，可準(zhǔn)確預(yù)測帶冰翼型的最大升力系數(shù)與失速迎角。

（2）超臨界翼型受到SLD溢流結(jié)冰影響時，翼型表面壓差減小，流場提前發(fā)生分離，最大升力系數(shù)和失速迎角顯著下降；平尾翼型受溢流結(jié)冰影響較小，在大迎角下帶冰翼型流場產(chǎn)生的分離區(qū)小于干凈翼型，失速迎角提前。超臨界翼型和平尾翼型結(jié)冰后阻力系數(shù)顯著增加。結(jié)冰后超臨界翼型的力矩系數(shù)顯著變化，平尾翼型的力矩變化較小。

（3）本文所做SLD溢流結(jié)冰對超臨界翼型和平尾翼型氣動影響的分析，對大型客機(jī)的設(shè)計和適航審定具有一定指導(dǎo)意義。

［1］FAA.Occurrence investigation report：in－flight icing encounter and crash into the sea，Transasia airways flight 791［R］.ASC－AOR－05－04－001，2005.

［2］BRAGG M B，BROEREN A P，BLUMENTHAL L A.Icedairfoil aerodynamics［J］.Progress in Aerospace Sciences，2005，41（5）：323－362.doi：10.1016／j.paerosci.2005.07.001

［3］LYNCH F T，KHODADOUST A.Effects of ice accretions on aircraft aerodynamics［J］.Progress in Aerospace Sciences，2001，37（8）：669－767.doi：10.1016／S0376－0421（01）00018－5

［4］Ice Protection Harmonization working group（IPHWG）.Task 2 Working group report on supercooled large droplet rulemaking［R］.2005.

［5］FAA.Advisory Circular：aircraft ice protection［Z］.20－73A.2006.

［6］BROEREN A P，BRAGG M B，ADDY H E，et al.Effect of high－fidelity ice accretion simulations on the performance of a full－scale airfoil model［R］.AIAA 2008－434.doi：10.2514／6.2008－434

［7］EDWARD A W，BROEREN A P，BRAGG M B.Aerodynamics of scaled runback ice accretions［J］.Journal of Aircraft，2008，45（2）：591－603.doi：10.2514／1.32274

［8］LEE S.Effects of supercooled large droplet icing on airfoil aerodynamic［D］.Urbana：University of Illinois，2001.

［9］BROEREN A P，BRAGG M B.Effect of airfoil geometry on performance with simulated intercycle ice accretions［R］.AIAA 2003－0728.doi：10.2514／1.4734

［10］CHEN K，CAO Y H，AN K W，et al.Application of hybird grid to analyzing complex iced airfoil aerodynamic performance［J］.ACTA Aeronautica et Astronautica Sinica，2007，28（s）：87－91.（in Chinese）陳科，曹義華，安克文，等.應(yīng)用混合網(wǎng)格分析復(fù)雜積冰翼型氣動性能［J］.航空學(xué)報，2007，28（s）：87－91.

［11］MORTENSEN K.CFD simulations of an airfoil with leading edge ice accretion［D］.Copenhagen：Technical University of Denmark，2008.

［12］DUNN T A，LOTH E，BRAGG M B.Computational investigation of simulated large－droplet ice shape on airfoil aerodynamics［J］.Journal of Aircraft，1999，36（5）：836－843.doi：10.2514／2.2517

［13］BROEREN A P，EDWARD A W，BUSCH G T，et al.Aerodynamic simulation of runback ice accretion［J］.Journal of Aircraft，2010，47（3）：924－939.doi：10.2514／1.46475.

［14］LEE S，BRAGG M B.Effects of simulated－spanwise ice shapes on airfoils：experimental investigation［R］.AIAA Paper 99－0092，1999.

［15］LAUNDER B E，SPALDING D B.The numerical computation of turbulent flows［J］.Comput.Methods in Appl.Mech.Eng.，1974，3（2）：269－289.doi：10.1016／0045－7825（74）90029－2.