李 璇，蘇 進(jìn)

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

0 引言

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

針對粘彈性圓柱繞流問題，為實(shí)現(xiàn)粘彈流高效的結(jié)構(gòu)分析和動(dòng)態(tài)預(yù)測，引入一種基于稀疏促進(jìn)優(yōu)化的DMD方法(SP-DMD)，與傳統(tǒng)的DMD相比，該方法可以剔除粘彈性流體中的非關(guān)鍵模態(tài)。在粘彈性流場演化中，SP-DMD使粘彈性流場中關(guān)鍵模態(tài)的數(shù)量和降維近似殘差之間處于平衡狀態(tài)，從而更利于分析和發(fā)現(xiàn)那些能夠?qū)φ硰椥詧A柱繞流場發(fā)展起關(guān)鍵性作用的因素。此外，當(dāng)模態(tài)數(shù)量增加時(shí)，SP-DMD重構(gòu)的流場可以得到更多的局部流動(dòng)細(xì)節(jié)，并且與原始流場進(jìn)行對比，局部較小尺度的渦團(tuán)形態(tài)以及空間分布的特點(diǎn)都能夠很好地重現(xiàn)出來。

1 DMD算法

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

(x1x2...xm+1)。

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

從快照序列中構(gòu)造兩個(gè)快照矩陣：

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

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

線性算子A記錄著快照序列的動(dòng)態(tài)特性。在粘彈性流體問題中，矩陣A中包含大量數(shù)據(jù)，且一般情況下為復(fù)數(shù)。每個(gè)快照xi中分量的數(shù)量通常遠(yuǎn)大于快照的數(shù)量，即n>>m。

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

因此，求出線性算子A是DMD的主要目標(biāo)。

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

假設(shè)數(shù)據(jù)總維數(shù)是m，上述SVD通過選擇適當(dāng)?shù)慕財(cái)嘀祌來減少數(shù)據(jù)矩陣的維數(shù)，被消除的冗余項(xiàng)就用rem(redundancy term)表示，其維數(shù)是m-r。

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

(1)

動(dòng)力學(xué)方程(1)所生成的每一組數(shù)據(jù)為：

且每個(gè)αi是對應(yīng)DMD模態(tài)的振幅。同樣，在矩陣形式中有：

則有

X≈ΦDαVand。

(2)

(3)

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

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

(4)

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

2 稀疏化動(dòng)態(tài)模態(tài)分解(SP-DMD)方法

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

2.1 SP-DMD核心思想

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

具有非零幅值的DMD模態(tài)可以通過求解如下凸優(yōu)化問題而得到：

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

(5)

式(5)中，矩陣E反映了振幅向量的稀疏結(jié)構(gòu)信息。E的列均是單位向量，其非零元素對應(yīng)于具有零幅值的DMD模態(tài)。獲得具有非零幅值的DMD模態(tài)后，通過改變其幅值實(shí)現(xiàn)降維近似。

將優(yōu)化問題(3)進(jìn)行修改，通過給式(3)中的目標(biāo)函數(shù)J(α)增加一個(gè)正則化項(xiàng)card(α)來解決稀疏性的問題，這就懲罰了未知振幅向量Dα中非零元素的數(shù)量：

(6)

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

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

(7)

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

2.2 交替方向乘子法

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

minf(x)+g(z)，

s.t.Ax+Bz=c。

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

引入增廣拉格朗日函數(shù)：

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

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

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

從初始點(diǎn)開始進(jìn)行迭代，直到滿足如下優(yōu)化條件和停止準(zhǔn)則：

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

zk+1-zk+12≤εdual。

其中，εprim，εdual表示停止準(zhǔn)則的閾值。

3 數(shù)值算例

3.1 粘彈性圓柱繞流模型參數(shù)設(shè)置

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

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

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

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

表1 數(shù)值模擬參數(shù)

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

出口處按充分發(fā)展條件處理，出口處的水平速度的法向?qū)?shù)為0，其中n是單位法向量，即

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

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

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

對于本構(gòu)方程采用Oldroyd-B方程的格子Boltzmann模型，應(yīng)力張量的分布函數(shù)上下壁面的邊界條件應(yīng)用非平衡態(tài)外推方法，計(jì)算區(qū)域的上下邊界的應(yīng)力分量為

τyy=0。

入口和出口邊界的應(yīng)力分量采用充分發(fā)展邊界條件，即

3.2 粘彈性圓柱繞流數(shù)值結(jié)果

不同的Re(Reynolds)數(shù)對應(yīng)的Wi(Weissenberg)不同，且Wi是Re的函數(shù)。根據(jù)3.1節(jié)所述，模擬出Re=1000時(shí)的粘彈性圓柱繞流流場的數(shù)據(jù)，并且對流體彈性特征數(shù)Wi=1的數(shù)據(jù)進(jìn)行研究。每一時(shí)刻的流場數(shù)據(jù)構(gòu)成一個(gè)301×601的網(wǎng)格。

在研究粘彈性圓柱繞流問題的過程中，首先要將二維數(shù)據(jù)轉(zhuǎn)化成一維，由此可知每一時(shí)刻所對應(yīng)的空間維數(shù)高達(dá)180 901維。由于各種運(yùn)動(dòng)形態(tài)和不同結(jié)構(gòu)的尾流流態(tài)會(huì)表現(xiàn)出不同的能量和頻率，僅通過對瞬態(tài)流場的觀察，很難對粘彈性圓柱繞流流場的特性進(jìn)行深入理解。因此研究過程中，選擇了連續(xù)1 000個(gè)時(shí)刻的二維粘彈性流場數(shù)據(jù)進(jìn)行DMD算法處理，可以得到302個(gè)模態(tài)。得到的302個(gè)DMD模態(tài)的特征值在復(fù)數(shù)坐標(biāo)系內(nèi)的分布情況如圖2所示，這些特征值進(jìn)行對數(shù)化后的分布情況如圖3所示。

圖2 DMD未對數(shù)化的特征值分布 圖3 DMD對數(shù)化的特征值分布

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

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

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

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

幅值與頻率之間的關(guān)系可以定義為DMD頻譜，在研究粘彈性流場問題的過程中，DMD頻譜能夠?qū)⒘鲌鼋Y(jié)構(gòu)所包含的能量情況反映出來，因此又將其稱之為能譜。采用兩種方法所獲得的模態(tài)幅值與特征值虛部(頻率)的關(guān)系如圖5所示，對其進(jìn)行局部放大如圖6所示。利用兩種方法所得到的模態(tài)幅值與特值實(shí)部(衰減率)的關(guān)系如圖7所示，其局部放大如圖8所示。

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

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

圖8 采用 DMD 與 SP-DMD 得到的模態(tài)幅值與特征值實(shí)部(衰減率)的關(guān)系的局部放大

由圖5、6、7、8可知，部分具有較大衰減率的非穩(wěn)定的DMD模態(tài)所對應(yīng)的流場結(jié)構(gòu)承載著較多的能量，其對流場初期演化的影響較強(qiáng)，但在后期的演化過程中，影響力明顯下降。SP-DMD將能夠最大程度反映流場信息的穩(wěn)定模態(tài)保留下來，并對保留的模態(tài)幅值進(jìn)行調(diào)整，從而達(dá)到對原始粘彈性流場降維的效果。

利用SP-DMD所獲得的100個(gè)模態(tài)與30個(gè)模態(tài)對流場進(jìn)行重構(gòu)后的粘彈性圓柱繞流流場結(jié)構(gòu)的空間分布如圖9(c)和圖9(d)所示。為了驗(yàn)證SP-DMD有著較好的優(yōu)化降維效果，采用標(biāo)準(zhǔn)DMD所獲得的302個(gè)模態(tài)對流場進(jìn)行重構(gòu)后的粘彈性圓柱繞流流場結(jié)構(gòu)的空間分布如圖9(b)所示。為了進(jìn)行比較，原始的瞬態(tài)流場如圖9(a)所示。

圖9 原始流場與DMD和SP-DND重構(gòu)流場的對比

由圖9可知，用30個(gè)模態(tài)對圓柱繞流流場進(jìn)行重構(gòu)就可以將原始流場的整體流動(dòng)形態(tài)重現(xiàn)，并且能夠較好地刻畫出局部較大尺度的流場結(jié)構(gòu)。當(dāng)模態(tài)的數(shù)量增加到100時(shí)，可以看到重構(gòu)流場中顯示出更多的局部細(xì)節(jié)。然而，標(biāo)準(zhǔn)DMD算法由302個(gè)模態(tài)才能較好地刻畫出原始流場的整體形態(tài)。由此可見，SP-DMD能夠?qū)⒘鲌龅姆顷P(guān)鍵信息進(jìn)行剔除，從而保留并提取具有穩(wěn)定狀態(tài)的關(guān)鍵流場結(jié)構(gòu)。

4 結(jié)論

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