基于變分模態(tài)分解的壁湍流擬序結(jié)構(gòu)形態(tài)研究

王文康, 潘 翀, 王晉軍

(北京航空航天大學(xué) 教育部流體力學(xué)重點(diǎn)實(shí)驗(yàn)室, 北京 100191)

0 引 言

盡管譜分析在尺度分離中有廣泛的應(yīng)用，但它也有一定的局限性。首先，傅立葉變換隱含地假設(shè)被分析信號(hào)達(dá)到統(tǒng)計(jì)穩(wěn)態(tài)。對(duì)于通過(guò)單點(diǎn)測(cè)量(例如熱線風(fēng)速儀)獲得的速度脈動(dòng)的時(shí)間序列，可以通過(guò)增加采樣的持續(xù)時(shí)間來(lái)滿足這一條件。然而，對(duì)于像粒子圖像測(cè)速( PIV )這樣的場(chǎng)測(cè)量，視野范圍( FOV )通常受到測(cè)試儀器的限制不可能太大，因此瞬時(shí)脈動(dòng)速度在空間上的分布，如流向速度u(x)，往往不能滿足統(tǒng)計(jì)穩(wěn)態(tài)條件。這使得針對(duì)短信號(hào)的傅里葉分析無(wú)法準(zhǔn)確解析大尺度(低波數(shù))處的能譜。其次，在中低Re的情況下，尺度分離并不明顯[4, 7]，這使得選擇合適的分離尺度相當(dāng)困難。另外，人工給定一個(gè)大、小尺度之間的固定的截?cái)嚅撝担谖锢砩吓c湍流結(jié)構(gòu)的時(shí)變特性相矛盾。

經(jīng)驗(yàn)?zāi)B(tài)分解( Empirical Mode Decomposition, EMD )[8]及其變體是另外一種尺度分離選擇，其本身具有處理非穩(wěn)態(tài)信號(hào)的能力。EMD是通過(guò)一系列的篩選(sifting)過(guò)程將時(shí)間序列信號(hào)分解為有限數(shù)量的本征模態(tài)函數(shù)(Intrinsic Mode Functions, IMFs，在下文中也簡(jiǎn)稱(chēng)為模態(tài))，而不依賴(lài)于預(yù)先確定的模態(tài)形式, 即可以自適應(yīng)地捕捉到信號(hào)的局部尺度。在湍流尺度分析中得到了較為廣泛的應(yīng)用[9-10]，Cheng[11]等最近對(duì)通過(guò)DNS模擬的槽道湍流速度場(chǎng)進(jìn)行了二維EMD (B-EMD )分析, 得到了流場(chǎng)中的條帶結(jié)構(gòu)和流向渦結(jié)構(gòu)，并證明這些結(jié)構(gòu)的尺度沿法向高度自相似增長(zhǎng)的特性，這與傳統(tǒng)的附著渦假設(shè)所得到的結(jié)論是一致的。

由Dragomiretskiy和Zosso[12]提出的變分模態(tài)分解( VMD )是傳統(tǒng)EMD的最新升級(jí)。VMD方法的基本思想是使所分解出的模態(tài)應(yīng)該具有最緊致的帶寬。在此基礎(chǔ)上，VMD用嚴(yán)格的數(shù)學(xué)定義明確地構(gòu)造了一個(gè)凸優(yōu)化問(wèn)題。由于最小帶寬規(guī)則的限制，不同VMD模態(tài)在尺度空間中相互重疊最小，這避免了傳統(tǒng)EMD的模態(tài)摻混問(wèn)題。這種優(yōu)勢(shì)使得VMD在信號(hào)處理方面具有巨大的潛力。 在此基礎(chǔ)上，Wang[13-14]的工作對(duì)VMD方法在多維多分量信號(hào)分解上的應(yīng)用進(jìn)行了探索，提出了多分量準(zhǔn)二維VMD(MC-QBVMD)方法并用于湍流邊界層擬序結(jié)構(gòu)的分析，驗(yàn)證了這種方法在湍流多尺度分析中的有效性。然而這種自適應(yīng)尺度分解方法得到的模態(tài)所對(duì)應(yīng)的物理結(jié)構(gòu)的本質(zhì)仍不清楚。

本文將利用MC-QBVMD方法，對(duì)直接數(shù)值模擬(DNS)得到的湍流邊界層進(jìn)行分解。重點(diǎn)對(duì)分解得到的各階模態(tài)中的擬序結(jié)構(gòu)的尺度和幾何特性進(jìn)行分析，揭示其物理本質(zhì)。本文分為以下幾個(gè)部分，第1節(jié)將介紹變分模態(tài)分解方法及其衍生的MC-QBVMD方法；第2節(jié)將介紹所使用的DNS邊界層數(shù)據(jù)；第3節(jié)將對(duì)分解的到的結(jié)果進(jìn)行分析討論；第4節(jié)給出結(jié)論。

1 變分模態(tài)分解方法

1.1 一維變分模態(tài)分解

在這里我們簡(jiǎn)要的介紹一維VMD算法的核心思想和實(shí)現(xiàn)方法，以方便后面的討論。關(guān)于算法本身的更多細(xì)節(jié)可以參考文獻(xiàn)[12]。

VMD算法的核心是將一列時(shí)序信號(hào)f(t)分解成一系列本征模態(tài)函數(shù)(Intrinsic mode function, IMF)的線性疊加，即：

(1)

其中uk為IMF，其特征是極值點(diǎn)的數(shù)量與過(guò)零點(diǎn)的數(shù)量相等或僅相差一個(gè)。模態(tài)數(shù)量K需要提前指定。VMD算法的核心目的是要使得各階模態(tài)在尺度上充分分離，這也就引出了VMD算法的核心優(yōu)化函數(shù)，即：?jiǎn)蝹€(gè)模態(tài)的帶寬要最窄。這一最窄帶寬準(zhǔn)則可以表述成以下的優(yōu)化函數(shù)，即：

(2)

ak(t)=uk(t)+iHuk(t)

(3)

公式(3)中，H為希爾伯特變換(Hilbert Transform)。分析函數(shù)ak(t)的傅里葉變換是一個(gè)只在正頻率有幅值的函數(shù)。 e-iωkt的作用是將整個(gè)信號(hào)的頻域上移動(dòng)ωk， 其中ωk可以認(rèn)為是模態(tài)uk的中心頻率，或者載波頻率。將經(jīng)過(guò)平移的分析信號(hào)關(guān)于t求偏導(dǎo)數(shù)，再求L2范數(shù)，相當(dāng)于求ak(t)的頻譜關(guān)于ωk的轉(zhuǎn)動(dòng)慣量，因此可以反應(yīng)uk的帶寬。由此可見(jiàn)，公式(2)所表達(dá)的優(yōu)化函數(shù)的含義是：求一個(gè)原始信號(hào)f的模態(tài)的集合{uk}和中心頻率ωk的組合，使得所有模態(tài)帶寬相加之和最小。對(duì)該目標(biāo)函數(shù)的約束條件是，所得全部模態(tài)得線性疊加能夠完全還原原始信號(hào)。

求解公式(2)中得優(yōu)化問(wèn)題即可得到原始信號(hào)f的模態(tài)，下文稱(chēng)為VMD模態(tài)?？梢允褂米兿虺俗臃?Alternate direction method of multipliers, ADMM)來(lái)求解這一優(yōu)化問(wèn)題。為了簡(jiǎn)明起見(jiàn)不再給出推導(dǎo)求解的過(guò)程，僅僅給出Dragomiretskiy等[12]給出的迭代求解公式，如下所示：

(4)

(5)

1.2 多分量準(zhǔn)二維變分模態(tài)分解

瞬時(shí)湍流速度場(chǎng)是三維三分量的，但是湍流邊界層的流場(chǎng)本身的特殊性，比如脈動(dòng)速度沿法向的分布具有不均勻性、流場(chǎng)中的擬序結(jié)構(gòu)在不同方向上的尺度差別巨大等等導(dǎo)致同時(shí)在兩個(gè)以上的維度上進(jìn)行分解非常困難。比如二維FFT無(wú)法對(duì)x-y平面的速度場(chǎng)進(jìn)行分解，因?yàn)榉ㄏ虿粷M足周期性條件；二維EMD分解x-z平面流場(chǎng)時(shí)往往無(wú)法有效的提取其中的條帶狀結(jié)構(gòu)等等[14]。

針對(duì)這個(gè)問(wèn)題,Wang等[13]提出了多分量準(zhǔn)二維變分模態(tài)分解(MC-QBVMD)算法來(lái)分解二維多分量數(shù)據(jù)。這個(gè)方法的思路是，對(duì)于一個(gè)二維多分量的數(shù)據(jù)信號(hào)f(x,y)={fm(x,y),m=1∶M}，假設(shè)x方向?yàn)榉纸夥较?，y方向?yàn)檩o助方向，M為分量個(gè)數(shù)，則f(x,y)可以看成一系列x方向上的空間信號(hào)在y方向上的組合。假設(shè)由于擬序結(jié)構(gòu)的存在，同一階模態(tài)中主方向上的尺度沿輔助方向沒(méi)有變化。因此在其迭代求解的過(guò)程中，可以先在主方向上進(jìn)行尺度分解，然后在輔助方向上對(duì)不同位置處的中心尺度進(jìn)行重新平衡，使得每階模態(tài)的主方向上的尺度在輔助方向上保持一致。對(duì)于不同的速度分量，MC-QBVMD要求不同分量之間的尺度保持一致，這是因?yàn)閷?duì)于模態(tài)uk中某個(gè)擬序結(jié)構(gòu)引起的脈動(dòng)，其在不同分量上的投影的尺度應(yīng)該是一致的。

MC-QBVMD的目標(biāo)函數(shù)可以表示為：

(6)

上式假設(shè)x為主分解方向。由于f為二維模態(tài)，因此目標(biāo)函數(shù)中的L2范數(shù)被替換成Frobenius范數(shù)‖·‖F(xiàn)來(lái)估計(jì)二維模態(tài)的帶寬。對(duì)于一個(gè)二維矩陣uij，對(duì)比公式(2)與MC-QBVMD的目標(biāo)函數(shù)可以發(fā)現(xiàn)，MC-QBVMD的目標(biāo)函數(shù)相當(dāng)于不同分量和不同y處的一維空間信號(hào)的帶寬之和。

由于MC-QBVMD的目標(biāo)函數(shù)相當(dāng)于不同空間位置處的多分量一維信號(hào)的帶寬的線性疊加，因此其求解的步驟與VMD非常相似。使用ADMM求解式(6)，可以得到迭代解為：

(7)

(8)

公式(7)說(shuō)明在MC-QBVMD中，uk的更新步相當(dāng)于對(duì)f(x,y)在每個(gè)y位置和每個(gè)分量進(jìn)行一維的VMD分解，而ωk的更新步不但在各個(gè)分量之間進(jìn)行平均，同時(shí)在不同的y位置進(jìn)行平衡，使得每一階一維模態(tài)um,k在整個(gè)y域上具有相同的中心頻率。關(guān)于MC-QBVMD算法的更多細(xì)節(jié)可參照Wang等[13]的工作。

2 DNS湍流邊界層數(shù)據(jù)

我們將在后文中分析馬德里理工大學(xué)(Universidad Politécnica de Madrid) Javier Jimenez 教授課題組公開(kāi)的DNS湍流邊界層數(shù)據(jù)(https://torroja.dmt.upm.es/turbdata/blayers)。該數(shù)據(jù)為體發(fā)展邊界層流動(dòng)，雷諾數(shù)跨度Reτ=400～2000，模擬了光滑平板之上的一個(gè)六面體內(nèi)的速度場(chǎng)，展向采用的是周期性邊界條件，流向和法向?yàn)榉侵芷谛赃吔鐥l件。 在流向和法向采用交錯(cuò)網(wǎng)格三點(diǎn)緊致差分方法。展向上使用傅里葉譜進(jìn)行計(jì)算，并使用2/3律來(lái)消除混淆誤差。求解壓力使用的Poisson方程中采用的是二階有限差分格式。時(shí)間推進(jìn)使用的是半隱式三階Runge-Kutta方法。關(guān)于這兩套DNS數(shù)據(jù)的更多細(xì)節(jié)可以參考原始文獻(xiàn)[15]。

在當(dāng)前工作中，我們選取6個(gè)特征雷諾數(shù)下的流向-法向平面作為研究對(duì)象，流向上選取各個(gè)雷諾數(shù)對(duì)應(yīng)的位置上下游各6δ(δ為特征雷諾數(shù)當(dāng)?shù)氐倪吔鐚雍穸?的范圍，法向上選取y≤δ的區(qū)域，基本涵蓋了邊界層中所有流向尺度的結(jié)構(gòu)。在當(dāng)前選擇的流向范圍內(nèi)，雷諾數(shù)的變化在±3%以?xún)?nèi)，因此后面將采用流向范圍中心位置處的邊界層參數(shù)作為全局參數(shù)進(jìn)行無(wú)量綱。各個(gè)工況的流向平均速度型U和流向脈動(dòng)速度型σu如圖1所示。不同工況下的邊界層參數(shù)如表1所示。

(b)

3 變分模態(tài)分解結(jié)果

3.1 瞬時(shí)流場(chǎng)結(jié)構(gòu)分解

使用MC-QBVMD對(duì)湍流邊界層中流向-法向(x-y)平面內(nèi)的多分量瞬時(shí)速度場(chǎng)進(jìn)行尺度分離。

表1 邊界層參數(shù)Table 1 The parameters of turbulent boundary layers

需要說(shuō)明的是，湍流邊界層中脈動(dòng)速度的統(tǒng)計(jì)量在平行壁面平面內(nèi)是準(zhǔn)均勻的，在垂直壁面方向是非均勻的。這種情況適合應(yīng)用準(zhǔn)二維分解方法進(jìn)行分解，分解的主方向是流向，輔助方向是法向。而傳統(tǒng)的二維分解方法，例如二維FFT、或者二維EMD，由于在兩個(gè)方向上均進(jìn)行分解，所以并不適用于流向-法向平面內(nèi)的二維數(shù)據(jù)集。參考Wang的工作[13]，模態(tài)的階數(shù)定為K=2, 罰系數(shù)定為α=100，λ=0。

圖2展示了湍流邊界層(Reτ=1750)某單幀速度場(chǎng)以及MC-VMD分解得到的兩階模態(tài)。這里我們僅僅展示原視野范圍的一小部分，以便更加清晰地顯示模態(tài)中的流動(dòng)結(jié)構(gòu)。 第一階模態(tài)(圖2(b))很好的捕捉到了原始流場(chǎng)中的大尺度結(jié)構(gòu)，即沿流向排布的交替排列的高低速流體，其流向尺度達(dá)到了邊界層厚度的量級(jí)。這種結(jié)構(gòu)與前人觀察到的大尺度結(jié)構(gòu)在x-y平面的是一致的[16-17]，其特征都是向流向傾斜，一端與壁面接觸，另一端向?qū)?shù)區(qū)延伸。由于u、v兩個(gè)分量是同時(shí)分解的，因此可以看到這些附著于壁面的結(jié)構(gòu)與大尺度的掃掠噴射事件密切相關(guān)：即高速區(qū)對(duì)應(yīng)大尺度的Q4事件，低速區(qū)對(duì)應(yīng)大尺度的Q2事件。

圖2 湍流邊界層MC-QBVMD瞬時(shí)模態(tài)(Reτ=1750)。 (a) 原始速度場(chǎng); (b) 一階MC-QBVMD模態(tài); (c) 二階MC-QBVMD模態(tài)。矢量場(chǎng)代表(u,v,w)，背景云圖為u。Fig.2 Instantaneous MC-QBVMD modes of turbulent boundary layer. (a) original field; (b) 1st mode; (c) 2nd mode. The vectors represent (u,v,w) and the color contours shows magnitude of u.

第二階模態(tài) (圖2(c))顯示了瞬時(shí)的小尺度模態(tài)，其特征尺度比第一階模態(tài)小一個(gè)量級(jí)，可用內(nèi)尺度進(jìn)行標(biāo)度(λx+～300)。小尺度模態(tài)中可以觀察到近壁面(y+<50)頻繁發(fā)生的小尺度噴射和掃掠事件。這些事件在流向交替排列，并且與大尺度結(jié)構(gòu)一樣與流向方向存在傾角。這些小尺度結(jié)構(gòu)可能與湍流邊界層內(nèi)區(qū)的自循環(huán)機(jī)制有關(guān)[18]。

3.2 大小尺度模態(tài)結(jié)構(gòu)特征

在觀察MC-QBVMD模態(tài)的瞬態(tài)結(jié)構(gòu)(如圖2)時(shí)，我們發(fā)現(xiàn)兩階模態(tài)均表征出相類(lèi)似的結(jié)構(gòu)，即：從壁面傾斜延伸向外的、沿流向交替排列的高低速條狀結(jié)構(gòu)。為了進(jìn)一步揭示這種結(jié)構(gòu)的形態(tài)學(xué)特征，我們統(tǒng)計(jì)了大/小尺度模態(tài)中不同速度分量的傾斜結(jié)構(gòu)的結(jié)構(gòu)傾角隨高度的變化。計(jì)算結(jié)構(gòu)傾角的方法如下，首先計(jì)算各個(gè)速度分量在不同法向高度處的兩點(diǎn)相關(guān)，即：

Ruu(rx,ry;y)=

(9)

(b) v分量脈動(dòng)

(c) w分量脈動(dòng)

對(duì)圖4顯示出的結(jié)果可以做以下三點(diǎn)更加深入的討論。第一，仔細(xì)分析結(jié)構(gòu)傾角的雷諾數(shù)影響可以發(fā)現(xiàn)，原始速度場(chǎng)的結(jié)構(gòu)傾角呈現(xiàn)出較為明顯的雷諾數(shù)依賴(lài)性，但在將原始速度場(chǎng)分解后，大/小尺度模態(tài)的結(jié)構(gòu)傾角的雷諾數(shù)依賴(lài)性都有不同程度的減弱，尤其是在大尺度u和w分量(圖4(a，c))上。這似乎暗示，對(duì)于平行壁面的速度分量，其傾斜的大尺度結(jié)構(gòu)的幾何特征在所研究的雷諾數(shù)范圍內(nèi)具有“普適性”，該推測(cè)需要在更高的雷諾數(shù)下去做進(jìn)一步的驗(yàn)證。第二，從象限分析的角度出發(fā)，v分量的結(jié)構(gòu)傾角達(dá)到90°這一現(xiàn)象，可以理解為Q2/Q4事件對(duì)v分量湍流脈動(dòng)的貢獻(xiàn)達(dá)到均衡。圖4(b)顯示，相比于大尺度而言，小尺度v分量的結(jié)構(gòu)傾角都在一個(gè)更低的流層上收斂至90°，且不同雷諾數(shù)之間的差異很小，這與小尺度的內(nèi)區(qū)自維持循環(huán)只限于近壁流區(qū)、且具有雷諾數(shù)無(wú)關(guān)性的傳統(tǒng)觀點(diǎn)相吻合。另一方面，大尺度v分量的結(jié)構(gòu)傾角收斂至90°的法向高度，似乎隨著雷諾數(shù)的增大而增大。這也許可以說(shuō)明，以Q2噴射事件為特征的大尺度低速結(jié)構(gòu)是更為常見(jiàn)的近壁大尺度結(jié)構(gòu)；或者可以說(shuō)明，壁面條件對(duì)大尺度v分量的影響延伸至很高的流區(qū)，且與雷諾數(shù)相關(guān)。第三，附著渦模型(attached eddy model)假設(shè)附著渦在平行于壁面的速度分量(u和w)上的流向尺度正比于附著渦所處的法向高度[16]，即：lx～y，這一假設(shè)暗示了附著渦在不同的流層上具有相同的結(jié)構(gòu)傾角。圖4(a, c)顯示，原始流場(chǎng)u、w分量的結(jié)構(gòu)傾角確實(shí)在一定的y+以上不再改變，但大/小兩種尺度的u、w分量的結(jié)構(gòu)傾角并未呈現(xiàn)出類(lèi)似的收斂趨勢(shì)。因此，不能簡(jiǎn)單認(rèn)為在大尺度模態(tài)中發(fā)現(xiàn)的從壁面延伸向外區(qū)的條狀結(jié)構(gòu)就是附著渦的物理表征。

4 結(jié) 論

本工作應(yīng)用MC-QBVMD對(duì)湍流邊界層x-y平面脈動(dòng)速度場(chǎng)進(jìn)行了分解。MC-QBVMD的大尺度模態(tài)對(duì)應(yīng)的是與LSM相似的結(jié)構(gòu)。小尺度模態(tài)則比較復(fù)雜，在近壁區(qū)傾斜的Q2和Q4事件比較相似，而在對(duì)數(shù)區(qū)則近似于垂直的渦結(jié)構(gòu)。MC-QBVMD小尺度模態(tài)的在流向尺度和形態(tài)特征上與流向渦的尺度非常接近，且兩者在不同雷諾數(shù)的工況下均有較好的對(duì)應(yīng)關(guān)系，表明x-y平面的MC-QBVMD小尺度模態(tài)應(yīng)該對(duì)應(yīng)于原始流場(chǎng)中的小尺度流向渦結(jié)構(gòu)在x-y平面的投影。

在當(dāng)前工作中我們只將原始流場(chǎng)分解為大/小兩種尺度的模態(tài)，這樣做雖然便于表征和分析，但可能并不能夠最好的反映湍流的多尺度特性。更多關(guān)于模態(tài)數(shù)量以及更高維度的分解的工作需要在后續(xù)的工作中開(kāi)展。