串列雙圓柱繞流的模態特征及流場重構

杜曉慶，林芝強，吳葛菲

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

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

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

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

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

1 數值方法

1.1 基本控制方程

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

(1)

(2)

(3)

式中：ρ為流體密度，p為流體壓力，μ為黏性系數，u和v分別x、y方向的速度分量，Fx和Fy分別為流體在x、y方向的體力。

1.2 計算模型

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

圖1 雙圓柱計算模型

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

圖2 計算域網格方案

1.3 動力學模態分解

假設1個線性動力系統來映射當前流場到后續流場，即假設流場xi+1可以通過流場xi的線性映射表示為

xi+1=Axi

(4)

利用1到m時刻的流場快照，構建快照矩陣：

(5)

根據(4)式的假設，可得：

(6)

任意時刻流場可進行重構或預測：

2 模型驗證

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

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

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

圖3 串列雙圓柱平均阻力系數

圖4 串列雙圓柱脈動升力系數

圖5 下游圓柱升力的斯托羅哈數

3 計算結果及分析

3.1 繞流場特性

圖6為P/D=1.5、3、4時串列雙圓柱的平均風壓系數、平均流線、平均風速比圖，其中平均風速比為流場內局部平均風速U與來流風速U0的比。由圖6可見，在P/D=1.5、3時，上游圓柱的尾流區沒有完整的渦脫落產生；而在P/D=4時，上游圓柱尾流區產生了穩定的渦脫落。結合圖3、4可知，氣動力在P/D=3～4之間發生了突變，說明該現象與上下游圓柱尾流發生變化直接相關。

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

圖6 串列雙圓柱的平均風壓系數、流線、風速比

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

3.2 模態分解特征值

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

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

圖7 串列雙圓柱的模態特征值

表2 DMD主模態的增長率和頻率

3.3 模態分析

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

圖8 DMD模態振幅與頻率的關系

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

圖9 渦量場的DMD前4階主導模態

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

3.4 流場重構與誤差分析

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

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

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

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

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

為進一步分析所用模態數對重構誤差的影響，本文對單圓柱與串列雙圓柱3種流態的重構渦量場的平均誤差進行分析，如圖11所示?？梢?，隨著重構渦量場所用模態數的增加，各工況平均誤差均大幅減小，當所用模態數達到6～7時，平均誤差基本接近于0。

圖11 重構渦量場的平均誤差

4 結 論

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

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

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

3)主導模態均具有明確的物理意義。頻率為0的模態體現了串列雙圓柱流動剪切層的分離在3種流態間的轉變過程；與圓柱渦脫主頻相對應的模態中存在相似的渦結構，該渦結構由下游圓柱尾流逐漸平移至上游圓柱尾流中，表征了3種流態的內在演變機制。

4)流場重構結果表明，基于DMD的降階模型可以較好地重構串列雙圓柱渦量場。當子模態數量相同時，單一鈍體與剪切層再附流態的重構精度高于單圓柱，而雙渦脫流態的重構精度則低于單圓柱;流場重構的誤差主要集中在尾流區，由旋渦的交替脫落導致，具有較強非線性。