多段翼型縫翼前緣結冰大迎角分離流動數值模擬

張恒， 李杰， 龔志斌

西北工業大學 航空學院， 西安 710072

多段翼型縫翼前緣結冰大迎角分離流動數值模擬

張恒， 李杰*， 龔志斌

西北工業大學 航空學院， 西安 710072

應用基于SST(Shear-Stress-Transport)湍流模型的IDDES(Improved Delayed Detached Eddy Simulation)方法，對大迎角狀態下多段翼型縫翼前緣典型角狀冰引起的復雜分離流動進行了數值模擬研究。采用后臺階流動標準算例和干凈無冰多段翼型分離流動算例對數值方法的可靠性和適用性進行了驗證。縫翼結冰狀態下的數值模擬結果表明：來流迎角較大時，前緣角狀冰將會導致結構相對穩定的流動分離泡產生，分離泡的非定常尾跡會對主翼前緣附近流場產生較為強烈的干擾，抑制了縫道流動的加速效應，使得縫翼增升效率降低。在失速點附近，由于分離泡回流強度隨來流迎角而增長，同時脫落旋渦的輸運方向逐漸向遠離壁面方向偏移，使得尾跡影響區域范圍和強度均有所增加。

IDDES方法； 多段翼型； 結冰； 分離流動； 湍流； 旋渦； 數值模擬

現代大型客機在起飛和著陸過程中通常使用多段翼形式的增升裝置，此時飛機飛行速度和高度相對較低，在過冷水滴條件下增升裝置各部件前緣容易產生結冰現象。其中縫翼前緣結角狀冰將導致多段翼最大升力系數顯著下降，失速迎角大幅提前，對飛行安全造成嚴重威脅[1]。目前關于該現象的流動機理認識仍然較為缺乏，當地流動特征變化與氣動性能損失之間的關聯性也并不明確。

對于多段翼型而言，部件流場間的相互耦合作用較為緊密，而結冰導致的分離流動被不同尺度旋渦的非定常運動過程和相互作用所支配[2-4]。因此，在大迎角狀態下縫翼結冰引起的分離流動和旋渦輸運過程將對多段翼型的整體繞流情況產生復雜的影響，相關數值模擬研究工作對湍流模擬方法的精度提出了較高的要求。

目前，雷諾平均Navier-Stokes(RANS)方法常用于翼型和部件結冰后流動問題的分析研究。該方法對湍流脈動信息以時均方式加以完全模化，對于完全附著或小分離流動能夠得到較為滿意的宏觀數值結果，但難以對復雜分離流動的細節特征及流動干擾現象進行合理描述。大渦模擬(Large Eddy Simulation, LES)方法通過引入某種過濾尺度，直接求解大尺度旋渦，利用亞格子應力模型模擬小尺度渦對大渦的影響；能夠有效描述分離區域的湍流瞬時脈動情況。但該方法對計算資源的要求相對較高，并且對于近壁附著流動而言，邊界層內部由小尺度、高頻率的各向同性旋渦結構主導，一般需要構造合理的壁面模型，因此該方法的應用廣泛程度仍較為有限。

DES[5](Detached Eddy Simulation)類方法結合了RANS方法和LES方法的特點，一定程度上能夠兼顧分離流動的計算精度和求解效率，近年來得到了快速發展。現有的研究工作[6-8]表明DES類方法能夠較為有效地描述單段翼型結冰后流場的分離特性；但針對多段翼型結冰后流場特征和氣動特性變化的分析研究目前相對較少。

Travin等[9]通過將DDES[10](Delayed DES)方法應用于LES壁面模型(Wall-modelling in LES, WMLES)，構造了IDDES(Improved DDES)方法。該方法不僅能夠解決DES方法直接應用于WMLES時產生的對數層不連續(Log-Layer Mismatch, LLM)問題[11]；并且就數值模擬效果而言，該方法有利于分離區域湍流結構的充分解析，同時在流動過渡區域也能夠取得更為滿意的結果，因此適宜分析附著和分離流動并存的結冰翼型分離流動問題。Xiao等[12-13]的工作表明該方法對不同類型分離流動分析的適用性較強，具備推廣到多段翼型結冰后復雜流動分析的潛力。

本文基于IDDES方法，就大迎角狀態下縫翼前緣角狀冰導致的流動分離問題進行了數值模擬研究。采用后臺階流動標準算例和干凈無冰多段翼型算例，驗證了數值方法模擬分離泡流動的可靠性和適用性。針對縫翼結冰分離流動算例，探討了大迎角狀態下角狀冰分離泡尾跡對多段翼型前緣繞流的影響，并分析了失速點附近迎角變化時分離泡特征及其尾跡流動的演化規律。

1 數值模擬方法

1.1 控制方程及其離散

在有限體積法基礎上，對三維可壓縮非定常Navier-Stokes方程進行求解。無黏通量項離散采用Roe-MUSCL(Roe-Monotone Upstream-centred Schemes for Conservation Law)[14]三階迎風通量差分分裂格式，黏性通量項離散采用二階中心差分格式。時間推進采用二階隱式近似因子分解方法。

1.2 湍流模擬方法

在Menterk-ω剪切應力輸運(SST)兩方程湍流模型[15]的基礎上，根據文獻[9]對IDDES方法進行構造，實現湍流數值模擬。該方法建立在原準DDES方法的基礎上，主要改進內容包括以下2方面。

1) 亞格子尺度重新定義

由于在亞格子尺度相同的情況下，模擬壁面附近自由剪切湍流所需的最優模型常數一般應小于分離區域各向同性湍流，這等價于在同一最優模型常數下對自由剪切湍流取較小的亞格子尺度。基于以上思路，文獻[9]中通過引入當地網格參數，對亞格子尺度做如下定義

Δ=min[max(Cwdw,Cwhmax,hwn),hmax]

(1)

式中:dw為網格單元與壁面距離;Cw=0.15為由LES解得到的經驗常數;hmax為網格單元三向最大尺度;hwn為當地網格單元壁面法向高度。在近壁面區域，對于長寬比較大的RANS薄層網格，一般總有max(Cwdw,Cwhmax,hwn)hmax。因此該定義能夠減小自由剪切湍流區域亞格子尺度，同時維持各向同性湍流區域亞格子尺度不變。

2) RANS/LES 混合長度構造

該混合長度由RANS長度尺度和LES長度尺度2部分構成，對于本文應用的k-ωSST湍流模型，混合長度函數形式為

(2)

式中:lRANS為RANS長度尺度，對SST湍流模型有lRANS=k1/2/(Cμω)，k為湍動能，Cμ為經驗常數，ω為比耗散率;CDES為亞格子應力模型常數。

當該混合函數在WMLES模式下工作時，函數fhyb能夠加快分離區域RANS方法到LES方法的轉換，同時函數frestore能夠防止RANS和LES區域交界面附近的雷諾應力損失。

混合函數fhyb包含了DDES分支和WMLES分支，其構造形式為

fhyb=max(1-fd,fstep)

(3)

式中：fd為DDES方法中的延遲函數，

fd=1-tanh[(8rd)3]

(4)

式中：rd為湍流延遲函數，

(5)

其中：νt為湍流運動黏性系數。

fstep是WMLES分支函數，構造形式為

fstep=min{2exp(-9α)2,1.0}

(6)

式中：函數α=0.25-dw/hmax，fstep的構造能夠使得在0.5

函數frestore在WMLES分支下啟動，其構造形式為

frestore=max{fhill-1,0}famp

(7)

(8)

famp=1.0-max{ft,fl}

(9)

(10)

(11)

(12)

式中:νl為層流運動黏性系數；cl和ct為湍流模型相關常數，對于本文應用的k-ωSST湍流模型，cl取5.00，ct取1.87。

2 后臺階流動驗證算例

由于翼型前緣結角狀冰后將產生典型分離泡結構，從流動相似性角度出發，選取后臺階流動標準算例對數值方法的可靠性進行考核驗證。

該算例是Travin等[9]關于IDDES方法的測試算例之一；本文計算域幾何尺寸與其保持一致，并沿用了網格拓撲形式和結點分布，總網格量約為1.5×106。圖1給出了計算網格空間截面，其中展向長度取臺階高度h的2倍。計算域入口給定速度型邊界條件，物面采用絕熱、無滑移和法向零壓力梯度條件，出口給定無反射邊界條件，展向設置周期性邊界條件。入口來流中心速度U=11.3 m/s，基于臺階高度的雷諾數為Re=2.8×104。

基于k-ωSST湍流模型進行非定常RANS計算，獲得充分發展的初始流場，在此初場基礎上進行后續IDDES計算。無量綱時間步長Δt*=UΔt/h=0.01，Δt為物理時間。計算至非定常流場基本穩定后進行時間平均；在時均流場的基礎上進行展向空間平均，得到時空平均流場。

圖2給出了時均摩擦系數Cf計算結果與試驗值的對比情況，其中RANS結果是在上述非定常初場基礎上作時空平均得到的。IDDES方法較好地描述了因流動分離形成的負摩阻區域。圖3 給出了分離泡內部垂直壁面方向速度型u/U時均計算結果，d為壁面距離，其中再附點x/h=7.03處的速度分布與試驗值吻合良好，表明IDDES方法能夠較為準確地描述分離泡的時均再附特征。圖4給出了計算所得瞬態旋渦結構表征參數Q等值面[16]分布，表明IDDES方法能夠解析分離區域的主要湍流結構，并刻畫剪切層失穩后內部脫落旋渦發生的滾轉和變形現象，但在出口位置附近未能較好地給出當地流動情況。

圖1 計算網格空間截面
Fig.1 Cross-section of computational grid

圖2 時均摩擦系數計算結果與試驗值對比
Fig.2 Comparison of time-averaged friction coefficients with test data

圖3 不同位置速度型時均計算結果與試驗值對比
Fig.3 Comparison of time-averaged velocity profile distributions with test data at different locations

圖4 瞬態Q等值面分布(Q=0.000 5)
Fig.4 Instantaneous Q iso-surface distribution
(Q=0.000 5)

該算例表明，IDDES方法能夠較好地描述典型后臺階流動分離泡的基本特征及其尾跡的發展變化過程，所構造的數值方法適用于分析存在分離泡的流動問題。

3 干凈無冰翼型驗證算例

由于現有多段翼型帶冰試驗結果提供的流場細節參數比較有限，且縫翼冰型不具備角狀特征[17]，故本文選取ONERA RA16SC1翼型算例[18-19]，以對數值方法與多段翼型相關的分離流動模擬效果進行校驗，同時為縫翼結冰多段翼型計算分析提供基準網格方案。

3.1 幾何模型與計算網格

該翼型是ONERA為歐洲EUROPIV-2流動顯示研究項目設計的標準翼型，也是DES類混合方法研究項目DESider的測試算例[18]。模型幾何尺寸與風洞試驗[20]一致，展向長度取為干凈翼型弦長c的0.1倍，圍繞模型建立15c×15c的矩形計算域。計算域遠場給定無反射邊界條件，物面采用絕熱、無滑移和法向零壓力梯度條件，展向設置周期性邊界條件。對于縫翼凹腔內的大渦區域以及后緣襟翼附近流動分離區域，網格單元三向尺度取為2%c。近壁面法向首層網格高度取為10-5c，滿足近壁面y+<1。計算域內網格結點總數約為1.4×107。縫翼凹腔附近計算網格對比情況如圖5所示。

圖5 關注區域計算網格空間截面
Fig.5 Grid cross-section of focus regions

3.2 計算結果及分析

采用IDDES方法進行數值模擬，計算條件按照風洞試驗[20]給定：來流馬赫數Ma=0.160 6，基于干凈翼型弦長的雷諾數為Re=1.7×106，來流迎角α=12°。由于洞壁效應的影響，數值模擬過程中通常要對來流迎角進行修正，根據文獻[19]所得的相關結論，計算迎角設為9°。

在非定常RANS初場的基礎上進行后續IDDES計算，得到時空平均流場。無量綱時間步長Δt*=0.004。以下將IDDES時空平均數值結果與風洞試驗結果[20]及其他數值方法結果進行對比，其中RANS結果是在非定常初場計算基礎上作時空平均得到的。

圖6給出了翼型表面時均壓力分布Cp計算結果與試驗結果對比情況，不同方法計算得到的壓力分布形態與試驗結果總體吻合。IDDES方法預測的襟翼前緣吸力峰值略高，襟翼后緣出現壓力平臺，與試驗得到的流動分離現象相對應。

圖6 時均壓力系數分布計算結果與試驗值對比
Fig.6 Comparison of time-averaged pressurecoefficients distribution with test data

圖7分別給出了縫翼凹腔附近時均流場速度分布U/U∞的試驗和計算結果。由圖可知，IDDES方法計算得到的凹腔流動分離形態與試驗較為相似，渦核位置比較準確，且縫道附近的流動再附位置與試驗基本相同。

圖8分別給出了后緣襟翼附近時均流場速度分布U/U∞及流線的試驗PIV(Particle Image Velocimetry)測量結果和計算結果。試驗結果顯示襟翼上表面后緣位置出現了低強度的流動分離現象。文獻[18]采用的ZDES(Zonal DES)方法難以反映分離特征，而文獻[19]基于ILES(Implicit LES)方法得到的分離流動強度偏高。本文計算得到的速度分布情況與試驗結果基本一致，反映IDDES方法具備良好的近壁分離流動預測能力。

圖7 縫翼附近時均流場速度云圖及流線圖對比
Fig.7 Comparison of time-averaged flow field velocity contour and streamlines around slat

圖8 襟翼附近時均流場速度云圖及流線圖對比
Fig.8 Comparison of time-averaged flow field velocity contour and streamlines around flap

圖9給出了壁面附近的6個速度型監測點位置。圖10給出了上述位置壁面法向速度型U/U∞時均計算結果與試驗結果的對比情況。在0號監測點所處的大渦區域，IDDES方法能夠較RANS方法更好地捕捉速度峰值。在1號到4號監測點，IDDES方法得到的計算結果較RANS方法與試驗值吻合程度更為良好。在5號監測點處，RANS方法獲得的速度分布與試驗結果存在較大差異，無法描述壁面附近的回流趨勢；IDDES方法獲得的速度型形態與試驗值更為接近，能夠捕捉到近壁流動的分離現象。

圖9 速度型監測點位置分布圖
Fig.9 Locations of measuring points of velocity profiles

圖10 不同監測點速度型時均結果
Fig.10 Time-averaged velocity profile distributions at different locations

該算例表明，就干凈多段翼型分離流動問題而言，IDDES方法對于能夠取得與風洞試驗吻合程度良好的計算結果，可以應用于結冰條件下多段翼型分離流場的分析研究。

4 結冰翼型幾何模型與計算網格

選取結冰風洞試驗常用的LB606b三段翼型進行數值模擬研究。Miller等[21]利用NASA Glenn研究中心的Lewis風洞就該翼型進行了多組結冰試驗，研究了不同成冰條件下各個部件的結冰情況，獲得了一系列有代表性的冰形。本文選取上述結冰試驗中成冰時間為6 min的角狀冰作為縫翼前緣結冰典型冰形。多段翼型縫翼前緣帶冰幾何形狀如圖11所示。

根據文獻[22]的觀點，多段翼型計算域展向長度應至少與關注區域內的大尺度旋渦尺寸相當；此處考慮到縫翼上表面區域可能形成分離泡，故令計算域展向長度與縫翼弦長相等。流場入口距縫翼前緣為15c，出口距襟翼后緣為15c，上下邊界距翼型表面各為12c。

沿用ONERA RA16SC1多段翼型算例中的網格基本拓撲形式，在此基礎上對角狀冰后方分離區域網格進行重點設計。對于此類流場擾動源比較明確的問題，參考Spalart[23]的觀點，分離區域網格單元尺度Δ0的選取原則是若要滿足對波長為λ的旋渦進行LES解析的要求，則Δ0應取為λ的1/5。對于本算例而言，由于縫翼前緣冰角高度約為10-2c，認為縫翼上表面及主翼前緣附近關注區域內主要旋渦結構的擾動波長與此相當。根據上述觀點，在關注區域內Δ0分別取2×10-3c、10-3c和0.5×10-3c，構造粗、中、細3套計算網格，以對網格密度進行驗證。近壁面法向首層網格高度維持10-5c，滿足近壁面y+<1。三套網格計算域內網格結點總數約為0.5×107、2.0×107、9.2×107，流動關注區域附近網格對比情況如圖12所示。

圖11 縫翼前緣帶角狀冰LB606b多段翼型
Fig.11 Multi-element airfoil LB606b with horn ice on leading edge of slat

圖12 流動關注區域附近計算網格對比
Fig.12 Comparison of computational grids near focus regions

5 結冰翼型計算結果及分析

結冰翼型相關各算例計算狀態參照結霜冰翼型氣動力測量風洞試驗[17]給定：來流馬赫數Ma∞=0.2，基于翼型弦長的雷諾數Re=9×106，采用IDDES方法進行數值模擬，非定常計算無量綱時間步長Δt*=0.003。其余計算手段與第3節干凈無冰翼型算例維持一致。

以下首先給出網格無關性驗證結果；在此基礎上對干凈無冰與縫翼結冰結果進行了對比分析；并研究了不同迎角下分離泡及其尾跡流動的演化過程與影響效應。

5.1 網格無關性驗證

選取迎角α=20° 計算狀態進行網格無關性驗證，對應多段翼型干凈無冰條件下的最大升力系數狀態。圖13以時均流場速度分布U/UINF的形式給出了不同網格基礎上計算相同步數后所得時均分離泡的大致形狀。由圖可知基準網格和密網格基礎上獲得的分離泡形狀大體相同，再附位置和法向高度基本一致；密網格獲得的流線邊界更加光滑，渦核位置略微靠后；兩套網格在相同計算步數下的分離泡大致收斂于同一結果。稀網格獲得的分離泡形狀與之存在較大差異，且渦核位置過于靠近縫翼后緣。由于本文主要研究分離泡對周圍流動的影響，參考單段結冰翼型數值模擬相關文獻[6]的結論，認為分離泡幾何形態基本收斂后獲得的計算結果即可為流場分析提供支持，綜合考慮計算精度和效率，選取中等密度網格進行后續長時間平均計算。

圖13 不同密度網格時均分離泡計算結果對比
Fig.13 Comparison of time-averaged separation bubble of different grid density

5.2 干凈無冰與縫翼結冰結果

在上述中等密度網格計算結果基礎上進行長時間平均續算。以下從時均和瞬態角度對縫翼結冰翼型計算結果與干凈翼型同一計算條件下的相應結果進行對比分析，以期揭示縫翼結冰對多段翼型大迎角流場與氣動特性的影響機制。

圖14給出了計算穩定階段所得氣動力系數的時間序列，圖中Iter表示迭代步數，離散點表示各氣動力系數的瞬時值，實線表示時間平均值。由時均結果可知縫翼結冰后多段翼型升力系數CL為3.25左右，相對于干凈無冰翼型相同迎角風洞試驗結果[17](CL=4.354 8)的下降幅值約為25%。

圖14 氣動力系數時間序列
Fig.14 Time history of aerodynamic force coefficients

圖15給出了結冰前后多段翼型時均壓力分布Cp變化情況，其中干凈翼型壓力分布數據取自風洞試驗[17]。由圖可知，干凈翼型呈現典型多段翼型大迎角壓力分布，縫翼上表面未出現明顯的流動分離現象，數值模擬結果與試驗值吻合情況良好。結冰條件下縫翼前緣負壓峰值消失，壓力分布呈現平臺特征。主翼上表面產生了較為明顯的全局負壓損失，但其分布特征基本與干凈翼型相似，這與文獻[17]采用霜冰模型開展風洞試驗所得結論定性一致，表明主翼后緣并未出現顯著的流動分離，流動仍能維持附著。除縫翼下表面由于不規則積冰的存在而產生鋸齒狀壓力分布外，整個翼型下表面壓力分布情況基本與干凈無冰狀態相同。由以上計算結果可知，結冰后多段翼型總升力下降的直接原因是主翼升力貢獻的降低，縫翼升力的絕對下降量則比較有限。

圖16給出了結冰前后縫翼附近時均流場速度分布U/UINF對比情況。由圖可知，在干凈無冰條件下，縫翼上表面流動處于完全附著狀態；當縫翼前緣存在積冰時，形成了尺度與縫翼弦長相當的時均分離泡，使外部繞流速度大小有所增長，但主翼前緣附近卻表現出了較為明顯的速度損失。此外，由于縫翼下表面存在積冰，改變了凹腔位置的當地來流方向，令回流渦大小有一定程度地增加，但結構基本保持不變。

在主翼上表面設置如圖17(a)所示的6個監測點P1～P6，以分析結冰前后主翼前緣附近流動特征的變化情況。圖17(b)統計了相同監測點附近法向同一高度(z-zsurf=0.015)處時均速度u/U改變量，z為法向距離；zsurf為物面高度。圖中顯示結冰后主翼前緣附近的速度損失約為15%左右，在7%弦長位置速度損失達到峰值。

圖15 時均流場壓力系數分布對比
Fig.15 Comparison of time-averaged pressurecoefficient distribution

圖16 縫翼附近時均流場對比
Fig.16 Comparison of time-averaged flow field near slat

圖17 主翼前緣附近時均速度變化情況對比
Fig.17 Comparison of time-averaged velocity changes near leading edge of main wing

由圖18給出的壁面法向速度型u/U時均計算結果可知，起始監測點1附近結冰前后速度分布情況大體相同，表明分離流動對縫道出口流動影響較小。從監測點2開始出現了明顯的速度損失。結冰前后各監測點處流動速度量值沿弦向的變化趨勢比較相似，速度型斜率基本相同，并且壁面附近速度分布基本維持不變。這表明結冰導致的速度損失在弦向和法向空間分布上相對均勻，與當地來流狀況的變化直接相關。

圖19給出了縫翼結冰前后主翼前緣附近的瞬態空間流線分布對比情況。由圖可知，結冰條件下，縫道出口位置處分離泡尾跡流動與縫道流動間存在較為復雜的相互摻混作用，較干凈無冰狀態顯著改變了主翼前緣附近的流動情況，從而抑制了縫道流動的加速效應，使得縫翼對主翼所起的增升作用大幅下降。

圖18 主翼前緣附近時均速度分布對比
Fig.18 Comparison of time-averaged velocity distribution near leading edge of main element

圖19 瞬態空間流線分布對比
Fig.19 Comparison of instantaneous streamlines distribution

圖20給出了瞬態流場無量綱展向渦量z-vorticity 分布情況，由圖可知，冰角后方的脫落剪切層具備相對較高的強度，在縫翼后緣點附近呈現出較為明顯的失穩和破碎現象；由于縫翼本體幾何尺寸較小，其上表面流動完全為剪切層的失穩脫落過程所支配，所形成分離泡的結構和大小隨時間的變化情況基本較為穩定，可認為流動處于一種“臨界再附”狀態。圖中顯示分離泡的非定常尾跡由剪切層失穩后旋渦脫落產生的一系列渦串結構組成，具備較為規則的分布形式，對主翼前緣繞流情況產生了強烈的干擾；但旋渦輸運方向基本與來流方向一致，并未直接與主翼前緣作用。隨著分離泡尾跡的耗散，尾跡渦的強度不斷減弱，對主翼上表面流動的影響也逐步降低，這與圖18給出的速度分布情況相吻合。

圖21給出了計算所得瞬態Q等值面分布，由圖可知分離泡尾跡流場中的主要三維旋渦基本能夠得到解析，反映了剪切層沿弦向逐步失穩破碎，并產生更為復雜的尾跡渦結構這一過程。圖中顯示尾跡后半段的湍流結構較縫道出口附近反而更為復雜，這表明縫道流動對分離泡尾跡也具備一定程度的吹除和抑制作用。

圖20 t=20時瞬態流場展向渦量分布
Fig.20 Instantaneous spanwise vortices distribution at t=20

圖21 t=20時瞬態Q等值面分布(Q=50)
Fig.21 Instantaneous Q iso-surface distribution at t=20 (Q=50)

5.3 縫翼結冰不同迎角結果

選取α=20° 附近兩個狀態α=18° 和α=22° 進行計算，以期對大迎角條件下分離泡的演化過程及其尾跡流動對主翼前緣附近流場的影響效應進行分析。

表1給出了時均升力系數對比情況，其中干凈翼型升力系數取自風洞試驗[17]數據。結果表明在α=18° 時，多段翼型就已進入失速狀態，失速迎角較干凈翼型有所提前，但升力系數變化情況同樣較為和緩。圖22給出了結冰翼型不同迎角下的時均壓力分布變化情況。由圖可知，隨著迎角增加，縫翼上表面負壓值呈平臺狀逐漸上升；主翼后緣、襟翼前緣負壓區縮小，但主翼前緣壓力峰值差異不大，表明分離泡尾跡流動可能抑制了負壓峰值的增加。

圖23給出了不同迎角下縫翼附近時均流場速度分布U/UINF的對比情況。在本文涉及的計算狀態下，縫翼上表面都出現了大尺度時均分離泡結構；隨著迎角增加，流動剪切效應增強，分離泡法向高度有所增長，渦核位置不斷后移，但再附位置變化并不明顯；分離泡尾跡附近加速區域范圍逐漸縮小，速度損失量值增加，表明分離泡強度的增加可能成為影響加速區域的主導效應。

圖24給出了不同迎角下主翼前緣3個監測點的時均速度分布u/U變化情況，為闡明分離泡尾跡流動對遠離壁面區域流場的影響，法向監測距離較圖18有所增加，其中速度型峰值體現了來流速度的損失程度。由圖可知，來流速度損失隨迎角而增加，表明尾跡流動的影響在不斷增強。同時法向減速區域沿主翼前緣逐漸擴大，表明尾跡的空間擴散范圍也在增加；但近壁速度有所恢復，尾跡對壁面附近流動的影響可能逐步減弱。

表1 時均升力系數對比Table 1 Comparison of time-averaged lift coefficients

圖22 不同迎角時均流場壓力分布對比
Fig.22 Comparison of time-averaged pressure distribution at different angles of attack

圖23 不同迎角縫翼附近時均流場對比
Fig.23 Comparison of time-averaged flow field near slat at different angles of attack

圖25給出了不同迎角下瞬態流場無量綱展向渦量z-vorticity的分布情況。對比各圖可知，在所研究的迎角變化范圍內，剪切層脫落渦能夠在縫翼后緣附近向壁面方向輸運，分離泡結構始終能得到維持。不同迎角下分離泡瞬態尾跡流場的基本結構較為類似，仍由一系列較為規則的渦串組成；且尾跡旋渦運動特征并不存在定性變化。隨著迎角增加，剪切層脫落角相應增加，弦向失穩位置提前；失穩后析出的旋渦尺度存在增長趨勢，強度有所提高；尾跡渦渦核與主翼壁面距離增加，尾渦列具備向外部自由來流流動區域移動的趨勢，這與圖24體現的近壁速度分布相吻合。

圖26給出了不同迎角下計算所得瞬態Q等值面對比情況。由圖可知，隨著迎角增加，剪切層起始失穩區域逐步前移，尾跡區域釋放出的三維湍流結構更加豐富，拉伸效應增強，但旋渦基本空間結構大體相同。

圖24 不同迎角主翼前緣附近時均速度分布對比
Fig.24 Comparison of time-averaged velocity distribution near leading edge ofmain wing at different angles of attack

圖25 不同迎角瞬態流場展向渦量分布
Fig.25 Instantaneous spanwise vortices distribution at different angles of attack

圖26 不同迎角瞬態Q等值面分布(Q=50)
Fig.26 Instantaneous Q iso-surface distribution at different angles of attack (Q=50)

6 結 論

1) 就典型后臺階流動算例及多段翼型后緣分離流動算例而言，IDDES方法能夠較好地描述分離泡的物理特征。

2) 大迎角狀態下，縫翼前緣結冰后多段翼型總體升力下降的直接原因是主翼升力貢獻的降低，縫翼本體升力的絕對下降量比較有限。

3) 縫翼前緣結冰后將形成大尺度分離泡結構，其非定常尾跡將對主翼前緣流動產生強烈干擾，導致縫道流動對主翼增升效應的下降。

4) 失速點附近，分離泡強度隨迎角而增長，尾跡旋渦輸運方向逐漸偏離壁面，尾跡影響范圍和強度均有所增加，但流場結構不存在定性變化。

[1] LYNCH F T, KHODADOUST A. Effects of ice accretions on aircraft aerodynamics[J]. Progress in Aerospace Sciences, 2001, 37(8): 669-767.

[2] GURBACKI H M, BRAGG M B. Unsteady flowfield about an iced airfoil：AIAA-2004-0562[R]. Reston: AIAA, 2004.

[3] ANSELL P J, BRAGG M B. Measurement of unsteady flow reattachment on an airfoil with a leading-edge horn-ice shape: AIAA-2012-2797[R]. Reston: AIAA, 2012.

[4] ANSELL P J, BRAGG M B. Characterization of ice-induced low-frequency flowfield oscillations and their effect on airfoil performance: AIAA-2013-2673[R]. Reston: AIAA, 2013.

[5] SPALART P R, JOU W H, STRELETS M, et al. Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach[M]. Advances in DNS/LES, Columbus: Greydon Press, 1997.

[6] THOMPSON D S, MOGILI P. Detached-eddy simulations of separated flow around wings with ice accretions: year one report: CR-2004-213379[R]. Washington, D.C.: NASA, 2004.

[7] MOGILI P, THOMPSON D S, CHOO Y, et al. RANS and DES computations for a wing with ice accretion: AIAA-2005-1372[R]. Reston: AIAA, 2005.

[8] LORENZO A, VALERO E, DE-PABLO V. DES/DDES post-stall study with iced airfoil: AIAA-2011-1103[R]. Reston: AIAA, 2011.

[9] TRAVIN A K, SHUR M L, SPALART P R, et al. Improvement of delayed detached-eddy simulation for LES with wall modeling[C]//European Conference on Computational Fluid Dynamics, 2006.

[10] SPALART P R, DECK S, SHUR M, et al. A new version of detached-eddy simulation, resistant to ambiguous grid densities[J]. Theoretical and Computational Fluid Dynamics, 2006, 20(3): 181-195.

[11] NIKITIN N V, NICOUD F, WASISTHO B, et al. An approach to wall modeling in large-eddy simulations[J]. Physic of Fluids, 2005, 12(7): 1629-1632.

[12] XIAO Z X, LIU J, LUO K Y, et al. Numerical investigation of massively separated flows past rudimentary landing gear using advanced DES approaches[J]. AIAA Journal, 2013, 51(1): 107-125.

[13] HUANG J B, XIAO Z X, LIU J, et al. Simulation of shock wave buffet and its suppression on an OAT15A supercritical airfoil by IDDES[J]. Science China Physics, Mechanics & Astronomy, 2012, 55(2): 260-271.

[14] VAN LEER B. Towards the ultimate conservative difference scheme V: A second order sequel to Godunov’s method[J]. Journal of Computational Physics, 1979, 32: 101-136.

[15] MENTER F R. Two-equation eddy viscosity turbulence models for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605.

[16] HUNT J, WRAY A, MOIN P. Eddies, streams and convergence zones in turbulent flows[C]//Proceedings of the 1988 summer program. Stanford: Center for Turbulence Research, 1988.

[17] KHODADOUST A, DOMINIK C. Effect of in-flight ice accretion on the performance of a multi-element airfoil： TM-112174[R]. Washington, D.C.: NASA, 1995.

[18] DECK S. Zonal-Detached-Eddy simulation of the flow around a high-lift configuration[J]. AIAA Journal， 2005, 43(11): 2372-2384.

[19] ZHONG B, SCHEURICH F, TITAREV V, et al. Turbulent flow simulations around a multi-element airfoil using URANS, DES and ILES approaches： AIAA-2009-3799[R]. Reston: AIAA, 2009.

[20] ARNOTT A D, SCHNEIDER G, NEITZKE K P, et al. Detailed characterization using PIV of the flow around an airfoil in high-lift configuration[M]//Particle Image Velocimetry: Recent Improvements. Berlin Heidelberg: Springer, 2004.

[21] MILLER D, SHIN J W, SHELDON D， et al. Further investigations of icing effects on an advanced high-lift multi-element airfoil：TM-106947[R]. Washington, D.C.: NASA, 1995.

[22] PENG S H, NEBENFUHR B, DAVIDSON LARS. Lessons learned from hybrid RANS-LES computations of a three-element airfoil flow： AIAA-2013-2741[R]. Reston: AIAA, 2013.

[23] SPALART P R. Young-person’s guide to detached-eddy simulation grids: CR-2001-211032[R]. Washington, D.C.: NASA, 2001.

(責任編輯： 李明敏)

URL：www.cnki.net/kcms/detail/11.1929.V.20161121.1439.002.html

Numericalsimulationofseparatedflowaroundamulti-elementairfoilathighangleofattackwithicedslat

ZHANGHeng,LIJie*,GONGZhibin

SchoolofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China

Theimproveddelayeddetachededdysimulation(IDDES)basedontheshear-stress-transport(SST)turbulentmodelisappliedinthenumericalsimulationofcomplexseparatedflowcausedbyatypicalhorn-likeiceontheslatleadingedgeofamulti-elementairfoilunderalargeangleofattack.Thereliabilityandapplicabilityofthenumericalmethodisverifiedbasedontheanalysisofthestandardexampleoftheseparationflowandtheexampleofthecleanmulti-elementairfoil.Resultsofnumericalsimulationoftheslatinicingconditionshowthatthehorniceontheleadingedgewillleadtoformationofalargescaleseparationbubblewithrelativelystablestructureathighangleofattack.Theunsteadywakeoftheseparationbubblewillgeneraterelativelystronginterferencewiththeflowfieldaroundtheleadingregionofthemainwing,resultinginthedeclineofaccelerationeffectoftheflowfromthegapandthedecreaseofslatliftaugmentationefficiency.Nearthestallpoint,asthebackflowintensityofseparationbubbleincreaseswiththeangleofattackoftheincomingflowandthewakevortextransportdirectiongraduallydeviatesfromthewallsurface,therangeandintensityoftheinfluenceareaofthewakeareincreased.

improveddelayeddetachededdysimulation(IDDES)method;multi-elementairfoil;icing;separationflow;turbulence;vortex;numericalsimulation

2016-08-31；Revised2016-09-26；Accepted2016-11-03；Publishedonline2016-11-211439

s：NationalBasicResearchProgramofChina(2015CB755800);NationalNaturalScienceFoundationofChina(11172240);AeronauticalScienceFoundationofChina(2014ZA53002)

.E-maillijieruihao@163.com

2016-08-31；退修日期2016-09-26；錄用日期2016-11-03； < class="emphasis_bold">網絡出版時間

時間：2016-11-211439

www.cnki.net/kcms/detail/11.1929.V.20161121.1439.002.html

國家“973”計劃 (2015CB755800)； 國家自然科學基金 (11172240)； 航空科學基金 (2014ZA53002)

.E-maillijieruihao@163.com

張恒， 李杰， 龔志斌． 多段翼型縫翼前緣結冰大迎角分離流動數值模擬J． 航空學報,2017,38(2):520733.ZHANGH,LIJ,GONGZB.Numericalsimulationofseparatedflowaroundamulti-elementairfoilathighangleofattackwithicedslatJ.ActaAeronauticaetAstronauticaSinica,2017,38(2):520733.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0285

V211.3

A

1000-6893(2017)02-520733-14