冰脊對翼型氣動特性影響的數值模擬研究

周莉，徐浩軍，楊哲，蔡軍

（1.空軍工程大學工程學院，陜西西安 710038；2.中國人民解放軍駐陜飛公司軍代表室，陜西漢中 723213）

引言

飛機穿越含有過冷水滴的云層時，機體碰撞過冷水滴后會在其部件表面出現結冰。對于機翼來說，水滴易撞擊在翼剖面的前緣，使得機翼前緣結冰現象較為普遍，因此，出現了大量關于機翼前緣結冰過程及其對飛機性能影響的研究文獻［1-3］。但當過冷水滴的尺寸超出FAR25附錄C所給的2～50 μm結冰包線范圍時，在機翼上表面防冰區域后部會形成突起的展向冰脊，造成飛機氣動性能惡化，嚴重時甚至會引發飛行事故。

近幾年來，國外對大尺度過冷水滴的碰撞特性以及冰脊的形成進行了一系列試驗與數值模擬研究。Tan等［4］通過風洞試驗研究了大尺度過冷水滴的碰撞、噴濺和分裂特性，為帶冰脊翼型氣動特性的分析奠定了基礎。Raimund等［5］采用歐拉法對大尺度過冷水滴的碰撞特性進行了數值模擬，所得收集效率的分布與試驗結果吻合較好。此外，Lee等［6］和Pan等［7］分別運用試驗手段以及雷諾數平均的N-S方程研究了展向冰脊對翼型氣動特性的影響，使人們對冰脊危害的認識不斷深入。相比之下，國內關于冰脊的研究工作開展得較少，僅文獻［8］對冰脊高度和弦向位置對氣動特性的影響進行了仿真，但所計算的迎角范圍較小，沒有充分給出帶冰脊翼型的失速特性，也沒有考慮冰脊幾何形狀以及雷諾數對翼型性能的影響。針對冰脊對翼型氣動的危害以及前期研究的不足，本文采用數值模擬方法，系統研究了不同冰脊幾何形狀、弦向位置、冰脊高度以及不同雷諾數對翼型氣動性能的影響。

1 數值計算方法

1.1 網格劃分

本文計算采用O型結構網格，選取具有代表性的前向1/4圓、后向1/4圓、半圓和前向斜坡來模擬不同形狀的冰脊，其中斜坡面的基底長度與高度比為3。考慮到國外針對NACA23012翼型已開展了大量結冰方面的研究工作，有較為豐富的試驗數據，本文選用該翼型進行數值模擬。圖1為NACA23012翼型在未結冰以及分別加入以上各種冰脊后的翼型網格。

圖1 計算網格

1.2 方法驗證

本文對未結冰的NACA23012翼型進行數值模擬，并將計算結果與風洞試驗結果對比，以此來驗證本文數值計算方法的可行性和準確性。

采用無反射壓力遠場作為遠場邊界條件，計算時給定無窮遠處來流馬赫數及其方向。翼型表面為固壁邊界，流動無滑移、無穿透。針對SSTk-ω在近壁面區有良好的精度和算法穩定性的特點，選用SSTk-ω湍流模型進行計算［9］。

計算驗證狀態為:Ma=0.2，Re=1.8 ×106，α=5°。經過流場計算，表1中給出了本文所得升力系數、阻力系數的計算結果。與文獻［6］的試驗值進行比較，結果表明本文升力系數及阻力系數的計算結果與試驗值較為接近。由圖2中給出的NACA23012翼型數值計算與風洞試驗所得壓力分布的對比情況同樣可以看出，本文計算結果與試驗值吻合較好。因此，運用本文的計算方法模擬翼型外流較為可靠。

表1 計算氣動力系數與試驗值的比較

圖2 NACA23012翼型計算與試驗壓力系數對比

2 帶冰脊翼型的數值模擬

在不同的氣象條件以及飛行條件下，不同弦向位置將出現不同形狀、不同高度的冰脊，對飛機的氣動特性會產生不同程度的影響。通過數值模擬，得到不同條件下的翼型氣動系數。計算基準條件為:Ma=0.12，Re=10.5 ×106，冰脊高度為k/c=1.39%，冰脊弦向位置為x/c=0.1。其中，k為冰脊高度，x為冰脊在弦向的坐標，c為翼型弦長。

2.1 冰脊幾何形狀的影響

在未結冰以及冰脊幾何形狀分別為前向1/4圓、后向1/4圓、半圓以及前向斜坡四種類型時，通過數值計算得到不同迎角下的升力系數、阻力系數如圖3所示。

由圖3（a）可以看出，與未結冰情況相比，結冰脊后最大升力系數下降、失速迎角減小，且不同冰脊形狀對升力系數的影響程度并不相同。冰脊為前向1/4圓和前向斜坡時對升力系數的影響較為接近、且最為嚴重；冰脊為半圓時對升力系數的影響最小；而冰脊為后向1/4圓時對升力系數的影響程度介于半圓和前向1/4圓之間。對比具有相同鈍迎風面形狀的前向1/4圓和前向斜坡冰脊的升力系數可知，雖然下游形狀差別較大，但兩者的升力系數曲線幾乎一致，可見鈍迎風面冰脊的下游形狀對性能的影響并不是很重要；而通過具有相同流線型迎風面形狀的后向1/4圓和半圓冰脊的升力系數比較、以及具有相同下游形狀的前向1/4圓和半圓冰脊的升力系數比較可知，流線型迎風面冰脊的下游形狀以及冰脊迎風面的形狀對性能的影響較大，且冰脊迎風面形狀的影響程度更大。

圖3（b）所示為不同冰脊形狀下的阻力系數變化。可以看出，結冰脊后阻力系數均增大，且與不同形狀的冰脊對升力系數的影響趨勢一致。

圖3 冰脊幾何形狀對氣動系數的影響

2.2 冰脊弦向位置的影響

冰脊選用前向1/4圓，冰脊弦向位置分別為x/c=0.02，0.10和0.20時，不同迎角下的升力系數、阻力系數曲線如圖4所示。

由圖4（a）的計算結果可以看出，x/c=0.10時的最大升力系數最小，并且最早失速，可見此處對冰脊最為敏感，因此x/c=0.10處被認為是 NACA23012翼型最危險的冰脊位置。圖4（b）中，阻力系數隨迎角的變化曲線同樣也說明在冰脊高度一定的情況下，x/c=0.10處的氣動特性惡化最為嚴重。

圖4 冰脊弦向位置對氣動系數的影響

圖5為α=0°時未結冰翼型上表面的壓力云圖。可以發現，最危險的冰脊位置與此處流速（壓力）大小及逆壓梯度大小有關。這是因為在流速越大（壓力越小）處結冰會消耗掉更多的邊界層能量，同時，當結冰區越靠近最大逆壓梯度區域，則冰脊引起的下游氣流分離越嚴重，所產生的分離氣泡越大，導致翼型氣動性能損失越大。而x/c=0.10恰好同時處于NACA23012翼型的最小壓力與最大逆壓梯度區域，因此冰脊在此處對翼型氣動特性的影響最為明顯，即此處為冰脊的危險位置。

圖5 未結冰翼型上表面的壓力云圖

圖6給出了迎角α=0°時冰脊分別在x/c=0.02，0.10和0.20處的壓力系數分布曲線。可以看出，計算所得Cp分布與文獻［6］的試驗值基本吻合，但是并沒有充分捕捉到冰脊后的分離區。對比圖6可得出以下結論:冰脊在x/c=0.02時對Cp分布的影響最小，這是因為x/c=0.02處于順壓梯度范圍，所以冰脊所引起的分離泡能夠較容易地重新附著到翼型表面，導致分離泡較小；與前面壓力云圖的分析結果一致，x/c=0.10時冰脊上的Cp變化量最大，即冰脊影響最為嚴重；冰脊在x/c=0.20時，Cp分布變化程度介于x/c=0.02和x/c=0.10之間，這主要是因為雖然冰脊處于逆壓梯度區導致分離泡較大，但是x/c=0.20并不在最低壓力范圍，導致此處冰脊的影響并不是最嚴重。

圖6 冰脊在不同弦向位置處的壓力系數分布

2.3 冰脊高度的影響

冰脊選用前向1/4圓，當冰脊高度分別為k/c=0（未結冰），0.83%和1.39%時，升力系數、阻力系數曲線如圖7所示。

圖7（a）中，與未結冰情況相比，帶有冰脊時的升力系數曲線斜率減小，且最大升力系數及失速迎角明顯降低。k/c=1.39%時，最大升力系數由未結冰時的1.69降為0.23，失速迎角由17.7°減小到3.93°。k/c=0.83%時，最大升力系數為0.41，略大于k/c=1.39%時的0.23，失速迎角為5.98°，也高于k/c=1.39%時的3.93°。可見，冰脊越高時最大升力系數及失速迎角越小。

由圖7（b）可以看出，與未結冰時的阻力系數相比，k/c=0.83%和1.39%時阻力系數均增大，并且當冰脊越高時，阻力系數的增加幅度越大。

圖7 冰脊高度對氣動系數的影響

綜上分析可知，冰脊使翼型的氣動性能下降明顯，且冰脊越高，影響越大，即翼型氣動特性越差。

圖8給出了迎角α=0°、冰脊高度分別為k/c=0，0.83%和1.39%時的壓力系數分布曲線。

圖8 不同冰脊高度的壓力系數分布

可以看出，計算所得Cp分布與文獻［6］的試驗值基本吻合，與圖6相似，并沒有很好地捕捉到冰脊后的分離區，但是也充分反映出了不同冰脊高度下的翼型壓力系數分布:冰脊的存在使得翼型前緣上表面的壓力增加，即前緣吸力峰被拉低，導致升力減小；同時，隨著冰脊高度的增加，下表面的壓力會減小，從而造成升力的進一步損失。

2.4 雷諾數的影響

冰脊選用前向1/4圓，雷諾數分別為Re=3.5×106，7.5 ×106，10.5 ×106時，計算所得升力系數、阻力系數曲線如圖9所示。

圖9 雷諾數對氣動系數的影響

從圖9可以看出，三種不同的雷諾數下升力系數及阻力系數變化并不明顯，可見雷諾數對帶冰脊翼型的氣動特性影響較小。雖然對于未結冰翼型來說，當雷諾數增大時轉捩點前移，導致紊流區增大，從而延遲分離，使最大升力系數和失速迎角增大。但是當翼型帶有冰脊后，在冰脊處會發生氣流分離，此時對失速起決定性作用的已不再是附面層的氣流特性，因此改變雷諾數對氣動特性影響不大。

3 結論

本文針對冰脊對翼型氣動特性的嚴重影響，運用數值模擬方法系統研究了不同冰脊特性條件下翼型氣動性能的變化。計算結果表明:

（1）飛機在結冰脊后，最大升力系數下降、失速迎角減小，且阻力系數增大。

（2）冰脊迎風面的形狀對翼型性能的影響較大，并且當迎風面不同時，冰脊下游形狀對翼型性能的影響程度也不同:流線型迎風面冰脊的下游形狀影響翼型性能，而鈍迎風面冰脊的下游形狀對性能的影響并不明顯。

（3）不同冰脊弦向位置對氣動特性的影響不同，其中，x/c=0.10為最危險冰脊位置。

（4）冰脊高度也會影響翼型氣動特性，在一定冰脊高度范圍內，氣動特性的惡化程度隨冰脊高度的增加而增大。

（5）雷諾數對帶冰脊翼型氣動性能的影響并不明顯。

［1］孫志國，朱程香，付斌，等.二維翼型結冰數值計算［J］.航空動力學報，2010，25（7）:1485-1490.

［2］李林，王立新，彭小東.結冰對民機飛行性能的影響研究［J］.飛行力學，2004，22（3）:12-16.

［3］鐘長生，杜亮，洪冠新.飛機結冰引起的飛行動力學問題探討［J］.飛行力學，2004，22（3）:64-68.

［4］Tan S C，Papadakis M，Miller D，et al.Experimental study of large droplet splashing and breakup［R］.AIAA-2007-904，2007.

［5］Raimund H，Wagdi G H.Eulerian modeling of in-flight icing due to supercooled large droplets［J］.Journal of Aircraft，2008，45（4）:1290-1296.

［6］Lee S，Bragg M B.Experimental investigation of simulated large-droplet ice shapes on airfoil aerodynamics［J］.Journal of Aircraft，1999，36（5）:844-850.

［7］Pan J P，Loth E.Reynolds-averaged Navier-Stokes simulations of airfoils and wings with ice shapes［J］.Journal of Aircraft，2004，41（4）:879-891.

［8］易賢.飛機結冰的數值計算與結冰試驗相似準則研究［D］.綿陽:中國空氣動力研究與發展中心，2007.

［9］王福軍.計算流體動力學分析:CFD軟件原理與應用［M］.北京:清華大學出版社，2004.