低雷諾數下串聯雙圓柱渦激振動機理的數值研究

郭曉玲， 唐國強,2， 劉名名， 呂 林,3， 滕 斌

(1．大連理工大學 海岸和近海工程國家重點實驗室，大連 116024；2．韓國海洋科學技術研究所，大田 305-343；3．大連理工大學 深海工程研究中心，大連 116024)

在海洋工程中，當水流經過海洋平臺的樁柱、支撐結構、立管及海底管線等這些非流線型結構時，會在尾流區內產生交替的旋渦脫落，使作用在結構上的流體作用力表現出明顯的周期性特征，進而誘發渦激振動(Vortex-Induced Vibration，VIV)。當渦旋的脫落頻率與結構的固有頻率相接近時，將發生“鎖定”現象。此時，結構物會以較大的振幅振動，導致結構物的壽命急劇降低。在過去的幾十年中，均勻流中單個圓柱的渦激振動問題得到廣泛的關注[1-6]。人們通過對高雷諾數實驗研究發現，在低質量-阻尼比(m*ξ)下，圓柱的位移響應可以分為三個階段[1]：初始分支(Initial branch)、上分支(Upper branch)及下分支(Lower branch)，而在高質量-阻尼比(m*ξ)情況下，圓柱的位移響應只存在初始分支和下分支[2]。通過強迫振動的實驗研究發現[3]：對于僅發生橫流向受迫振動的圓柱結構，其尾渦脫落模式通常可劃分為2P、2S和P+S三種形式，從Khalak等[1]中可以看到，位移響應的初始分支對應的渦脫落模式為2S，下分支對應的為2P模式。

對于海洋工程結構中近距離放置的多個圓柱的情況，由于圓柱之間存在明顯的水動力干涉作用，相應的渦激振動響應與單個圓柱的情況存在顯著的差別[7-11]。針對上游圓柱固定不動，下游圓柱發生渦激振動的實驗研究表明[7]，由于受到上游圓柱尾跡區的影響，下游圓柱開始發生鎖定時所對應的約化速度大于單個圓柱的情況。并且，下游圓柱的鎖定帶寬與兩圓柱的間距有關。例如，Brika等[7]的結果表明，當兩圓柱間的相對圓心距L/D=10.0時，以約化速度表征的下游圓柱的鎖定區間約為單圓柱情況的2倍。Carmo等[8]分別計算了二維(Re = 150)和三維(Re = 300)條件下，彈性支撐圓柱在上游固定圓柱尾流干涉作用下的橫流向單自由度振動問題。計算條件為圓心間距分別為LX/D=1.5、3.0、5.0和8.0，質量比為m*=2.0，阻尼比為ξ=0.007。數值結果表明，在Re=150條件下，當發生共振時，下游圓柱的振幅較單圓柱情況增大近50%，而對于三維情況，在不同的圓心間距下，下游圓柱的振動幅值隨約化速度的變化趨勢與對應的二維情況類似。在以往的研究工作中，考慮到圓柱發生渦激振動時順流向位移的振幅比較小，因此，順流向的振動往往被忽略[8-9]。而Vandiver等[10]的研究工作表明，順流向和橫流向存在明顯的相互作用。陳文曲等[11]也對這一方面開展了數值分析工作。但是，人們對低雷諾數流動中，處于固定圓柱尾跡干涉作用下，同時具有順流向和橫流向雙自由度的圓柱渦激振動問題仍缺乏深刻的認識，特別是圓柱間距及質量比對圓柱受力和渦激振動響應的影響還需要更加細致的研究。

以目前現有的計算能力，對大尺度(幾千米水深)范圍內的實際立管渦激振動開展全場三維數值分析還面臨很大的挑戰，計算時間過長，不能滿足實際工程設計的需要。但有關低雷諾數下，多圓柱渦激振動問題的數值研究對于發展可靠的數值模型以及建立快捷的近似工程分析預報方法是有借鑒意義的。因此本文的工作主要在于加深對多圓柱渦激振動現象物理本質的認識，更側重于基礎理論方面的研究。本文采用三步有限元數值方法對二維Navier-Stokes方程進行了求解，并結合任意拉格朗日-歐拉(ALE)動網格方法，建立和開發了相應的數值模型和計算程序，以Re=150為代表，重點研究了兩個等直徑串聯圓柱的渦激振動問題，其中上游圓柱固定不動，下游圓柱在彈性支撐和阻尼作用下允許同時發生順流向和橫流向振動。通過數值計算，來研究圓柱間隙比、質量比以及約化速度等因素對下游圓柱渦激振動特性的影響。

1 數值模型

1.1 流動控制方程

二維不可壓縮均勻粘性牛頓流體運動的基本控制方程為連續性方程和Navier-Stokes方程，在ALE參考坐標系下，可以表示成如下的無量綱形式：

(1)

(2)

其中，xi表示笛卡爾坐標(二維情況下i= 1、2，分別對應x和y方向)，ui為xi方向的流速分量，p為壓力，Re=UD/υ為雷諾數，U為均勻來流流速，D為圓柱直徑，υ為流體的運動學粘性系數,cj為對流速度，并有:

(3)

應用三步有限元方法[12-13]，可得到如下離散形式的動量方程：

(4)

(5)

(6)

(7)

對于時間推進，本文根據CFL(Courant-Friedrichs-Lewy)條件，采用如下的動態時間步長：

式中，Se為網格的面積，ue為網格中心點的流速，min(*)表示在計算域內取最小值，Cs為安全系數，文中取為Cs= 0.2。

由于上述離散格式具有高階Tayler-Galerkin性質，同時速度與壓力通過投影(Projection)過程進行了解耦，因此式(4)～式(6)以及式(7)均可采用統一的標準Galerkin有限元方法進行空間離散，相關的離散方法可參見文獻[12～14]。

當獲得流場和壓力場后，圓柱所受到的流體作用力可通過對圓柱表面的壓力及粘性剪切力進行表面積分求得，無量綱化的拖曳力系數CD和升力系數CL分別為：

(8)

(9)

式中，p,u,v分別為圓柱表面上某點的壓力、x方向速度、y方向速度，θ為該點與圓心的連線與x軸正向之間的夾角。

1.2 圓柱的運動方程

(10)

(11)

當圓柱的位移確定后，本文通過微分網格變換方法進行網格更新[15]，進而獲得新的網格坐標以及網格運動速度，用于下一時間步的流場求解。

1.3 計算模型及邊界條件

計算模型及邊界條件如圖1所示。兩等直徑圓柱串聯放置，無量綱的圓柱直徑為d=1。設坐標原點位于下游圓柱的初始圓心，兩圓柱圓心間距為LX，相應的無量綱參數為LX/D。計算域的順流向長度為55D，橫流向寬度為50D。在入口處指定無量綱的速度u=1，v=0；側壁采用了對稱邊界條件?u/?y=0，v= 0；出口處的速度邊界條件為?ui/?t+c?ui/?xi=0，其中c為局部平均流速；圓柱表面施加不可滑移邊界條件u=dx/dt，v=dy/dt。計算中，在出口處指定相對壓力p= 0，在其他邊界采用?p/?n=0的壓力邊界條件，n為指出流體域的外法向單位矢量。在初始時刻，流場中的速度及相對壓力分布均設為零(即初始速度場滿足連續方程)。

圖1 計算域及邊界條件示意圖

2 模型驗證

2.1 網格收斂性驗證

表1 網格收斂性驗證結果

2.2 圓柱橫流向渦激振動的數值驗證

為進一步驗證本文數值模型的可靠性，同時考慮到在之前的工作中已對固定單圓柱繞流問題進行了詳細的驗證[13]，以下僅對單個圓柱(Isolated)橫流向自激振動以及上游圓柱固定，下游圓柱(直徑同于上游圓柱)的橫流向自激振動問題分別進行計算。在計算的過程中，均采用了相同的雷諾數Re=150、質量比m*=2.0和阻尼比ξ=0.007。圖2給出了單圓柱橫流向最大振幅(Amax/D)隨約化速度(Ur)的變化及其與Carmo等[8]的對比結果。對于串聯雙圓柱的情況，考慮了LX/D= 1.5和3.0兩種條件下，下游圓柱的橫流向最大振幅對約化速度的依賴關系。從圖2中可以看出，本文的數值計算結果(present)與已有數值結果吻合良好，表明本文的數值模型具有良好的數值精度。

3 計算結果及分析

在Carmo等[8]的數值分析工作中，考慮了上游圓柱固定不動的情況下，下游圓柱只在橫流向發生渦激振動的問題。本文將對下游圓柱同時發生橫流向和順流向雙自由度渦激振動的問題進行研究。為便于與文獻[8]的結果進行對比分析，本文選擇與該文相同的雷諾數Re=150、阻尼比ξ=0.007以及圓心間距比LX/D=3.0、5.0和8.0，在約化速度Ur= 3.0 ～ 12.0的范圍內對下游圓柱的兩自由度渦激振動響應特性開展研究。同時，本文進一步以m*= 5.0、10.0、20.0為例，研究質量比對下游圓柱發生渦激振動時動力響應的影響。

圖2 單圓柱及下游圓柱橫流向最大振幅Amax /D隨約化速度Ur的變化(與文獻[8]的對比，Re=150，m* = 2.0，ξ = 0.007)

3.1 圓心間距比的影響

根據固定串聯圓柱的不同渦脫落模態[16-17]，本文選擇了三個代表性的圓心間距比，LX/D= 3.0、5.0和8.0，并與相同雷諾數Re=150下的單圓柱(Isolated)渦激振動情況進行了比較，分析了圓柱的位移響應、振動頻率、受力特性、尾渦脫落模式及運動軌跡的情況。

3.1.1 下游圓柱的位移及振動頻率

圖3(a)以質量比m*=10.0、阻尼比ξ=0.007為例，給出不同圓心間距下，下游圓柱橫流向位移均方根YRMS/D隨Ur的變化。從圖中可以看出，當Lx/D=3.0時，其橫流向位移響應曲線同單圓柱類似，只有一個峰，而Lx/D=5.0和8.0的YRMS/D曲線均出現了兩個峰值。在Lx/D=3.0時，下游圓柱的最大位移均方根比單圓柱的情況增大近50%，而對于Lx/D=5.0和8.0，最大位移均方根約等于單圓柱的位移。從圖3(a)中還可以看出，在Ur> 7.0的高約化速度范圍內，與單圓柱的位移響應相比，下游圓柱的橫流向位移依然維持在比較高的水平，并且Lx/D=5.0和8.0的位移均方根在幅值和變化趨勢上都很接近，二者均比Lx/D=3.0所對應的位移明顯偏大。

圖3 不同圓心間距比LX / D下，下游圓柱渦激振動響應特性隨約化速度Ur的變化關系(Re = 150，m* = 10.0，ξ = 0.007)

由圖3(b)中下游圓柱的順流向平均位移XM/D隨約化速度Ur的變化曲線可以看到，隨著Ur的增大，下游圓柱的順流向平均位移也逐漸增大。同時，由于受到上游固定圓柱遮蔽效應的影響，下游圓柱的順流向平均位移在各圓心間距下都小于單圓柱的情況。對于Lx/D=3.0的情況，在Ur≤5.5時，下游圓柱的XM/D出現負值，即下游圓柱向上游偏移，當發生鎖定時，XM/D明顯增大，在Ur>7.0后，XM/D基本隨Ur呈線性增加。對于Lx/D=5.0和8.0的情況，順流向平均位移與橫流向位移均方根的變化趨勢類似，在發生鎖定前各出現一個次峰，鎖定后，二者接近重合，并且隨著Ur的增大基本呈線性增加，在量值上與Lx/D=3.0的情況接近。

圖4給出了不同圓心間距比下，下游圓柱相對振動頻率f/fn隨Ur的變化關系及其與單圓柱情況的比較。對于圓柱結構的渦激振動問題，當圓柱的振動頻率與固有頻率接近時(即f/fn趨近于1.0，以下以0.95≤f/fn≤1.05為界)，一般可認為發生鎖定[18]。圖4表明，單圓柱的鎖定區間約為5.0 ≤Ur≤6.5，而對于串聯雙圓柱的情況，在Lx/D=3.0時，下游圓柱的鎖定區間為6.0≤Ur≤7.0，當Lx/D=5.0和8.0時，在本文考慮的約化速度范圍內，鎖定區間分別為 6.0≤Ur≤12.0和6.5≤Ur≤12.0。圖4關于鎖定區間的結果與圖3 (a)的位移響應呈現出良好的對應關系，即圓柱振動的鎖定區間對應著橫流向位移比較大的區間。從圖4中同時也可以發現，隨著Lx/D的增大，開始發生鎖定的臨界約化速度也呈現出增大的趨勢，并且對于Lx/D=5.0和8.0的情況，鎖定區間明顯大于單圓柱及Lx/D=3.0的情況。值得注意的是，在Lx/D=5.0和8.0時，橫流向位移的第一個峰值出現在鎖定發生之前。

圖4 不同圓心間距比LX / D下，下游圓柱橫流向相對振動頻率f/fn隨約化速度Ur的變化關系(Re = 150，m* = 10.0，ξ = 0.007)

3.1.2 下游圓柱的受力特性

圖5 不同圓心間距比LX / D下，下游圓柱受力特性隨約化速度Ur的變化關系(Re=150，m*=10.0，ξ= 0.007)

3.1.3 尾渦脫落模式

圖6分析了Lx/D及Ur對尾渦脫落模式的影響作用。對于Lx/D=3.0的串聯雙圓柱，數值結果表明，在Ur<6.0的范圍內，上游固定圓柱后方沒有渦旋脫落，圓柱系統尾流區內的旋渦以圖6(a1)和6(a2)的模式進行整體脫落；當Ur≥6.0，即開始發生鎖定后，上游圓柱后方開始出現渦脫落，如圖6(a3)所示；隨著Ur的進一步增大，下游圓柱不再發生鎖定時，其尾渦脫落模式(參見圖6(a4))相對鎖定區內的圖6(a3)未見明顯變化。當Lx/D=5.0時，在Ur<5.5的范圍內，典型尾渦脫落模式如圖6(b1)所示，上游固定圓柱后方發生交替的渦旋脫落，而下游圓柱的尾渦則整齊地排成兩排；當Ur=5.5時，尾渦模式如圖6(b2)所示，整體上與圖6(a1)及6(a2)的情況類似；當Ur≥6.0后，尾渦模式如圖6(b3)和6(b4)所示，上、下游圓柱尾流區內分別可見旋渦的脫落。從前文的圖3和圖5中看到，在Lx/D=5.0，Ur=5.5的條件下，其位移和受力會突然變小，這正是由于當Ur=5.5時，尾渦脫落模式的突然改變所致。對于Lx/D=8.0，由于圓柱間距比較大，隨著Ur的變化，上下游圓柱的尾流區內始終同時存在旋渦的脫落，參見圖6(c1)-(c4)，上游圓柱脫落的渦旋經充分發展后，與下游圓柱的尾渦發生更為明顯的干涉作用。

3.1.4 下游圓柱的運動軌跡

為說明下游圓柱在x和y兩個自由度方向運動的相互影響作用，圖7給出了不同圓心間距比和約化速度下的運動軌跡。其中，圖7(a1)、圖7(a2)、圖7(b1)和圖7(c1)中右下角分別給出了放大后的軌跡曲線。從圖7(a1)-(a4)中可以看到，在相同的圓心間距比Lx/D=3.0下，隨著Ur的增大，下游圓柱的運動軌跡始終是“8”字形的封閉曲線。而對于Lx/D=5.0和8.0，當Ur≤5.0時，下游圓柱的運動軌跡類似于Lx/D=3.0的情況，也為“8”字形。當Ur>5.0后，圓柱的軌跡則變為不封閉的且不規則的曲線，如圖7(b3)-(b4)及7(c3)-(c4)所示。這說明圓心間距比和約化速度對下游圓柱的運動有明顯的影響。當Lx/D比較小時，下游圓柱在兩自由度方向上運動的相關性并不隨約化速度有明顯的變化。而當Lx/D較大時，Ur> 5.0后，圓柱在x和y方向運動的相互影響作用則由簡單變得復雜。

圖6 不同圓心間距比Lx /D及約化速度Ur下的尾渦脫落模式(Re = 150，m*=10.0，ξ=0.007)

圖7 不同圓心間距比Lx / D及約化速度Ur下下游圓柱的運動軌跡(Re=150，m* = 10.0，ξ = 0.007)

3.2 質量比對下游圓柱渦激振動的影響

質量比是圓柱結構渦激振動響應的重要影響因素之一，但對于多圓柱問題，相對于前述的間隙比和約化速度，人們對質量比的影響作用還缺乏細致的研究和足夠的了解。在此，本文將以Lx/D=5.0為例，考慮三種不同的質量比m*= 5.0、10.0和20.0，對這一問題開展進一步的數值計算，分析質量比對圓柱的位移響應、振動頻率、鎖定區間、受力特性、尾渦脫落模式及運動軌跡的影響。

3.2.1 下游圓柱的位移及振動頻率

圖8首先給出了下游圓柱橫流向相對振動頻率f/fn隨約化速度Ur的變化情況。從圖中可以看出，在本文所考慮的參數范圍內，質量比對鎖定區間的影響作用不大，三種質量比下，鎖定區間所對應的約化速度范圍均為6.0≤Ur≤12.0。

圖9給出了不同質量比下，下游圓柱的振動響應隨約化速度的變化情況。圖9(a)的結果表明，在鎖定區間內，下游圓柱的橫流向位移均方根總體上隨約化速度的增大而逐漸減小。同時，在相同Ur下，m*越小，YRMS/D越大。而在未發生鎖定的低約化速度范圍內，在m*=20.0的條件下，YRMS/D隨著Ur的增大持續增大，直到發生鎖定時達到極值；對于相對較小的質量比m*=5.0和10.0，二者的YRMS/D在發生鎖定前均出現跳躍。進一步的流場分析結果表明，m*=5.0，Ur=4.5和m*=10.0，Ur=5.5兩個跳躍點均由尾渦脫落模式的改變所致，這說明質量比的變化也會對渦激振動過程中的流動模式和位移響應產生影響作用。

圖9(b)考察了下游圓柱順流向平均位移XM/D隨Ur的變化情況。從圖中可以看出，順流向平均位移總體上均隨著Ur的增大呈上升趨勢，質量比越大，圓柱向下游偏移的平均位移越小。對于m*=5.0和10.0，XM/D也分別在Ur=4.5和5.5處發生跳躍，這與圖9(a)的情況類似。

值得說明的是，在m*=5.0和20.0的情況下，相對于m*=10.0，下游圓柱在x,y兩個方向的聯合位移形態沒有本質性的改變，與圖7(b1)-(b4)所給出的運動形式在總體上基本一致。這表明，在本文所考查的參數范圍內，質量比對下游圓柱縱向-橫向運動的相干性的影響作用有限。

圖8 不同質量比m*下，下游圓柱橫流向相對振動頻率f / fn隨約化速度Ur的變化關系(Re = 150，Lx / D = 5.0，ξ = 0.007)

圖9 不同質量比m*下，下游圓柱渦激振動響應特性隨約化速度Ur的變化關系(Re =150，Lx / D = 5.0，ξ = 0.007)

3.2.2 下游圓柱的受力特性

圖10 不同質量比m*下，下游圓柱受力特性隨約化速度Ur的變化關系(Re=150，Lx/D=5.0，ξ=0.007)

4 結 論

本文利用迎風有限元數值方法對不可壓縮粘性流體的Navier-Stokes方程進行了求解，并結合ALE動網格方法，對Re=150條件下，串聯雙圓柱的渦激振動問題開展了二維數值分析。其中，上游圓柱固定不動，下游圓柱在彈簧和阻尼約束下允許同時發生順流向和橫流向的運動。在對數值模型的可靠性進行了驗證的基礎上，重點研究了兩等直徑圓柱圓心間距、質量比以及約化速度對下游圓柱渦激振動響應的影響作用，包括位移響應、振動頻率特征、鎖定區間、流體作用力、尾渦模式以及運動軌跡等方面。數值結果表明，在本文的計算范圍內，隨著兩圓柱間相對圓心間距Lx/D的增大，下游圓柱開始發生鎖定時對應的約化速度增大，在較大的間距比下(Lx/D=5.0和8.0)，鎖定區間較單圓柱及小間距情況(Lx/D=3.0)明顯變寬。在Lx/D=3.0的小間距下，下游圓柱的橫流向振動位移可達到單圓柱情況下的1.5倍，同時下游圓柱的順流平均位移出現負值，即下游圓柱的平衡位置向上游圓柱移動。數值結果也表明，隨著約化速度和間距比的變化，串聯圓柱的尾渦模式會發生顯著改變，進而深刻影響圓柱的受力和運動響應。本文還以Lx/D= 5.0為典型間距，研究了質量比對下游圓柱渦激振動響應和受力特性的影響。結果表明，質量比對鎖定區間的影響不明顯，但當鎖定發生時，在相同約化速度下，隨著質量比的增大，下游圓柱的橫流向運動位移響應會有所減小。進一步的流場分析表明，質量比的變化也會導致尾渦脫落模式的改變，進而影響到圓柱的受力和運動特性，這在進入鎖定區間前非常明顯。

目前關于兩自由度渦激振動問題的研究，通常固定雷諾數，考慮運動響應等對約化速度的依賴關系。但對實際問題而言，往往是圓柱結構形式和流體物性一定，來流流速發生改變，此時雷諾數和約化速度同時發生變化。因此，固定雷諾數而僅改變約化速度的研究方法不能滿足全面了解渦激振動特性的需要，這是在以后的研究工作中需要加以關注的問題。另一方面，隨著雷諾數的不斷增加，二維模型將不再適用，需要開展三維數值計算。本文的計算結果也表明，即便是在簡單的二維低雷諾數條件下，約化速度的微小變化也會引起渦激振動響應特性的劇烈改變，這就要求采用很小的約化速度增量來開展細致的數值分析工作，這些都會導致計算量的激增，對現有的計算條件提出了較大的挑戰。此外，在高雷諾數條件下，直接數值模擬(Direct Numerical Simulation)在理論上是可靠的，但計算量巨大。采用湍流模型雖可緩解計算量的壓力，但湍流模式本身的適用性問題至今仍是全人類共同面對科學難題之一。

參 考 文 獻

[1]Khalak A，Williamson C H K.Motions,forces and mode transitions in vortex-induced vibrations at low mass-damping [J].Journal of Fluids and Structures，1999，13(7-8)：813-851.

[2]Feng C C.The measurement of vortex induced effects in flow past stationary and oscillating circular and D-section cylinders [D].Vancouver：University of British Columbia，1968.

[3]Williamson C H K，Roshko A.Vortex formation in the wake of an oscillating cylinder [J].Journal of Fluids and Structures，1988，2(4)：355-381.

[4]何長江，段忠東.二維圓柱渦激振動的數值模擬 [J].海洋工程，2008，26(1)：57-63.

HE Chang-jiang，DUAN Zhong-dong.Numerical simulation of vortex-induced vibration on 2D circular cylinders [J].The Ocean Engineering，2008，26(1)：57-63.

[5]方平治，顧 明.圓柱兩自由度渦激振動的數值模擬研究 [J].同濟大學學報(自然科學版)，2008，36(3)：295-298.

FANG Ping-zhi，GU Ming.Numerical simulation for vortex-induced vibration of circular cylinder with two degree of freedoms [J].Journal of Tongji University (Nature Science)，2008，36(3)：295-298.

[6]唐國強.立管渦激振動數值模擬方法及物理模型實驗 [D].大連：大連理工大學，2011.

[7]Brika D，Laneville A.The flow interaction between a stationary cylinder and a downstream flexible cylinder [J].Journal of Fluids and Structures，1999，13(5)：579-606.

[8]Carmo B S，Sherwin S J，Bearman P W，et al.Flow-induced vibration of a circular cylinder subjected to wake interference at low Reynolds number [J].Journal of Fluids and Structures，2011，27(4)：503-522.

[9]Bao Y，Zhou D，Tu J H.Flow interference between a stationary cylinder and an elastically mounted cylinder arranged in proximity [J].Journal of Fluids and Structures.2011，27(8)：1425-1446.

[10]Vandiver J K，Jong J Y.The relationship between in-line and cross-flow vortex-induced vibration of cylinders [J].Journal of Fluids and structures，1987，1(4)：381-399.

[11]陳文曲，任安祿，李廣望.串聯雙圓柱繞流下游圓柱兩自由度渦致振動研究 [J].力學學報，2004，36(6)：732-738.

CHEN Wen-qu，REN An-lu，LI Guang-wang.The Numerical study of two-degree-of -freedom vortex-induced vibration of the downstream cylinder in tandem arrangement [J].Chinese Journal of Theoretical and Applied Mechanics，2004，36(6)：732-738.

[12]Jiang C B，Kawahara M.The analysis of unsteady incompressible flows by a three-step finite element method [J].International Journal for Numerical Methods in Fluids，1993，16(9)：793-811.

[13]Lu L，Qin J M，Teng B，et al.Numerical investigations of life suppression by feedback rotary oscillation of circular cylinder at low Reynolds number [J].Physics of Fluids，2011，23(3)：033601.

[14]Lu L，Cheng L， Teng B，et al.Numerical investigation of fluid resonance in two narrow gaps of three identical rectangular structures [J].Applied Ocean Research，2010，32(2)：177-190.

[15]呂 林.海洋工程中小尺度物體的相關水動力數值計算 [D].大連：大連理工大學，2006.

[16]Zdravkovich M M.The effects of interference between circular cylinders in cross flow [J].Journal of Fluids and Structures，1987，1(2)：239-261.

[17]Carmo B S，Meneghini J R，Sherwin S J.Secondary instabilities in the flow around two circular cylinders in tandem [J].Journal of Fluid Mechanics，2010，644：395-431.

[18]Leontini J S，Stewart B E，Thompson M C，et al.Predicting vortex-induced vibration from driven oscillation results [J].Applied Mathematical Modelling，2006，30(10)：1096-1102.