近失速形態下冰脊分離非定常流的IDDES和模態分析

譚 雪,張 辰,徐文浩,王福新，文敏華

(上海交通大學 a.航空航天學院；b.船舶海洋與建筑工程學院；c.高性能計算中心，上海 200240)

溢流冰脊通常形成于機翼防除冰裝置后方[1]，其引起的流動大尺度分離容易導致飛機控制面失效，是過冷大水滴環境影響飛機飛行安全的典型特征之一[2].冰脊雖然結冰量較少，對氣動性能產生的影響卻十分嚴重[3]——冰脊前發生流動分離，并誘導剪切層振蕩，伴隨著復雜的渦脫落、破碎和再附著現象，是研究人員關注的熱點問題[4-5].美國聯邦航空管理局(FAA)在最新發布的結冰適航修正案[6]和咨詢通告中多次強調：溢流冰比傳統認知的臨界冰型更加危險，其造成的壓力低頻脈動可能是導致某些控制面失效，并最終引發飛行事故的根本原因.

在早期研究中，冰脊所造成的氣動損失被重點關注，對分離流的探討主要集中于最大升力系數和失速攻角的損失，其中的代表性成果是Bragg團隊基于伊利諾伊大學低湍流度風洞試驗的相關工作.文獻[7]發現NACA23012翼型的最大升力系數從1.47下降至0.7附近，下降了52%.文獻[8]結合壓力系數分布，進一步指出雷諾數對冰脊翼型近失速狀態下的升阻力影響并不大，這主要歸因于冰脊翼型的前緣發生大尺度分離，減少了壁面流動對雷諾數的依賴.文獻[9]后續進一步研究了帶舵面翼型發生溢流冰現象的氣動損失，證明了這類冰型的破壞性.

為了解析這一物理現象，人們嘗試引入計算流體動力學(CFD)方法對非定常流動過程作進一步研究.由于雷諾平均Navier-Stoks(RANS)方程的時均處理方式難以描述分離流的不穩定性，文獻[10]把RANS方程和大渦模擬(LES)的混合算法引入冰脊翼型的非定常流計算，最早實現了對冰脊剪切層發展、振蕩、脫落和附著過程的模擬.研究表明，冰脊翼型存在瞬態氣動力大幅波動的現象.文獻[11]提出了改進延遲脫體渦模擬(IDDES)的模擬方案，用于改進翼面壓力分布的模擬.文獻[12]提出采用壁面模化大渦模擬(WMLES)改善網格計算不足時的模化應力導致的計算“灰區”問題.文獻[13-14]分別提出5階加權基本無振蕩 (WENO)低耗散格式、自適應長度尺度等方法進一步改善剪切層失穩破碎的模擬能力.這些研究均為冰脊剪切流條件下的高分辨率解析提供了重要的技術支撐.

在最新的試驗研究中，文獻[15]發現冰脊流場中存在一種區別于常規模態的低頻模態，表現為剪切層揮舞形態，與壓力波動存在關聯.剪切區的譜分析結果表明，這種振蕩特性的Strouhal數在0.02附近[16]，與低雷諾數下后臺階流的流動特性比較相似.結合后臺階流理論[17]不難發現，剪切層尾跡中的渦破裂，引發流體向外卷吸速度和向內的回注速度同時發生變化，這種“流量夾帶”和“卷吸流量回注”之間的瞬時不平衡造成了剪切層發展的不穩定，剪切層在這種作用下會伴隨尾流發生低頻自維持振蕩，然而這種現象的本質并沒有完全得到解釋.

本文的研究目的在于采用IDDES方法[18]，對溢流冰脊誘導剪切層振蕩的非定常流特性進行高分辨率模擬，期望通過結合正交分解方法，進一步獲取近失速形態下，造成冰脊分離流壓力波動的主導模態，深入探討冰脊分離流動不穩定性的時域和頻域特征，探索挖掘大尺度相干結構與升力波動的關聯性.本文的研究揭示了近失速形態下冰脊分離流導致升力劇烈振蕩的物理機制.

1 研究方法

1.1 計算條件

采用數值模擬方法對近失速形態下的冰脊翼型分離流模態特征進行計算分析.來流馬赫數Ma設置為0.21，來流迎角α為5°，弦長參考雷諾數Rec為2.0×106.翼型模型采用NACA23012，冰脊在x方向位置與翼型弦長c之比x/c為0.1，冰脊高度k與翼型弦長c之比為 0.013 9.該模型幾何外形規則無毛刺，便于生成高質量的貼體網絡，通過合適的網格分布設計，可以保證冰脊附近的網格正交性.由于文獻[7]在早期的NACA23012縮比模型實驗中已經證明：冰脊采用后向1/4圓時會產生最不利氣動危害，并以此開展了詳細的流場特征研究，得到一系列氣動數據，所以本研究選擇該冰脊外形特征進行對比.計算條件與文獻[10]分離渦模擬(DES)的模擬結果(Ma=0.21,Re=2.0×106)相同.

1.2 計算方法

計算模型的幾何外形如圖1所示，其中：y/c為y方向位置與弦長之比.采用基于CFL3D框架開發的高分辨率流場數值求解器[19]實現三維冰脊流場的計算模擬.基于連續性方程、動量方程和理想氣體方程實現ρ、p、u、v、w、T這6個未知量的計算封閉，其中：ρ為密度；p為壓力；u,v,w分別為x,y,z方向的速度；T為溫度.忽略體積力和外部熱源，其守恒形式可寫成矢量形式如下：

圖1 NACA23012-QR溢流冰脊翼型的幾何形狀Fig.1 Geometry of NACA23012-QR ridge ice airfoil

(1)

(2)

i=1,2,3

式中：H為狀態變量；Fi為i方向的對流矢量通量；Fv,i為i方向的黏性矢量通量；i=1,2,3 分別為笛卡爾坐標系下的x,y,z方向；xi為笛卡爾坐標系下的不同方向；t為時間；E為能量；ui為笛卡爾坐標系下x,y,z方向的速度分量；δij為克羅內克爾符號，i≠j,j=1,2,3;Wi為i方向黏性應力所做的功.按照Stokes假設，黏性應力τij及其所做的功Wi可以表示為

(3)

(4)

式中：xj為笛卡爾坐標系下的不同方向,i≠j,j=1,2,3;uk為笛卡爾坐標系下x,y,z方向的速度分量且k=1,2,3,i≠j≠k；uj為在笛卡爾坐標系下x,y,z方向的速度分量,i≠j；μ為動力黏度；κ為導熱系數.對于理想氣體，根據狀態方程，對單位體積氣體總能量定義的數學表達式做恒等變換后可以得到：

(5)

式中：γ為比熱比.動力黏度μ與溫度T有關，可通過下式計算，則有：

(6)

式中：T0=273.16 K，為參考溫度；對于空氣，有μ0=17.161 μPa·s，為參考溫度對應的動力黏度；Ts=124 K，為與空氣對應的常數.

湍流渦黏性采用重新定義的長度尺度和LES/RANS特征自適應混合函數對Spalart-Allmaras(S-A)模型進行模化.LES/RANS特征自適應混合函數lhyb可由下式定義：

(7)

(8)

lWM=fB(1+fe)lRANS+(1-fB)lLES

(9)

此時，WMLES被激活，并快速從RANS模式切換到LES模式，確保分離區除近壁面區域外的大部分湍流流動能夠計算準確.

當進口條件不包含湍流時，亞格子長度尺度變為DDES的長度尺度，即lhyb=lDDES，

(10)

這時，亞格子長度尺度與延遲脫體渦(DDES)方法保持一致，保留了原始DDES方法模擬分離流場的能力.

綜上所述，新的亞格子尺度可以用下式表達：

Δ=

min{max{Cwdw,Cwhmax,hwn},hmax}

(11)

式中：hwn為沿壁面的法向單元網格尺度；Cw為壁面常數；hmax為網格在各方向上的最大尺度；dw為網格距離壁面的長度.由式(11)可以看出，新定義的長度尺度中既考慮了局部網格的影響，又考慮了壁面距離的影響.近壁面的亞格子長度尺度減小，導致近壁面渦黏系數急劇下降.這種效果增強了長度尺度對近壁面和自由湍流模擬的適應性，更加適合用于模擬復雜的湍流分離流動.

數值離散方面，空間上對流項計算采用基于5階WENO和6階中心差分的混合格式進行離散[19]，黏性項采用4階中心格式；時間上采用雙時間步進策略[20]和隱式近似因子分解內迭代法以2階精度推進求解.物理時間導數項采用2階精度的三點向后差分格式離散，偽時間導數項采用1階向后差分格式逼近.

1.3 計算網格

空間離散采用結構化網格.展向計算域長度為10%翼型弦長，展向面采用周期邊界條件進行關聯.網格劃分的分區策略如圖2所示.由圖2可知，翼型上方冰脊附近和尾跡區域被重點加密，以保證對冰脊誘導的非定常分離流動具備足夠的分辨率；這兩個區域的網格點被精心布置，以保證其具有絕對正交性；且通過局部空間的網格加密可以有效改善DES方法受灰區的影響范圍；近壁區域網格按照無量綱壁面距離y+準則生成(第1層網格到壁面的距離約為3×10-6)，以保證壁面壓力分布的預測要求.

圖2 冰脊翼型網格策略示意圖Fig.2 Schematic diagram of mesh strategy for ridge ice airfoil

表1 無關性驗證的網格參數Tab.1 Irrelevance verification of mesh parameter

1.4 計算驗證

圖3 冰脊翼型的Cp及驗證Fig.3 Validation with Cp and of ridge ice airfoil

2 結果與討論

2.1 流場分析

基于Q準則的NACA23012-QR冰脊翼型瞬時流場特征(Rec=2.0×106，Ma=0.21，α=5° )，如圖4所示.其中：Q=0.01，為流場中速度梯度張量的第2矩陣不變；渦結構采用馬赫數進行著色.從圖4中可以清晰地看出冰脊分離流的空間演化特征，相比于傳統的角狀冰[12]，冰脊在前緣提前誘導出流動分離，在冰脊前方產生駐渦并誘導出大速度的剪切層流動.這種流動呈現出剪切層Kelvin-Helmholtz不穩定性演化出的管狀發卡渦結構，向空間各向扭曲，最終演化為湍流漩渦，與角狀冰的渦脫落和破碎現象有一定相似之處.由于前方駐渦的影響，冰脊流的滯止線與翼型切面的夾角增大，使得分離流相比角冰更難以完全再附到翼型表面，這很可能是其在較小攻角下發生尾流劇烈不穩定振蕩的主要原因.

圖4 基于馬赫數著色的溢流冰脊Q準則瞬時流場(0≤Ma≤0.32)Fig.4 Instantaneous flow of ridge ice Q criterion colored by Mach number (0≤Ma≤0.32)

在α=5° 情況下，冰脊分離流平均馬赫數云圖如圖5所示.由圖5可以看出，冰脊和下翼面尾緣同時誘導出兩個沿流向的剪切層.其中，冰脊誘導的剪切層并未再附到壁面，與下翼面的剪切流動在遠場匯聚.兩個剪切流之間形成大尺度的低能態結構，導致整個上翼面幾乎不存在附面層流動.

為了研究其不穩定性，進一步分析其湍動能Ek.如圖6所示.由圖6可知，在靠近冰脊的回流區的湍動能較弱，這說明剪切層在x/c=0.1～0.3區域的流量夾帶很少影響回流區.在x/c=0.4以后，剪切層充分發展，翼型上方的湍動能明顯增強.湍動能的最大值出現在翼型上方的后緣處，這可能是由于冰脊和下翼面尾緣的剪切層發生相互干擾，從而對整個流場的不穩定性產生影響.關于這種影響的量化描述，在下文中結合主導模態做更詳細的分析.

升力系數的歷史曲線圖如圖7(a)所示，其中：Φ為升力系數功率密度譜(PSD).由圖7(a)可知，升力波動的范圍在0.2～0.6之間，與Loth采用DES方法計算的升力系數范圍大致吻合.利用Fourier分析獲得的Φ如圖7(b)所示,其中：頻率用斯坦頓數St進行歸一化.由圖7(b)可知，冰脊翼型的瞬態升力系數雖然出現了多個頻率，但主導頻率出現在0.02～0.025之間.

圖7 升力系數歷史曲線及其PSDFig.7 History of lift coefficient and its PSD

壓力脈動p′的云圖如圖8所示.由圖8可知，壓力脈動在x/c=0.4達到最大，這可能與此處剪切層渦對流形正在逐漸演化為大尺度渦結構相關.在x/c=0.2，0.3，0.4，0.5，0.8，1.0處(L1～L6)提取的壓力系數脈動Cp,rms特征分布如圖9所示，其中：rms為均方根值(見圖8中L1～L6).由圖9可知，x/c=0.2～0.3處的壓力脈動受剪切層約束，壓力脈動空間主要位于近壁區；x/c=0.4后方，由于在尾跡中的大尺度結構相互干擾，壓力脈動分布趨于平緩.

圖8 壓力脈動云圖Fig.8 Contour of pressure pulsation

2.2 頻譜分析

t0，t1，t2時刻的瞬時壓力脈動云圖如圖10所示.每個時刻間隔1個時間步長，從t0～t2分離渦逐漸向下游移動，前序分離渦在發展中逐漸變大.圖11提取了P1，P2，P3，P4這4個位置的壓力脈動能譜.P1點的壓力脈動能譜存在St=0.132 5 及St=0.025兩個峰值；文獻[12]在研究角冰時將其歸因于剪切層的渦配對和脫落，這在本研究中同樣成立.P2和P3點位置只存在St=0.025的峰值，說明此時剪切層已演化為渦脫落；P4點位于剪切層尾跡中，只能檢測到寬頻信號，此時剪切層已經演化為低能態大尺度結構.

圖10 壓力脈動的時序演化Fig.10 Temporal evolution of pressure pulsation

圖11 P1～P4的壓力脈動功率密度譜圖Fig.11 PSD of pressure fluctuations at locations of P1 to P4

2.3 模態分析

為了進一步分析冰脊剪切流的不穩定性特征，采用快照正交分解(POD)方法[21]對流場的壓力脈動時域信息進行處理.POD的基本思想是將時間和空間上連續的物理量分解.給定M個時刻的流場快照，任意脈動量可以寫作M維POD基和對應模態系數的乘積，則有：

(12)

式中：p′(x,y,tl)為第tl個時刻流場的壓力脈動值；φn(x,y)為第n個POD模態基函數；an(tl)為第tl個時刻中第n個POD基的模態系數；l,n=1,2,…,M，分別為原始流場的第l個時刻以及POD處理后的第n個POD模態.要計算POD基，首先需要定義相關性矩陣R，具體表示為

R=VTV

(13)

式中：V=[p′(x,y,t1)p′(x,y,t2) …p′(x,y,tM)]，為減掉平均流場的脈動量組成的序列向量.

由于R為Gram矩陣，其特征值必大于等于0，特征值對應的特征向量正交.在POD方法中，相關性矩陣R的特征向量決定了POD模態的形式，而R的特征值大小表征對應模態的能量.對于矩陣RN×N,可計算獲得M個特征值和特征向量，即得到M個POD模態.由于特征值越大，對應于權重更高的POD模態，只需保留少數高能量的模態，就能有效捕捉原始數據中的主要特征.對相關性矩陣作奇異值分解，可以得到：

RAn=λnAn

(14)

式中：An=[an(t1)an(t2) …an(tM)]T，為模態系數矩陣；λn為特征值.POD的基表達式如下：

(15)

由式(6)～(9)可知，POD方法能夠將流場脈動分解為3個特征量：① POD特征值；② POD基模態；③ 各模態對應的時間系數.

基于奈奎斯特取樣定理的樣本數驗證如圖12(a)所示.通過對比200，400，800張樣本量的模態能量占比λn/∑λn，基本排除了快照樣本數對分析結果的影響，后續研究采用800×800的取樣標準.圖12(b)為POD前20階模態能量占比以及前100階模態累計能量∑λj如圖12(b)所示.由圖12(b)可以看出，低階模態的特征值較大，是流場不穩定性的主要因素.前20階模態包含了流場中80%的壓力脈動能力，前200階模態占據了99%，因此本研究主要關注其前200階的POD模態特征.

圖12 POD模態歸一化能量及其累計能量Fig.12 Normalized energy and its cumulative distribution of POD mode

在α=5° 情況下，流場脈動的POD基模態如圖13所示.模態1和模態2的能量占比分別為18.2%和13.3%，其空間分布呈現正負間隔的能量脈動序列，主要分布在翼型后緣處,表現為明顯的大尺度相干結構，影響區域涉及后緣的近壁區，對后緣壁面壓力系數造成影響.模態3和模態4的空間分布與前兩階模態特點相似，尺度略小但更靠近壁面.模態49和模態50的能量占比為0.21%，且其影響區間范圍與近壁面的相關性較小.模態199和200的能量占比僅為0.09%，其壓力脈動范圍幾乎不涉及壁面，可推斷其對壁面壓力波動幾乎不產生影響.

圖13 冰脊分離流的典型POD基模態Fig.13 Typical POD modes of ice-ridge separation flow

圖14對比了升力系數頻譜與1/2階POD模態、3/4階POD模態、49/50階POD模態以及199/200階POD模態的時間系數功率密度譜，其中：an(t)為第n個POD基的模態時間系數.由圖14可知，升力脈動的頻譜特征幾乎與前4階模態的POD時間系數重合，這充分說明大尺度相干結構對升力脈動的影響.相比較而言，高階模態的時間系數功率譜與升力系數頻譜并不吻合，這暗示了在近失速形態下，冰脊誘導的剪切層振蕩結構對升力波動的影響并不大.

圖14 典型模態時間系數與升力系數能譜圖的比較Fig.14 Comparison of typical mode temporal coefficients and lift coefficients

3 結論

本文采用改進延遲脫體渦模擬方法，對近失速形態下，溢流冰脊誘導剪切層振蕩的非定常流特性進行高分辨率模擬.結合湍流統計、譜分析和正交分解等方法，本研究進一步分析了冰脊分離流的非定常特性，提取了剪切區的壓力振蕩頻譜特征和全域流場壓力脈動的主導模態，探索挖掘了POD基模態與升力波動的關聯性.本研究主要獲得以下結論.

(1)近失速形態下，溢流冰脊和下翼面尾緣同時誘導出兩個剪切流動.冰脊誘導的剪切層與下翼面上洗流動相互干擾，形成大尺度低能態結構，從而導致整個上翼面幾乎不存在附面層流動.

(2)冰脊誘導的剪切層中存在著兩種典型的脈動頻率與Kelvin-Helmholtz不穩定性相關，具體表現為剪切層誘導的渦配對和渦脫落；剪切層在演化為大尺度結構位置x/c=0.4的壓力脈動最為劇烈；x/c=0.4之前受剪切層約束局限于近壁區，x/c=0.4之后受大尺度結構干擾趨于平緩.

(3)基于正交分解提取獲得的主導模態空間形態為冰脊剪切流尾緣區相互作用的大尺度結構，其時間系數譜與升力脈動的功率密度譜的峰值高度相似，這表明真正造成冰脊翼型升力波動的主要因素來自于尾緣中相互作用的大尺度相干結構.