杜曉慶，林芝強(qiáng)，吳葛菲

(1.上海大學(xué) 土木工程系， 上海 200444; 2.上海大學(xué) 風(fēng)工程和氣動控制研究中心， 上海 200444)

圓柱型結(jié)構(gòu)在土木工程中應(yīng)用廣泛，且常常以柱群的形式出現(xiàn)，如纜索承重橋并列索、煙囪群、冷卻塔群、海洋立管、多分裂導(dǎo)線等[1-2]。當(dāng)流體流經(jīng)鈍體時，出現(xiàn)剪切層分離、失穩(wěn)及轉(zhuǎn)捩等現(xiàn)象，并在尾流中形成交替脫落的旋渦[3]。為了進(jìn)一步了解復(fù)雜的流動結(jié)構(gòu)及潛在的流態(tài)轉(zhuǎn)變機(jī)制，從高維非線性流場數(shù)據(jù)中提取低維空間下的流動特征，是當(dāng)前研究的熱點[4]。

近些年來，隨著高精度數(shù)值模擬技術(shù)的迅猛發(fā)展，低維降階技術(shù)在復(fù)雜流動分析中的應(yīng)用日益廣泛[4-5]。其中，本征正交分解(POD)、動力學(xué)模態(tài)分解(DMD)是較典型的兩種方法。POD將流場分解為1組完全正交的模態(tài)，并根據(jù)其模態(tài)能量大小進(jìn)行排序。該方法由Berkooz等[6]首先提出，用來識別湍流中的相干結(jié)構(gòu)，捕捉流場中的主導(dǎo)模態(tài)。DMD方法是庫普曼算子的一種近似，Schmid[7]首先提出DMD，用來分析試驗或數(shù)值模擬得到的高維、大規(guī)模的流場數(shù)據(jù)，提取的低維子空間的特征值和特征向量包含流場的主要動力學(xué)特征，相比于POD，DMD得到的模態(tài)對應(yīng)單一的頻率與增長率?？芗覒c等[8]對DMD的研究進(jìn)展進(jìn)行總結(jié)，提出DMD未來的發(fā)展方向。Sakai等[9]、葉坤等[10]、Bai等[11]將POD與DMD分別應(yīng)用于單圓柱與并列雙圓柱，分析其穩(wěn)定性以及尾流流態(tài)的流動特性。

目前，對串列雙圓柱繞流的研究主要集中在不同雷諾數(shù)下氣動性能分析及臨界間距比的確定[12-14]。Zdravkovich[1]對串列雙圓柱的繞流場特性進(jìn)行試驗研究，并根據(jù)間距比將柱間流態(tài)分成：單一鈍體、交替再附、定常再附、間歇性再附和雙渦脫流態(tài)。Sumner[2]則將串列雙圓柱流態(tài)分成3類：單一鈍體流態(tài)、剪切層再附流態(tài)、雙渦脫流態(tài)。然而，目前的研究僅通過試驗或CFD對流場的流動特性進(jìn)行分析，這意味著得到的數(shù)據(jù)通常包含多種信息成分，并且很難確定主導(dǎo)流態(tài)變化的流場特征。

為了澄清串列雙圓柱的流態(tài)轉(zhuǎn)變機(jī)制，分析各流態(tài)的模態(tài)特征，探究DMD流場重構(gòu)的精度及誤差。本文在雷諾數(shù)Re=100下，通過數(shù)值模擬對間距為1.1D～5D的串列雙圓柱進(jìn)行了研究，將數(shù)值模擬得到的氣動力系數(shù)與文獻(xiàn)結(jié)果進(jìn)行了比較；并采用DMD方法對不同間距比下的串列雙圓柱繞流場的動力學(xué)特性進(jìn)行分析，通過對比分析3種流態(tài)下的各個主導(dǎo)模態(tài)，闡明3種流態(tài)的變化規(guī)律；最后根據(jù)DMD的主導(dǎo)模態(tài)建立降階模型，從而對雙圓柱繞流場進(jìn)行了流場重構(gòu)，并對其重構(gòu)的精度與誤差進(jìn)行了分析。

1 數(shù)值方法

1.1 基本控制方程

本文采用Fluent中的層流模型，在雷諾數(shù)Re=100下，對串列雙圓柱繞流進(jìn)行了模擬。對于此二維流動問題，流體視為黏性不可壓流體，其在直角坐標(biāo)系oxy下的連續(xù)方程和動量方程(N-S方程)分別為：

(1)

(2)

(3)

式中：ρ為流體密度，p為流體壓力，μ為黏性系數(shù)，u和v分別x、y方向的速度分量，F(xiàn)x和Fy分別為流體在x、y方向的體力。

1.2 計算模型

圖1為串列雙圓柱的計算模型。兩圓柱的直徑均為D，中心間距為P，來流風(fēng)速為U0，雷諾數(shù)為Re=100。圖中的CD1和CL1分別表示上游圓柱的平均阻力系數(shù)和平均升力系數(shù)，C′D1和C′L1分別代表上游圓柱的脈動阻力系數(shù)與脈動升力系數(shù)。其中，阻力系數(shù)以順流方向為正， 升力系數(shù)以向上為正， 下游圓柱同理。本文共計算了6種不同的中心間距：P/D=1.1、1.5、2、3、4、5。

圖1 雙圓柱計算模型

數(shù)值模擬采用O型計算域，計算域直徑為46D，阻塞率為2%，入流面采用速度入口邊界條件；出流面采用自由出流邊界條件，圓柱壁面為無滑移圓柱表面采用無滑移壁面條件。計算模型采用結(jié)構(gòu)化網(wǎng)格，圓柱周向400個單元，徑向180層單元；近壁面最小網(wǎng)格為0.000 1D，近壁面y+≈1；無量綱時間步Δt*為0.01(Δt*=ΔtU0/D，其中Δt為計算時間步，U0為來流風(fēng)速)。計算模型網(wǎng)格總單元數(shù)從17萬至20萬不等，網(wǎng)格方案見圖2。

圖2 計算域網(wǎng)格方案

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

假設(shè)1個線性動力系統(tǒng)來映射當(dāng)前流場到后續(xù)流場，即假設(shè)流場xi+1可以通過流場xi的線性映射表示為

xi+1=Axi

(4)

利用1到m時刻的流場快照，構(gòu)建快照矩陣：

(5)

根據(jù)(4)式的假設(shè)，可得：

(6)

任意時刻流場可進(jìn)行重構(gòu)或預(yù)測：

2 模型驗證

為了驗證本文所采用的計算方法和計算參數(shù)的合理性，本文首先以雷諾數(shù)為Re=100的固定單圓柱為研究對象進(jìn)行計算模型的結(jié)果驗證。 表1總結(jié)了本文固定單圓柱的網(wǎng)格驗證結(jié)果以及文獻(xiàn)[13,16-18]已發(fā)表的試驗和數(shù)值模擬結(jié)果， 考慮周向網(wǎng)格數(shù)量、無量綱時間步長、阻塞率(B=D/H，H為計算域橫向?qū)挾?對計算結(jié)果的影響?？梢园l(fā)現(xiàn)，周向網(wǎng)格數(shù)量從256增加至400， 模擬結(jié)果基本一致。考慮到雙圓柱周圍流場的復(fù)雜性，選擇周向網(wǎng)格數(shù)量為400的工況進(jìn)一步確認(rèn)無量綱時間步長與阻塞率對模擬結(jié)果的影響?？梢钥闯?，在無量綱時間步長為0.01、 阻塞率為2%時， 模擬結(jié)果更接近文獻(xiàn)結(jié)果， 因此最終采用A3工況的網(wǎng)格參數(shù)進(jìn)行后續(xù)模擬計算。

表1 固定單圓柱模型的網(wǎng)格方案和結(jié)果驗證(Re=100)

為了進(jìn)一步驗證本文計算模型的可靠性，圖3～圖5給出了間距比P/D=1.1～5串列雙圓柱的平均阻力系數(shù)，脈動升力系數(shù)以及下游圓柱的St數(shù)，并與文獻(xiàn)的結(jié)果進(jìn)行對比。通過與文獻(xiàn)[12]中相同雷諾數(shù)的結(jié)果比較可知，本文上、下游圓柱的平均阻力系數(shù)、脈動升力系數(shù)以及St數(shù)與文獻(xiàn)[12]的誤差均在5%以內(nèi)。對比文獻(xiàn)[19-20]相近雷諾數(shù)的結(jié)果可發(fā)現(xiàn)，隨著間距比的變化，本文的氣動力系數(shù)與St數(shù)的變化趨勢與文獻(xiàn)結(jié)果一致。

圖3 串列雙圓柱平均阻力系數(shù)

圖4 串列雙圓柱脈動升力系數(shù)

圖5 下游圓柱升力的斯托羅哈數(shù)

3 計算結(jié)果及分析

3.1 繞流場特性

圖6為P/D=1.5、3、4時串列雙圓柱的平均風(fēng)壓系數(shù)、平均流線、平均風(fēng)速比圖，其中平均風(fēng)速比為流場內(nèi)局部平均風(fēng)速U與來流風(fēng)速U0的比。由圖6可見，在P/D=1.5、3時，上游圓柱的尾流區(qū)沒有完整的渦脫落產(chǎn)生；而在P/D=4時，上游圓柱尾流區(qū)產(chǎn)生了穩(wěn)定的渦脫落。結(jié)合圖3、4可知，氣動力在P/D=3～4之間發(fā)生了突變，說明該現(xiàn)象與上下游圓柱尾流發(fā)生變化直接相關(guān)。

由圖6(a)可見，兩圓柱間距較近時，上游圓柱的剪切層包住下游圓柱，兩圓柱間隙存在兩個回流區(qū)，回流區(qū)會經(jīng)流下游圓柱的迎風(fēng)面，從而導(dǎo)致下游圓柱迎風(fēng)面的風(fēng)壓為負(fù)值，由上文分析可知，此時下游圓柱的平均阻力為負(fù)值；圖6(b)與圖6(a)類似，此時兩圓柱間距比增大，兩圓柱間的回流區(qū)進(jìn)一步發(fā)展，回流的規(guī)律性逐漸被打破，導(dǎo)致間隙間的負(fù)壓絕對值不斷減小，從而使下游圓柱受到的負(fù)阻力絕對值相應(yīng)減??；當(dāng)P/D=4時，由圖6(c)可見，此時兩圓柱間不存在回流區(qū)，上游圓柱的剪切層不再再附到下游圓柱表面，而是在上游圓柱的尾流中形成充分發(fā)展的漩渦，導(dǎo)致上游圓柱的背風(fēng)面受到一定的負(fù)壓，而迎風(fēng)面存在較強(qiáng)的正壓，這就使得上游圓柱所受阻力更大，也就導(dǎo)致氣動力在P/D=3～4發(fā)生了明顯的突變。

圖6 串列雙圓柱的平均風(fēng)壓系數(shù)、流線、風(fēng)速比

因此，通過如上分析可知，本文成功捕捉了串列雙圓柱的3種流態(tài)，即當(dāng)P/D≤2時為單一鈍體流態(tài)，剪切層在下游圓柱后方形成充分發(fā)展的漩渦；當(dāng)P/D=3時為剪切層再附流態(tài)，此時上游圓柱分離的剪切層再附到下游圓柱表面；當(dāng)P/D=4、5時為雙渦脫流態(tài)，此時上游圓柱尾流區(qū)產(chǎn)生旋渦脫落。

3.2 模態(tài)分解特征值

為了進(jìn)一步闡明3種流態(tài)的變化機(jī)制，本文采用DMD對3種流態(tài)的渦量場進(jìn)行分析。分別對P/D=1.5、3、4時的串列雙圓柱穩(wěn)定階段提取了6個周期的非定常渦量場數(shù)據(jù)，每個周期采樣32次，選取前4階主要模態(tài)進(jìn)行分析。這3個間距比分別代表單一鈍體流態(tài)、剪切層再附流態(tài)以及雙渦脫流態(tài)。

圖7為3個間距比下DMD特征值的實部與虛部在單位圓上的分布。橫坐標(biāo)為特征值的實部，對應(yīng)模態(tài)的增長率。橫坐標(biāo)為正，表明該模態(tài)系數(shù)發(fā)散；橫坐標(biāo)為負(fù)，表明該模態(tài)系數(shù)收斂?？v坐標(biāo)為特征值的虛部，對應(yīng)模態(tài)的頻率。由于本文所考慮的階段為穩(wěn)定極限階段，因此3個間距比下串列雙圓柱模態(tài)特征值均在單位圓內(nèi)或單位圓上，對應(yīng)周期性模態(tài)或穩(wěn)定模態(tài)[10]。圖7給出了所選取的前4階模態(tài)，模態(tài)特征值均處于單位圓上，表明本文所有模態(tài)均為周期性模態(tài)或穩(wěn)定模態(tài)。同時，表2給出了提取的前4階模態(tài)所對應(yīng)的頻率與增長率，可見，前4階模態(tài)增長率都為0或接近于0，說明選取的前4階模態(tài)均屬于穩(wěn)定模態(tài)。

圖7 串列雙圓柱的模態(tài)特征值

表2 DMD主模態(tài)的增長率和頻率

3.3 模態(tài)分析

圖8為3個間距比下DMD模態(tài)振幅與頻率的關(guān)系，其中縱坐標(biāo)振幅的大小即為標(biāo)準(zhǔn)DMD的排序標(biāo)準(zhǔn)。圖8中每個點都代表流場的1個空間模態(tài)。本文通過對渦量場快照使用動力學(xué)模態(tài)分解，采用了文獻(xiàn)[15]發(fā)展的排序準(zhǔn)則，依據(jù)模態(tài)對流場的貢獻(xiàn)進(jìn)行排序，并且提取了前4階模態(tài)進(jìn)行分析，這些模態(tài)包含大多數(shù)的動態(tài)信息。在P/D=1.5和3時，提取的模態(tài)沒有完全按照振幅的大小進(jìn)行排序，而是按照該模態(tài)對流場的貢獻(xiàn)進(jìn)行排序。模態(tài)1與模態(tài)2對流場的貢獻(xiàn)最大，處于主導(dǎo)地位，而其他大部分模態(tài)對流場貢獻(xiàn)相對較小。

圖8 DMD模態(tài)振幅與頻率的關(guān)系

圖9為P/D=1.5、3、4時DMD捕捉到的前4階動力學(xué)模態(tài)。可見，隨著模態(tài)階數(shù)的增加，渦的尺度逐漸減小。當(dāng)P/D=1.5時，模態(tài)分解的結(jié)果與單圓柱類似，區(qū)別在于模態(tài)1的剪切層相對于單圓柱延伸的更長且尺度更大，而對于模態(tài)2～4來說，單圓柱的模態(tài)中渦結(jié)構(gòu)更接近圓柱表面；P/D=3時，其各階模態(tài)沿流向的尾流長度均長于P/D=1.5時的串列雙圓柱以及單圓柱；而在P/D=4時，其各階模態(tài)上游圓柱的尾流與P/D=1.5、2以及單圓柱類似，而下游圓柱后方的近尾流區(qū)則有很大的不同。

圖9 渦量場的DMD前4階主導(dǎo)模態(tài)

進(jìn)一步分析可知，模態(tài)1所對應(yīng)的頻率為0，近似于平均流場，對應(yīng)于采樣段內(nèi)流場的平均值。當(dāng)P/D=1.5、3時，上游圓柱的剪切層在下游圓柱后方分離；當(dāng)P/D=4時，上游圓柱的剪切層在上游圓柱后方分離，這表明模態(tài)1體現(xiàn)了串列雙圓柱流動剪切層的分離在3種流態(tài)間的轉(zhuǎn)變過程。模態(tài)2所對應(yīng)的頻率與圓柱的渦脫主頻相同，表明該模態(tài)與圓柱漩渦脫落的動力學(xué)特性有關(guān)；此外，3個間距比下串列雙圓柱尾流中存在相似的渦結(jié)構(gòu)，如圖9中紅色虛線框所示，該渦結(jié)構(gòu)由下游圓柱尾流逐漸平移至上游圓柱尾流中，表征了3種流態(tài)下旋渦脫落的內(nèi)在演變機(jī)制。模態(tài)3則為2倍頻模態(tài)，對應(yīng)的頻率為圓柱漩渦脫落頻率的2倍，且呈反對稱結(jié)構(gòu)，該模態(tài)表明了圓柱后方交替脫落的旋渦在尾流中的發(fā)展。模態(tài)4所對應(yīng)的頻率為圓柱旋渦脫落頻率的3倍，且呈正對稱結(jié)構(gòu)，為流場的高階諧波模態(tài)，表征了流場中的高頻行為，是對尾流沿流向延伸這一特征的補(bǔ)充。

3.4 流場重構(gòu)與誤差分析

為了進(jìn)一步評估DMD對串列雙圓柱渦量場特征的提取效果，利用得到的DMD模態(tài)進(jìn)行渦量場重構(gòu)，選擇圓柱采樣段內(nèi)升力峰值時刻的渦量場進(jìn)行重構(gòu)。同時，采用均方根誤差對DMD模態(tài)重構(gòu)的渦量場誤差進(jìn)行全面評估，均方根誤差定義為

式中Np為重構(gòu)快照的數(shù)目，xCFD(i)和xDMD(i)分別為CFD計算與DMD模態(tài)重構(gòu)的第i個渦量場快照。

圖10為實際渦量場與DMD重構(gòu)渦量場的對比以及均方根誤差云圖。首先，將本文的3種流態(tài)與文獻(xiàn)[2]中Re=1.2×104的流態(tài)對比可知，均依次呈現(xiàn)單一鈍體流態(tài)、剪切層再附流態(tài)、雙渦脫流態(tài)。對于DMD重構(gòu)渦量場，單圓柱與P/D=1.5、3的串列雙圓柱使用前4階模態(tài)，P/D=4時的串列雙圓柱使用前6階模態(tài)進(jìn)行流場重構(gòu)，從而使均方根誤差達(dá)到相同量級。

圖10 重構(gòu)渦量場與真實渦量場的對比及均方根誤差云圖

由圖10可見，4種工況的重構(gòu)渦量場均能捕捉到主要結(jié)構(gòu)，達(dá)到準(zhǔn)確的重構(gòu)效果。P/D=4的渦量場重構(gòu)需要更多的模態(tài)來達(dá)到與其他工況相同的誤差，原因在于其上下游圓柱尾流均發(fā)生旋渦脫落，僅用前4階模態(tài)重構(gòu)的渦量場誤差較大。重構(gòu)誤差均集中在圓柱尾流處，也就是漩渦脫落發(fā)展的區(qū)域，這說明沒有參與重構(gòu)的高階模態(tài)主要對圓柱尾流中旋渦發(fā)展的區(qū)域有貢獻(xiàn)。同時可以發(fā)現(xiàn)，DMD可以較精確的捕捉低頻大尺度的旋渦，但在捕捉高頻小尺度渦上存在一定誤差。

為進(jìn)一步分析所用模態(tài)數(shù)對重構(gòu)誤差的影響，本文對單圓柱與串列雙圓柱3種流態(tài)的重構(gòu)渦量場的平均誤差進(jìn)行分析，如圖11所示?？梢姡S著重構(gòu)渦量場所用模態(tài)數(shù)的增加，各工況平均誤差均大幅減小，當(dāng)所用模態(tài)數(shù)達(dá)到6～7時，平均誤差基本接近于0。

圖11 重構(gòu)渦量場的平均誤差

4 結(jié) 論

本文利用動力學(xué)模態(tài)分解(DMD)方法，對串列雙圓柱的渦量場進(jìn)行模態(tài)分解，在低維空間分析了不同模態(tài)的流場特征;并探討了不同間距串列雙圓柱流態(tài)轉(zhuǎn)變的內(nèi)在機(jī)制;最后基于分解的主模態(tài)建立了降階模型，并對串列雙圓柱流場重構(gòu)精度以及誤差進(jìn)行分析，主要結(jié)論如下：

1)各階模態(tài)的流場特征存在顯著差異。第1階模態(tài)表示平均流場，頻率為0；第2階模態(tài)所對應(yīng)的頻率與圓柱的渦脫主頻相等，該模態(tài)與圓柱旋渦脫落的動力學(xué)特性有關(guān)；第3階模態(tài)則為2倍頻模態(tài)，該模態(tài)表明了圓柱后方交替脫落的旋渦在尾流中的發(fā)展。第4階模態(tài)為3倍頻模態(tài)，是流場的高階諧波模態(tài)，反應(yīng)了尾流中更小尺度的旋渦在尾流中沿橫流向的延伸。此外，隨著模態(tài)階數(shù)的增加，渦的尺度逐漸減小。

2)3種流態(tài)的模態(tài)特征也有顯著區(qū)別。單一鈍體流態(tài)(P/D=1.1～2)時，各階模態(tài)特征與單圓柱類似；剪切層再附流態(tài)(P/D=3)時，各階模態(tài)沿流向的尾流長度均長于單一鈍體流態(tài)以及單圓柱；而在雙渦脫流態(tài)(P/D=4～5)時，各階模態(tài)上游圓柱的尾流與前2種流態(tài)以及單圓柱類似，而下游圓柱后方的近尾流區(qū)則有很大的不同。

3)主導(dǎo)模態(tài)均具有明確的物理意義。頻率為0的模態(tài)體現(xiàn)了串列雙圓柱流動剪切層的分離在3種流態(tài)間的轉(zhuǎn)變過程；與圓柱渦脫主頻相對應(yīng)的模態(tài)中存在相似的渦結(jié)構(gòu)，該渦結(jié)構(gòu)由下游圓柱尾流逐漸平移至上游圓柱尾流中，表征了3種流態(tài)的內(nèi)在演變機(jī)制。

4)流場重構(gòu)結(jié)果表明，基于DMD的降階模型可以較好地重構(gòu)串列雙圓柱渦量場。當(dāng)子模態(tài)數(shù)量相同時，單一鈍體與剪切層再附流態(tài)的重構(gòu)精度高于單圓柱，而雙渦脫流態(tài)的重構(gòu)精度則低于單圓柱;流場重構(gòu)的誤差主要集中在尾流區(qū)，由旋渦的交替脫落導(dǎo)致，具有較強(qiáng)非線性。