翼型結冰狀態復雜分離流動數值模擬綜述

趙賓賓，張恒，李杰

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

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

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

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

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

翼型前緣結冰導致的失速分離形態及機制與干凈翼型大迎角狀態存在本質差異。干凈翼型失速分離特征一般與流向逆壓梯度直接相關。而結冰翼型的分離生成通常是由于積冰當地幾何特征發生突變，使得剪切層直接失穩脫落，剪切層渦系結構的發展促進了回流區域和主流區域的流動摻混效應，最終在下游與壁面相互作用、發生再附現象，形成湍流分離泡。因此，驅動分離流動的實質因素是幾何間斷觸發Kelvin-Helmholtz（K-H）不穩定性和迎角效應產生壓力梯度兩者的綜合效應，主導失速分離特征的核心是分離泡結構的生長、擴展和膨脹現象，這無疑大大增加了分離流動數值模擬分析的復雜度。

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

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

1 結冰分離流動模擬概述

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

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

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

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

2 RANS方法

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

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

圖2 明冰冰形附近二維結構網格［14］Fig.2 Two-dimensional structured grid near glaze ice shape［14］

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

圖4 明冰翼型流向速度分布對比（α=6°，上：S-A模型計算結果；下：試驗結果）［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］就 簡 化 冰 形 對 不 同 翼 型 氣 動力造成的影響進行了較為系統的分析，歸納了來流馬赫數和雷諾數效應，分析了冰形弦向位置與壓力分布形態間的關聯，同樣表明小迎角狀態下取得的數值結果基本可用，但同時也指出RANS在失速點附近無法準確計算氣動力，特別對于幾何尺寸較大的冰形而言，由于分離強度顯著增加，數值結果與試驗值間的差異不容忽視。

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

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

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

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

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

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

3 LES方法

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

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

圖8 三維冰形表面非結構網格［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］在 非 常 靠 近 壁 面 的 邊 界 層 內 體 現RANS特征，而在其余區域則體現經過修正的LES方法基本性質。Xiao M C等［32］結合中心-迎風混合格式，基于WMLES系統分析了GLC305及NLF0414翼型結冰狀態下的分離演化規律，表明數值方法能夠準確刻畫自由剪切層的K-H不穩定性，在此基礎上就分離流場的非定常頻域特征進行了深入探討。圖11［32］給出了冰角后方剪切層失穩后二維渦管卷起生成三維大尺度類發卡渦的演化過程。但是，由于結冰分離流場中存在典型再附特征；就小尺度、高頻率湍流結構主導的附著邊界層流動而言，LES在近壁模型的構造和完善方面尚存在一些缺陷，附著區域的合理界定和判斷、?；?解析區域的快速轉換仍然是制約LES廣泛應用的固有難題。因此，LES在翼型結冰分離流動分析中的推廣仍有相當多的前置工作需要開展。

圖11 NLF0414結冰翼型剪切層渦系結構［32］Fig.11 Vortex structures of shear layer of NLF0414 iced airfoil［32］

4 RANS/LES混合方法

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

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

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

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

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

圖15 冰形附近混合網格［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 時均壓力分布結果對比（α=8°）［44］Fig.19 Comparison of time-averaged pressure distribu?tions （α=8°）［44］

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

圖20 時均壓力分布結果對比［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類方法對結冰翼型失速后強分離流場的預測精度仍存在一定提升空間。其中的首要問題是冰角后方RANS/LES區域的延遲轉換特征與剪切層失穩過程解析需求的高度不匹配，由此生成的強亞格子黏性嚴重制約了K-H不穩定性的合理釋放，預測得到的三維旋渦結構生成位置通常較試驗結果更為接近下游，使得剪切層擺動-失穩、二維渦管釋放-扭曲、發卡渦結構發展變化等典型剪切層相關精細湍流特征難以得到充分辨識，制約了分離流動機理特別是非定常特征的深入剖析。

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

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

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

圖23 IDDES-SLA刻畫的剪切層渦系結構生成演化規律［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混合方法。該方法根據分離流場的基本結構，在基準DES方法［28］的基礎上，以區域交界面的形式將流場人為劃分為RANS區域和DES區域，實現湍流模擬方法的混合［63］。Duclercq等［64］基 于 圖25所 示 的 多 塊 結 構 化 網 格，利 用ZDES方法對展向冰脊誘導產生流場的非定常脈動特征進行了較為細致的研究。由于該方法能夠在冰角后方人為強制切換到LES模式，因而具備相對較好的K-H不穩定性預測能力。由圖26［64］可知數值方法能夠較為準確地預測分離泡的渦核位置和時均速度分布情況。通過在冰脊后方剪切層脫落處設置如圖27（a）［64］所示的2個監測點，給出了如圖27（b）［64］所示的鄰近冰角當地速度脈動功率譜密度（PSD）。圖中給出的特征頻率f0表征了剪切層內部旋渦的脫落頻率；而特征頻率f1可視為剪切層垂直運動的基頻。圖28［64］以Q等值面的形式顯示了冰脊后方的瞬時流場，表明ZDES對分離區域小尺度三維湍流結構具備較強的解析能力。

圖25 冰脊附近多塊結構化空間網格［64］Fig.25 Muti-block structured grid near ice ridge［64］

圖26 時均流向速度分布對比（上：PIV試驗結果；下：ZDES計算結果）［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］通過對長度尺度項進行修正，實現了ZDES的 改 進。Zhang等［65］基 于 該 方 法 分 析 了GLC305和NACA23012兩類翼型不同結冰條件下的分離流動，指出計算域展向長度的選取將影響尾跡區域湍流結構的精確解析。由圖29［65-66］可知，在展向長度增加到與分離泡長度接近的情況下，尾跡區域獲得的大尺度旋渦基本形狀和小尺度湍流結構精細程度明顯高于展向長度較小的算例。Costes等［66-67］應用嵌套網格對類似的結冰分離流動算例進行分析時，認為ZDES表現出了較強的網格敏感性，影響流動再附位置的判定和旋渦脫落過程的解析；這與RANS區域和DES區域的空間劃分位置、網格密度差異都有著密切聯系。

圖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］基于附加應力項加權混合的思路發展了一種動態RANS/LES混合模型（Dy?namic Hybrid RANS-LES model， DHRL），Alam等［47-48］基于該方法在非結構網格的基礎上針對GLC305結冰翼型臨界失速算例進行數值模擬，從近壁面速度分布和湍流脈動強度等流場細節角度驗證了數值方法的可靠性。由圖30［48］和圖31［48］給出的粗細2套網格計算結果對比情況可知，由于該方法對RANS方法和LES方法的混合方式與當地網格參數并不直接相關，而是根據當地湍流信息直接判定，因此該混合模型對網格布置的依賴性相對較小，能夠在網格整體較為稀疏的情況下，給出較DDES更為良好的近壁流動細節預測結果。圖32［48］給出的流向速度等值面表明該方法對瞬時流場結構的描述也是較為精細的，能夠較DDES更早地預測到因剪切層失穩導致的旋渦析出和運動現象。但該方法的計算效率目前尚不十分滿意，相同推進步數下耗費的計算時間較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］

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

5 結束語

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

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

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

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

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