多段翼型縫翼前緣結(jié)冰大迎角分離流動數(shù)值模擬

張恒， 李杰， 龔志斌

西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072

多段翼型縫翼前緣結(jié)冰大迎角分離流動數(shù)值模擬

張恒， 李杰*， 龔志斌

西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072

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

IDDES方法； 多段翼型； 結(jié)冰； 分離流動； 湍流； 旋渦； 數(shù)值模擬

現(xiàn)代大型客機在起飛和著陸過程中通常使用多段翼形式的增升裝置，此時飛機飛行速度和高度相對較低，在過冷水滴條件下增升裝置各部件前緣容易產(chǎn)生結(jié)冰現(xiàn)象。其中縫翼前緣結(jié)角狀冰將導(dǎo)致多段翼最大升力系數(shù)顯著下降，失速迎角大幅提前，對飛行安全造成嚴(yán)重威脅[1]。目前關(guān)于該現(xiàn)象的流動機理認(rèn)識仍然較為缺乏，當(dāng)?shù)亓鲃犹卣髯兓c氣動性能損失之間的關(guān)聯(lián)性也并不明確。

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

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

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

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

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

1 數(shù)值模擬方法

1.1 控制方程及其離散

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

1.2 湍流模擬方法

在Menterk-ω剪切應(yīng)力輸運(SST)兩方程湍流模型[15]的基礎(chǔ)上，根據(jù)文獻(xiàn)[9]對IDDES方法進(jìn)行構(gòu)造，實現(xiàn)湍流數(shù)值模擬。該方法建立在原準(zhǔn)DDES方法的基礎(chǔ)上，主要改進(jìn)內(nèi)容包括以下2方面。

1) 亞格子尺度重新定義

由于在亞格子尺度相同的情況下，模擬壁面附近自由剪切湍流所需的最優(yōu)模型常數(shù)一般應(yīng)小于分離區(qū)域各向同性湍流，這等價于在同一最優(yōu)模型常數(shù)下對自由剪切湍流取較小的亞格子尺度。基于以上思路，文獻(xiàn)[9]中通過引入當(dāng)?shù)鼐W(wǎng)格參數(shù)，對亞格子尺度做如下定義

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

(1)

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

2) RANS/LES 混合長度構(gòu)造

該混合長度由RANS長度尺度和LES長度尺度2部分構(gòu)成，對于本文應(yīng)用的k-ωSST湍流模型，混合長度函數(shù)形式為

(2)

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

當(dāng)該混合函數(shù)在WMLES模式下工作時，函數(shù)fhyb能夠加快分離區(qū)域RANS方法到LES方法的轉(zhuǎn)換，同時函數(shù)frestore能夠防止RANS和LES區(qū)域交界面附近的雷諾應(yīng)力損失。

混合函數(shù)fhyb包含了DDES分支和WMLES分支，其構(gòu)造形式為

fhyb=max(1-fd,fstep)

(3)

式中：fd為DDES方法中的延遲函數(shù)，

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

(4)

式中：rd為湍流延遲函數(shù)，

(5)

其中：νt為湍流運動黏性系數(shù)。

fstep是WMLES分支函數(shù)，構(gòu)造形式為

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

(6)

式中：函數(shù)α=0.25-dw/hmax，fstep的構(gòu)造能夠使得在0.5

函數(shù)frestore在WMLES分支下啟動，其構(gòu)造形式為

frestore=max{fhill-1,0}famp

(7)

(8)

famp=1.0-max{ft,fl}

(9)

(10)

(11)

(12)

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

2 后臺階流動驗證算例

由于翼型前緣結(jié)角狀冰后將產(chǎn)生典型分離泡結(jié)構(gòu)，從流動相似性角度出發(fā)，選取后臺階流動標(biāo)準(zhǔn)算例對數(shù)值方法的可靠性進(jìn)行考核驗證。

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

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

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

圖1 計算網(wǎng)格空間截面
Fig.1 Cross-section of computational grid

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

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

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

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

3 干凈無冰翼型驗證算例

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

3.1 幾何模型與計算網(wǎng)格

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

圖5 關(guān)注區(qū)域計算網(wǎng)格空間截面
Fig.5 Grid cross-section of focus regions

3.2 計算結(jié)果及分析

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

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

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

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

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

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

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

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

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

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

4 結(jié)冰翼型幾何模型與計算網(wǎng)格

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

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

沿用ONERA RA16SC1多段翼型算例中的網(wǎng)格基本拓?fù)湫问剑诖嘶A(chǔ)上對角狀冰后方分離區(qū)域網(wǎng)格進(jìn)行重點設(shè)計。對于此類流場擾動源比較明確的問題，參考Spalart[23]的觀點，分離區(qū)域網(wǎng)格單元尺度Δ0的選取原則是若要滿足對波長為λ的旋渦進(jìn)行LES解析的要求，則Δ0應(yīng)取為λ的1/5。對于本算例而言，由于縫翼前緣冰角高度約為10-2c，認(rèn)為縫翼上表面及主翼前緣附近關(guān)注區(qū)域內(nèi)主要旋渦結(jié)構(gòu)的擾動波長與此相當(dāng)。根據(jù)上述觀點，在關(guān)注區(qū)域內(nèi)Δ0分別取2×10-3c、10-3c和0.5×10-3c，構(gòu)造粗、中、細(xì)3套計算網(wǎng)格，以對網(wǎng)格密度進(jìn)行驗證。近壁面法向首層網(wǎng)格高度維持10-5c，滿足近壁面y+<1。三套網(wǎng)格計算域內(nèi)網(wǎng)格結(jié)點總數(shù)約為0.5×107、2.0×107、9.2×107，流動關(guān)注區(qū)域附近網(wǎng)格對比情況如圖12所示。

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

圖12 流動關(guān)注區(qū)域附近計算網(wǎng)格對比
Fig.12 Comparison of computational grids near focus regions

5 結(jié)冰翼型計算結(jié)果及分析

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

以下首先給出網(wǎng)格無關(guān)性驗證結(jié)果；在此基礎(chǔ)上對干凈無冰與縫翼結(jié)冰結(jié)果進(jìn)行了對比分析；并研究了不同迎角下分離泡及其尾跡流動的演化過程與影響效應(yīng)。

5.1 網(wǎng)格無關(guān)性驗證

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

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

5.2 干凈無冰與縫翼結(jié)冰結(jié)果

在上述中等密度網(wǎng)格計算結(jié)果基礎(chǔ)上進(jìn)行長時間平均續(xù)算。以下從時均和瞬態(tài)角度對縫翼結(jié)冰翼型計算結(jié)果與干凈翼型同一計算條件下的相應(yīng)結(jié)果進(jìn)行對比分析，以期揭示縫翼結(jié)冰對多段翼型大迎角流場與氣動特性的影響機制。

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

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

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

圖16給出了結(jié)冰前后縫翼附近時均流場速度分布U/UINF對比情況。由圖可知，在干凈無冰條件下，縫翼上表面流動處于完全附著狀態(tài)；當(dāng)縫翼前緣存在積冰時，形成了尺度與縫翼弦長相當(dāng)?shù)臅r均分離泡，使外部繞流速度大小有所增長，但主翼前緣附近卻表現(xiàn)出了較為明顯的速度損失。此外，由于縫翼下表面存在積冰，改變了凹腔位置的當(dāng)?shù)貋砹鞣较颍罨亓鳒u大小有一定程度地增加，但結(jié)構(gòu)基本保持不變。

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

圖15 時均流場壓力系數(shù)分布對比
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時均計算結(jié)果可知，起始監(jiān)測點1附近結(jié)冰前后速度分布情況大體相同，表明分離流動對縫道出口流動影響較小。從監(jiān)測點2開始出現(xiàn)了明顯的速度損失。結(jié)冰前后各監(jiān)測點處流動速度量值沿弦向的變化趨勢比較相似，速度型斜率基本相同，并且壁面附近速度分布基本維持不變。這表明結(jié)冰導(dǎo)致的速度損失在弦向和法向空間分布上相對均勻，與當(dāng)?shù)貋砹鳡顩r的變化直接相關(guān)。

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

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

圖19 瞬態(tài)空間流線分布對比
Fig.19 Comparison of instantaneous streamlines distribution

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

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

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

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

5.3 縫翼結(jié)冰不同迎角結(jié)果

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

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

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

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

表1 時均升力系數(shù)對比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給出了不同迎角下瞬態(tài)流場無量綱展向渦量z-vorticity的分布情況。對比各圖可知，在所研究的迎角變化范圍內(nèi)，剪切層脫落渦能夠在縫翼后緣附近向壁面方向輸運，分離泡結(jié)構(gòu)始終能得到維持。不同迎角下分離泡瞬態(tài)尾跡流場的基本結(jié)構(gòu)較為類似，仍由一系列較為規(guī)則的渦串組成；且尾跡旋渦運動特征并不存在定性變化。隨著迎角增加，剪切層脫落角相應(yīng)增加，弦向失穩(wěn)位置提前；失穩(wěn)后析出的旋渦尺度存在增長趨勢，強度有所提高；尾跡渦渦核與主翼壁面距離增加，尾渦列具備向外部自由來流流動區(qū)域移動的趨勢，這與圖24體現(xiàn)的近壁速度分布相吻合。

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

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

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

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

6 結(jié) 論

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

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

3) 縫翼前緣結(jié)冰后將形成大尺度分離泡結(jié)構(gòu)，其非定常尾跡將對主翼前緣流動產(chǎn)生強烈干擾，導(dǎo)致縫道流動對主翼增升效應(yīng)的下降。

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

[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.

(責(zé)任編輯： 李明敏)

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">網(wǎng)絡(luò)出版時間

時間：2016-11-211439

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

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

.E-maillijieruihao@163.com

張恒， 李杰， 龔志斌． 多段翼型縫翼前緣結(jié)冰大迎角分離流動數(shù)值模擬J． 航空學(xué)報,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