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

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

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

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

為充分了解多鈍體的流致振動規(guī)律，本文采用拓?fù)渚W(wǎng)格變形技術(shù)，結(jié)合耦合界面，實現(xiàn)流體與多個運動鈍體之間的耦合計算。分別對二維雙圓柱和三圓柱、三維雙圓柱模型進行了數(shù)值計算，得到了圓柱的振幅和頻率響應(yīng)，觀察圓柱尾流旋渦結(jié)構(gòu)，并與實驗進行比較分析。

1 物理模型

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

圖1 物理模型Fig.1 Physical model

2 數(shù)值方法

2.1 流體控制方程

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

(1)

(2)

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

(3)

τij=2vTSij

(4)

式中湍動黏度vT的相關(guān)定義為

(5)

(6)

(7)

(8)

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

2.2 流致振動控制方程

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

(9)

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

(10)

(11)

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

(12)

2.3 網(wǎng)格模型

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

為避免柱體運動造成網(wǎng)格拉伸或擠壓變形，保證在數(shù)值計算過程中維持較高的網(wǎng)格質(zhì)量，通過網(wǎng)格拓?fù)渥冃渭夹g(shù)實現(xiàn)運動子塊上下網(wǎng)格的自動生成或塌縮，避免網(wǎng)格單元體積的劇烈變化。如圖2所示，振動柱體所在子塊的上下側(cè)為動態(tài)層，h為動態(tài)層與鄰近網(wǎng)格(n)的距離。當(dāng)流致振動發(fā)生時，子塊內(nèi)網(wǎng)格整體運動，過程如下：

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

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

其中hmin和hmax分別是根據(jù)整體網(wǎng)格質(zhì)量而設(shè)定的最小和最大網(wǎng)格厚度。拓?fù)渥冃渭夹g(shù)能自動對網(wǎng)格進行更新，有效解決旋轉(zhuǎn)、拉伸和擠壓造成的網(wǎng)格變形問題，有效減小網(wǎng)格變形造成的計算誤差。

圖2 拓?fù)渚W(wǎng)格及耦合界面Fig.2 Topological mesh and coupling interface

3 結(jié)果與討論

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

3.1 二維雙圓柱流致振動

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

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

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

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

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

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

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

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

圖6 d=2D時串列雙圓柱尾渦結(jié)構(gòu)Fig.6 Vortex structure of two cylinders in tandem with d=2D

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

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

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

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

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

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

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

3.2 二維三圓柱流致振動

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

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

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

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

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

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

3.3 三維多柱體流致振動

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

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