趙賓賓，張恒，李杰

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

2．中國商用飛機有限責任公司 上海飛機設(shè)計研究院，上海 201210

3．清華大學 航天航空學院，北京 100084

結(jié)冰是民機飛行安全的主要威脅之一。翼面結(jié)冰狀態(tài)下的氣動特性評估是適航取證驗證、結(jié)冰防護研究、容冰安全飛行的重要依據(jù)［1-2］。前緣復(fù)雜分離是翼型結(jié)冰狀態(tài)氣動特性影響研究的核心問題［3］。此時分離流場由剪切層失穩(wěn)和多尺度旋渦運動支配，其本質(zhì)是幾何間斷-壓力梯度共同作用下的分離-再附過程，具備較強的非線性/非定常特征，將造成最大升力系數(shù)損失、失速迎角提前、抗擾動能力下降、進入失速狀態(tài)的可能性大幅提高。相關(guān)流動機理的認識和總結(jié)仍然不盡完備，流動特性改變與氣動力變化之間的關(guān)聯(lián)性也并不明確。很大程度上對結(jié)冰致災(zāi)機理的深入研究構(gòu)成了嚴重制約。

目前翼型結(jié)冰后流場和氣動力的相關(guān)研究手段涵蓋飛行試驗、風洞試驗和數(shù)值模擬。飛行試驗可信程度較高，是確定結(jié)冰前后氣動特性變化情況的重要依據(jù)；但存在耗費巨大、危險性高、流場信息相對匱乏的問題。風洞試驗可在降低研究成本和風險的前提下，得到較為貼近實際的氣動數(shù)據(jù)；但存在測量參數(shù)有限、數(shù)據(jù)相關(guān)性需要修正的問題。相對而言，數(shù)值模擬能夠獲取較為全面的流場信息，實現(xiàn)相對簡捷和經(jīng)濟，能夠與氣動力設(shè)計分析流程緊密結(jié)合，是評估結(jié)冰影響、剖析流動機理的有效方法。

翼型前緣結(jié)冰導(dǎo)致的失速分離形態(tài)及機制與干凈翼型大迎角狀態(tài)存在本質(zhì)差異。干凈翼型失速分離特征一般與流向逆壓梯度直接相關(guān)。而結(jié)冰翼型的分離生成通常是由于積冰當?shù)貛缀翁卣靼l(fā)生突變，使得剪切層直接失穩(wěn)脫落，剪切層渦系結(jié)構(gòu)的發(fā)展促進了回流區(qū)域和主流區(qū)域的流動摻混效應(yīng)，最終在下游與壁面相互作用、發(fā)生再附現(xiàn)象，形成湍流分離泡。因此，驅(qū)動分離流動的實質(zhì)因素是幾何間斷觸發(fā)Kelvin-Helmholtz（K-H）不穩(wěn)定性和迎角效應(yīng)產(chǎn)生壓力梯度兩者的綜合效應(yīng)，主導(dǎo)失速分離特征的核心是分離泡結(jié)構(gòu)的生長、擴展和膨脹現(xiàn)象，這無疑大大增加了分離流動數(shù)值模擬分析的復(fù)雜度。

近年來，隨著計算流體力學（Computational Fluid Dynamics， CFD）方法的不斷進步，一系列新型數(shù)值模擬方法得到了發(fā)展，獲得了許多有價值的結(jié)論和數(shù)據(jù)，為結(jié)冰狀態(tài)下的飛行安全與防護研究提供了重要的理論儲備與技術(shù)支撐。然而在Bragg［3］和Lynch［4］等知名學者發(fā)表的有關(guān)綜述中，相對注重結(jié)冰后分離流動問題的物理機制探討，而缺乏從數(shù)值模擬方法角度出發(fā)的總結(jié)性分析。雖然Stebbins等［5］于2019年業(yè)已完成了關(guān)于結(jié)冰翼面氣動特性計算評估手段的綜述性文章，但未針對湍流模擬方法的應(yīng)用效果進行詳細分析和評論。鑒于在描述分離流場結(jié)構(gòu)及其演化特征層面的重要地位，仍然有必要進一步從湍流模擬方法的角度，針對結(jié)冰后氣動特性預(yù)測及分離流場特征描述等方面的沿革進展進行評述；不僅為結(jié)冰影響分析手段的選取及應(yīng)用提供參考依據(jù)，也為模擬方法的修正及完善提供評判準則。

本文首先對翼型結(jié)冰后數(shù)值模擬研究的基本框架進行介紹；在此基礎(chǔ)上從典型湍流模擬方法的實際應(yīng)用層面，回顧近年來數(shù)值方法在結(jié)冰翼型失速特性預(yù)測、分離流場特征刻畫等方面的進展；最后對研究手段和研究內(nèi)容的相關(guān)發(fā)展趨勢進行總結(jié)和展望。

1 結(jié)冰分離流動模擬概述

由于冰生成過程相對非定常流場演化特征而言具備較大的時間尺度，在目前的翼型結(jié)冰分離流場數(shù)值模擬工作中，通常將結(jié)冰過程研究和流場分析研究解耦，采用結(jié)冰試驗或數(shù)值模擬研究中某時刻獲得的冰形建立幾何模型。伴隨幾何重構(gòu)技術(shù)的發(fā)展，大致經(jīng)歷了二維簡化冰形—二維真實冰形—考慮三維效應(yīng)的展向拉伸冰形—具備完全三維特征的真實冰形這一發(fā)展過程。得益于多塊結(jié)構(gòu)化網(wǎng)格、非結(jié)構(gòu)混合網(wǎng)格和嵌套網(wǎng)格等技術(shù)的不斷完善，目前已能夠針對復(fù)雜結(jié)冰外形生成大規(guī)模貼體計算網(wǎng)格，為精細化數(shù)值模擬方法特別是新型湍流模擬方法的成功應(yīng)用建立了良好的基礎(chǔ)。

鑒于分離流動的高度復(fù)雜性和對氣動特性的強烈影響，數(shù)值模擬中通常選取角狀冰和冰脊作為典型冰形。圖1［3］給出了翼型前緣角狀冰誘導(dǎo)產(chǎn)生的瞬時流場結(jié)構(gòu)示意圖（圖中α為迎角，U∞為來流速度）。圖中顯示在K-H不穩(wěn)定性的驅(qū)動下，冰角頂端的剪切層失穩(wěn)伴隨復(fù)雜的多尺度旋渦脫落和輸運過程，促進了外部高速流動與近壁面回流流動之間的混合，生成幾何尺度較大的分離泡結(jié)構(gòu)。在壓力梯度的作用下，這種復(fù)雜流動結(jié)構(gòu)通常在翼型失速點附近表現(xiàn)出較強的非定常特征，不僅使得分離泡的幾何特征參數(shù)隨時間變化，并且關(guān)于來流擾動變化的拉伸-膨脹效應(yīng)及演化機制高度復(fù)雜，進而導(dǎo)致宏觀氣動力的高度不穩(wěn)定，是當前翼型結(jié)冰狀態(tài)分離流動數(shù)值模擬研究關(guān)注的重點內(nèi)容。

圖1 角狀冰后方瞬時流場結(jié)構(gòu)示意圖［3］Fig.1 Schematic diagram of transient flow field struc?ture behind horn ice［3］

就氣動特性分析而言，數(shù)值模擬研究重點關(guān)注結(jié)冰翼型失速迎角的確定、最大升力系數(shù)及失速形態(tài)的變化情況。顯而易見地，翼型失速特性的準確預(yù)測依賴于把握上述前緣分離泡的生成演化過程，這對湍流模擬方法的精度提出了較高要求。隨著雷諾平均（Reynolds Averaged Navier-Stokes， RANS）方法、大渦模擬（Large Eddy Simulation， LES）方 法 和RANS/LES混合方法的相繼建立和完善，數(shù)值模擬對湍流流場的描述更加準確，獲得的流動特征量更加豐富，能夠更好地反映分離流動的非定常本質(zhì)。以下從不同類型湍流模擬方法的應(yīng)用角度，對近年來數(shù)值模擬研究的進展情況進行簡要評述。

2 RANS方法

RANS方法將瞬時流動變量分解為平均量與脈動量兩部分，使得方程中出現(xiàn)了額外的雷諾應(yīng)力張量，從而需要引入某種假設(shè)對該變量進行估計。依靠對湍流內(nèi)在演化機制的描述，能夠?qū)字Z應(yīng)力做出各種假設(shè)，建立附加方程模型從而刻畫湍流平均量。對于翼型結(jié)冰后流動的計算分析，通常采用的湍流模型包括一方程Spalart-Allmaras（S-A）模型［6］、兩方程k-ε模型［7］、k-ω模型［8］和Shear-Stress-Transport（SST）模型［9］等。采用的幾何模型涵蓋二維到三維，允許使用單元尺度較大的計算網(wǎng)格，相應(yīng)的空間離散格式精度要求也較低，可基于定常或非定常方式進行流場求解，對計算資源的要求較易滿足，因此在數(shù)值模擬研究的早期階段及現(xiàn)階段工程領(lǐng)域得到了廣泛應(yīng)用。

應(yīng)用RANS分析翼型結(jié)冰后氣動特性變化的工作始于20世紀80年代。Potapczuk和Gerhart［10］首先應(yīng)用薄層RANS方法進行了一些基本的結(jié)冰 翼 型 氣 動 特 性 計 算。Kwon和Sankar［11-12］對 翼型結(jié)冰引起的局部流動分離現(xiàn)象進行了數(shù)值模擬。Shim等［13］初步研究了冰形幾何形狀同翼型氣動特性變化間的相關(guān)性，就不同湍流模型的效果進行了分析。Chi等［14-15］基于圖2所示的二維結(jié)構(gòu)網(wǎng)格，分別針對霜冰和明冰條件下的GLC305翼型進行了宏觀氣動力計算，就湍流模型對計算結(jié)果的影響進行了對比。對于圖3（a）［14］給出的霜冰算例（圖中CL為升力系數(shù)），不同湍流模型在升力線性段附近獲得的結(jié)果基本類似，與試驗吻合程度良好；但最大升力系數(shù)和失速迎角均與試驗值存在一定差距。而對于圖3（b）［14］給出的明冰翼型，湍流模型的影響擴展到升力線性段，失速點附近的計算結(jié)果差別較大，均無法較好地反映升力變化的基本特征。就該算例而言，SST等復(fù)雜湍流模型的效果并不明顯優(yōu)于S-A模型。圖4［15］給出的明冰翼型流向速度（u/U∞）分布對比情況表明，S-A模型未能在失速臨界迎角附近準確預(yù)測分離泡的基本形態(tài)，獲得的再附位置較試驗結(jié)果明顯偏后，反映了常用湍流模型在刻畫結(jié)冰誘導(dǎo)分離流動特征層面的固有缺陷。

圖2 明冰冰形附近二維結(jié)構(gòu)網(wǎng)格［14］Fig.2 Two-dimensional structured grid near glaze ice shape［14］

圖3 翼型升力系數(shù)計算結(jié)果對比［14］Fig.3 Comparison of computed results of lift coeffi?cients for airfoil［14］

圖4 明冰翼型流向速度分布對比（α=6°，上：S-A模型計算結(jié)果；下：試驗結(jié)果）［15］Fig.4 Comparison of streamwise velocity contours of airfoil with glaze ice（α=6°，top： result of S-A model；bot?tom： result of test）［15］

Pan和Loth［16］就 簡 化 冰 形 對 不 同 翼 型 氣 動力造成的影響進行了較為系統(tǒng)的分析，歸納了來流馬赫數(shù)和雷諾數(shù)效應(yīng)，分析了冰形弦向位置與壓力分布形態(tài)間的關(guān)聯(lián)，同樣表明小迎角狀態(tài)下取得的數(shù)值結(jié)果基本可用，但同時也指出RANS在失速點附近無法準確計算氣動力，特別對于幾何尺寸較大的冰形而言，由于分離強度顯著增加，數(shù)值結(jié)果與試驗值間的差異不容忽視。

Marongiu等［17］針對層流翼型NLF0414結(jié)典型雙角狀明冰后的氣動特性改變問題，在不同求解器的基礎(chǔ)上對比了S-A模型和SST模型的模擬 效 果。圖5（a）［17］給 出 的 結(jié) 果 表 明RANS在 小迎角下獲得的壓力分布情況均比較滿意，但由于失速點附近的壓力恢復(fù)過程和旋渦空間輸運過程直接相關(guān)，圖5（b）［17］所示的計算結(jié)果與試驗值之間出現(xiàn)了較大偏差。該算例中SST模型在分離流動較強時取得的氣動力預(yù)測結(jié)果優(yōu)于S-A模型。此外，文獻［17］還指出大迎角下即使采用非定常RANS也難以有效描述流動的分離特性，因此借助RANS/LES混合方法或LES方法等更為精細的湍流模擬方法進行流場求解是必要的。

圖5 壓力系數(shù)Cp分布對比（左：S-A模型；右：SST模型）［17］Fig.5 Comparison of pressure coefficientCpdistributions （left： S-A model； right： SST model）［17］

Jun等［18］通 過 激 光 掃 描 技 術(shù) 建 立 了NACA 23012翼型前緣真實冰形的 三 維 數(shù) 模。圖6［18］體現(xiàn)了非結(jié)構(gòu)網(wǎng)格較強的幾何外形適應(yīng)性?；赗ANS方法獲得了物面壓力分布和升力特性曲線等基本氣動數(shù)據(jù)。圖7［18］表明由于幾何模型復(fù)雜度的提高，即使在小迎角下，壓力分布的計算精度也并不十分滿意。Mirzaei等［19］基于兩方程k-ε模型分析了NLF0414翼型角狀冰影響下的宏觀分離流場演化過程，在中小迎角下獲得了與試驗結(jié)果吻合良好的再附位置/湍流強度極值等關(guān)鍵特征量。

圖7 α=0°時模型中面壓力分布對比［18］Fig.7 Comparison of midspan pressure distributions of model atα=0°［18］

近年來國內(nèi)研究者也相繼開展了一些有代表性的研究工作。陳科等［20-21］基于近壁面非結(jié)構(gòu)網(wǎng)格-遠場結(jié)構(gòu)化網(wǎng)格的思路，發(fā)展了一種混合網(wǎng)格生成方法，避免了復(fù)雜結(jié)冰外形難以劃分結(jié)構(gòu)化網(wǎng)格的問題；針對一類具備不規(guī)則前緣冰形的翼型，分別基于S-A模型和SST模型進行了氣動特性計算分析，較好地預(yù)測了翼型失速點之前的升力特性變化情況。李焱鑫等［22］結(jié)合雷諾應(yīng)力模型（Reynolds Stress Model， RSM）［23］，針對大型客機典型超臨界翼型及平尾翼型帶溢流冰條件下的氣動力進行模擬分析，表明該模型可較為準確地預(yù)測結(jié)冰翼型的最大升力系數(shù)與失速迎角。Li等［24-25］在考慮湍流非平衡特性的前提下發(fā)展了2種分離修正湍流模型，有效提升了翼型結(jié)冰狀態(tài)失速分離流場的預(yù)測精度。黃冉冉等［26］基于S-A模型量化分析了具備不同微觀特征參數(shù)的冰形對翼型失速特性的影響規(guī)律。

以上結(jié)果表明，基于簡化的結(jié)冰幾何構(gòu)型，在分離流動強度有限的小迎角及臨界迎角條件下，RANS能夠取得較好的宏觀流場和氣動力數(shù)值結(jié)果。因而在機理研究的初始階段及工程應(yīng)用領(lǐng)域，不失為一種能夠有效描述分離流場概貌、高效獲取氣動數(shù)據(jù)的方法［27］。但是，由于結(jié)冰翼型失速點附近存在豐富的旋渦生成演化及相互作用，與湍流脈動相關(guān)的非定常效應(yīng)占據(jù)主導(dǎo)地位；而RANS對脈動信息采用完全?；奶幚矸绞?，難以合理描述復(fù)雜流動細節(jié)，能夠得到的流動特征量比較有限；并且對于不同計算條件和結(jié)冰外形，同一湍流模型的氣動力預(yù)測準確程度也存在差別。因此，在結(jié)冰誘導(dǎo)分離流動的機理研究工作中，仍然有必要應(yīng)用更為符合物理事實、精度更高的湍流模擬方法。

3 LES方法

LES認為湍動能由大尺度運動主導(dǎo)，難以基于普適性模型描述；耗散則由小尺度運動主導(dǎo)，且各向同性效應(yīng)顯著。因此對大尺度運動進行直接求解，而小尺度運動則采用?；幚?，尺度劃分通過流動變量在空間上的加權(quán)積分完成。這種處理方式使得方程中出現(xiàn)了亞格子應(yīng)力張量，需要進行?；幚?。在遠離壁面的分離區(qū)域，LES能夠有效給出湍流瞬時脈動信息，較為精確地描述流動細節(jié)特征。對于翼型結(jié)冰分離流動問題而言，該方法應(yīng)用的計算模型一般為具備一定展向長度的準三維模型，要求計算網(wǎng)格在分離區(qū)域具備較高的各向同性性質(zhì)且足夠精細。同時還需要配合高精度空間離散格式，并以較小時間步長進行非定常推進，因而對計算資源的要求總體較高。

目前應(yīng)用LES方法分析翼型結(jié)冰后分離流動的工作還比較有限。Brown等［28］應(yīng)用Implicit LES（ILES）方法［29］，針對一類結(jié)冰旋翼翼型開展了數(shù)值模擬研究。在三維冰形的構(gòu)建過程中應(yīng)用 了 高 精 度CT掃 描 技 術(shù)，圖8［28］顯 示 該 冰 形 具備結(jié)構(gòu)復(fù)雜、隨機性較強和不規(guī)則度較高的特點。圖9［28］以Q等值面［30］的形式給出了翼型不同迎角下的瞬時流場結(jié)構(gòu)，表明ILES較為精細地描述了剪切層失穩(wěn)和旋渦脫落過程隨迎角的變化情況，有效地刻畫了空間三維小尺度湍流結(jié)構(gòu)。但圖10［28］顯示計算所得的翼型失速特性與試驗值之間仍存在差距。表明僅僅增加湍流模擬方法復(fù)雜度并不能有效改善結(jié)冰翼型失速點附近氣動力的預(yù)測結(jié)果。

圖8 三維冰形表面非結(jié)構(gòu)網(wǎng)格［28］Fig.8 Unstructured surface mesh of three-dimensional ice shape［28］

圖9 不同迎角下瞬時流場Q等值面［28］Fig.9 Iso-surface ofQ-criterion of instantaneous flow structure of different angles of attack［28］

圖10 升力特性曲線對比［28］Fig.10 Comparison of lift characteristic curves［28］

LES壁 面 模 型（Wall-Modelling in LES，WMLES）［31］在 非 常 靠 近 壁 面 的 邊 界 層 內(nèi) 體 現(xiàn)RANS特征，而在其余區(qū)域則體現(xiàn)經(jīng)過修正的LES方法基本性質(zhì)。Xiao M C等［32］結(jié)合中心-迎風混合格式，基于WMLES系統(tǒng)分析了GLC305及NLF0414翼型結(jié)冰狀態(tài)下的分離演化規(guī)律，表明數(shù)值方法能夠準確刻畫自由剪切層的K-H不穩(wěn)定性，在此基礎(chǔ)上就分離流場的非定常頻域特征進行了深入探討。圖11［32］給出了冰角后方剪切層失穩(wěn)后二維渦管卷起生成三維大尺度類發(fā)卡渦的演化過程。但是，由于結(jié)冰分離流場中存在典型再附特征；就小尺度、高頻率湍流結(jié)構(gòu)主導(dǎo)的附著邊界層流動而言，LES在近壁模型的構(gòu)造和完善方面尚存在一些缺陷，附著區(qū)域的合理界定和判斷、?；?解析區(qū)域的快速轉(zhuǎn)換仍然是制約LES廣泛應(yīng)用的固有難題。因此，LES在翼型結(jié)冰分離流動分析中的推廣仍有相當多的前置工作需要開展。

圖11 NLF0414結(jié)冰翼型剪切層渦系結(jié)構(gòu)［32］Fig.11 Vortex structures of shear layer of NLF0414 iced airfoil［32］

4 RANS/LES混合方法

由于RANS方法和LES方法對湍流附加應(yīng)力項的描述能夠以相似的形式給出；因此可以采取較為簡潔的手段對兩者加以混合。其基本思想是在壁面附近流動區(qū)域?；飨蛲缘男〕叨韧牧鹘Y(jié)構(gòu)，體現(xiàn)RANS性質(zhì)，在降低計算資源需求的同時，維持壁面附近流動信息的求解準確程度，避免雷諾應(yīng)力損失；在遠離壁面的分離區(qū)域則對大尺度旋渦進行解析，體現(xiàn)LES性質(zhì)，保證分離區(qū)域湍流運動的模擬精度。該方法要求分離區(qū)域網(wǎng)格三向同性并且劃分足夠細密，但在近壁面可使用長寬比較大的薄層網(wǎng)格［33］。一般采用多塊結(jié)構(gòu)化網(wǎng)格或結(jié)構(gòu)/非結(jié)構(gòu)混合網(wǎng)格，以滿足不同性質(zhì)湍流模擬方法對網(wǎng)格布置的需求。對空間離散格式和非定常計算時間步長的要求較LES略低，所需的計算資源總體高于非定常RANS但低于LES［34］，并且易于程序?qū)崿F(xiàn)，因而成為了近年來的研究熱點，開始逐步應(yīng)用于結(jié)冰后復(fù)雜分離流動的計算分析工作中。

一般根據(jù)構(gòu)造形式將RANS/LES混合方法歸納為兩大類，即湍流量混合法和混合界面法。湍流量混合法的基本思想是采用加權(quán)思想對RANS雷諾應(yīng)力項和LES亞格子應(yīng)力項進行混合；其重點是如何構(gòu)建合理的加權(quán)混合函數(shù)?；旌辖缑娣▽⒘鲌龇譃镽ANS區(qū)域和LES區(qū)域分別進行計算。令湍流模擬方法在部分流場區(qū)域（例如近壁面區(qū)域和遠場區(qū)域）內(nèi)體現(xiàn)RANS性質(zhì)，而在其余流場區(qū)域（例如流動分離區(qū)域）內(nèi)則體現(xiàn)LES性質(zhì)。采用這種方法對流場進行數(shù)值模擬，則在2個流動區(qū)域之間必然存在分界面。該方法的重點是對分界面進行合理界定，并保持分界面兩側(cè)流動變量的光滑過渡［35］。

混合界面法中具備代表性的一類方法是Spalart等［36］于1997年 提 出 的 分 離 渦 模 擬 方 法（Detached Eddy Simulation， DES/DES97）。該方法將S-A湍流模型和LES亞格子應(yīng)力模型結(jié)合為統(tǒng)一模型，通過比較當?shù)鼐W(wǎng)格尺度與湍流混合長度進行方法轉(zhuǎn)換。在近壁面區(qū)域表現(xiàn)為傳統(tǒng)RANS方法；在以大渦輸運為主要特征的區(qū)域快速降低當?shù)乩字Z應(yīng)力水平，體現(xiàn)LES方法的特點，解析求解大尺度空間湍流結(jié)構(gòu)。但采用該方法會出現(xiàn)模化應(yīng)力損失（Modeled Stress Depletion，MSD）［37］問題，引起網(wǎng)格誘導(dǎo)分離（Grid Induced Separation， GIS）。造成這種現(xiàn)象的主要原因是壁面附近區(qū)域的網(wǎng)格單元各向尺度大致相當時，DES97方法會將此區(qū)域標定為LES區(qū)域，此時壁面邊界層內(nèi)雷諾應(yīng)力模化不足，計算過程中邊界層內(nèi)渦黏性將會過度減少。為了解決MSD問題，Menter和Kuntz［38］通 過 引入SST湍流模型中 的 過渡函數(shù)，實現(xiàn)了壁面邊界層到主流區(qū)域的漸進轉(zhuǎn)換。Spalart等［37］基于相似的思想，采用“延遲LES函數(shù)”使LES亞格子應(yīng)力模型推遲進入壁面邊界層區(qū)域，此類方法稱為延遲DES方法（Delayed DES， DDES）。通 過 將DDES與WMLES相 結(jié)合，Shur等［39］提出了IDDES方法。IDDES對亞格子尺度進行了重新定義，并且在同時考慮RANS和LES長度尺度影響的前提下構(gòu)造了RANS/LES混合長度，一定程度上解決了DES類方法直接應(yīng)用于WMLES時產(chǎn)生的“對數(shù)層不匹配”（Log-Layer Mismatch，LLM）問題。

Pan和Loth［40-41］較 早 地 利 用DES方 法 對 帶簡化冰脊模型的NACA23012翼型繞流進行了分析，冰形附近的多塊結(jié)構(gòu)網(wǎng)格如圖12［41］所示。結(jié)果表明在失速點附近，DES獲得的宏觀氣動力結(jié)果較RANS方法有所提升。圖13［41］顯示DES計算所得的壓力恢復(fù)過程明顯優(yōu)于RANS；但由于分離區(qū)域計算網(wǎng)格量較為有限，只獲得了圖14［41］所示的旋渦結(jié)構(gòu)概貌，但初步顯示了DES解析空間分離區(qū)域大尺度湍流特征的能力。

圖12 簡化冰形附近多塊結(jié)構(gòu)網(wǎng)格［41］Fig.12 Multi-block structured grid near simplified ice shape［41］

圖13 壓力分布對比（α=5°）［41］Fig.13 Comparison of predicted pressure distributions（α=5°）［41］

圖14 瞬時流場展向渦量分布云圖［41］Fig.14 Contours of spanwise vorticity of instantaneous flow field［41］

Thompson［42］和Mogili［43］等 基 于 非 結(jié) 構(gòu) 混 合網(wǎng)格，利用DES對GLC305翼型結(jié)角狀冰后不同迎角下的流場和氣動力變化問題做了較為全面的分析研究。提出了將冰角高度與關(guān)注區(qū)域單元尺度相互關(guān)聯(lián)的網(wǎng)格設(shè)計思路，就數(shù)值方法的網(wǎng)格敏感性開展了系統(tǒng)的對比分析，基準網(wǎng)格和加密網(wǎng) 格 空 間 分 布 對 比 情 況 如 圖15［42］所 示。圖16［42］和圖17［42］給出的失速點附近壓力和速度分布計算結(jié)果表明，由于DES有效描述了剪切層渦系演化過程的時間累積效應(yīng)，因此不僅能夠較為準確地描述壓力分布的平臺特征及其恢復(fù)過程，同時能夠更好地預(yù)測時均分離泡的再附位置和其幾何形態(tài)。圖18［42］給出的湍流強度分布情況也與風洞試驗結(jié)果基本類似，表明DES在分離區(qū)域能夠有效切換到LES模式，充分解析湍流脈動信息。

圖15 冰形附近混合網(wǎng)格［42］Fig.15 Hybrid mesh near ice shape［42］

圖16 壓力分布對比（α=6°）［42］Fig.16 Comparison of predicted pressure distributions（α=6°）［42］

圖17 時均流向速度分布對比（α=6°）［42］Fig.17 Comparison of time-averaged streamwise ve?locity distributions（α=6°）［42］

圖18 流向速度脈動均方根對比（α=6°）［42］Fig.18 Comparison of root-mean-square of fluctuations in streamwise velocity component （α=6°）［42］

圖19 時均壓力分布結(jié)果對比（α=8°）［44］Fig.19 Comparison of time-averaged pressure distribu?tions （α=8°）［44］

Lorenzo等［44］分 別 采 取DES和DDES就 平尾翼型M5-6結(jié)冰失速前后的氣動力變化情況進行了分析，認為失速點之前的宏觀氣動力計算結(jié)果是比較滿意的。圖19［44］表明在臨界失速條件下，DDES較DES能夠更為準確地反映結(jié)冰翼型的時均壓力分布平臺特征，集中體現(xiàn)于描述了更為和緩而非相對急促的壓力恢復(fù)過程，對應(yīng)更加符合物理事實、更為靠近下游的旋渦結(jié)構(gòu)生成及摻混現(xiàn)象；其主要原因是DDES通過引入延遲函數(shù)，避免了LES亞格子尺度過度侵入RANS區(qū)，從而減緩了原始DES存在的?；瘧?yīng)力衰減、K-H不穩(wěn)定性過度釋放問題。Lakshmipathy和To?giti［45］對 比 了PANS （Partially Averaged Navier-Stokes）［46］方 法 和DDES方 法 對NACA23012翼型結(jié)冰后強分離流動的模擬效果，圖20［45］和圖21［45］表明DDES能夠提供更好的時均壓力分布結(jié)果，同時能夠更為充分地解析尾跡區(qū)域的三維湍 流 結(jié)構(gòu)。Alam等［47-48］在GLC305翼型 臨 界失速狀態(tài)算例中也在一定程度上肯定了DDES的流場細節(jié)刻畫能力。Molina等［49］采用DDES和非結(jié)構(gòu)網(wǎng)格開展了基于SU2程序的GLC305結(jié)冰翼型臨界失速分離流場分析研究，肯定了引入自適應(yīng)亞格子尺度對預(yù)測結(jié)果的改善作用。

圖20 時均壓力分布結(jié)果對比［45］Fig.20 Comparison of time-averaged pressure distributions［45］

圖21 瞬時流場Q等值面對比［45］Fig.21 Comparison of iso-surface ofQ-criterion of in?stantaneous flow field［45］

但是，上述研究工作中同時也指出DES類方法對結(jié)冰翼型失速后強分離流場的預(yù)測精度仍存在一定提升空間。其中的首要問題是冰角后方RANS/LES區(qū)域的延遲轉(zhuǎn)換特征與剪切層失穩(wěn)過程解析需求的高度不匹配，由此生成的強亞格子黏性嚴重制約了K-H不穩(wěn)定性的合理釋放，預(yù)測得到的三維旋渦結(jié)構(gòu)生成位置通常較試驗結(jié)果更為接近下游，使得剪切層擺動-失穩(wěn)、二維渦管釋放-扭曲、發(fā)卡渦結(jié)構(gòu)發(fā)展變化等典型剪切層相關(guān)精細湍流特征難以得到充分辨識，制約了分離流動機理特別是非定常特征的深入剖析。

近年來，IDDES方法逐漸成為了結(jié)冰分離特征精細預(yù)測的重要手段。國內(nèi)Xiao Z X等［50-52］率先基于該方法開展了一系列卓有成效的研究工作，表明該方法在強分離區(qū)域具備與DDES相當?shù)耐牧髂M精度，同時在附著和分離并存的過渡區(qū)域能夠取得更為滿意的計算結(jié)果，因此適用于分析剪切層分離-再附效應(yīng)主導(dǎo)的結(jié)冰失速分離問題。在該方法的基礎(chǔ)上，翼型結(jié)冰后分離流場的研究范疇得到了進一步拓展，不僅包含傳統(tǒng)意義上的臨界失速分離泡模擬分析，并且涵蓋了失速點后分離泡的膨脹機制、剪切層失穩(wěn)相關(guān)的非定常模態(tài)特征等領(lǐng)域，國內(nèi)在該領(lǐng)域積累了較為豐富的研究經(jīng)驗：張恒等［53-54］采用該方法分析了GLC305結(jié)冰翼型失速階段分離泡的拉伸-膨脹歷程，指出發(fā)卡渦抬升導(dǎo)致剪切層渦系摻混融合作用降低是驅(qū)動結(jié)冰翼型失速分離流場演化的關(guān)鍵機制；Hu等［55-56］厘清了角冰/冰脊誘導(dǎo)剪切層的失穩(wěn)及旋渦脫落模式，初步建立了流場關(guān)鍵模態(tài)與不同尺度旋渦之間的關(guān)聯(lián)；Bao等［57］將該方法延拓到結(jié)冰翼型氣動噪聲特性的分析研究工作當中；譚雪等［58-59］基于該方法進一步識別了冰脊觸發(fā)分離的各階模態(tài)特征，表明大尺度旋渦結(jié)構(gòu)是導(dǎo)致翼型臨界失速狀態(tài)升力波動的根本原因。

由 于 傳 統(tǒng)DES方 法 通 常 采 用Δ=max(Δx，Δy，Δz)（Δx、Δy、Δz分別為x、y、z向的計算網(wǎng)格間距），即Δmax作為亞格子模型。當計算網(wǎng)格嚴格各向同性時，該亞格子尺度等同于經(jīng)典LES尺度，即1/3網(wǎng)格單元體積的立方根。但對于以近壁區(qū)域為代表的典型各向異性網(wǎng)格而言，平行物面方向單元尺度一般遠大于法向單元尺度，顯然上述定義將會導(dǎo)致當?shù)亻L度尺度過大，造成當?shù)販u黏特征過度預(yù)測，降低了DES類方法對湍流精細結(jié)構(gòu)的解析能力。針對上述問題，Shur等［60］提出了一種綜合考慮計算網(wǎng)格當?shù)貛缀翁卣髋c湍流流場信息的亞格子模型，利用湍流流場信息指示準二維流動特征，能夠較為合理地降低各向異性網(wǎng)格/準二維分離區(qū)域的亞格子黏性，同時在各向同性網(wǎng)格/充分發(fā)展的三維分離區(qū)域恢復(fù)傳統(tǒng)LES亞格子模型的特征。國內(nèi)Xiao M C和Zhang［61］結(jié) 合 不 同 亞 格 子 應(yīng) 力 模 型［62］分 析 了GLC305及NLF0414翼型角冰/冰脊影響下的翼型失速分離流場，基于IDDES-SLA［60］方法精細刻畫了剪切層渦結(jié)構(gòu)的生成演化模式及空間壓力脈動規(guī)律，如圖22~圖24所示。由于SLA（Shear Layer Adapted）模型的基本思路為在各向異性較大的網(wǎng)格處采用各個網(wǎng)格節(jié)點位置加權(quán)處理的亞格子尺度與相關(guān)流場信息相乘的新亞格子尺度，因此具備在各向異性較強、剪切作用顯著的流動區(qū)域有效降低當?shù)販u黏性、提高K-H不穩(wěn)定性解析能力的特點，從而增強了冰角后方剪切層失穩(wěn)過程的刻畫精度。因此，亞格子模型或亞格子常數(shù)的合理修正及完善同樣是結(jié)冰翼型分離流場模擬方法改進階段值得關(guān)注的重點研究內(nèi)容之一。

圖22 IDDES-SLA方法捕獲的多尺度旋渦結(jié)構(gòu)［61］Fig.22 Muti-scale vortices structure captured by IDDESSLA［61］

圖23 IDDES-SLA刻畫的剪切層渦系結(jié)構(gòu)生成演化規(guī)律［61］Fig.23 Generation and evolution of shear layer vortex structure described by IDDES-SLA［61］

圖24 不同亞格子模型捕獲的空間壓力脈動特征［61］Fig.24 Characteristics of spatial pressure fluctuations captured by different subgrid models［61］

Zonal DES（ZDES）方法是另一類有代表性的RANS/LES混合方法。該方法根據(jù)分離流場的基本結(jié)構(gòu)，在基準DES方法［28］的基礎(chǔ)上，以區(qū)域交界面的形式將流場人為劃分為RANS區(qū)域和DES區(qū)域，實現(xiàn)湍流模擬方法的混合［63］。Duclercq等［64］基 于 圖25所 示 的 多 塊 結(jié) 構(gòu) 化 網(wǎng) 格，利 用ZDES方法對展向冰脊誘導(dǎo)產(chǎn)生流場的非定常脈動特征進行了較為細致的研究。由于該方法能夠在冰角后方人為強制切換到LES模式，因而具備相對較好的K-H不穩(wěn)定性預(yù)測能力。由圖26［64］可知數(shù)值方法能夠較為準確地預(yù)測分離泡的渦核位置和時均速度分布情況。通過在冰脊后方剪切層脫落處設(shè)置如圖27（a）［64］所示的2個監(jiān)測點，給出了如圖27（b）［64］所示的鄰近冰角當?shù)厮俣让}動功率譜密度（PSD）。圖中給出的特征頻率f0表征了剪切層內(nèi)部旋渦的脫落頻率；而特征頻率f1可視為剪切層垂直運動的基頻。圖28［64］以Q等值面的形式顯示了冰脊后方的瞬時流場，表明ZDES對分離區(qū)域小尺度三維湍流結(jié)構(gòu)具備較強的解析能力。

圖25 冰脊附近多塊結(jié)構(gòu)化空間網(wǎng)格［64］Fig.25 Muti-block structured grid near ice ridge［64］

圖26 時均流向速度分布對比（上：PIV試驗結(jié)果；下：ZDES計算結(jié)果）［64］Fig.26 Comparison of time-averaged streamwise ve?locity contours （top： result of PIV； bottom： re?sult of ZDES）［64］

圖27 冰脊后方流動速度脈動功率譜密度［64］Fig.27 Power spectral densities of velocity fluctuations behind ice ridge［64］

Deck［63］通過對長度尺度項進行修正，實現(xiàn)了ZDES的 改 進。Zhang等［65］基 于 該 方 法 分 析 了GLC305和NACA23012兩類翼型不同結(jié)冰條件下的分離流動，指出計算域展向長度的選取將影響尾跡區(qū)域湍流結(jié)構(gòu)的精確解析。由圖29［65-66］可知，在展向長度增加到與分離泡長度接近的情況下，尾跡區(qū)域獲得的大尺度旋渦基本形狀和小尺度湍流結(jié)構(gòu)精細程度明顯高于展向長度較小的算例。Costes等［66-67］應(yīng)用嵌套網(wǎng)格對類似的結(jié)冰分離流動算例進行分析時，認為ZDES表現(xiàn)出了較強的網(wǎng)格敏感性，影響流動再附位置的判定和旋渦脫落過程的解析；這與RANS區(qū)域和DES區(qū)域的空間劃分位置、網(wǎng)格密度差異都有著密切聯(lián)系。

圖28 Q等值面顯示的冰脊后方瞬時流場［64］Fig.28 Iso-surface ofQ-criterion of instantaneous flow field behind ice ridge［64］

圖29 ZDES方法獲得的瞬時流場Fig.29 Instantaneous wake flow field of ZDES

Bhushan等［68-69］基于附加應(yīng)力項加權(quán)混合的思路發(fā)展了一種動態(tài)RANS/LES混合模型（Dy?namic Hybrid RANS-LES model， DHRL），Alam等［47-48］基于該方法在非結(jié)構(gòu)網(wǎng)格的基礎(chǔ)上針對GLC305結(jié)冰翼型臨界失速算例進行數(shù)值模擬，從近壁面速度分布和湍流脈動強度等流場細節(jié)角度驗證了數(shù)值方法的可靠性。由圖30［48］和圖31［48］給出的粗細2套網(wǎng)格計算結(jié)果對比情況可知，由于該方法對RANS方法和LES方法的混合方式與當?shù)鼐W(wǎng)格參數(shù)并不直接相關(guān)，而是根據(jù)當?shù)赝牧餍畔⒅苯优卸?，因此該混合模型對網(wǎng)格布置的依賴性相對較小，能夠在網(wǎng)格整體較為稀疏的情況下，給出較DDES更為良好的近壁流動細節(jié)預(yù)測結(jié)果。圖32［48］給出的流向速度等值面表明該方法對瞬時流場結(jié)構(gòu)的描述也是較為精細的，能夠較DDES更早地預(yù)測到因剪切層失穩(wěn)導(dǎo)致的旋渦析出和運動現(xiàn)象。但該方法的計算效率目前尚不十分滿意，相同推進步數(shù)下耗費的計算時間較DDES方法高出15%左右。

圖30 壁面附近平均流向速度型對比［48］Fig.30 Comparison of mean streamwise velocity pro?files near wall［48］

圖31 壁面附近流向速度脈動均方根對比［48］Fig.31 Comparison of root-mean-square of fluctuations in streamwise velocity component near wall［48］

圖32 瞬時流向速度等值面［48］Fig.32 Iso-surface of instantaneous streamwise velocity［48］

綜上所述，就結(jié)冰后翼型分離流動的精細化數(shù)值模擬而言，現(xiàn)階段RANS/LES混合方法能夠較為準確地預(yù)測失速分離流動的基本形態(tài)和多尺度結(jié)構(gòu)的演化過程，并對流場速度分布細節(jié)和湍流脈動情況進行合理描述。然而，目前此類方法在結(jié)冰分離流動問題研究中所涉及的結(jié)冰幾何外形和流動狀態(tài)還比較有限，不同混合模型的適用性仍有待驗證，特別是國內(nèi)外近年來提出的一些新 型RANS/LES混 合 方 法，例 如PANS［70］、SAS［71］、CLES［72］、窗 口 嵌 入 式RANS/LES［73］等仍有待在該領(lǐng)域加以推廣應(yīng)用。此外，冰角后方剪切層區(qū)域LES模式的快速合理轉(zhuǎn)換仍然是需要解決的首要問題；從分離-再附流動的結(jié)構(gòu)完整性出發(fā)，目前開展的大多數(shù)工作主要關(guān)注數(shù)值方法在分離起始階段的表現(xiàn)，但對于與再附過程密切關(guān)聯(lián)的多尺度旋渦結(jié)構(gòu)-近壁面流動相互干擾現(xiàn)象預(yù)測效能仍然有待針對性改進，從而進一步提升再附位置的精準捕捉及壓力恢復(fù)過程的合理描述能力，同時為剪切層-再附點擺動運動模式的深入分析提供更為可靠的方法手段［74-75］。

5 結(jié)束語

綜上所述，現(xiàn)階段關(guān)于翼型結(jié)冰后復(fù)雜分離流動的數(shù)值模擬研究取得了長足進步，前緣分離泡演化過程的刻畫精度不斷提高，給出的流場特征量更加準確和全面；所關(guān)注的問題從定常狀態(tài)下的基本氣動力和流場結(jié)構(gòu)等宏觀問題延拓到非定常渦脫落和湍流脈動機制等細節(jié)特征。但是，由于結(jié)冰誘導(dǎo)分離流動與生俱來的影響因素高度復(fù)雜、物理問題尺度差異巨大、多學科深度交叉融合等特征，相關(guān)工作仍然有待進一步開展，尤其在復(fù)雜冰形構(gòu)造、湍流模擬方法的發(fā)展及應(yīng)用、流動非定常特性等方面還需要深入研究。從而更為完備和清晰地認識翼型結(jié)冰狀態(tài)復(fù)雜分離流動演化蘊含的深層次物理機制，更為準確和快捷地提供氣動性能特別是失速特性的預(yù)測結(jié)果，以期更好地為民機適航取證驗證、飛機結(jié)冰致災(zāi)機理及飛行安全防護研究提供支撐，具體包括以下幾個方面：

1）考慮隨機性的高精度冰形模型構(gòu)造研究。由于現(xiàn)階段冰形測量、提取及重構(gòu)手段的限制，當前研究通常采用二維幾何輪廓沿展向拉伸的冰形構(gòu)造方式。但是，實際飛行過程自然生成的翼面結(jié)冰并不滿足上述連續(xù)光滑過渡條件，而是具備顯著的展向非規(guī)則、不連續(xù)及粗糙度等特征。上述因素導(dǎo)致的剪切層不穩(wěn)定性改變/類射流渦系結(jié)構(gòu)生成等特殊現(xiàn)象無疑增加了數(shù)值模擬分析研究的復(fù)雜度。如何在數(shù)值模擬分析研究中如實反映上述精細特征，不僅需要進一步完善物面復(fù)雜幾何邊界計算網(wǎng)格的生成策略，并且對湍流模擬方法特別是RANS/LES混合方法及WMLES方法的近壁面模型構(gòu)造提出了更高要求，期望長度尺度具備更為靈敏的切換能力以及更為良好的復(fù)雜幾何邊界適應(yīng)性。

2）新型湍流模擬方法的構(gòu)造和應(yīng)用研究。對于結(jié)冰翼型失速點附近分離流動演化過程的精準模擬和分析而言，目前所應(yīng)用的湍流模擬方法仍然是不盡完備的，關(guān)于壓力梯度驅(qū)動下K-H不穩(wěn)定性的捕捉精度及旋渦-壁面相互作用現(xiàn)象的辨識能力仍然在很大程度上與冰形-翼型幾何特征直接相關(guān)，同時依賴于冰角后方計算網(wǎng)格密度的提升及拓撲結(jié)構(gòu)的人為設(shè)計。因此，就RANS/LES混合方法而言，期望K-H不穩(wěn)定性的解析降低對于冰角后方區(qū)域當?shù)赜嬎憔W(wǎng)格設(shè)計及布置的嚴格要求、RANS/LES模式的轉(zhuǎn)換更為全面地考慮當?shù)亓鲌鑫锢頇C制及湍流信息而非僅注重邊界幾何特征因素；從而更為精確地識別冰角觸發(fā)剪切層多尺度渦系結(jié)構(gòu)的失穩(wěn)—卷起—壯大—輸運這一完整物理過程，同時更為合理地描述再附過程完成之后的邊界層特征。此外，RANS/LES方法與自適應(yīng)網(wǎng)格的有機結(jié)合也是可供參考借鑒的解決途徑之一。

3）結(jié)冰后分離流動的深層次非定常特性研究。雖然現(xiàn)階段數(shù)值模擬研究能夠獲得不同時刻瞬態(tài)流場的豐富信息，但對功率譜密度/模態(tài)特征等統(tǒng)計學方法量化的時域/頻域量研究涉及不多。需要更為深入地結(jié)合本征正交分解（Proper Orthogonal Decomposition， POD）或動力學模態(tài)分解（Dynamic Mode Decomposition， DMD）等方法，從時間維度上更為精細地辨識和提取深層次非定常因素，以期更為全面地剖析結(jié)冰狀態(tài)下多尺度渦系結(jié)構(gòu)生成演化過程的時空相關(guān)規(guī)律，更好地發(fā)揮RANS/LES混合方法和LES方法對分離區(qū)域湍流脈動描述較為準確的優(yōu)勢。

4）結(jié)冰過程-流動特性的實時耦合分析研究。由于真實飛行條件下結(jié)冰對翼型流動的影響必然屬于時變過程，因而基于實時耦合的分析研究更能體現(xiàn)結(jié)冰生成的隨機性質(zhì)，能夠反映由于冰層累積和增長導(dǎo)致的非定常特性，數(shù)值模擬結(jié)果具備更高的科學意義及實用價值。然而，由于結(jié)冰過程與分離流動兩者涉及的研究領(lǐng)域跨度相對較大，屬于典型的多學科交叉耦合問題，目前較為全面和系統(tǒng)的研究工作尚不多見。