基于拓?fù)渚W(wǎng)格方法的多鈍體流致振動(dòng)分析

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

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

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

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

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

1 物理模型

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

圖1 物理模型Fig.1 Physical model

2 數(shù)值方法

2.1 流體控制方程

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

(1)

(2)

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

(3)

τij=2vTSij

(4)

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

(5)

(6)

(7)

(8)

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

2.2 流致振動(dòng)控制方程

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

(9)

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

(10)

(11)

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

(12)

2.3 網(wǎng)格模型

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

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

(a)子塊向上運(yùn)動(dòng)時(shí)，h變小，當(dāng)h

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

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

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

3 結(jié)果與討論

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

3.1 二維雙圓柱流致振動(dòng)

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

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

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

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

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

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

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

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

圖6 d=2D時(shí)串列雙圓柱尾渦結(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時(shí)，上下游雙圓柱所在網(wǎng)格子區(qū)域之間存在固定網(wǎng)格流場(chǎng)區(qū)域，包含兩組耦合界面，這與d=2D的工況有所不同。如圖7所示，在計(jì)算中拓?fù)渚W(wǎng)格成功實(shí)現(xiàn)了d=2.5D時(shí)雙圓柱流致振動(dòng)網(wǎng)格變化。

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

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

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

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

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

3.2 二維三圓柱流致振動(dòng)

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

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

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

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

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

圖12 串列雙圓柱(d=2D,d=2.5D)、 三圓柱(d=2.5D)流致振動(dòng)響應(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 三維多柱體流致振動(dòng)

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

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