杜孟華，田琳琳，朱春玲

(南京航空航天大學(xué) 航空學(xué)院，南京 210016)

0 引言

飛機(jī)在起飛或降落過(guò)程中，會(huì)遭遇過(guò)冷云層，過(guò)冷云層中的過(guò)冷水滴撞向飛機(jī)表面產(chǎn)生積冰，而機(jī)翼結(jié)冰后其前緣形狀被破壞，氣動(dòng)特性發(fā)生改變，對(duì)航空安全造成嚴(yán)重威脅。其中，角狀積冰對(duì)飛機(jī)飛行安全的影響最大。角狀積冰后的流場(chǎng)結(jié)構(gòu)比較復(fù)雜，旋渦相互作用，冰角后方脫出的自由剪切層出現(xiàn)失穩(wěn)[1]，非定常效應(yīng)占主導(dǎo)地位[2]。另一方面，對(duì)帶冰機(jī)翼流動(dòng)分離的預(yù)測(cè)和分析，直接影響氣動(dòng)力和飛機(jī)操縱穩(wěn)定性分析，進(jìn)一步影響飛機(jī)在結(jié)冰狀態(tài)下的安全性評(píng)估。因此對(duì)帶冰機(jī)翼進(jìn)行數(shù)值模擬和氣動(dòng)分析成為研究熱點(diǎn)和難點(diǎn)。

計(jì)算流體力學(xué)（CFD）方法是研究帶冰翼型氣動(dòng)特性的有力工具。早期研究中，使用雷諾平均（Reynolds averaged Navier-Stokes,RANS）方法進(jìn)行數(shù)值模擬，但此方法采用時(shí)均化處理，難以捕捉到帶冰翼型的精細(xì)化流場(chǎng)結(jié)構(gòu)[3]。而后，為了提高結(jié)冰引起的流動(dòng)分離預(yù)測(cè)精度，張恒等[4]使用改進(jìn)延遲脫體渦（improve delayed detached-eddy simulation,IDDES）方法對(duì)帶冰翼型進(jìn)行模擬，結(jié)果表明IDDES 方法有效解析了分離區(qū)小尺度的湍流結(jié)構(gòu)，且較為準(zhǔn)確地預(yù)測(cè)了大尺度分離泡的再附位置。Hu 等[5]采用IDDES 方法對(duì)脊冰和角冰進(jìn)行了數(shù)值模擬，成功捕捉到多尺度渦及相關(guān)流動(dòng)結(jié)構(gòu)。為進(jìn)一步精細(xì)化研究帶冰翼型問(wèn)題，清華大學(xué)Xiao 等[6]通過(guò)IDDES方法結(jié)合渦度相關(guān)的亞格子長(zhǎng)度尺度和剪切層的亞格子長(zhǎng)度尺度，對(duì)角冰問(wèn)題進(jìn)行研究，新的長(zhǎng)度尺度準(zhǔn)確地捕捉到初始自由剪切層渦的卷曲和配對(duì)模式，同時(shí)正確地預(yù)測(cè)了再附區(qū)渦脫落的頻率。Xiao 等在另一項(xiàng)研究中[7]采用壁模型大渦模擬對(duì)三段翼結(jié)冰流動(dòng)分離進(jìn)行了模擬研究，捕捉其非定常流動(dòng)特性。

在非定常復(fù)雜流場(chǎng)研究中，構(gòu)建流場(chǎng)的降階模型可得到復(fù)雜流動(dòng)現(xiàn)象的主要特征，有助于進(jìn)一步分析其流動(dòng)機(jī)理。基于流場(chǎng)特征提取，主要存在兩種典型的降階模型：本征正交分解（proper orthogonal decomposition,POD）[8]方法和動(dòng)態(tài)模態(tài)分解（dynamic mode decomposition,DMD）[9]方法。

前人采用這兩種模態(tài)分解方法對(duì)機(jī)翼繞流后的流場(chǎng)進(jìn)行了一系列的分析研究。例如，Ohmichi 等[10]采用DMD 方法對(duì)三維后掠翼跨聲速抖振的主要流動(dòng)結(jié)構(gòu)和非定常特性進(jìn)行了研究。Ricciardi 等[11]采用不同的正交分解技術(shù)研究了NACA 0012 翼型中的湍流相干結(jié)構(gòu)，主要使用了濾波POD 和頻域POD 描述相干結(jié)構(gòu)的動(dòng)力學(xué)特征，并將它們與噪聲聯(lián)系起來(lái)。國(guó)內(nèi)，王晉軍團(tuán)隊(duì)[12]使用DMD 方法對(duì)實(shí)驗(yàn)中NACA 0015翼型加裝Gurney 襟翼后的尾流速度場(chǎng)進(jìn)行了分析。葉正寅團(tuán)隊(duì)[13]應(yīng)用DMD 方法提取了NACA 0012 翼型在大迎角流動(dòng)中的主頻和高頻倍率及對(duì)應(yīng)的流場(chǎng)模態(tài)結(jié)構(gòu)。張偉偉團(tuán)隊(duì)[14]對(duì)機(jī)翼跨聲速抖振進(jìn)行數(shù)值模擬和模態(tài)分析，通過(guò)DMD 方法研究了激波展向失穩(wěn)特性，結(jié)果顯示機(jī)翼的展長(zhǎng)和后掠是主要的誘因。張辰等[15-16]使用POD 方法對(duì)NACA 23012 帶脊冰的模型提取特征模態(tài)，研究脊冰后的三維湍流尾跡流動(dòng)特征。綜上，采用模態(tài)分析對(duì)機(jī)翼繞流作研究，主要圍繞尾跡結(jié)構(gòu)[17]、幾何形狀[18]、翼尖渦[19]、氣動(dòng)聲學(xué)[20]和抖振[21]等方面，而利用模態(tài)分析方法開展結(jié)冰問(wèn)題穩(wěn)定性分析的研究較少。

在對(duì)非定常流場(chǎng)進(jìn)行模態(tài)分解后，往往會(huì)使用分解后的模態(tài)對(duì)流場(chǎng)進(jìn)行重構(gòu)，用較少的模態(tài)即可對(duì)流場(chǎng)還原，達(dá)到數(shù)據(jù)壓縮的效果。對(duì)于流場(chǎng)重構(gòu)，Rajmohan 等[22]采用POD 方法對(duì)艦面非定常流場(chǎng)進(jìn)行重構(gòu)，以解決艦面流場(chǎng)數(shù)據(jù)量過(guò)大的問(wèn)題。Zhan 等[23]基于POD 方法和機(jī)器學(xué)習(xí)方法構(gòu)建降階模型以快速重構(gòu)飛機(jī)結(jié)冰現(xiàn)象、支持飛機(jī)結(jié)冰的適航認(rèn)定，從而提高航空安全。Sun 等[24]通過(guò)POD 方法建立了一套風(fēng)場(chǎng)快速重構(gòu)的方法，可以提供精確的空間分布律，以克服傳統(tǒng)計(jì)算方法計(jì)算量大的問(wèn)題。

從前述可知，一方面，當(dāng)前針對(duì)翼型結(jié)冰問(wèn)題的研究較少，且大多數(shù)的研究都集中在結(jié)冰后的穩(wěn)態(tài)影響上，不能描述結(jié)冰后的不穩(wěn)定性；另一方面，采用模態(tài)分析對(duì)帶冰翼型流場(chǎng)分析的研究更少，特別是對(duì)角狀積冰的研究。前人有關(guān)結(jié)冰研究的結(jié)果表明，分離泡的發(fā)展和誘導(dǎo)穩(wěn)定性是理解角冰流場(chǎng)的關(guān)鍵；模態(tài)分析方法對(duì)分析主要模態(tài)和運(yùn)動(dòng)頻率存在優(yōu)勢(shì)。同時(shí)，當(dāng)前航空領(lǐng)域?qū)︼w機(jī)防除冰問(wèn)題精細(xì)化模擬和對(duì)帶冰機(jī)翼的全局穩(wěn)定性分析有著迫切需求。因此有必要對(duì)帶冰機(jī)翼開展數(shù)值模擬和模態(tài)分析研究。

本文采用IDDES 方法開展帶角狀積冰翼型繞流非定常模擬研究。在此基礎(chǔ)上分別采用POD 和DMD方法對(duì)非定常流場(chǎng)數(shù)據(jù)進(jìn)行模態(tài)分解，分析帶冰翼型在近失速條件下的流動(dòng)變化特征。此外，在多步長(zhǎng)結(jié)冰數(shù)值模擬中，因角狀積冰而產(chǎn)生的非定常流場(chǎng)數(shù)據(jù)量大，計(jì)算效率低。因此對(duì)某一時(shí)刻的流場(chǎng)分別進(jìn)行POD 和DMD 重構(gòu)，以驗(yàn)證兩種方法對(duì)流場(chǎng)特征提取的效果，同時(shí)也為在結(jié)冰模擬過(guò)程中使用模態(tài)分析方法進(jìn)行流場(chǎng)重構(gòu)數(shù)據(jù)壓縮提供參考。

1 計(jì)算方法

1.1 非定常流場(chǎng)求解方法

考慮控制體流體微團(tuán)Ω，三維積分形式的N-S 方程可以表達(dá)如下：

其中：Ω為控制體，W為守恒變量，F(xiàn)c和Fυ分別為對(duì)流矢通量和黏性矢通量。

使用有限體積法對(duì)方程進(jìn)行離散。其中，對(duì)流通量采用HLLC 格式，時(shí)間推進(jìn)采用隱式LU-SGS 方法，求解非定常流場(chǎng)采用雙時(shí)間步推進(jìn)。選用基于SSTk-ω兩方程模型的IDDES 方法需要構(gòu)造新的湍流混合長(zhǎng)度，具體公式如下：

式中，fhyb為混合函數(shù)，fres為升降函數(shù)，lRANS為RANS 的長(zhǎng)度尺度，CDES為亞格子模型常數(shù)，Δ為亞格子尺度，定義如下：

式中，Cw為經(jīng)驗(yàn)常數(shù)，dw為網(wǎng)格單元到壁面的距離，hmax為三個(gè)方向上網(wǎng)格最大步長(zhǎng)，hw,n為壁面法向網(wǎng)格單元尺度。

1.2 POD 方法簡(jiǎn)介

POD 方法是一種模態(tài)分解技術(shù)，它本質(zhì)上是將流場(chǎng)分解為若干正交基，選擇流動(dòng)的主要模態(tài)對(duì)流場(chǎng)進(jìn)行近似描述。POD 給出具有最小基向量的正交集，在構(gòu)建流場(chǎng)降階模型時(shí)非常有用。POD 方法的應(yīng)用包括流動(dòng)的基本分析、降階建模、數(shù)據(jù)壓縮重建、流量控制和空氣動(dòng)力學(xué)優(yōu)化設(shè)計(jì)。

POD 方法將流場(chǎng)分解為時(shí)均量和脈動(dòng)量，且任意脈動(dòng)量定義為模態(tài)基與時(shí)間系數(shù)的乘積總和，第i個(gè)POD 特征模態(tài)表示為：

式中，M為流場(chǎng)快照個(gè)數(shù)，特征值λ代表模態(tài)能量，ai(tj) 為模態(tài)系數(shù)矩陣，p′(x,tj)為第j個(gè)流場(chǎng)快照脈動(dòng)量，POD 模態(tài)通常按照模態(tài)能量大小進(jìn)行排序。具體推導(dǎo)見參考文獻(xiàn)[8]。

1.3 DMD 方法簡(jiǎn)介

DMD 方法提供了一種將時(shí)間分辨數(shù)據(jù)分解為模態(tài)的方法，且每一個(gè)模態(tài)都有單一的特征頻率和增長(zhǎng)率。DMD 方法可以被視為結(jié)合了POD 和傅里葉變換的有利方面，從數(shù)據(jù)中識(shí)別時(shí)空相干結(jié)構(gòu)。DMD方法起源于流體力學(xué)，由于是一種純粹的數(shù)據(jù)驅(qū)動(dòng)算法，不需要控制方程，因此在流體動(dòng)力學(xué)中廣泛應(yīng)用。

首先假設(shè)線性映射矩陣A將流場(chǎng)vi映射到vi+1，即：

其中，系統(tǒng)矩陣A刻畫了流場(chǎng)的時(shí)間演化特性，但A的數(shù)據(jù)量極大，因此使用矩陣A的近似值來(lái)估計(jì)其特征值。DMD 的模態(tài)可表示為：

具體推導(dǎo)過(guò)程詳見參考文獻(xiàn)[9]。

2 計(jì)算工況及流場(chǎng)計(jì)算結(jié)果

2.1 計(jì)算網(wǎng)格與計(jì)算工況

本文選取GLC305 翼型帶944 號(hào)冰形[25]進(jìn)行三維流動(dòng)模擬研究。其中，944 號(hào)冰形是該翼型在結(jié)明冰的溫度條件下暴露22.5 min 產(chǎn)生的典型雙角狀積冰。翼型幾何如圖1 所示，弦長(zhǎng)為0.914 4 m。網(wǎng)格前后計(jì)算域均為20 倍弦長(zhǎng)。采用混合網(wǎng)格策略，表面生成附面層，第一層網(wǎng)格到壁面的距離為5.0×10-6c（c為翼型弦長(zhǎng)），網(wǎng)格增長(zhǎng)率為1.1，其余用非結(jié)構(gòu)網(wǎng)格填充。計(jì)算域遠(yuǎn)場(chǎng)分別設(shè)置為速度入口和壓力出口邊界條件，物面采用絕熱無(wú)滑移邊界條件，展向?yàn)橹芷谶吔鐥l件。

計(jì)算工況：來(lái)流馬赫數(shù)Ma=0.12，此時(shí)基于翼型弦長(zhǎng)的雷諾數(shù)Re=1.8×106，迎角為α=6°。此迎角對(duì)應(yīng)風(fēng)洞失速點(diǎn)附近狀態(tài)。設(shè)置非定常時(shí)間步長(zhǎng)Δt=0.001 s，每個(gè)時(shí)間步取30 步內(nèi)迭代。

首先開展網(wǎng)格無(wú)關(guān)性驗(yàn)證研究。分別采用細(xì)、中、粗三套網(wǎng)格，具體信息如表1 所示，二維網(wǎng)格如圖2 所示。保持第一層網(wǎng)格高度不變，改變節(jié)點(diǎn)數(shù)量和網(wǎng)格增長(zhǎng)率，總網(wǎng)格量分別達(dá)到6.0×106、1.0×107、1.4×107。計(jì)算結(jié)果表明，基于中網(wǎng)格和細(xì)網(wǎng)格得到的升力系數(shù)趨于一致，具體結(jié)果見表1，且與實(shí)驗(yàn)測(cè)量結(jié)果較為相符，因此本文選用1.0×107網(wǎng)格開展后續(xù)的計(jì)算分析研究。

表1 網(wǎng)格無(wú)關(guān)性驗(yàn)證Table 1 Grid convergence study

圖2 帶冰翼型網(wǎng)格劃分示意圖（網(wǎng)格無(wú)關(guān)性驗(yàn)證）Fig.2 The generated computational meshes of the iced airfoil for grid convergence study

2.2 計(jì)算結(jié)果

首先從時(shí)均流場(chǎng)結(jié)果進(jìn)行分析，圖3 為帶冰翼型表面壓力系數(shù)分布?？梢钥闯?，在帶冰翼型的下表面，基于SSTk-ω湍流模型的計(jì)算結(jié)果（標(biāo)記為SSTRANS）和IDDES 方法的計(jì)算結(jié)果均與實(shí)驗(yàn)值吻合。而對(duì)于翼型上表面，IDDES 方法計(jì)算結(jié)果明顯優(yōu)于SST-RANS 結(jié)果，與實(shí)驗(yàn)值更加貼合。具體而言，在x/c=0～0.25 處，IDDES 方法準(zhǔn)確預(yù)測(cè)了冰角下游的壓力平臺(tái)（由于分離泡存在，導(dǎo)致此處的壓力系數(shù)恒定），但SST-RANS 的結(jié)果并沒(méi)有很好體現(xiàn)這一現(xiàn)象；在x/c＞ 0.25 的壓力恢復(fù)區(qū)域，IDDES 方法計(jì)算的壓力恢復(fù)過(guò)程相對(duì)于SST-RANS 結(jié)果向后移動(dòng)，與實(shí)驗(yàn)值更為相符。

圖3 帶冰翼型表面時(shí)均壓力系數(shù)對(duì)比Fig.3 Comparisons of the obtained time-averaged distributions of the pressure coefficient

表2 給出了時(shí)均氣動(dòng)力的對(duì)比（取相對(duì)穩(wěn)定階段的流場(chǎng)做時(shí)均）。由表可知，IDDES 方法對(duì)于宏觀氣動(dòng)力的預(yù)測(cè)明顯優(yōu)于SST-RANS 方法，尤其是升力系數(shù)的預(yù)測(cè)，與實(shí)驗(yàn)的誤差僅為5.1%。上述結(jié)果表明IDDES 方法能夠準(zhǔn)確地預(yù)測(cè)翼型表面的壓力分布和氣動(dòng)力大小，可為后續(xù)基于數(shù)據(jù)驅(qū)動(dòng)的模態(tài)分析提供準(zhǔn)確的數(shù)據(jù)樣本。

表2 時(shí)均氣動(dòng)力與實(shí)驗(yàn)值的對(duì)比Table 2 Comparisons of the time-averaged aerodynamic forces that obtained with different models

以下針對(duì)瞬態(tài)結(jié)果進(jìn)行分析。圖4 為基于Q準(zhǔn)則的帶冰翼型瞬時(shí)流場(chǎng)特征，渦結(jié)構(gòu)采用湍動(dòng)能染色，清晰地反映了帶冰翼型上翼面冰角后方分離流的空間演化特征。具體而言，在翼型的上表面，自由剪切層由冰角的后緣脫出，由于Kelvin-Helmholtz 不穩(wěn)定性形成展向的渦結(jié)構(gòu)；緊接著，由于展向渦結(jié)構(gòu)無(wú)法自我維持，卷起形成了不同方向的管狀幾何結(jié)構(gòu)。管狀渦在向下游流動(dòng)的過(guò)程中，發(fā)生破碎、合并，使得整個(gè)渦結(jié)構(gòu)更加混亂、更加復(fù)雜。繼續(xù)向下游發(fā)展時(shí)，剪切層中由于強(qiáng)烈的逆壓梯度而卷起旋渦，內(nèi)部的低能分離流和外部的高能主流發(fā)生摻混，使得旋渦重新附著在翼型上表面，形成流動(dòng)再附現(xiàn)象，如圖5（a～c）所示。同時(shí)在冰角后形成分離泡，即翼型表面和分離剪切層之間的一個(gè)循環(huán)流體。

圖4 基于湍動(dòng)能著色的Q 準(zhǔn)則帶冰翼型瞬時(shí)流場(chǎng)Fig.4 Characteristic coherent structures depicted by Q-criterion and colored by turbulent kinetic energy behind the 3-D horn ice

圖5 瞬態(tài)流場(chǎng)展向渦量Fig.5 Evolutions and distributions of the instantaneous spanwise vortices

為了進(jìn)一步研究其不穩(wěn)定性，對(duì)升力系數(shù)歷史曲線進(jìn)行了研究，如圖6 所示。升力系數(shù)在一定范圍內(nèi)波動(dòng)，表現(xiàn)了IDDES 對(duì)小尺度渦精確捕捉的能力，升力系數(shù)振蕩具有不規(guī)則性。同時(shí)觀察到大部分時(shí)間內(nèi)獲得的升力系數(shù)高于平均值0.693，表明這段時(shí)間氣流基本附著在翼型表面。

圖6 升力系數(shù)隨時(shí)間的變化趨勢(shì)Fig.6 Evolutions of the lift coefficient with time

利用快速傅里葉變換（FFT）得到升力系數(shù)功率譜密度（PSD）如圖7，頻率用無(wú)量綱數(shù)St進(jìn)行歸一化。使用翼型弦長(zhǎng)作為特征長(zhǎng)度標(biāo)度，自由來(lái)流速度作為特征速度標(biāo)度進(jìn)行計(jì)算。通常，低頻模式在機(jī)翼上產(chǎn)生的非定常力較大，可在升力系數(shù)中識(shí)別出低頻模式。由圖可知，帶冰翼型升力系數(shù)出現(xiàn)了多個(gè)頻率，主導(dǎo)頻率為St=0.025。因此可以認(rèn)為此帶冰翼型在近失速的工況下出現(xiàn)了低頻流動(dòng)振蕩現(xiàn)象。這種振蕩表明了流動(dòng)在失速和非失速狀態(tài)下的準(zhǔn)周期切換。對(duì)于不同翼型以及不同迎角，其低頻模式的頻率是不相同的[26]。為了進(jìn)一步明晰帶冰翼型流場(chǎng)中的頻率特征及主要流動(dòng)特征，下面結(jié)合模態(tài)分析進(jìn)行研究。

圖7 升力系數(shù)能譜圖Fig.7 PSDs of the iced airfoil's lift coefficients

3 模態(tài)分析

針對(duì)上述非定常計(jì)算結(jié)果，分別采用POD 和DMD 方法對(duì)角狀冰形下游的流動(dòng)不穩(wěn)定性進(jìn)行研究。在流場(chǎng)充分發(fā)展后的穩(wěn)定階段，選取400 個(gè)流場(chǎng)快照，每個(gè)快照的時(shí)間間隔Δt=0.01。考慮到展向流動(dòng)比其他方向流動(dòng)弱得多，因此將非定常三維流場(chǎng)沿流向切片進(jìn)行特性分析。

3.1 POD 模態(tài)分解結(jié)果

POD 方法可提取主要模態(tài)，通過(guò)獲取主要模態(tài)的特征，分析對(duì)流動(dòng)分離產(chǎn)生影響的主要因素。對(duì)POD 分解后的模態(tài)按照能量大小進(jìn)行排序。前100階模態(tài)能量占比和前100 階累積模態(tài)能量占比如圖8所示。第1 階和第2 階模態(tài)能量占比分別為16.16%和14.7%，為流場(chǎng)的主要模態(tài)。前10 階模態(tài)能量基本包含了整個(gè)流場(chǎng)80%的能量，前30 階累積能量超過(guò)整個(gè)流場(chǎng)90%的能量。因此，前30 階的模態(tài)可代表此狀態(tài)下流動(dòng)不穩(wěn)定的最主要特征。POD 分析表明，僅用少量的模態(tài)即可很好地捕捉流場(chǎng)中的波動(dòng)，這有利于對(duì)流場(chǎng)數(shù)據(jù)進(jìn)行顯著壓縮，同時(shí)也利于進(jìn)行流動(dòng)分析。下面對(duì)前幾階主要模態(tài)進(jìn)行分析。

圖8 各階模態(tài)能量占比和累積模態(tài)能量Fig.8 Percentage of the energy and cumulative energy contributions with different POD modes

圖9 給出了前6 階POD 模態(tài)壓力分布，可以清楚地看到POD 方法清晰地提取了流場(chǎng)相干結(jié)構(gòu)。翼型上表面呈現(xiàn)交替相間的能量序列，在前6 階POD 模態(tài)中，已清楚地捕捉到沿剪切層發(fā)展并向后緣傳播的大波動(dòng)。模態(tài)1 和模態(tài)2 由低頻大尺度流動(dòng)結(jié)構(gòu)組成，顯示在結(jié)冰翼型前緣有一個(gè)較大的渦結(jié)構(gòu)，對(duì)前緣的壁面壓力系數(shù)有顯著的影響。3～6 階模態(tài)均包含較小的模式結(jié)構(gòu)?？臻g相干結(jié)構(gòu)從x=0.2處逐漸出現(xiàn)，這種結(jié)構(gòu)可能與剪切層（因Kelvin-Helmholtz 不穩(wěn)定性而失穩(wěn)）開始脫落有關(guān)。

圖9 前6 階POD 模態(tài)壓力分布Fig.9 The first six POD modes of pressure field

圖10 為快速傅里葉變換得到的前4 階POD 模態(tài)時(shí)間系數(shù)。從圖中可以看出，主要頻率在1～50 Hz范圍內(nèi)，且存在從低頻到高頻多個(gè)主導(dǎo)頻率。這意味著此類流場(chǎng)中，旋渦具有多樣性。圖10（a、b）給出了模態(tài)1 和模態(tài)2 的時(shí)間系數(shù)經(jīng)過(guò)快速傅里葉變換得到的結(jié)果，據(jù)此可得高幅低頻的主要成分。該模態(tài)下主要頻率為f=1.5 Hz 和f=2.25 Hz，與文獻(xiàn)[26]的剪切層拍打頻率范圍一致。該頻率由剪切層失穩(wěn)而產(chǎn)生，且與冰角后的分離泡密切相關(guān)。在圖10（c、d）中，模態(tài)3 和模態(tài)4 對(duì)應(yīng)的主導(dǎo)頻率分別為29 Hz 和45 Hz，這說(shuō)明與剪切層失穩(wěn)的相關(guān)頻率主要在1 階和2 階主導(dǎo)模態(tài)中，進(jìn)一步證明了剪切層的拍打模式主導(dǎo)了此類非定常流動(dòng)。

圖10 通過(guò)快速傅里葉變換得到的前4 階POD 模態(tài)時(shí)間系數(shù)Fig.10 Corresponding temporal coefficient for the first four POD mode via fast Fourier transformation

圖10（c、d）中低頻信號(hào)開始轉(zhuǎn)變?yōu)楦哳l信號(hào)，主導(dǎo)頻率變?yōu)?5 Hz。用分離泡的平均長(zhǎng)度和自由來(lái)流速度對(duì)頻率進(jìn)行無(wú)量綱化，得到非定常無(wú)量綱數(shù)St=0.68。這與Ansell 和Bragg[26]根據(jù)各類研究總結(jié)到的規(guī)則模式的St取值范圍（0.5＜St＜0.8）一致。這種規(guī)則模式源于流動(dòng)再附位置附近的壓力波動(dòng)（主要由剪切層內(nèi)的旋渦運(yùn)動(dòng)和渦脫落引起的）。圖9 的模態(tài)3 和模態(tài)4 結(jié)果可以清晰地看到規(guī)則的旋渦脫落，更加證實(shí)了規(guī)則模式與旋渦運(yùn)動(dòng)有關(guān)，且規(guī)則模式的不穩(wěn)定特征為分離流動(dòng)的次要模式。

3.2 DMD 的模態(tài)分解結(jié)果

上述已采用POD 方法對(duì)流場(chǎng)進(jìn)行了分析。此外，相較于POD 方法分解后的模態(tài)對(duì)應(yīng)多個(gè)頻率，DMD 分解后的模態(tài)對(duì)應(yīng)單一頻率和增長(zhǎng)/放大率，據(jù)此可對(duì)每個(gè)模態(tài)的穩(wěn)定性進(jìn)行分析。因此，下文將對(duì)流場(chǎng)樣本進(jìn)行DMD 模態(tài)分析，并將DMD 和POD 方法進(jìn)行對(duì)比研究。

首先進(jìn)行DMD 分解，然后將DMD 模態(tài)進(jìn)行排序。當(dāng)前存在兩種不同的排序方式[27]：1）直接使用模態(tài)振幅進(jìn)行排序，這種排序的優(yōu)勢(shì)是計(jì)算簡(jiǎn)單、計(jì)算量小，缺點(diǎn)是容易受到噪聲的影響而忽略初始幅值小但增長(zhǎng)很快的信號(hào)；2）采用模態(tài)范數(shù)（表示各個(gè)模態(tài)對(duì)流場(chǎng)的貢獻(xiàn)和影響）的方法進(jìn)行排序，優(yōu)點(diǎn)是容易衡量各模態(tài)對(duì)流場(chǎng)的貢獻(xiàn)值，缺點(diǎn)是計(jì)算量比較大?？紤]到本文樣本量少，計(jì)算量較小，因此采用模態(tài)范數(shù)進(jìn)行排序。

DMD 的模態(tài)特征值以共軛復(fù)數(shù)的形式呈現(xiàn)，將每一對(duì)特征值看作一個(gè)DMD 模態(tài)。圖11 為DMD模態(tài)特征值的實(shí)部和虛部在復(fù)平面上的分布情況。其中，紅色實(shí)心點(diǎn)表示在單位圓外，對(duì)應(yīng)不穩(wěn)定狀態(tài)；空心點(diǎn)表示在單位圓上或單位圓內(nèi)，對(duì)應(yīng)周期模態(tài)或穩(wěn)定模態(tài)。結(jié)果表明，絕大部分的特征值都在單位圓上。從圖11（b）的局部特征值分布可以看出，紅色實(shí)心點(diǎn)雖然在單位圓外，但十分靠近單位圓，說(shuō)明不穩(wěn)定模態(tài)處于弱發(fā)散狀態(tài)。

圖11 DMD 特征值的實(shí)部和虛部在單位圓上的分布Fig.11 The real and imaginary parts of the DMD eigenvalues on the unit circle

圖12 為放大率和模態(tài)幅值之間的關(guān)系，其中紅色圓點(diǎn)表示放大率大于0，黑色圓點(diǎn)表示放大率小于0 或等于0。放大率大于0 對(duì)應(yīng)模態(tài)系數(shù)發(fā)散，反之模態(tài)系數(shù)為收斂狀態(tài)。帶冰翼型在此迎角狀態(tài)下不穩(wěn)定模態(tài)數(shù)占總模態(tài)數(shù)的32.25%，對(duì)應(yīng)流動(dòng)結(jié)構(gòu)隨著時(shí)間不受控制增長(zhǎng)，從而導(dǎo)致流動(dòng)分離的發(fā)生。

圖12 放大率與模態(tài)幅值之間的關(guān)系Fig.12 Relationship between the growth rate and the amplitude

圖13 為前6 階DMD 模態(tài)壓力分布。首先DMD的第1 階模態(tài)為主導(dǎo)模態(tài)，對(duì)應(yīng)靜態(tài)模態(tài)，代表時(shí)均的壓力分布，與所求流場(chǎng)的時(shí)均分布非常相似，因此可以作為數(shù)據(jù)集統(tǒng)計(jì)平均穩(wěn)定性的簡(jiǎn)單驗(yàn)證。如果流場(chǎng)快照不具有收斂性，則1 階模態(tài)與時(shí)均壓力分布不符，所以證實(shí)了快照的數(shù)量足以進(jìn)行當(dāng)前分析。模態(tài)2 和模態(tài)3 為除去平均模態(tài)最主要的DMD 模態(tài)，主要反映翼型中段交替出現(xiàn)的渦結(jié)構(gòu)。與POD 模態(tài)類似，渦結(jié)構(gòu)開始出現(xiàn)的位置與剪切層開始失穩(wěn)的位置相同。但是與POD 模式不同的是所識(shí)別的結(jié)構(gòu)分解成了小的不規(guī)則結(jié)構(gòu)，同時(shí)結(jié)構(gòu)之間似乎表現(xiàn)出復(fù)雜的相互作用。

圖13 DMD 模態(tài)壓力分布Fig.13 The first six DMD mode of pressure field

3.3 流場(chǎng)重構(gòu)

在結(jié)冰數(shù)值模擬多步長(zhǎng)計(jì)算中，帶角狀積冰翼型的非定常流場(chǎng)數(shù)據(jù)量大。對(duì)此，本部分對(duì)角狀積冰分解后的模態(tài)進(jìn)行流場(chǎng)重構(gòu)，為解決多步長(zhǎng)結(jié)冰中非定常流場(chǎng)數(shù)據(jù)量大的問(wèn)題提供參考；同時(shí)，還進(jìn)一步校核兩種模態(tài)分解方式對(duì)流場(chǎng)的提取效果。

計(jì)算中選取任一時(shí)刻的流場(chǎng)（本文選取t=10.74），分別使用POD 和DMD 分解后的模態(tài)對(duì)流場(chǎng)進(jìn)行重構(gòu)。其中，對(duì)于POD 方法，選擇前10 階模態(tài)和前20 階模態(tài)加上平均流場(chǎng)進(jìn)行重構(gòu)。圖14（a）為原始流場(chǎng)結(jié)果，圖14（b、c）為POD 重構(gòu)后的結(jié)果，可以看出前10 階模態(tài)重構(gòu)基本與原始流場(chǎng)一致，低壓區(qū)的位置以及大小基本吻合。前20 階重構(gòu)結(jié)果與前10 階相差較小，因此前10 階模態(tài)即可還原流場(chǎng)特征。

圖14 POD 與DMD 流場(chǎng)重構(gòu)結(jié)果(t=10.74)Fig.14 Flow field reconstructions based on the POD and DMD methods，respectively (t=10.74)

對(duì)于DMD 方法，前20 階模態(tài)重構(gòu)的結(jié)果與原始流場(chǎng)差距很大，直到前300 階模態(tài)重構(gòu)才與原始流場(chǎng)比較接近。圖15 顯示了模態(tài)數(shù)和損失函數(shù)之間的關(guān)系。損失函數(shù)[27]可衡量重構(gòu)流場(chǎng)和原始流場(chǎng)之間的性能損失，用于選擇適當(dāng)?shù)哪B(tài)數(shù)量?？梢钥吹剑?50 階模態(tài)的損失函數(shù)小于20%，前300 階模態(tài)的損失函數(shù)小于10%。因此，DMD 方法需要比POD 方法更多的模態(tài)數(shù)才能準(zhǔn)確地重構(gòu)流場(chǎng)。

圖15 模態(tài)數(shù)和損失函數(shù)之間的關(guān)系Fig.15 Relationship between the extracted DMD modes number and the loss function

4 結(jié)論

本文采用IDDES 方法對(duì)帶角狀冰形的翼型流動(dòng)進(jìn)行了非定常數(shù)值模擬?；诖耍謩e采用POD 和DMD 兩種典型的模態(tài)分解方法對(duì)數(shù)值模擬結(jié)果進(jìn)行模態(tài)分析，同時(shí)使用分解后的模態(tài)對(duì)流場(chǎng)進(jìn)行重構(gòu)。主要結(jié)論如下：

1）POD 方法識(shí)別出與剪切層拍打相關(guān)的低頻高幅的頻率，且此頻率出現(xiàn)在主導(dǎo)模態(tài)中。證明剪切層拍打模式主導(dǎo)了此類非定常流動(dòng)。

2）DMD 方法通過(guò)分析模態(tài)頻率和放大率，發(fā)現(xiàn)32.25%的模態(tài)處于發(fā)散狀態(tài)，對(duì)應(yīng)的流動(dòng)結(jié)構(gòu)隨時(shí)間不受控制增長(zhǎng)，導(dǎo)致流動(dòng)分離的發(fā)生。

3）POD 方法還原流場(chǎng)特征所需要的模態(tài)階數(shù)遠(yuǎn)小于DMD 方法，因此在今后的結(jié)冰模擬中，可采用POD 方法重構(gòu)部分流場(chǎng)，使用更少的數(shù)據(jù)達(dá)到數(shù)據(jù)壓縮的效果。