基于動力學(xué)模態(tài)分解法的繞水翼非定常空化流場演化分析1)

謝慶墨 陳 亮 張桂勇 孫鐵志

(大連理工大學(xué)船舶工程學(xué)院,大連 116024)

(哈爾濱工業(yè)大學(xué)機電工程學(xué)院,哈爾濱 150001)

(工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室,大連 116024)

引言

當(dāng)局部壓力降低至飽和蒸氣壓時，液體發(fā)生由液相轉(zhuǎn)變?yōu)闅庀嗟默F(xiàn)象稱為空化現(xiàn)象[1].空化現(xiàn)象常常出現(xiàn)在高速船舶的舵、螺旋槳、水翼等高速運動的流體機械表面以及高速出入水的航行體表面[2-4].空泡的周期性脫落和潰滅會產(chǎn)生壓力脈沖和輻射噪聲[5].一方面，這會使流體機械表面出現(xiàn)明顯的振動和材料剝蝕[6]，另一方面，當(dāng)壓力脈動頻率與槳自身的固有頻率接近時，還會引發(fā)共振，對流體機械的水動力性能帶來很大的影響[7].

空化的產(chǎn)生及其流體動力學(xué)特性一直是研究的重點[8-10].Knapp 等[11]首先對片云空化進行了廣泛的研究，觀察了片云空化從片云空化腔中剝離/脫落的過程.Leroux 等[12]通過壓力傳感器試驗研究了NACA66 水翼周圍空化的不穩(wěn)定性，并進行了數(shù)值模擬，結(jié)果表明回射流是導(dǎo)致空泡脫落的主要原因.王雅赟等[13]采用不同的湍流模型模擬二維水翼空化流動，發(fā)現(xiàn)空泡尾部大尺度旋渦的形成導(dǎo)致了空泡的斷裂，為后人進行數(shù)值模擬時應(yīng)該選用哪種湍流模型提供了有效的參考.除了理想流體，流體的可壓縮性與熱力學(xué)效應(yīng)也會對流場產(chǎn)生影響[14].王暢暢等[15]模擬了可壓縮流體的空化現(xiàn)象，對空泡潰滅時的激波進行了重點研究，詳細說明了大尺度云空泡的潰滅過程.孫鐵志等[16]研究了不同空化模型在空化流場數(shù)值模擬中的影響，并考慮了熱力學(xué)效應(yīng)對空化流動帶來的影響，為低溫流體的空化現(xiàn)象提供參考.隨著計算機運算能力的發(fā)展，更加復(fù)雜的泄漏渦空化也開始成為研究的熱點[17].彭凱等[18]研究了軸流泵頂隙附近的渦流結(jié)構(gòu)，采用滑移網(wǎng)格技術(shù)模擬了漿葉的旋轉(zhuǎn)，并對旋渦與旋渦、旋渦與壁面之間的相互作用進行了分析.Cheng 等[19]使用大渦模擬(LES)捕獲不穩(wěn)定的尖端泄漏空化流，成功模擬出了尖端泄漏渦、分離渦、誘導(dǎo)渦這些復(fù)雜的渦結(jié)構(gòu)，并將尖端泄漏空化流的空間演化分為3 個階段,揭示了空化對渦度和湍流的影響機理.

由于空化具有強烈的非定常性，并伴隨有復(fù)雜的旋渦結(jié)構(gòu)[20].空泡與壁面、空泡與旋渦之間具有強烈的相互作用，是空泡水動力學(xué)研究的熱點[21-23].空化流場這一復(fù)雜的流動結(jié)構(gòu)具有多尺度特性.目前的研究大多都在對整體流場進行分析，若可以將空化流場分解成若干個不同尺度的流動結(jié)構(gòu)，將對進一步分析空化流動機理提供參考.近年來，動力學(xué)模態(tài)分解(dynamic mode decomposition，DMD)已逐漸應(yīng)用于計算流體動力學(xué)，被認為是評估非線性波動數(shù)據(jù)的重要分析工具[24].DMD 方法可以有效地將流場分解成若干個模態(tài)，并識別出主導(dǎo)模態(tài)進行分析.且具有模態(tài)頻率和增長率單一的特點，在對周期性流動的分析中有較大優(yōu)勢.目前形式的DMD 算法最初由Schmid[25]提出，Rowley 等[26]進一步對其進行擴展，他們將DMD 解釋為Koopman 分析的一種形式.寇家慶和張偉偉[27]首次對現(xiàn)有的DMD 方法與類似的本征模態(tài)分解(POD)進行了比較，闡述了DMD 方法在流體力學(xué)分析領(lǐng)域的優(yōu)勢.并總結(jié)了現(xiàn)有的DMD 方法以及應(yīng)用領(lǐng)域，對DMD方法的現(xiàn)狀與前景作了深入討論.對于空化現(xiàn)象的DMD 分解，韓亞東和韓磊[28]使用DMD 方法研究了文丘里管的空化流動，討論了各階模態(tài)的選取原則以及模態(tài)的特征.Liu 等[29]使用DMD 方法研究了高雷諾數(shù)下的水翼空化流動特征.

目前的研究工作來看，DMD 方法用于空化流場研究的取得了令人鼓舞的進展，但基于DMD 方法探究非定常空化流場特性仍需要進一步研究.本文采用RANS 方法，結(jié)合Schnerr-Sauer 空化模型求解非定常的水翼空化流場，獲取流場數(shù)據(jù)，并利用動力學(xué)模態(tài)分解方法，研究了水翼空化流場各階模態(tài)的主要特征.同時對前四階模態(tài)進行重構(gòu)，得到前四階模態(tài)的對應(yīng)流場，研究了各階模態(tài)的流動特征.最后對不同空化數(shù)下模態(tài)2 的對應(yīng)流場進行了比較，研究了空化數(shù)對空泡脫落產(chǎn)生旋渦對流場的影響機理.

1 計算模型

1.1 基本控制方程

均相流模型是空化流場數(shù)值模擬使用的常規(guī)模型，其質(zhì)量守恒和動量守恒控制方程為

其中，p是壓力，ui是沿第i方向的速度.混合密度ρ 為

其中，αl是水體積分數(shù)，αv是水蒸氣體積分數(shù).

1.2 空化模型

空化模型通常由輸運方程表征，本文采用Schnerr-Sauer 模型[30]，其輸運方程為

其中，αv是水蒸氣體積分數(shù)，源項m+和m?表示為

其中，m+為蒸發(fā)源項，m?為凝結(jié)源項，pv為當(dāng)?shù)仫柡驼魵鈮?Rb為空泡半徑，表示為

其中Nb=1.0×1013.

1.3 湍流模型

本文擬采用SSTk-ω 湍流模型，SSTk-ω 模型由Menter 提出[31]，具體方程如下

其中，α，β為經(jīng)驗系數(shù)，α=0.44，β=0.082 8；k代表湍流的動能，ω表示比耗散率，Pk表示湍動能的產(chǎn)生項，Dω為擴散項，μt為渦黏系數(shù).它們分別定義為

其中，F(xiàn)1，F(xiàn)2為混合函數(shù)；S為剪應(yīng)力張量的常數(shù)項；σω,2為經(jīng)驗系數(shù).

1.4 動力學(xué)模態(tài)分解

采用DMD 方法進行模態(tài)分解時，首先需要將非定常流場的時間序列以快照序列矩陣的方式呈現(xiàn)，快照數(shù)據(jù)來源為物理試驗或數(shù)值仿真，且任意兩個快照之間的時間間隔均為?t，即

線性動力學(xué)矩陣A被定義為非線性映射的近似，該矩陣控制系統(tǒng)從一個時刻到下一個時刻的演變

由于A的維度很高，難以直接計算，因此需要通過降維的方式從數(shù)據(jù)序列中計算出A.因此，通過將最終快照表示為前幾張快照的線性組合，建立一個伴隨矩陣S作為所需的近似

通過友矩陣變換，S的形式為

由于S中的未知量僅為a向量，即可以計算使殘差r最小的a來求得S

當(dāng)殘差r很小時，S的特征值就近似于A的特征值.S矩陣可以看作是A矩陣的低維形式，因此，S的特征值能夠代表A矩陣中的主要特征值.對S矩陣進行特征值分解

然后給出DMD 模態(tài)?i

則任意時刻的流場快照可以通過前m個流場快照表示

其中，m為DMD 模態(tài)數(shù)，μj為模態(tài)D的特征值，?j為模態(tài)?i的列向量.

在流體力學(xué)的分析過程中，往往只能夠獲得數(shù)據(jù)的實數(shù)部分，如速度，壓力等，對于實數(shù)數(shù)據(jù)進行DMD 分解，所得模態(tài)的特征值是共軛的，他們所對應(yīng)的特征向量也是共軛的.流場數(shù)據(jù)的虛部，即流場變化的相位信息，在標(biāo)準(zhǔn)DMD 中并沒有體現(xiàn)，因此會造成一定誤差[25].因此，在進行標(biāo)準(zhǔn)DMD 方法前，使用了希爾伯特--黃變換[32]對原始數(shù)據(jù)進行預(yù)處理.得到原始數(shù)據(jù)波動的相位信息，與原始數(shù)據(jù)疊加得到復(fù)數(shù)形式的數(shù)據(jù).通過希爾伯特--黃變換可以得到更精確的原始數(shù)據(jù)，使DMD 分解結(jié)果更加準(zhǔn)確.其中數(shù)據(jù)的實部代表流場的物理信息，例如速度、壓力，虛部代表這些物理信息的瞬時相位

其中，xi(t)為第i個探測點的時間進程,IMFi(t)為Huang 經(jīng)驗?zāi)B(tài)分解法中的本征模態(tài)函數(shù).

1.5 模型設(shè)置

本文模擬的模型為NACA66 水翼，該翼型的最大厚度位于距前緣45%處，厚度與弦長的比值為12%.最大拱度比位于距前緣50%處，弧度與弦長的比值為2%.來流攻角α=6?，弦長C=0.150 m.計算域與網(wǎng)格如圖1 所示，設(shè)水翼弦長為C，進流段長度為2C，出流段長度為5C，上下壁面相距2C.采用C 型網(wǎng)格，為了提高計算精度，在水翼周圍及尾流段進行加密，網(wǎng)格數(shù)為38 750.邊界條件設(shè)置為：入口處為速度進口，出口處為壓力出口，在水翼表面與上下壁面均施加無滑移條件.流動速度為V=5.33 m/s，其對應(yīng)的雷諾數(shù)為Re=8.93×105.出口壓力由空化數(shù)σ 來確定，初始溫度設(shè)置為298 K，時間步長為1.0×10?4s.

圖1 計算域及網(wǎng)格劃分Fig.1 Computing domains and mesh generation around the hydrofoil

2 數(shù)值計算方法驗證

為驗證數(shù)值模擬的有效性，將空化數(shù)σ=1.00的結(jié)果與試驗[12]進行對比，所有條件設(shè)置與試驗一致.圖2 為一個脫落周期內(nèi)水翼空泡形態(tài)的試驗與數(shù)值模擬對比.可以看出，數(shù)值模擬可以有效模擬出空泡周期性演變中發(fā)展、脫落、潰滅的特性，且數(shù)值模擬的形態(tài)與試驗結(jié)果吻合較好.

圖3 為翼型吸力側(cè)中部x/C=0.6，x/C=0.7 兩個壓力監(jiān)測點壓力值的數(shù)值預(yù)報與試驗數(shù)據(jù)對比.從圖中可以看出，各觀測點的壓力存在周期性波動，這與空泡周期性生成和脫落有關(guān)，其中壓力值較低的區(qū)域處在空化區(qū)域，壓力值高主要是空化潰滅引起的.但每個周期壓力波動的峰值都不相同，體現(xiàn)了空化現(xiàn)象的強非定常性.對比試驗和數(shù)值結(jié)果可以看出，試驗和數(shù)值結(jié)果在振蕩頻率和壓力波動的大小上都有很好的一致性.綜上所述，本文的空化流場數(shù)據(jù)與試驗結(jié)果吻合較好，從而驗證數(shù)值方法的有效性.

圖2 數(shù)值預(yù)報的空泡形態(tài)周期演化與試驗結(jié)果對比Fig.2 Comparison of the cavity evolution between numerical and experimental results

圖3 數(shù)值預(yù)測的水翼吸力面中部壓力值與試驗結(jié)果[6]對比Fig.3 Comparison of the pressure fluctuation on the suction surface of hydrofoil between numerical and experimental results

3 結(jié)果與分析

3.1 水翼空化流場動力學(xué)模態(tài)分解

本文對空化數(shù)σ=1.00 的空化流場中的速度場進行DMD 分解，圖4 表示特征值在復(fù)域上的分布情況.與標(biāo)準(zhǔn)DMD 有所不同，改進后的DMD 方法的不以共軛的形式存在，而是單獨存在.黑色圓圈為單位圓，特征值點與單位圓的關(guān)系揭示了對應(yīng)模態(tài)的非穩(wěn)態(tài)特征.其中特征值位于單位圓外代表該模態(tài)是發(fā)散的，位于單位圓內(nèi)代表該模態(tài)是收斂的，位于單位圓上代表該模態(tài)穩(wěn)定.由于原始數(shù)據(jù)經(jīng)過了希爾伯特--黃變換，還原出的原始數(shù)據(jù)虛部為負值，因此流場數(shù)據(jù)的有效信息大部分集中在虛部為負值的區(qū)域.虛部為正值的區(qū)域為對應(yīng)頻率大，振幅小，這是由于希爾伯特--黃變換存在一定的誤差，這個誤差存在于數(shù)據(jù)的開始部分，所占的比例很小，在此可以忽略[32].結(jié)果表明，在虛部為負值的區(qū)域，大多數(shù)特征值點都在單位圓的范圍內(nèi)，這說明該流動條件下的空化流動是周期性的.圖4 給出了選取的前四階模態(tài)，均在單位圓上，表1 給出了前四階模態(tài)對應(yīng)的頻率和放大率，其中模態(tài)的頻率和放大率可以由下式求出

其中，gj為模態(tài)的增長率，ωj為模態(tài)的對應(yīng)頻率.由于前四階模態(tài)的放大率接近于0，因此選取的前四階模態(tài)都是穩(wěn)定的.前四階模態(tài)的選取方法將在下文進行詳細論述.

表1 前四階模態(tài)對應(yīng)頻率和放大率Table 1 Frequency and magnification of the first fourth modes

圖4 模態(tài)特征值的分布Fig.4 Distribution of mode eigenvalues

動力學(xué)模態(tài)中包含的能量高低是評價模態(tài)在流場中主導(dǎo)作用的重要依據(jù).模態(tài)中的能量高低可由模態(tài)?i的范數(shù)表示[26]

圖5 動力學(xué)模態(tài)的能量分布Fig.5 Energy distribution of dynamic mode

圖6 表示DMD 分解得到的前四階動力學(xué)模態(tài).模態(tài)值只取實部，虛部與實部主要體現(xiàn)在相位上的差別，在分析結(jié)果時一般可忽略.模態(tài)1 對應(yīng)于頻率為零的平均流場.模態(tài)在水翼表面以及水翼尾部的上方，體現(xiàn)了片狀空泡與云狀脫落空泡產(chǎn)生的區(qū)域.模態(tài)2 的特征主要集中在水翼的前端，結(jié)構(gòu)與數(shù)值模擬中的空泡演化形態(tài)相似，體現(xiàn)了空泡脫落的特征結(jié)構(gòu).同時也表明片狀空泡的周期演化是空化現(xiàn)象呈現(xiàn)周期性的主要原因.因此模態(tài)2 可以描述出空泡的發(fā)展趨勢以及周期性演化特征.模態(tài)3 的結(jié)構(gòu)具有相反的正負值，且呈現(xiàn)交替分布的特征.這表明流場中存在明顯的相干結(jié)構(gòu)，圖中正值區(qū)域為空泡的脫落點，負值區(qū)域描述了空泡的潰滅階段的流動特性，這意味著當(dāng)一個結(jié)構(gòu)由于空化而生長，另一個結(jié)構(gòu)由于傳輸和凝結(jié)而減弱.模態(tài)4 的相干結(jié)構(gòu)尺度更小，這說明模態(tài)4 捕捉了能量更少，占比較小的結(jié)構(gòu)，對應(yīng)的頻率為空泡脫落頻率的3 倍，體現(xiàn)了空泡云脫落現(xiàn)象中的高頻特征.隨著頻率增加，模態(tài)的振幅減小，相干結(jié)構(gòu)的尺寸也變小.

圖6 DMD 前四階動力學(xué)模態(tài)Fig.6 The first 4 dynamic modes of DMD

3.2 基于模態(tài)分解的空化流場分析

為了進一步揭示DMD 分解中各階模態(tài)表達的物理意義，對空化數(shù)σ=1.00 的工況進行DMD 分解，取前四階模態(tài)特征值使用重構(gòu)公式，還原出單階模態(tài)對應(yīng)的流場.其中模態(tài)1 為平均流場，對應(yīng)流場不隨時間變化而變化.圖7 表示模態(tài)2,3,4 對應(yīng)的流場，由于不同模態(tài)對應(yīng)的特征頻率不同，因此對應(yīng)的周期表示為T1,T2,T3.對于模態(tài)2、模態(tài)3、模態(tài)4，相對應(yīng)的，流場描述了空化在生長的過程中水翼后緣的大尺度旋渦與回射流的形成，旋渦的周期脫落，以及水翼尾緣處形成的旋渦.觀察流線圖，在模態(tài)2 中，旋渦Ⅰ先從水翼吸力面頭部開始形成，在向尾部移動的過程中不斷生長，移動到接近尾緣時，尾緣處開始產(chǎn)生新的旋渦Ⅱ.旋渦Ⅰ在脫離水翼后停止移動，并開始向下游水平分裂出新的旋渦Ⅲ，在旋渦Ⅲ完全形成后，旋渦Ⅰ退化成節(jié)點，旋渦Ⅲ繼續(xù)向下游移動.旋渦Ⅱ形成后不斷向下游運動.在模態(tài)3 中，水翼吸力面形成的旋渦Ⅰ在脫離水翼后，與水翼尾緣形成的旋渦Ⅱ之間存在融合行為.融合的過程中，旋渦Ⅰ逐漸縮小，旋渦Ⅱ在下游移動的同時向上移動，形成新的旋渦Ⅱ繼續(xù)向下游移動.模態(tài)4 中同時存在多個旋渦，其旋渦行為與模態(tài)2 和模態(tài)3 基本相似，這揭示了兩階模態(tài)的旋渦在運動的過程中存在高頻行為.

圖7 模態(tài)2，3，4 的對應(yīng)流場流線圖Fig.7 Flow field streamline of mode 2,3 and 4

3.3 不同空化數(shù)下的重構(gòu)流場分析

為了對不同空化數(shù)下的空化流場進行比較，分別將DMD 方法應(yīng)用于σ=1.25，σ=0.75，σ=0.50的流場數(shù)據(jù).圖8 表示3 種不同空化數(shù)條件下模態(tài)2 的對應(yīng)流場.對比不同空化數(shù)流場可以看出，在模態(tài)2 中，不同空化數(shù)下的旋渦運動模式基本一致.都包含主旋渦向下游運動，主旋渦分裂形成新的小旋渦.不同的是，隨著空化數(shù)的減小，脫落的主旋渦的尺寸增大，這說明空化數(shù)的減小使流場變化更加劇烈，對水翼的影響更強.

圖8 σ=1.25，σ=0.75，σ=0.50 工況下模態(tài)2 的對應(yīng)流場流線圖Fig.8 Flow field streamline of mode 2 for σ=1.25，σ=0.75，σ=0.50

表2 表示各階模態(tài)對應(yīng)的頻率值，對于不同空化數(shù)下的流場，模態(tài)1 的頻率均為0 Hz，這表明平均流場均為流場中的主要結(jié)構(gòu)，占據(jù)著流場中的大部分能量，從模態(tài)2 的頻率變化可以看出，頻率隨著空化數(shù)的減小而減小.這主要是由于隨著空化數(shù)的降低，空化現(xiàn)象逐漸劇烈，空穴長度逐漸增長，導(dǎo)致空化的生長、脫落時間變長.模態(tài)3、模態(tài)4分別為模態(tài)2 的2,3 倍頻，這是對流場中的主要信息進行補充描述.

表2 各階模態(tài)頻率Table 2 Frequency of each mode

4 結(jié)論

本文針對水翼開展非定常空化流動數(shù)值預(yù)報，并利用動力學(xué)模態(tài)分解方法對空化流場進行分析，得到如下結(jié)論：

(1)利用DMD 方法對速度場進行分解，得到各階流場模態(tài).第一階模態(tài)表示平均流場，其頻率為0 Hz，是流場的主要結(jié)構(gòu).第二階模態(tài)與空泡周期性脫落頻率相同，與空泡脫落的動力特性有關(guān).第3,4 階模態(tài)頻率分別為第二階模態(tài)的2,3 倍頻，這表示流場中的高頻行為.

(2)通過觀察各階模態(tài)對應(yīng)流場一個周期內(nèi)的流場演化，結(jié)果表明，第二階模態(tài)描述了旋渦1在水翼吸力面生長，脫離水翼后形成新的旋渦并向下游移動.旋渦2 在水翼尾緣形成并向下游移動.第三階模態(tài)揭示了旋渦1 與旋渦2 在向下游移動的過程中存在融合行為.第四階模態(tài)中存在更多的旋渦，旋渦運動較為復(fù)雜，包含旋渦分裂、融合、退化.

(3)各模態(tài)對應(yīng)流場中的旋渦結(jié)構(gòu)與DMD 方法得到的模態(tài)特征值呈現(xiàn)的相干結(jié)構(gòu)相對應(yīng).與空泡脫落頻率相同的第二階模態(tài)反應(yīng)了流場的主要旋渦結(jié)構(gòu).旋渦結(jié)構(gòu)的尺寸隨著模態(tài)階數(shù)的增大而減小.

(4)隨著空化數(shù)的降低，空泡脫落頻率逐漸減小.模態(tài)2 的主旋渦尺寸隨著空化數(shù)的降低而增大，說明空化數(shù)的降低加劇了流場的演化，形成的旋渦更大.

本文通過采用動力學(xué)模態(tài)分解方法對繞水翼空化流動進行了分解、分析以及重構(gòu)，重點通過DMD 方法的運用，進一步獲得了空化流場流動結(jié)構(gòu)特性.雖然取得了一定的研究進步，但是在結(jié)合空泡形態(tài)演化和動力學(xué)模態(tài)方法之間的內(nèi)在關(guān)聯(lián)機制仍略顯不足.此外，在低空化數(shù)條件下，RANS 方法存在過度預(yù)報湍流黏性問題，導(dǎo)致對空化瞬變過程的預(yù)報存在一定的不足.在今后的研究中將提升對非定常空化流動的預(yù)報，從而進一步深入開展相關(guān)研究.