粘彈性圓柱繞流場的稀疏化動態模態分解方法

李 璇，蘇 進

(西安工程大學 理學院，陜西 西安 710048)

0 引言

粘彈性流體是一種具有粘性和彈性的特殊非時變性非牛頓流體，在生物工程等多種領域有著廣泛應用。針對粘彈性圓柱繞流問題，不僅要考慮圓柱壁面(邊界)的流動特點，還需考慮粘彈性圓柱繞流所具有的減阻特性。因而，分析彈性誘導的渦團結構特征對粘彈性流體的工程應用具有重要價值[1-2]。近年來，學者們對流體流動結構進行了較多研究：Thomases等[3]針對四輥軋機幾何中的粘彈性流體，利用本征正交分解(proper orthogonal decomposition，POD)提取了彈性模態，并根據其能量貢獻分析了流體動力學特征；Henshaw等[4]采用POD量化流體動力學系統的模型結構；Shakeri等[5]利用POD提取了聚合物的彈性湍流結構并發現反向旋轉渦對流動的總動能有顯著貢獻。最近在流體動力學領域內，動態模態分解[6](dynamic mode decomposition，DMD)廣泛用于分析非定常流場[7]的流動特征和構建低階的流場動力學模型[8-9]。與POD形成的完全基于空間相關性和能量含量的模態層次結構相比，DMD方法分解后的每個模態由空間的相關結構組成，在時間上具有相同的線性行為[10]。DMD可以看作是將奇異值分解(singular value decomposition，SVD)在空間降維方面的優點與快速傅里葉變換(fast fourier transform，FFT)在時間頻率識別方面的優點相結合。然而，在高速繞流條件下，粘彈性圓柱繞流場流體粘性減小阻力較低，隨著速度的增加減阻效果逐漸減弱，這導致粘彈性圓柱繞流流場結構中有很多復雜的模態特征。DMD雖然在牛頓流體流場數據降維方面起到了一定的作用[11-12]，但是對于粘彈性圓柱繞流數據的減阻現象往往不能起到很好的降維效果。

針對粘彈性圓柱繞流問題，為實現粘彈流高效的結構分析和動態預測，引入一種基于稀疏促進優化的DMD方法(SP-DMD)，與傳統的DMD相比，該方法可以剔除粘彈性流體中的非關鍵模態。在粘彈性流場演化中，SP-DMD使粘彈性流場中關鍵模態的數量和降維近似殘差之間處于平衡狀態，從而更利于分析和發現那些能夠對粘彈性圓柱繞流場發展起關鍵性作用的因素。此外，當模態數量增加時，SP-DMD重構的流場可以得到更多的局部流動細節，并且與原始流場進行對比，局部較小尺度的渦團形態以及空間分布的特點都能夠很好地重現出來。

1 DMD算法

從數學角度看，DMD是對Koopman算子模態的近似，從諧波的角度反映流場的特性。DMD方法的運用需要從數值模擬或物理實驗中收集一系列快照x1，x2，...，xm+1，并形成一個數據矩陣：

(x1x2...xm+1)。

其中矩陣的每一列表示每一時刻的空間狀態。此外，假定數據在時間上是等距采樣的，時間步長為Δt。通常情況下，每一個快照xi=x(iΔt)是一個有n分量的復向量，即xi∈n。

從快照序列中構造兩個快照矩陣：

其中，X∈n×m，X′∈n×m。并假設快照是由離散時間線性時變系統生成的，該系統表示如下：

xk+1=Axk，k=1，2，...，m。

線性算子A記錄著快照序列的動態特性。在粘彈性流體問題中，矩陣A中包含大量數據，且一般情況下為復數。每個快照xi中分量的數量通常遠大于快照的數量，即n>>m。

根據兩個連續時間步長之間的線性關系xk+1=Axk，k=1，2，...，m，將數據矩陣X和X′通過線性算子A連接起來，則矩陣X′有如下表示：

因此，求出線性算子A是DMD的主要目標。

(1)對矩陣X∈n×m進行奇異值分解(SVD)：

假設數據總維數是m，上述SVD通過選擇適當的截斷值r來減少數據矩陣的維數，被消除的冗余項就用rem(redundancy term)表示，其維數是m-r。

其中，矩陣Q的Frobenius范數由下式來決定：

(1)

動力學方程(1)所生成的每一組數據為：

且每個αi是對應DMD模態的振幅。同樣，在矩陣形式中有：

則有

X≈ΦDαVand。

(2)

(3)

J(α)可以等價表示為：

J(α)=α*pα-q*α-α*q+s，

(4)

(7)通過最小化關于α的二次函數(4)可以得到該優化問題(3)的DMD振幅的最優矢量為：

2 稀疏化動態模態分解(SP-DMD)方法

DMD算法雖然能夠求出振幅最優向量和提取出關鍵模態，但由于在維數較高的粘彈流中還存在部分非關鍵模態未剔除，因此引入稀疏化動態模態分解(SP-DMD)方法來剔除這些非關鍵模態。

2.1 SP-DMD核心思想

SP-DMD需要在所有DMD模態中提取出若干關鍵的具有非零幅值的模態，然后通過算法來調整其幅值，具體方法如下：

具有非零幅值的DMD模態可以通過求解如下凸優化問題而得到：

minJ(α)，
s.t.ETα=0。

(5)

式(5)中，矩陣E反映了振幅向量的稀疏結構信息。E的列均是單位向量，其非零元素對應于具有零幅值的DMD模態。獲得具有非零幅值的DMD模態后，通過改變其幅值實現降維近似。

將優化問題(3)進行修改，通過給式(3)中的目標函數J(α)增加一個正則化項card(α)來解決稀疏性的問題，這就懲罰了未知振幅向量Dα中非零元素的數量：

(6)

在修改后的優化問題(6)中，γ是正則化參數，它反映了“稀疏化”振幅向量Dα的權重性。γ越大，說明越重視向量Dα中非零元素的數量。因此，鼓勵用(6)這個更稀疏的解決方案。

一般來說，找到問題(6)的解決方案相當于組合搜索，對于任何問題，組合搜索會變得很棘手。為了繞過這個問題，引入問題(6)的更寬泛的表述，通過用向量Dα的L1范數來代替基函數：

(7)

SP-DMD的核心思想就是找到式(7)的解，這是一個凸優化問題，解決該問題使用了交替方向乘子法(alternating direction method of multipliers，ADMM)，在實驗或數值快照降維近似質量和DMD模態數量之間達到預期的平衡后，從而確定最終優化后的稀疏化振幅向量Dα。

2.2 交替方向乘子法

交替方向乘子法(ADMM)是一種用于求解優化問題的計算框架，適用于求解分布式凸優化問題。ADMM算法將大的全局問題分解為多個較小且較容易求解的局部子問題，并通過協調子問題的解而得到大的全局問題的解。ADMM通常用于求解兩個優化變量且只含等式約束條件的優化問題，其一般表示形式為：

minf(x)+g(z)，

s.t.Ax+Bz=c。

其中，x和z是優化變量，f(·)和g(·)都是凸函數。

引入增廣拉格朗日函數：

其中：λ為拉格朗日乘子；ρ為懲罰系數且ρ>0。

由于ADMM算法是基于增廣拉格朗日函數最小化的迭代算法，迭代方式如下所示：

λk+1=λk+ρ(Axk+1+Bzk+1-c)。

從初始點開始進行迭代，直到滿足如下優化條件和停止準則：

Axk+1+Bzk+1-c2≤εprim，

zk+1-zk+12≤εdual。

其中，εprim，εdual表示停止準則的閾值。

3 數值算例

3.1 粘彈性圓柱繞流模型參數設置

為了模擬粘彈性圓柱繞流在高Re數下的湍流減阻現象，設置流場中的平均流速為U，圓柱直徑的特征長度尺寸為d，幾何特征停留時間為d/U。流體彈性特征的參數為Wi=λU/d，表征流體黏性特征的參數為Re=ρUd/η。

選取如圖1所示的幾何示意圖進行二維粘彈性圓柱繞流湍流減阻的模擬。

圖1 粘彈性圓柱繞流的幾何模型

采用雙分布 LBM-IBM 耦合算法模擬粘彈圓柱繞流，設置Re=1000，粘度比β=ηs/η=0.01。數值模擬參數如表1所示。

表1 數值模擬參數

宏觀方程的處理采用IBM耦合D2Q9格子模型，D2Q9格子模型中碰撞遷移過程的邊界處理為非平衡外推格式，區域上下邊界采用無滑移條件，即

出口處按充分發展條件處理，出口處的水平速度的法向導數為0，其中n是單位法向量，即

這里入口處的初始水平速度u設置成拋物狀，初始垂直速度v=0，入口的初始速度

u(0，y)=3Umaxy(10-y)/50，

其中，Umax是通道水平中線處的最大速度且U=2Umax/3。

對于本構方程采用Oldroyd-B方程的格子Boltzmann模型，應力張量的分布函數上下壁面的邊界條件應用非平衡態外推方法，計算區域的上下邊界的應力分量為

τyy=0。

入口和出口邊界的應力分量采用充分發展邊界條件，即

3.2 粘彈性圓柱繞流數值結果

不同的Re(Reynolds)數對應的Wi(Weissenberg)不同，且Wi是Re的函數。根據3.1節所述，模擬出Re=1000時的粘彈性圓柱繞流流場的數據，并且對流體彈性特征數Wi=1的數據進行研究。每一時刻的流場數據構成一個301×601的網格。

在研究粘彈性圓柱繞流問題的過程中，首先要將二維數據轉化成一維，由此可知每一時刻所對應的空間維數高達180 901維。由于各種運動形態和不同結構的尾流流態會表現出不同的能量和頻率，僅通過對瞬態流場的觀察，很難對粘彈性圓柱繞流流場的特性進行深入理解。因此研究過程中，選擇了連續1 000個時刻的二維粘彈性流場數據進行DMD算法處理，可以得到302個模態。得到的302個DMD模態的特征值在復數坐標系內的分布情況如圖2所示，這些特征值進行對數化后的分布情況如圖3所示。

圖2 DMD未對數化的特征值分布 圖3 DMD對數化的特征值分布

由圖2可知，大部分的DMD模態捕捉到的流場結構強度近似處于穩定的狀態。此外，還可以發現有個別特征值落在單位圓內部，并且距離單位圓較遠，此現象表示這部分特征值所對應的流場結構強度的衰減率比較大，流場結構的狀態不穩定，將在隨后的圓柱繞流流場演化過程中逐漸消散。由圖3可知，每個特征值的實部(即衰減率)都基本接近于0，這更直觀地表明了大部分DMD模態所對應的流場結構強度是近似穩定的。

采用SP-DMD和標準DMD算法得到的粘彈性流場非對數化特征值的對比如圖4(a)所示。由于該圖不能清晰地看出特征值的分布情況，將其局部放大，如圖4(b)所示。

圖4 DMD和SP-DMD未對數化的特征值分布

在采用SP-DMD時，分別保留了30個模態和100個模態進行對比。SP-DMD關注的是DMD模態在整個粘彈性圓柱繞流流場的演化過程中對流場的貢獻，通過SP-DMD能夠更加方便地識別并提取在整個流場的演化過程中起到比較重要影響作用的模態信息。由圖4可知，利用標準DMD算法所得到的具有較大衰減率的模態在SP-DMD中都被有效剔除了，從而保留了有穩定狀態的DMD模態。

幅值與頻率之間的關系可以定義為DMD頻譜，在研究粘彈性流場問題的過程中，DMD頻譜能夠將流場結構所包含的能量情況反映出來，因此又將其稱之為能譜。采用兩種方法所獲得的模態幅值與特征值虛部(頻率)的關系如圖5所示，對其進行局部放大如圖6所示。利用兩種方法所得到的模態幅值與特值實部(衰減率)的關系如圖7所示，其局部放大如圖8所示。

圖5 采用 DMD 與 SP-DMD 得到的模態幅值與特征值虛部(頻率)的關系

圖6 采用 DMD 與 SP-DMD 得到的模態幅值與特征值虛部(頻率)的關系的局部放大

圖8 采用 DMD 與 SP-DMD 得到的模態幅值與特征值實部(衰減率)的關系的局部放大

由圖5、6、7、8可知，部分具有較大衰減率的非穩定的DMD模態所對應的流場結構承載著較多的能量，其對流場初期演化的影響較強，但在后期的演化過程中，影響力明顯下降。SP-DMD將能夠最大程度反映流場信息的穩定模態保留下來，并對保留的模態幅值進行調整，從而達到對原始粘彈性流場降維的效果。

利用SP-DMD所獲得的100個模態與30個模態對流場進行重構后的粘彈性圓柱繞流流場結構的空間分布如圖9(c)和圖9(d)所示。為了驗證SP-DMD有著較好的優化降維效果，采用標準DMD所獲得的302個模態對流場進行重構后的粘彈性圓柱繞流流場結構的空間分布如圖9(b)所示。為了進行比較，原始的瞬態流場如圖9(a)所示。

圖9 原始流場與DMD和SP-DND重構流場的對比

由圖9可知，用30個模態對圓柱繞流流場進行重構就可以將原始流場的整體流動形態重現，并且能夠較好地刻畫出局部較大尺度的流場結構。當模態的數量增加到100時，可以看到重構流場中顯示出更多的局部細節。然而，標準DMD算法由302個模態才能較好地刻畫出原始流場的整體形態。由此可見，SP-DMD能夠對將流場的非關鍵信息進行剔除，從而保留并提取具有穩定狀態的關鍵流場結構。

4 結論

針對粘彈性圓柱繞流渦團結構的研究，引入了DMD與稀疏促進相結合的特征分析方法。該方法在提取粘彈性流場關鍵模態時，不只考慮幅值，而且通過計算模態對流場貢獻率來剔除非關鍵模態，從而留下產生影響較大的關鍵模態，使重構后的粘彈性流場仍然可以保留原始粘彈流流場中的重要動態信息。因此，利用SP-DMD可以更好地解析粘彈性圓柱繞流流場的特性，為粘彈性圓柱繞流流場的特性分析研究提供了新的視角。在未來的發展中，能否將包含控制輸入的模型預測控制與SP-DMD相結合來解決實際問題，還有待進一步測試和改進。