吳 浩， 孫大鵬， 呂 林

(大連理工大學(xué) 海岸和近海工程國家重點實驗室， 遼寧 大連 116023)

低雷諾數(shù)下多根圓柱橫向受迫振動的數(shù)值模擬研究

吳 浩， 孫大鵬， 呂 林

(大連理工大學(xué) 海岸和近海工程國家重點實驗室， 遼寧 大連 116023)

采用三步有限元方法直接求解ALE描述的N-S方程組，結(jié)合ALE動網(wǎng)格技術(shù)，建立了2D橫向受迫振動的數(shù)值模型。利用該模型計算了1個主圓柱體與4個附屬小圓柱體的橫向受迫振動問題，Re=185，振幅比Ae/D=0.2，頻率比fe/fs=0.86～1.04，間隔0.01。通過分析總流體力及其頻譜隨頻率比的變化，發(fā)現(xiàn)總流體力出現(xiàn)了突變，且流體力的振動模態(tài)可以分為4種不同的類型，其中包括基本鎖定狀態(tài)。最后分析了升力突變前后流體-結(jié)構(gòu)能量傳遞關(guān)系、結(jié)構(gòu)位移和升力的相位關(guān)系以及渦脫落模式的變化。

多根圓柱；橫向受迫振動；三步有限元方法

0 引言

海洋立管的渦激振動問題是近年來學(xué)術(shù)界高度關(guān)注的熱點，立管渦激振動的簡化模型即流體力學(xué)中經(jīng)典的圓柱繞流問題。圓柱體在穩(wěn)定流場中的運動可分為自激振動和受迫振動兩種類型，圍繞著這兩種運動類型，學(xué)術(shù)界開展了廣泛而深入的研究，對研究成果比較全面的回顧可見Bearman[1]、Williamson[2]和Sarpkaya[3]的綜述文獻。自激振動的研究表明，隨著縮減速度的變化，圓柱體的振幅、振動頻率和渦脫落模態(tài)都迅速發(fā)生變化，因此，要確定這些影響因素之間的相互聯(lián)系是非常困難的。而在受迫振動的研究中，圓柱體的振幅、振動頻率可以精確的控制。

對于在流場中發(fā)生受迫振動的圓柱體，若振動頻率fe接近圓柱體固定不動時的渦脫落頻率fs，則圓柱體所受流體力和渦脫落模態(tài)都會發(fā)生顯著變化，這是受迫振動最主要的水動力特性。Gu等[4]通過粒徑追蹤和粒子圖像測速方法將受迫振動圓柱的近尾跡可視化，并獲得了流線圖和速度場分布。對于給定的雷諾數(shù)和振幅，隨著fe/fs的增加，從圓柱一側(cè)最初形成的尾渦逐漸靠近圓柱底部，最終達到極限位置，然后尾渦的生成位置和符號都發(fā)生突變，這種現(xiàn)象在低雷諾數(shù)和高雷諾數(shù)時都會發(fā)生。其后，Lu和Dalton[5]、Guilmineau和Queutey[6]分別采用不同的數(shù)值計算方法對Gu的實驗結(jié)果進行了驗證，給出了圓柱體上的流體力系數(shù)。Blackburn和Henderson[7]通過數(shù)值模擬方法研究了單圓柱受迫振動問題，在計算中固定雷諾數(shù)和振幅，關(guān)注fe/fs=1附近振動頻率的變化對水動力特性和尾跡結(jié)構(gòu)的影響，通過分析流體力系數(shù)和流體-結(jié)構(gòu)能量傳遞關(guān)系，將受迫振動時的流動狀態(tài)劃分為4種類型：周期狀態(tài)(Periodic)、準周期狀態(tài)(Quasi-periodic)、弱混沌狀態(tài)(Weaklychaotic)和混沌狀態(tài)(Chaotic)，并發(fā)現(xiàn)在周期狀態(tài)下存在不對稱尾跡結(jié)構(gòu)和非零的平均升力系數(shù)。

國內(nèi)方面，潘志遠[8]總結(jié)了國內(nèi)外圓柱體自激振動和受迫振動的主要研究成果，并討論了兩種運動類型的區(qū)別和聯(lián)系。樊娟娟等[9]在N-S方程和k-ω湍流模型的基礎(chǔ)上，結(jié)合Fluent動網(wǎng)格技術(shù)，研究了較高雷諾數(shù)下大振幅比圓柱橫向受迫振動問題。趙靜等[10]在N-S方程和k-ω湍流模型的基礎(chǔ)上，利用流線迎風(fēng)有限元方法結(jié)合ALE動網(wǎng)格技術(shù)，研究了亞臨界雷諾數(shù)下圓柱體的受迫振動問題，分析了流體力與振動頻率的聯(lián)系并給出了鎖定區(qū)間。

實際海洋工程中的立管周圍都安裝了多根附屬管線，如節(jié)流壓井管線、BOP液控管線和水龍帶等。從局部看，當發(fā)生渦激振動時，主管與附屬管線作為一個整體進行運動，附屬管線將對主管周圍流場的水動力特性造成明顯影響?；谶@一工程背景，Wu等[11]開展水池拖曳實驗，研究了長細比為1 750的柔性立管上安裝4根附屬管線的渦激振動問題，通過改變間距和附屬管線沿管長的覆蓋率，詳細分析了附屬管線對立管位移振幅、主控頻率和振動模態(tài)等結(jié)構(gòu)響應(yīng)特性的影響。利用三步有限元方法直接求解ALE描述的N-S方程組，結(jié)合ALE動網(wǎng)格技術(shù)，建立了2D橫向受迫振動數(shù)值模型。首先求解經(jīng)典的單圓柱橫向受迫振動問題，對數(shù)值模型進行驗證，然后進一步將該數(shù)值模型應(yīng)用于計算1個主圓柱體與4個附屬小圓柱體的橫向受迫振動問題，在計算中給定雷諾數(shù)和振幅，關(guān)注fe/fs對圓柱體的總流體力系數(shù)和渦脫落模態(tài)的影響。

1 數(shù)值模型

對于密度為常數(shù)、忽略質(zhì)量力、層流態(tài)的二維不可壓縮N-S方程組簡化形式，選取圓柱體直徑D(對于多根圓柱體，指主圓柱體直徑)、來流流速U和時間尺度D/U作為無量綱化尺度參數(shù)，并采用ALE方法處理圓柱體運動邊界，可以得到ALE描述的無量綱形式N-S方程組：

式中：u′=u/U、x′=x/D、t′=tU/D、p′=p/ρU2、Re=UD/υ，υ為運動學(xué)粘性系數(shù)；u′m=um/U為無量綱ALE網(wǎng)格運動速度分量；i,j=1, 2為坐標分量。為方便計算，將U和D均取為常數(shù)1。

采用高階時間格式和Galerkin有限元空間離散結(jié)合的三步有限元方法求解ALE描述的無量綱形式N-S方程組，該方法的具體離散求解過程、網(wǎng)格依賴性和計算穩(wěn)定性分析可參考呂林[12]和唐國強[13]的工作。

在受迫振動的計算中，將圓柱體周圍的有限范圍指定為ALE網(wǎng)格更新區(qū)域，在該區(qū)域內(nèi)將網(wǎng)格看作彈性模量為常數(shù)的彈性體，這樣由于邊界條件改變而引起的網(wǎng)格坐標變化可以用Laplace方程表示。在每一時間步，預(yù)先給定結(jié)構(gòu)的運動位置，即網(wǎng)格的運動邊界條件，然后在ALE區(qū)域內(nèi)求解Laplace方程自動生成新的網(wǎng)格坐標，從而得到網(wǎng)格運動速度u′m，ALE區(qū)域外的網(wǎng)格則不用更新，有效地減小了計算量。測試表明，通過選擇合適的ALE區(qū)域大小，網(wǎng)格在運動中不會發(fā)生畸變，且保持了良好的正交性。

2 無量綱參數(shù)

2.1 流體力系數(shù)

通過求解N-S方程得到流場數(shù)值解，即可根據(jù)圓柱表面的壓力分布以及速度梯度計算圓柱受力，圓柱單位長度上的正向力Fx和橫向力Fy分別為：

2.2 流體-結(jié)構(gòu)能量傳遞系數(shù)

結(jié)構(gòu)在發(fā)生橫向受迫振動時的瞬時位移為y(t)，其無量綱形式為y′(t)=y(t)/D，在每一個振動周期內(nèi)，流體向結(jié)構(gòu)傳遞的機械能可以寫成如下的無量綱形式：

式中：T為結(jié)構(gòu)受迫振動的周期。

E被稱為流體-結(jié)構(gòu)能量傳遞系數(shù)，E為正值則流體對結(jié)構(gòu)做功，E為負值則結(jié)構(gòu)對流體做功。結(jié)構(gòu)位移和升力的基本諧波之間的相位角φ也常用于衡量流體-結(jié)構(gòu)能量傳遞，當E為正值時，0°<φ<180°。但是相位角僅僅反映了流體力的基本諧振動與結(jié)構(gòu)之間的能量傳遞信息，當流場處于非周期狀態(tài)時，升力的基本諧波可能是不存在的，因此在結(jié)果分析中采用能量傳遞系數(shù)更加準確直觀。

Blackburn和Henderson[7]指出，單根圓柱體受迫振動時，尾渦生成和圓柱體運動之間的相位關(guān)系變化伴隨著E符號的變化。在單根圓柱體受迫振動的基本鎖定區(qū)，尾跡處于周期性的渦脫落狀態(tài)，升力為穩(wěn)定的周期振動，此時結(jié)構(gòu)位移和升力的相平面圖為一閉合環(huán)路，且E的符號與閉合環(huán)路的方向有關(guān)，對于正的能量傳遞，環(huán)路方向為順時針，反之對于負的能量傳遞，環(huán)路方向為逆時針。

3 計算結(jié)果

3.1 數(shù)值模型的驗證

圖1 單圓柱橫向受迫振動的流體力系數(shù)，Re=185，Ae/D=0.2

圖2 多根圓柱體計算模型

CDCLrmsSt該文1．3180．4470．193Lu和Dalton[5]1．3100．4220．195Guilmineau和Queutey[6]1．2870．4430．195

3.2 多根圓柱體橫向受迫振動

多根圓柱體數(shù)值計算模型如圖2所示，主圓柱的無量綱直徑D=1，以主圓柱中心為原點，建立笛卡爾直角坐標系，x軸方向為順流方向(IL)，y軸方向為橫流方向(CF)。分別對應(yīng)IL和CF方向，對稱放置4個附屬小圓柱，其無量綱直徑d=0.25，主圓柱外表面到附屬小圓柱外表面的無量綱距離為0.375。

如圖3所示，計算域范圍為60D×50D，主圓柱中心距離上、下、左邊界25D，距離右邊界35D，在主圓柱周圍8D×8D范圍內(nèi)為ALE網(wǎng)格運動區(qū)域。整個計算域采用結(jié)構(gòu)化網(wǎng)格進行剖分，總單元數(shù)為37 470，總節(jié)點數(shù)為37 840，主圓柱表面網(wǎng)格尺寸約為0.02D，附屬小圓柱表面網(wǎng)格尺寸約為0.01D，圓柱體附近的計算網(wǎng)格剖分如圖4所示。計算時間步長Δt′為無量綱時間0.001。

圖3 計算域示意圖 圖4 圓柱體附近的網(wǎng)格剖分

圖5 多根固定圓柱體繞流的總力系數(shù)時程曲線和FFT頻譜，Re=185,Ae=0.2

圖6 多圓柱橫向受迫振動的總流體力系數(shù)，Re=185，Ae/D=0.2

在多根圓柱體橫向受迫振動的計算中，主圓柱和4個附屬小圓柱保持相對位置不變，整體沿橫流方向做簡諧運動，運動方程為y=Aesin(2πfet)。在計算中限定Re=185，Ae/D=0.2，則流體力和渦量場是唯一變量頻率比fe/fs的函數(shù)，因此在fe/fs=0.86～1.04范圍內(nèi)，將fe/fs的計算間隔精細到0.01，從而能夠準確地捕獲流體力和渦脫落模式發(fā)生轉(zhuǎn)變的時機。

選取總流體力系數(shù)CD和CL的振動穩(wěn)定段進行傅里葉變換，則得到了總流體力系數(shù)頻譜在頻率比fe/fs上的分布情況，如圖7所示。圖7中f′ -fe/fs平面的虛線從左至右依次為f′=1feD/U,2feD/U,3feD/U,4feD/U，這意味著在這些虛線上，無量綱頻率f′依次為結(jié)構(gòu)無量綱振動頻率feD/U的1～4倍。從圖7中可以看出，在各頻率比上，CD的主控?zé)o量綱頻率為2feD/U，CL的主控?zé)o量綱頻率為feD/U。這與單根圓柱體橫向受迫振動是相同的，說明了多根圓柱體在受迫振動時，總流體力仍然受結(jié)構(gòu)運動控制。與CDrms和CLrms相對應(yīng)，CD主控頻率的幅值在計算范圍內(nèi)隨fe/fs的增加而增加。CL主控頻率的幅值在fe/fs<0.94時隨fe/fs的增加而減小，在fe/fs>0.94時隨fe/fs的增加而增加。

CL頻譜的頻率組成在以下4個頻率比區(qū)間上具有明顯不同的特征：在fe/fs=0.86～0.89區(qū)間，除主控?zé)o量綱頻率feD/U以外，在3倍頻上也有單峰出現(xiàn)，但幅值非常小，同時在超低頻段和1、2倍頻附近都有寬頻帶小幅值分布，可見該區(qū)間的CL振動以feD/U的奇數(shù)倍頻諧振為主，并具有較弱的隨機性；在fe/fs=0.90～0.93區(qū)間，僅feD/U的1、3倍頻上有單峰出現(xiàn)，且3倍頻的幅值非常小，但不存在寬頻帶小幅值分布，可見該區(qū)間的CL振動是feD/U的周期性奇數(shù)倍頻諧振；在fe/fs=0.94～0.99區(qū)間，與前兩個區(qū)間相比，頻率組成發(fā)生了明顯的變化，主控?zé)o量綱頻率feD/U的幅值突然減小，同時在多個非倍頻位置也出現(xiàn)峰值，并且在整個頻段上都有寬頻帶小幅值分布，可見該區(qū)間的CL振動具有明顯的隨機性；在fe/fs=1.00～1.04區(qū)間，feD/U的1、3倍頻的幅值顯著增加，而非倍頻位置的幅值變化不大，整個頻段的寬頻特性突然消失了，可見該區(qū)間的CL振動是feD/U的奇數(shù)倍頻諧振和非倍頻諧振結(jié)合的形式，并具有一定的周期性。通過對CL頻譜的分析可知，從低到高4個頻率比區(qū)間上的CL振動分別對應(yīng)Blackburn和Henderson[7]描述的弱混沌狀態(tài)、周期狀態(tài)、混沌狀態(tài)和準周期狀態(tài)，將這4個區(qū)依次標記為I、II、III和IV區(qū)。CD頻譜隨fe/fs的變化與CL頻譜很相似，也可以根據(jù)頻率組成劃分為4個頻率比區(qū)間，不同的是CD頻譜除主控?zé)o量綱頻率2feD/U以外，在4倍頻上有單峰出現(xiàn)且幅值非常小，即CD振動中存在偶數(shù)倍頻諧振。

圖7 多圓柱橫向受迫振動的總流體力系數(shù)頻譜(虛線從左至右依次為f′=1, 2, 3, 4 feD/U)

圖9給出了流體力的振動穩(wěn)定段上能量傳遞系數(shù)的均方根Erms隨fe/fs的變化。從圖9中可以看出，只有在II區(qū)(周期狀態(tài))Erms接近于0，說明升力處于穩(wěn)定的周期振動狀態(tài)。III區(qū)(混沌狀態(tài))的Erms是最大的，說明升力振動的隨機性很強。I區(qū)(弱混沌狀態(tài))和IV區(qū)(準周期狀態(tài))的Erms則介于II和III區(qū)之間。

圖8 能量傳遞系數(shù)平均值隨fe/fs的變化 圖9 能量傳遞系數(shù)均方根隨fe/fs的變化

圖10給出了流體-結(jié)構(gòu)能量傳遞關(guān)系發(fā)生突變前后，兩個頻率比fe/fs=0.93和0.94上的總升力系數(shù)CL時間過程線。從圖10中可以看出，當fe/fs=0.93時，升力是穩(wěn)定的周期振動，且如前所述，因為存在較弱的奇數(shù)階諧振動，在每半個周期上時程線并不關(guān)于極值左右對稱。頻率比fe/fs增加了0.01，升力就變?yōu)殡S機振動。相應(yīng)的，在這兩個頻率比上能量傳遞系數(shù)E時間過程線以及結(jié)構(gòu)位移和升力的相平面圖也具有完全不同的特征。如圖11、圖12所示，當fe/fs=0.93時，E為負值且?guī)缀鹾愣ú蛔儯嗥矫鎴D為一閉合環(huán)路，環(huán)路方向為逆時針與E的符號對應(yīng)，而fe/fs=0.94時，E絕大多數(shù)時間為正值且波動很明顯，相平面圖不是閉合環(huán)路，其方向大體上為順時針仍然與E的符號對應(yīng)。因此，升力系數(shù)變化趨勢和振動模態(tài)的突變與結(jié)構(gòu)位移和升力的能量傳遞關(guān)系以及相位關(guān)系的突變有關(guān)。

圖10 總升力系數(shù)時間過程線，fe/fs = 0.93, 0.94

圖11 能量傳遞系數(shù)時間過程線 圖12 橫向位移和升力的相平面圖

圖13 一個運動周期內(nèi)的瞬時渦量場，fe/fs = 0.93, 0.94 (實線：渦量為正；虛線：渦量為負)

圖13給出了fe/fs=0.93和0.94時，一個運動周期內(nèi)的無量綱位移y′和總升力系數(shù)CL時間過程線，并給出四個代表性時刻：y′=0(向上運動)、y′=0.2、y′=0(向下運動)和y′=-0.2的瞬時渦量場。從圖13中可以看出，在流體-結(jié)構(gòu)能量傳遞關(guān)系發(fā)生突變前后，兩種頻率比上的渦脫落模式差別很大。當fe/fs=0.93，即發(fā)生基本鎖定時(圖13(a)～圖13(d))，渦脫落模式與單圓柱體卡門渦街非常相似，但是在近尾流處的渦結(jié)構(gòu)更加復(fù)雜。(a)與(c)，(b)與(d)的渦量場近似于上下對稱結(jié)構(gòu)，由主圓柱和各附屬小圓柱上形成的正負渦相互干涉，在近尾流區(qū)的上下兩側(cè)交替形成一個由3個正(負)渦圍繞1個負(正)渦組成的渦團，然后該渦團逐漸變形融合，形成一個由多個正(負)渦組成的不規(guī)則渦團并脫落。在一個運動周期內(nèi)，多根圓柱體的上下兩側(cè)對稱的各脫落一個渦團，從而造成總升力的周期性振動，且振動主控頻率等于結(jié)構(gòu)振動頻率。fe/fs=0.94時(圖13(e)～圖13(h))，渦脫落模式變得極不穩(wěn)定，近尾流區(qū)上下兩側(cè)的渦脫落形態(tài)完全是不對稱的，正負渦的干涉非常嚴重，渦團的結(jié)構(gòu)和脫落頻率在一個周期內(nèi)不斷地隨機變化，無明顯規(guī)律可循。結(jié)合之前的分析，fe/fs由0.93增加到0.94時，渦脫落模式由周期性的穩(wěn)定脫落變?yōu)殡S機的不穩(wěn)定脫落，能量傳遞關(guān)系由結(jié)構(gòu)向流體傳遞變?yōu)榱黧w向結(jié)構(gòu)傳遞，正是這些水動力特性的突變造成了總流體力的突變。

4 結(jié) 論

[ 1 ] Bearman P W. Vortex shedding from oscillating bluff bodies [J]. Annual Review of Fluid Mechanics, 1984, 16:195-222.

[ 2 ] Williamson C H K, Govardhan R. Vortex-induced vibrations [J]. Annual Review of Fluid Mechanics, 2004. 36:413-55.

[ 3 ] Sarpkaya T. A critical review of the intrinsic nature of vortex-induced vibrations [J]. Journal of Fluids and Structures, 2004, 19:389-447.

[ 4 ] Gu W, Chyu C, Rockwell D. Timing of vortex formation from an oscillating cylinder [J]. Physics of Fluids, 1994, 6:3677-3682.

[ 5 ] Lu X Y, Dalton C. Calculation of the timing of vortex formation from an oscillating cylinder [J]. Journal of Fluids and Structures, 1996, 10:527-541.

[ 6 ] Guilmineau E, Queutey P. A numerical simulation of vortex shedding from an oscillating circular cylinder [J]. Journal of Fluids and Structures, 2002, 16(6):773-794.

[ 7 ] Blackburn H M, Henderson R D. A study of two-dimensional flow past an oscillating cylinder [J]. Journal of Fluid Mechanics, 1999, 38(5):255-286.

[ 8 ] 潘志遠. 海洋立管渦激振動機理與預(yù)報方法研究 [D]. 上海: 上海交通大學(xué), 2006.

[ 9 ] 樊娟娟, 唐友剛, 張若瑜, 等. 高雷諾數(shù)下圓柱繞流與大振幅比受迫振動的數(shù)值模擬 [J]. 水動力學(xué)研究與進展, 2012, 27(1): 24-32.

[10] 趙靜, 呂林, 董國海, 等. 亞臨界雷諾數(shù)下圓柱受迫振動的數(shù)值研究[J]. 計算力學(xué)學(xué)報, 2012, 29(1): 74-79.

[11] Wu H, Sun D P, Lu L, et al. Experimental investigation on the suppression of vortex-induced vibration of long flexible riser by multiple control rods [J]. Journal of Fluids and Structures, 2012, 30:115-132.

[12] 呂林. 海洋工程中小尺度物體的相關(guān)水動力數(shù)值計算 [D]. 大連: 大連理工大學(xué), 2006.

[13] 唐國強. 立管渦激振動數(shù)值模擬方法及物理模型實驗 [D]. 大連: 大連理工大學(xué), 2011.

Numerical Simulation of Lateral Forced Vibration of Multiple Circular Cylinders at Low Reynolds Number

(State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Liaoning Dalian 116023, China)

Three steps finite element method is used to directly solve the N-S equations described by ALE method, and combining with ALE dynamic grid technique, a numerical model for 2D lateral forced vibration is developed. The lateral forced vibration of one main circular cylinder with additional four smaller circular cylinders is calculated, at Re=185, amplitude ratioAe/D=0.2andfrequencyratiofe/fs=0.86～1.04withintervalof0.01.Byanalyzingthetotalfluidforcesandtheirspectrumsasafunctionoffe/fs,anabruptchangeontotalliftforceisfound,andtheoscillatingmodeoffluidforcecanbedividedintofourtypesincludingtheprimarylock-inregime.Thenthefluid-structureenergytransfer,thephasebetweenstructuralmovementandliftforceandthevortexsheddingmodeareanalyzedtorealizethenatureonaccountforabruptchangeoffluidforce.

multiple circular cylinders; forced vibration; three steps finite element method

2014-07-01

國家自然科學(xué)基金項目(項目編號：51409041，51279027)。

吳 浩(1980-)，男，工程師。

1001-4500(2015)02-0048-09

P

A