基于拓撲網格方法的多鈍體流致振動分析

丁 林， 鄒 瑞， 張 力， 鄒群峰

(1. 重慶大學 低品位能源利用技術及系統教育部重點實驗室，重慶 400044; 2.重慶大學 能源與動力工程學院，重慶 400044)

鈍體是工程中一種常見的非流線型結構，在一定的流速下，流體繞流鈍體后會在鈍體兩側交替脫落旋渦。旋渦脫落會在結構表面形成周期性脈動力，當鈍體固定方式為彈性支撐或允許發生形變時，脈動力將引起鈍體振動，即流致振動[1]。在一切與流體有關的機械與工程中，流致振動都是一個涉及安全性的重大問題，在工程實踐中受到廣泛關注。多鈍體流致振動是一個較為復雜的流固耦合過程，流場與多個運動結構的耦合計算還需進一步優化，這不僅有利于分析多鈍體流致振動動力學響應的影響參數，而且可以更加深入地了解多鈍體流致振動特性及旋渦演變過程，具有非常重要的理論研究意義與工程應用價值。

圓柱是一種常見的鈍體結構，目前關于單圓柱繞流流致振動問題的研究已經取得了非常豐碩的成果[2-3]，孫麗萍等[4]采用了改進的k-ε湍流模型，研究了加速度對渦激振動響應的影響，優化了數值模擬方法，實現對單圓柱雙自由度渦激振動較精確的數值模擬；陳正壽等[5]試驗發現，低質量比圓柱對應的鎖振區范圍要廣于高質量比圓柱，低質量比圓柱橫向振幅較大，渦激振動現象更為顯著。相比單圓柱繞流流致振動，在實際工程中多圓柱繞流流致振動的情況更為普遍。雙圓柱串行排列是最為簡單的多鈍體系統，固定雙圓柱繞流的研究已經取得諸多成果[6-7]，但是雙圓柱流致振動的研究卻有待深入，目前主要集中于最簡單的尾流振動模型[8-9]。雙圓柱流致振動特性與圓柱間距密切相關[10]，Kim等[11]實驗發現當雙圓柱距離超過2.7倍圓柱直徑時，上游圓柱振動規律與單圓柱流致振動相似；關德寶等[12]試驗發現，下游圓柱振幅隨圓柱間距的增加先增大后減小。由于雙圓柱流致振動問題本身的復雜性，大部分研究都以實驗為主。在非線性流固耦合問題的數值研究中，關鍵在于流體和鈍體之間的移動邊界的處理，任意拉格朗日-歐拉方法是目前比較成功的用于解決移動邊界問題的方法[13]。但是，當渦致振動或馳振發生時，鈍體具有較高振幅，傳統的動網格方法已無法滿足計算需要。

為充分了解多鈍體的流致振動規律，本文采用拓撲網格變形技術，結合耦合界面，實現流體與多個運動鈍體之間的耦合計算。分別對二維雙圓柱和三圓柱、三維雙圓柱模型進行了數值計算，得到了圓柱的振幅和頻率響應，觀察圓柱尾流旋渦結構，并與實驗進行比較分析。

1 物理模型

當彈性支撐的鈍體產生流致振動時，橫向振幅遠大于流向振幅，因此本文僅考慮鈍體的橫向振動，針對兩個或兩個以上鈍體的流致振動展開研究。以串行排列雙圓柱物理模型為例，如圖1所示，每個振動系統的主要參數有：剛性圓柱的直徑為D、長度為L，支撐彈簧的彈性系數為K和摩擦等引起的系統阻尼為C。本文在圓柱表面增加附屬粗糙帶，以此加強柱體流致振動[14]。雙圓柱軸心距離為d，均為單自由度振動系統，振動方向垂直于來流速度和圓柱軸向。

圖1 物理模型Fig.1 Physical model

2 數值方法

2.1 流體控制方程

本文通過求解非穩態雷諾平均納維-斯托克斯方程組(Unsteady Reynolds-averaged Navier-Stokes)，獲得繞流彈性支撐多鈍體流場的數值解。不可壓縮流體的連續性方程和動量方程為

(1)

(2)

式中:Ui為平均流速；v為分子運動黏度；Sij為應變率張量

(3)

τij=2vTSij

(4)

式中湍動黏度vT的相關定義為

(5)

(6)

(7)

(8)

為了使時間和空間上的數值離散具有二階精度，發散項、梯度項和對流項采用二階高斯積分格式，時間積分采用二階向后歐拉積分法，動量方程和連續性方程的求解采用壓力隱式算子分裂算法。

2.2 流致振動控制方程

流場中的振動圓柱簡化為質量-彈簧-阻尼系統，其單自由度振動方程為

(9)

式中：m為振動系統總的慣性質量，由振動柱體質量和三分之一彈簧質量組成；f(t)為垂直流向的流體作用力，在計算過程中通過對圓柱表面的壓力和黏性力積分獲得。運動方程采用二階混合顯隱式時間積分法進行求解，過程為

(10)

(11)

yn+1=yn+Δt·vn+1

(12)

2.3 網格模型

由于彈性支撐多鈍體繞流的流固耦合問題本身的復雜性，多鈍體流致振動的數值模擬方法還需進一步優化。在彈性支撐多鈍體流致振動數值計算過程中，為了實現鈍體與流體介質之間作用力和運動關系的傳遞，需要不斷調整運動壁面區域的計算網格。為了減小因網格扭曲而引起的計算誤差，本文網格模型采用了耦合界面(Coupling Interface，CI)[16]和拓撲網格變形技術[17-18]。拓撲網格是一種底層靜態的網格，該網格本身并不隨運動子塊產生運動，不會因拉伸或擠壓產生網格形變，其特點在于子塊區域外圍網格會隨子塊運動而自動生成或崩塌消失。將每一個柱體附近計算區域劃分成方形網格子塊，柱體位于子塊中心。在流致振動發生時，子塊內網格隨鈍體壁面整體運動，子塊上下網格通過拓撲網格變形技術進行更新；子塊左右通過CI界面與上下游計算區域網格連接，實現多柱體流致振動的數值模擬。CI界面用來連接多個不連續網格區域，不要求相鄰界面的網格節點一一匹配。在具有旋轉部件的數值計算中要求所有界面網格相互匹配往往非常困難，CI界面是處理這類問題的有效途徑。對于界面網格不連續的區域，分別進行網格劃分，然后使用一個或者多個CI界面進行連接。與滑移網格不同，CI界面包含主界面和副界面，使用加權插值進行流率計算和傳遞。

為避免柱體運動造成網格拉伸或擠壓變形，保證在數值計算過程中維持較高的網格質量，通過網格拓撲變形技術實現運動子塊上下網格的自動生成或塌縮，避免網格單元體積的劇烈變化。如圖2所示，振動柱體所在子塊的上下側為動態層，h為動態層與鄰近網格(n)的距離。當流致振動發生時，子塊內網格整體運動，過程如下：

(a)子塊向上運動時，h變小，當h

(b)子塊向下運動時，h增大，當h>hmax時，將會自動生成網格n-1。

其中hmin和hmax分別是根據整體網格質量而設定的最小和最大網格厚度。拓撲變形技術能自動對網格進行更新，有效解決旋轉、拉伸和擠壓造成的網格變形問題，有效減小網格變形造成的計算誤差。

圖2 拓撲網格及耦合界面Fig.2 Topological mesh and coupling interface

3 結果與討論

為了研究彈性支撐多鈍體的流致振動響應特性，同時驗證本文提出的數值方法尤其是動網格模型的可靠性，本文分別對二維雙圓柱和三圓柱、三維雙圓柱模型進行了數值計算。

3.1 二維雙圓柱流致振動

二維雙圓柱計算區域如圖3所示，頂部和底部距離為9D，與實驗水深相同，速度入口距上游圓柱軸心為25D，下游圓柱軸心距離壓力出口為25D，兩種不同間距的串列雙圓柱模型的軸心間距d分別為2D和2.5D。頂部和底面為壁面邊界。為了強化圓柱的流致振動，在圓柱表面沿軸向安裝兩條附屬粗糙帶，粗糙帶覆蓋16°圓柱表面，厚度為0.847 mm，數值計算中采用壁面函數進行求解。本文計算采用結構化網格，均經過網格無關性驗證，并將計算結果與實驗進行比較。

圖3 計算區域及邊界條件Fig.3 Calculation region and boundary conditions

3.1.1 雙圓柱軸心距離d=2D

當串列雙圓柱軸心距離d=2D時，雙圓柱所在方形網格子區域之間存在共用耦合界面。當圓柱發生流致振動時，圓柱所在子區域整體運動，并在動態層經過區域成功實現網格自動崩塌與生成，如圖4所示。計算雷諾數(Re)范圍為3×104≤Re≤10×104，對應折減速度(U*=U/fnD)范圍為3.84≤U*≤12.81，fn為上游圓柱固有頻率。

圖4 串列雙圓柱流致振動網格變化(d=2D)Fig.4 Mesh changes of two cylinders in tandem(d=2D)

雙圓柱流致振動的振幅和頻率響應如圖5所示，數值結果與實驗值基本吻合。雙圓柱振動響應曲線明確顯示出流致振動的各個分支：渦致振動(Vortex Induced Vibration, VIV)初始分支(3×104≤Re<4×104)、上部分支(4×104≤Re<8×104)、過渡分支(8×104≤Re<9.5×104)和馳振(Re≥9.5×104)。當Re<3×104時，在實驗和數值模擬中都沒有觀察到圓柱發生流致振動。在渦致振動初始分支，上游圓柱振幅的數值結果較實驗值高，這是由于本文將圓柱振動系統阻尼簡化為線性阻尼，與實際的非線性阻尼存在差異[19]，而該差異在流速較低時表現尤為明顯，使計算中渦致振動初始分支的激發速度小于實驗中初始分支的激發速度。當振動處于渦致振動上部分支的初始段(4×104≤Re<5.5×104)時，實驗結果顯示下游圓柱振幅陡降，幾乎降至0值附近，這是由于下游圓柱的振動受到上游圓柱振動和自身旋渦脫落的共同作用，當Re≤4×104時，下游圓柱主要受上游圓柱振動的影響，表現出與上游圓柱相似的頻率和振幅響應；當Re=5×104時，下游圓柱振動頻率突降，低于上游圓柱振動頻率，下游圓柱振幅明顯下降；當Re≥5.5×104時，下游圓柱主要受自身旋渦脫落的影響，頻率和振幅逐漸上升。但是數值模擬中卻沒有捕捉到此現象，這是由于圓柱振動系統阻尼為線性阻尼，渦致振動初始分支更易激發，下游圓柱振動主要受自身旋渦脫落的影響，因而為表現出振幅和頻率的明顯波動。除了上部分支初始段的差異外，數值模擬和實驗結果吻合較好。當Re≥9×104圓柱發生馳振后，振動頻率比穩定在1.0左右，模擬值與實驗值結果保持一致。

圖5 串列雙圓柱(d=2D)流致振動響應Fig.5 The FIM response of two cylinders in tandem(d=2D)

圓柱的近尾跡渦量分布是反應流致振動特性的重要參數，渦量分布進一步形成不同的圓柱尾渦結構，數值模擬在描述和分析圓柱尾渦結構時具有極大優勢。由雙圓柱振幅頻率響應(見圖5)分析可知，Re=3×104時，上下游圓柱位于初始分支，從圖6(a)中可以清晰觀察到上游圓柱的旋渦脫落為2S模式。兩個旋渦隨流體向下游移動，在通過下游圓柱區域時，順時針旋轉的旋渦緊貼下游圓柱上側移動，逆時針旋轉的旋渦緊貼下游圓柱正下方移動。由于雙圓柱之間的特定距離，導致來自上游圓柱的旋渦吸收了下游圓柱分離的旋轉方向相同的旋渦，從而阻止下游圓柱大尺度Kármán旋渦的形成，進而抑制流致振動的產生，這與實驗現象相吻合。圖6(b)所示為Re=6×104時的近尾跡旋渦結構，上下游圓柱均為2P形態，即一個振動周期內脫落兩對旋渦，這與Khalak等[20]在渦致振動上部分支觀察到的圓柱尾渦形態相同。圖6(c)所示為Re=9.3×104時的近尾跡旋渦結構，此時雙圓柱流致振動由渦致振動向馳振過渡。對于上游圓柱，大部分振動周期為2P+2S模式，尾渦具有多種形態，但在某些周期內，上游圓柱在上行過程中會多脫落一個附加旋渦。由于受到來自上游圓柱旋渦的強烈干擾，下游圓柱的尾渦形態難以捕捉。圖6(d)所示為Re=10×104時的近尾跡旋渦結構，在一個振動周期內，上游圓柱產生的脫體旋渦高達8個，但是下游圓柱由于受到上游圓柱的嚴重影響，尾渦形態不易觀察。由以上分析可見，采用本文數值方法能有效描述雙圓柱流致振動。

圖6 d=2D時串列雙圓柱尾渦結構Fig.6 Vortex structure of two cylinders in tandem with d=2D

3.1.2 雙圓柱軸心距離d=2.5D

當串列雙圓柱軸心距離d=2.5D時，上下游雙圓柱所在網格子區域之間存在固定網格流場區域，包含兩組耦合界面，這與d=2D的工況有所不同。如圖7所示，在計算中拓撲網格成功實現了d=2.5D時雙圓柱流致振動網格變化。

圖7 串列圓柱流致振動網格變化(d=2.5D)Fig.7 Mesh changes of two cylinders in tandem(d=2.5D)

圖8所示為上下游圓柱的振幅和頻率響應曲線，可以觀察到渦激振動的初始分支(3×104≤Re<4×104)、上部分支(4×104≤Re<8×104)、過渡分支(8×104≤Re<9.5×104)和馳振(Re≥9.5×104)，在各個分支切換時，振幅和頻率都發生了明顯變化，與d=2D時雙圓柱振動相比，d=2.5D時下游圓柱振動雖然也受到上游旋渦抑制，但抑制作用減弱，下游圓柱流致振動也因上游圓柱的存在表現出不同的振動特性。在渦激振動的初始分支，上游圓柱振幅的模擬數值也略高于實驗值。當上游圓柱由渦致振動初始分支過渡到上部分支后，下游圓柱振幅隨著Re增加而降低，并在Re=5×104附近出現極小值。Re>9.5×104時，上游圓柱發生馳振，下游圓柱振幅隨Re增加先降低后增加，這是由于當Re≈9.5×104時，上游圓柱脫落的漩渦剛好附著于下游圓柱表面，抑制了下游圓柱漩渦的形成和脫落，從而使下游圓柱的振幅出現明顯下降；在數值計算中也捕捉到了這一現象，因而下游圓柱表現出明顯低于上游圓柱的振幅。在Re=5.5×104附近，實驗結果顯示下游圓柱頻率突然降低，之后隨Re增加緩慢增加。在渦致振動上部分支，實驗和模擬結果均表明，上下游圓柱振動頻率隨著Re增加出現小幅度增大；當Re≥9.5×104時，上下游圓柱頻率比均趨于穩定，最終穩定在1.0左右。

圖8 串列雙圓柱(d=2.5D)流致振動響應Fig.8 The FIM response of two cylinders in tandem(d=2.5D)

以Re=6×104時雙圓柱振動為例，上下游圓柱振動均位于渦致振動上部分支。圖9為d=2.5D時串列雙圓柱尾流旋渦云圖，從圖中可以看出雙圓柱的尾渦形態為2P，上游圓柱的尾渦形態較下游圓柱更為清晰。下游圓柱尾流旋渦在形成過程中受到來自上游圓柱的脫體旋渦影響，旋渦在演變過程中與來自上游的旋渦相互融合,進而使下游圓柱表現出明顯區別于上游圓柱的振動特性。

圖9 d=2.5D時串列雙圓柱尾渦結構(Re=6×104)Fig.9 Vortex structure of two cylinders in tandem with d=2.5D(Re=6×104)

3.2 二維三圓柱流致振動

為了進一步驗證本文數值方法的正確性，對二維三圓柱進行了數值模擬，相鄰圓柱之間軸心距離d=2.5D，其余參數與二維雙圓柱相同。圖10所示為流致振動發生時，串列三圓柱網格變化情況，在整個計算過程中，網格并不發生扭曲變形，當柱體振幅較大時，可有效減少由此帶來的計算誤差。

圖10 串列三圓柱流致振動網格變化(d=2.5D)Fig.10 Mesh changes of three cylinders in tandem(d=2.5D)

圖11為三圓柱振幅和頻率響應曲線，與二維雙圓柱的結果相似，可以觀察到流致振動依舊產生了四個分支，在各個分支切換時，振幅和頻率變化明顯。在渦致振動初始分支，上游圓柱和中間圓柱振幅的模擬值高于實驗值。上游圓柱的振動受下游其他圓柱影響較小；與二維雙圓柱相比，對應雷諾數下，上游圓柱的振幅變化也較小；上游圓柱最大振幅比A/D達到2.70。中間圓柱和下游圓柱則明顯受到其他圓柱的影響，在Re=5×104附近，實驗結果顯示中間圓柱和下游圓柱振幅和頻率突然降低，但振幅下降幅度明顯小于二維雙圓柱。當Re=9×104時，中間圓柱振幅出現明顯波動，當Re>9.5×104時，振幅出現小幅下降；當Re>9×104時，下游圓柱振幅比逐漸穩定在2.5左右，而上中下游圓柱頻率比均趨于穩定于1.0左右。

圖11 串列三圓柱(d=2.5D)流致振動響應Fig.11 The FIM response of three cylinders in tandem(d=2.5D)

為深入分析不同間距、數量對圓柱流致振動的影響，圖12對雙圓柱(d=2D,d=2.5D)和三圓柱(d=2.5D)流致振動三組工況的振幅和頻率響應進行了對比。三組工況中上游圓柱振幅的主要區別體現在渦致振動上部分支和過渡區域。如圖12(a)所示，在上部分支(如Re=5×104)，對于雙圓柱工況，d=2D上游柱體振幅明顯高于d=2.5D對應振幅，三圓柱工況中上游圓柱振幅最大；在過渡區域(8×104≤Re<9.5×104)，d=2.5D雙圓柱和三圓柱的上游圓柱振幅均高于d=2D雙圓柱；Re=1×105馳振發生時，三組工況上游圓柱振幅均達到2.5D以上。對于雙圓柱工況中的下游圓柱，在Re=3×104時柱體均無明顯振動；在Re=5×104時，由于上游圓柱的影響，d=2.5D對應下游圓柱振幅出現降低；當柱體發生馳振時，d=2D對應下游圓柱振動加強，最大振幅遠高于d=2.5D工況。值得注意的是，三圓柱的中間圓柱與d=2.5D雙圓柱的下游圓柱位置類似，除了Re=3×104時，振幅變化趨勢極為相似；受到第三柱體的影響，中間圓柱在Re=3×104時振動激發，振幅明顯高于雙圓柱下游圓柱。對于三圓柱的下游圓柱，由于上游兩圓柱的存在，振幅相對較低，振幅響應曲線與其他圓柱有明顯區別。對于三組工況的頻率響應特性，主要區別體現在雷諾數較小時(Re≤5×104)，如圖12(b)所示。當Re=4×104時，雙圓柱(d=2D)和三圓柱的上游圓柱的頻率高于d=2.5D時；隨著雷諾數增大，上下游圓柱頻率比趨勢吻合，最終穩定在1.0左右。由于三圓柱相互影響，當Re=8×104時，三圓柱頻率不匹配，上游圓柱頻率比最大，中間圓柱頻率比最小。

圖12 串列雙圓柱(d=2D,d=2.5D)、 三圓柱(d=2.5D)流致振動響應Fig.12 The FIM response of two cylinders (d=2D,d=2.5D) and three cylinders (d=2.5D) in tandem

3.3 三維多柱體流致振動

本文數值方法也可用于進行三維多柱體流致振動計算，與二維計算相比，三維模擬需要更高的硬件條件，而且需時更長。圖13為三維雙圓柱計算模型，計算區域尺寸24D×11D×6D，上游圓柱軸心距離入口7D，圓柱軸心間距d=2D，圓柱長L=4D。流道四周均為壁面邊界，入口為均勻來流，出口為壓力出口。

圖13 三維雙圓柱計算模型Fig.13 Calculation model of 3-dimensional two cylinders

圖14所示為三維雙圓柱發生流致振動時的網格變化情況。計算過程中，圓柱所在長方體區域隨圓柱振動發生整體移動，并伴隨網格自動生成與崩塌。圖15所示為三維雙圓柱流致振動響應曲線，上下游圓柱出現了明顯的振動分區。在渦激振動初始分支，上游圓柱振幅明顯降低，更加趨于實驗數值。在Re=5×104附近，下游圓柱振幅出現小幅下降，但仍然高于實驗數值；下游圓柱振動頻率與上游圓柱近似，并未捕捉到試驗中下游圓柱振動頻率突降的現象。當Re≥9×104時，上下游圓柱產生馳振，上下游圓柱振動頻率逐漸穩定，最終穩定在1.0附近，與實驗結果相吻合。結合圖5可知，三維計算中上下圓柱振幅和頻率響應曲線與二維計算相比更接近實驗結果。由此可見，三維計算結果更接近于實驗值，證明了本文數值方法對于計算三維多圓柱流致振動的可行性。

圖14 三維雙圓柱流致振動網格變化(d=2D)Fig.14 Mesh changes of 3-dimensional two cylinders in tandem(d=2D)

圖15 三維雙圓柱流致振動響應Fig.15 The FIM response of 3-dimensional two cylinders

4 結 論

本文針對彈性支撐的多鈍體流致振動進行數值分析，獲得多鈍體流致振動規律，并與實驗進行比較，得出以下結論：

(1)采用耦合界面結合拓撲網格技術，對彈性支撐多鈍體流致振動進行二維和三維計算，數值結果與實驗結果相吻合，證明該方法是處理高振幅多鈍體流致振動的有效方法。

(2)對于d=2D的串列雙圓柱流致振動，圓柱振幅和頻率響應的數值結果與實驗測試趨勢基本一致，能清晰觀察到渦致振動初始分支和上部分支，并且當Re>8×104時，振動由渦致振動向馳振過渡；Re≥9.5×104時，流致振動發展為馳振，上下游圓柱振幅降低，頻率比穩定在1.0左右。另外，圓柱尾渦形態隨流致振動分支切換發生變化，當馳振發生時，下游圓柱的尾渦形態受上游圓柱影響難以捕捉。

(3)對于d=2.5D的串列雙圓柱流致振動，數值模擬與實驗結果清晰顯示出流致振動各個分支。與d=2D相比，Re=3×104時下游圓柱受到上游圓柱的抑制作用減弱。Re=6×104時，雙圓柱尾渦結構與d=2D時相同，上下游圓柱均為2P模式。

(4)對于d=2.5D的串列三圓柱流致振動，依舊可以觀察到流致振動的各個分支。上游圓柱受到其下游圓柱的影響較小；在Re=5×104附近，中間圓柱和下游圓柱取得最小值，在Re=9×104附近，中間圓柱振幅出現明顯波動，下游圓柱振幅逐漸趨于穩定。

(5)與二維模擬相比，本文數值方法對三維高振幅多柱體流致振動仍能進行有效計算，所得結果更為接近實驗值，但三維計算需要更高的硬件條件，如何提高三維計算速度將是下一步研究工作的重點。