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

謝慶墨 陳 亮 張桂勇 孫鐵志

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

(哈爾濱工業大學機電工程學院,哈爾濱 150001)

(工業裝備結構分析國家重點實驗室,大連 116024)

引言

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

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

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

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

1 計算模型

1.1 基本控制方程

均相流模型是空化流場數值模擬使用的常規模型，其質量守恒和動量守恒控制方程為

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

其中，αl是水體積分數，αv是水蒸氣體積分數.

1.2 空化模型

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

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

其中，m+為蒸發源項，m?為凝結源項，pv為當地飽和蒸氣壓.Rb為空泡半徑，表示為

其中Nb=1.0×1013.

1.3 湍流模型

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

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

其中，F1，F2為混合函數；S為剪應力張量的常數項；σω,2為經驗系數.

1.4 動力學模態分解

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

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

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

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

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

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

然后給出DMD 模態?i

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

其中，m為DMD 模態數，μj為模態D的特征值，?j為模態?i的列向量.

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

其中，xi(t)為第i個探測點的時間進程,IMFi(t)為Huang 經驗模態分解法中的本征模態函數.

1.5 模型設置

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

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

2 數值計算方法驗證

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

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

圖2 數值預報的空泡形態周期演化與試驗結果對比Fig.2 Comparison of the cavity evolution between numerical and experimental results

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

3 結果與分析

3.1 水翼空化流場動力學模態分解

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

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

表1 前四階模態對應頻率和放大率Table 1 Frequency and magnification of the first fourth modes

圖4 模態特征值的分布Fig.4 Distribution of mode eigenvalues

動力學模態中包含的能量高低是評價模態在流場中主導作用的重要依據.模態中的能量高低可由模態?i的范數表示[26]

圖5 動力學模態的能量分布Fig.5 Energy distribution of dynamic mode

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

圖6 DMD 前四階動力學模態Fig.6 The first 4 dynamic modes of DMD

3.2 基于模態分解的空化流場分析

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

圖7 模態2，3，4 的對應流場流線圖Fig.7 Flow field streamline of mode 2,3 and 4

3.3 不同空化數下的重構流場分析

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

圖8 σ=1.25，σ=0.75，σ=0.50 工況下模態2 的對應流場流線圖Fig.8 Flow field streamline of mode 2 for σ=1.25，σ=0.75，σ=0.50

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

表2 各階模態頻率Table 2 Frequency of each mode

4 結論

本文針對水翼開展非定常空化流動數值預報，并利用動力學模態分解方法對空化流場進行分析，得到如下結論：

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

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

(3)各模態對應流場中的旋渦結構與DMD 方法得到的模態特征值呈現的相干結構相對應.與空泡脫落頻率相同的第二階模態反應了流場的主要旋渦結構.旋渦結構的尺寸隨著模態階數的增大而減小.

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

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