不同間隙比和邊界層厚度下的近壁面彈性圓柱渦激振動特性

劉旭菲， 花 陽， 衛昱含， 及春寧

(1. 浙江水利水電學院 水利與環境工程學院， 杭州 310018； 2. 上海勘測設計研究院有限公司， 上海 200335；3. 天津大學 水利工程仿真與安全國家重點實驗室， 天津 300072)

海上油氣的集疏運主要通過管道運輸實現，其破壞后，修復難，經濟成本和生態代價高[1]。造成其破壞的原因為懸空管道在近底水流作用下的渦激振動。當水流流經管道時，從管道的上下兩側交替泄放旋渦，產生周期性或擬周期性的升、阻力，導致管道發生彎曲振動[2-4]。在一定的流速范圍內，結構的振動將對旋渦脫落產生控制作用，使得旋渦的脫落頻率偏離Strauhal頻率，與管道振動的固有頻率相接近，發生“鎖定”，導致管道大幅振動，進而誘發疲勞破壞[5-6]。因此，對于近壁面海底管道的渦激振動現象及其誘發機理進行研究，具有十分重要的工程價值和理論意義。

近壁面圓柱由于受壁面和壁面邊界層流動的影響，其渦激振動響應與孤立圓柱明顯不同。Taneda[7]開展了試驗研究，探討了低雷諾數(Re=140)條件下，圓柱與壁面之間的間隙G對圓柱尾渦的影響。此后，Bearman等[8-9]試驗研究了高雷諾數(Re=2×104～2×105)內近壁圓柱的振動響應，當間隙比G/D小于臨界值0.3時，振動響應被明顯抑制；大于臨界值后，振動響應逐漸恢復。然而，Jacobsen等[10]的試驗研究結果表明，小間隙比工況下，即便圓柱泄渦被抑制，但圓柱仍發生微小幅值的振動。同樣，Wang等[11-12]發現即使間隙比很小，圓柱下部泄渦被抑制，但上部仍泄放單列渦街，圓柱發生微幅振動。Tsahalis[13]的試驗研究發現：在近壁影響下，順流向振幅接近橫流向振幅，圓柱的振動軌跡為橢圓形，這與無壁面圓柱渦激振動的“8”字形軌跡明顯不同。Fu等[14]發現在近壁面的影響下，圓柱渦激振動的振動頻率相較孤立圓柱工況有所增加，與速度同相位的升力值出現增大。Yang等[15-17]通過系列試驗發現，隨著間隙比的減小，渦激振動鎖定區間的寬度不斷增大，圓柱橫流向振幅也逐漸增加。Nielsen等[18]考察了圓柱初始下垂對渦激振動的影響，并探討了大長細比圓柱渦激振動的多模態特征。Li等[19-20]試驗研究了水流和波浪作用下彈性圓柱的渦激振動，同樣發現了高折合流速條件下的二階振動模態。婁敏等[21-22]進一步分析了管道內流流速、管內壓強等參數對渦激振動特性的影響。

在數值模擬研究方面，Zhao等[23]開展了間隙比為0.002和0.3條件下近壁面圓柱的二維數值模擬研究，重點討論了壁面碰撞反彈對振動響應的影響。Jiang等[24]發現近壁圓柱尾渦的穩定性隨間隙比減小而增強。Chern等[25-26]分別從壁面剪切力和圓柱升阻力等方面對近壁面圓柱渦激振動進行了數值模擬研究。

綜合以上研究可以發現，由于近壁圓柱渦激振動問題相對復雜，與無壁面工況相比，涉及間隙比和邊界層厚度的影響。另一方面，已有的近壁圓柱渦激振動數值模擬研究多針對彈性支撐的二維剛性圓柱，近壁彈性圓柱渦激振動的三維數值模擬研究開展得較少。本文對不同間隙比和邊界層厚度下的近壁面彈性圓柱渦激振動開展了三維數值模擬研究，以期對以下兩個問題作出解答。① 不同間隙比和邊界層厚度下流場旋渦結構有何變化。② 圓柱振動振幅、振動頻率、振動同步性與振動軌跡隨圓柱和邊界層相對位置改變發生怎樣的變化。

1 數值模擬方法

流固耦合的數值模擬方法采用浸入邊界法(immersed boundary method, IBM)。在IBM方法框架下，流體域通過固定的正交網格進行空間離散，而可移動的固體采用離散的拉格朗日點(浸入邊界點)表示。流固之間的耦合作用通過在動量守恒方程中引入額外的附加體積力項來實現。由于通常情況下正交網格與可移動物面并不重合，故需要應用插值函數進行流固邊界上的信息交換。基于IBM方法的流固耦合控制方程表示如下

(1)

?·u=0

(2)

式中：u為速度；t為時間；p為壓強；ρ為密度；ν為運動黏滯系數；?為梯度算子；f為基于浸入邊界法的附加體積力矢量，表示流固耦合邊界條件，具體可參考文獻[27]。

時間離散后的控制方程守恒形式如下

(3)

?·un+1=0

(4)

式中，h=?·[-uu+ν(?u+?uT)]由對流項與擴散項組成，T為矩陣轉置，附加體積力表示為

(5)

針對傳統浸入邊界法施加邊界條件精度不高的問題，Ji等提出了的嵌入式迭代的浸入邊界法，將浸入邊界法嵌入到壓強泊松方程的迭代求解中，利用壓強的中間解比初始值更接近真實值的特點，迭代修正附加體積力，在不顯著增加計算耗時的前提下，顯著提高了整個算法的求解精度。詳細求解過程參見Ji等的研究。

本文采用三維二結點梁單元模擬管道彎曲振動。每個結點有6個自由度，即3個位移分量和3個轉動分量。將管道進行空間離散，得到矩陣形式有限元控制方程為

(6)

式中：M，C和K分別為質量、阻尼和剛度矩陣；x(t)和F(t)分別為位移和流體力向量。阻尼矩陣采用瑞利阻尼形式，表示為

C=αM+βK

(7)

式中，α和β分別為與質量矩陣和剛度矩陣成比例的阻尼參數。

2 數值模型

2.1 問題描述

本文對近壁面圓柱渦激振動進行了三維直接數值模擬研究，彈性圓柱兩端為鉸接。數值模擬相關參數設置如下：圓柱直徑為D；根據實際工程中海底管道的最大懸跨長度，設置圓柱展向長度L=25D；圓柱質量比m*=m/mf=6，m為圓柱的質量，mf為圓柱排開水的質量；圓柱阻尼系數分別為α=2.4×10-3和β=2.7×10-4，泊松比μ=0.3；間隙比G/D=0.6和1.0；雷諾數Re=U∞D/v=350，U∞為來流流速。考慮到亞臨界雷諾數范圍內(300.0

2.2 模型設置

本文的計算域與邊界條件設置，如圖1所示。來流為均勻來流，壁面邊界層自由發展形成。通過調整圓柱距上游入口的距離，來生成不同厚度的邊界層(δ/D=0，0.7和1.6)。在無邊界層(δ/D=0)和薄邊界層(δ/D=0.7)工況中，計算域大小為51D(X)×40D(Y)×25D(Z)，圓柱距速度入口7D(X1)；在厚邊界層(δ/D=1.6)工況中，計算域大小為80D×40D×25D，圓柱距速度入口36D。圓柱周圍16D(X)×6D(Y)的范圍內采用正交均勻網格，網格尺寸為Δx=Δy=D/32，加密區以外采用正交漸變網格。計算域在展向上采用均勻網格，對于δ/D=0和0.7的工況，計算域沿展向共劃分為256層網格，對于δ/D=1.6的工況，計算域沿展向共劃分為384層網格。

圖1 計算域大小與邊界條件設置Fig.1 The computational domain size and boundary conditions

計算域的速度入口采用狄利克雷型邊界條件，u=U，v=w=0。速度出口采用諾依曼型邊界條件，?u/?x=?v/?x=?w/?x=0。上邊界取自由滑移邊界條件，?u/?y=0，v=0，?w/?y=0。對于無邊界層工況，壁面邊界設置為自由滑移邊界條件，?u/?y=0，v=0，?w/?y=0；對于有邊界層工況，壁面邊界設置為不可滑移邊界條件，u=v=w=0。展向邊界為周期邊界，圓柱表面設置為不可滑移邊界條件。速度出口壓強為零，上、下邊界和速度入口采用諾依曼型邊界條件，即壓強的法向梯度為零。

2.3 網格無關性和計算精度驗證

首先對網格無關性進行驗證。針對G/D=0.6，δ/D=0和Re=350的工況，選取兩種網格，即Δx=Δy=D/32(粗網格)和Δx=Δy=D/48(細網格)，進行數值模擬，結果如表1所示。可以看到兩組數值結果之間的差異可以忽略不計，最大誤差僅為0.62%。故后續模擬采用計算效率更高的粗網格。

表1 網格無關性驗證Tab.1 Grid independence verification

進一步驗證本文數值模擬方法的準確性。對Re=350的孤立圓柱繞流進行三維數值模擬，并于已有結果進行對比，如表2所示。可以看出本文的計算結果與Williamson等[28]的試驗結果和Jiang等[29]的數值結果吻合良好，證明了本文數值模擬方法的準確性。

表2 Re=350條件下孤立圓柱繞流水動力系數對比Tab.2 Comparison of hydrodynamic coefficients for flow around an isolated cylinder at Re = 350

3 結果和討論

為了全面考察圓柱和邊界層相對位置對流場特性和振動特性的影響，本文精心設計了5種工況，分別采用兩種間隙比(G/D=0.6和1.0)和3種邊界層厚度(δ/D=0，0.7和1.6)。本文采用無量綱參數δ/G定量表示圓柱和邊界層相對位置，其值越大，表明圓柱越浸沒入邊界層中。5種工況對應的δ/G，如表3所示。G/D=1.0，δ/D=0和G/D=1.0，δ/D=0.7工況的圓柱完全處于邊界層外；G/D=0.6，δ/D=0.7和G/D=1.0，δ/D=1.6工況的圓柱部分處于邊界層中；G/D=0.6，δ/D=1.6工況的圓柱完全處于邊界層中。

表3 圓柱和邊界層相對位置Tab.3 The relative postion of the cylinder in the boundary layer

3.1 流場特性

瞬時流場和旋渦結構的空間分布，如圖2所示。旋渦結構采用λ2等值面進行可視化。λ2為張量S2+Ω2的第二特征根，其中，S和Ω分別為速度梯度張量?u的對稱部分和不對稱部分。整體來看，流場中尾渦沿圓柱展向的分布表現出三維特性，即圓柱兩端泄渦窄，跨中泄渦寬，展向旋渦成“弓”形分布，這與圓柱兩端振幅小，跨中振幅大有關。這一點尤其在大間隙比和薄邊界層工況下(見圖2(c)和2(d))更為明顯。與之相應，流向渦也表現出類似的特征，并且，隨著邊界層厚度的增大，流向渦的數量減少，強度降低。

(a) G/D=0.6, δ/D=0.7

(b) G/D=0.6, δ/D=1.6

(c) G/D=1.0, δ/D=0

(d) G/D=1.0, δ/D=0.7

(e) G/D=1.0, δ/D=1.6圖2 不同間隙比和邊界層厚度下的流場渦核分布圖Fig.2 The vortex-shedding patterns with different gap ratio and boundary layer thickness

跨中位置(Z/D=12.5)截面的瞬時展向渦量圖，如圖3所示。整體來看，圓柱上下側泄渦不對稱，尾渦出現較明顯的向上偏斜。

(a) G/D=0.6, δ/D=0.7

(b) G/D=0.6, δ/D=1.6

(c) G/D=1.0, δ/D=0

(d) G/D=1.0, δ/D=0.7

(e) G/D=1.0, δ/D=1.6圖3 不同間隙比和邊界層厚度下的流場瞬時渦量圖Fig.3 The contours of instantaneous vorticity with different gap ratio and boundary layer thickness

當G/D=0.6時，壁面邊界層與圓柱上部剪切層的耦合作用明顯，壁面邊界層向上傾斜，直接與從圓柱上部脫落的旋渦合并。而下部旋渦的強度和尺寸均較小，泄渦受到抑制。在薄邊界層工況(δ/D=0.7)中，圓柱上下側交替泄渦，并略微向上偏斜，形成2S泄渦模式。而在厚邊界層工況(δ/D=1.6)中，圓柱下部剪切層脫落較小的的逆時針旋渦，并在向下游移動的過程中很快消散，圓柱上部泄渦與壁面平行，圓柱下游的壁面附近渦量較低。

當G/D=1.0時，圓柱泄渦寬度比G/D=0.6工況的更大，這與圓柱更大的振幅有關。在無邊界層工況(δ/D=0)中，一個周期內圓柱上下側各泄放一對旋渦，形成2P泄渦模式，上側的渦對傾斜向上，而下側的渦對由于受到壁面的限制與壁面平行。由于下側渦對中的順時針旋轉旋渦(深灰色)強度較低，并很快耗散，因此泄渦模式逐漸變為P+S模式。隨著邊界層厚度的增加，圓柱的下側泄渦受到抑制，圓柱上部泄渦與壁面邊界層相互融合，形成傾斜向上的2S泄渦模式。

3.2 振動響應特性

分別作圓柱順流向和橫流向振動位移時程圖，如圖4所示。由于本文彈性圓柱的長細比較小，各工況下，近壁圓柱的振動響應均表現為一階振型。圓柱跨中位置處的振幅最大，受圓柱端部支撐條件的限制，圓柱兩端的振幅為零。此外，由于近壁面的影響，順流向振動頻率等于橫流向振動頻率，這與無壁面孤立圓柱渦激振動的“順流向振動頻率是橫流向振動頻率兩倍”的結論不同。

(a) G/D=0.6, δ/D=0.7

(b) G/D=0.6, δ/D=1.6

(c) G/D=1.0, δ/D=0

(d) G/D=1.0, δ/D=0.7

(e) G/D=1.0, δ/D=1.6

圖4 不同間隙比和邊界層厚度下圓柱振動位移時程圖Fig.4 The time history of the vibration amplitude of the cylinder with different gap ratio and boundary layer thickness

跨中截面振幅的最大值和均方根值，如表4～表7所示。當近壁圓柱處于邊界層外時(δ/G<1.0)，在間隙比不變的條件下，增加邊界層厚度，可使圓柱順流向振幅增大，而橫流向振幅減小。例如：當G/D=1.0時，隨著δ/D從0增加到0.7，Ax,max由0.057增大到0.138，Ay,max由0.621減小到0.6。升阻力的均方根值也觀察到相似的變化。其原因為：一方面隨著邊界層厚度的增加，圓柱下部的泄渦受到抑制，因此圓柱所受的升力減小，橫流向振幅減小；另一方面，由于邊界層中順流向速度存在梯度，振動圓柱在較高位置處阻力較大，而較低位置處阻力較小。這個隨位置變化的阻力疊加上圓柱泄渦產生的周期變化的阻力，導致了圓柱阻力變幅增大，順流向振幅增大。Chen等[30]的二維近壁面兩自由度圓柱渦激振動結果中也觀察到了類似的現象。

表4 圓柱跨中位置順流向最大振幅Tab.4 The maximum in-line vibration amplitude of the cylinder at the mid-span position

表5 圓柱跨中位置橫流向最大振幅Tab.5 The maximum transverse vibration amplitude of the cylinder at the mid-span position

表6 圓柱跨中位置順流向均方根振幅Tab.6 The root mean square of in-line vibration amplitude of the cylinder at the mid-span position

表7 圓柱跨中位置橫流向均方根振幅Tab.7 The root mean square of transverse vibration amplitude of the cylinder at the mid-span position

然而，當圓柱由處于邊界層外過渡到部分處于邊界層中時(1.0<δ/G<2.0)，在間隙比不變的條件下，增加邊界層的厚度，圓柱的順流向振幅進一步增大，而橫流向振幅則小幅增加。例如：當G/D=1.0時，隨著δ/D從0.7增加到1.6，Ax,max由0.138增大到0.146，Ay,max由0.6小幅增大到0.606。產生這樣結果的原因與順流向振動對橫流向振動的促進作用有關。Williamson等[31]指出：當圓柱發生橫流向渦激振動時，無量綱最大振幅約為1.0；而圓柱發生兩自由度渦激振動時，由于順流向振動的促進作用，圓柱的橫流向無量綱振幅可達1.5。

最后，當圓柱由部分沒入邊界層到完全沒入邊界層時，在間隙比不變的條件下，圓柱的順流向和橫流向振幅隨著邊界層厚度增加均出現減小的趨勢。例如：當G/D=0.6時，隨著δ/D從0.7增加到1.6，Ax,max由0.152減小到0.046，Ay,max由0.539減小到0.396。這是由于隨著近壁圓柱完全浸沒入邊界層中(G/D=0.6和δ/D=1.6工況)，圓柱所在位置處的來流速度明顯降低，導致了圓柱所受升、阻力大幅減小。

可見，近壁圓柱的振動幅值與圓柱和邊界層之間的相對位置密切相關。當圓柱處于邊界層之外時，隨著邊界層厚度的增大，順流向振幅增大，而橫流向振幅減小；當圓柱由在邊界層外到部分沒入邊界層時，順流向振幅增大，而橫流向振幅小幅增大；當圓柱由部分處于邊界層中到全部淹沒在邊界層中時，圓柱的順流向和橫流向振幅均大幅減小。

圓柱順流向和橫流向振動能量譜，如圖5所示。整體來看，兩個方向的振動能量均集中跨中部分，且具有相同的頻率，表現出明顯的單頻特性。需要說明的是，G/D=1.0和δ/D=0工況的順流向振動在跨中部分出現了微幅的二階振動分量，這與此時圓柱的2P泄渦模式有關(見圖3)。

(a) G/D=0.6, δ/D=0.7

(b) G/D=0.6, δ/D=1.6

(c) G/D=1.0, δ/D=0

(d) G/D=1.0, δ/D=0.7

(e) G/D=1.0, δ/D=1.6圖5 不同間隙比和邊界層厚度下圓柱兩向振動頻譜圖Fig.5 The vibration frequency spectrum of cylinder with different gap ratio and boundary layer thickness

圓柱的順流向和橫流向振動頻率，如表8、表9所示。可見，兩個方向的振動頻率均保持一致。當圓柱處于邊界層外，或者部分處于邊界層中的工況中，圓柱兩向振動頻率保持不變，均為St=0.222。然而，當圓柱完全浸沒于邊界層內時(G/D=0.6和δ/D=1.6工況)，圓柱兩向振動頻率出現明顯下降，這與降低的來流流速有關。

表8 順流向振動頻率Tab.8 The in-line vibration frequency

表9 橫流向振動頻率Tab.9 The transverse vibration frequency

通過Hilbert變換求得圓柱順流向和橫流向的瞬時相位(即φx和φy)，進而通過式求得圓柱兩向振動的相位差φx,y。

φx,y=[pφx-qφy,mod 360°]

(8)

式中，p和q為整數，并滿足fx/fy=p/q，fx和fy分別為圓柱順流向和橫流向的振動頻率。在近壁面條件下，圓柱的兩向振動頻率基本保持一致，故取p=1，q=1。

不同間隙比和邊界層厚度下圓柱兩向振動相位差均值，如表10所示。在兩個方向同頻振動的條件下，圓柱振動軌跡的旋轉方向可以通過兩向振動的相位差表征，而振動軌跡的偏斜程度可以通過兩向振動的幅值和相位差共同表征。當δ/D=0，0.7時，圓柱兩向振動相位差小于180°，圓柱振動表現為逆時針軌跡。而當δ/D=1.6時，圓柱順流向和橫流向振動相位差超過了180°，圓柱振動表現為順時針軌跡。可以看出，隨著間隙比的減小和邊界層厚度的增加，即隨著圓柱逐漸沒入壁面邊界層中，兩向振動的相位差均值逐漸增大。這與壁面邊界層中的剪切流動有關。當圓柱靠近壁面時，圓柱周圍的流速較低；反之，當圓柱遠離壁面時，圓柱周圍的流速較高。這就造成了部分或全部浸沒的圓柱在遠離壁面過程中，同時向下游運動，即圓柱的振動軌跡為順時針方向，兩向振動的相位差較大。并且，圓柱沒入邊界層越深，兩向振動的相位差越大。

表10 圓柱順流向和橫流向振動相位差均值Tab.10 The mean value of the vibration phase difference between in-line and transverse vibration (°)

沿圓柱展向等間距截取10個截面，以振動平衡位置為中心，作不同間隙比和邊界層厚度下圓柱振動軌跡。圓柱的振動軌跡為雨滴形，這與無壁面孤立圓柱渦激振動的“8”字形軌跡明顯不同，如圖6所示。

圖6 不同間隙比和邊界層厚度下圓柱各截面振動軌跡圖Fig.6 The XY-trajectory of each section of cylinder with different gap ratio and boundary layer thickness

當G/D=0.6，δ/D=0.7時，由于兩向振動的相位差均值(171°)接近180°，振動軌跡表現為瘦長雨滴形，且下部較為尖銳，上部向來流方向偏斜。當G/D=0.6，δ/D=1.6時，相位差均值(290°)遠離180°，圓柱逆時針振動，且軌跡較為圓潤，上部略偏向下游。其他3個工況具有類似的結果，此處不再贅述。

4 結 論

本文應用嵌入式迭代浸入邊界法對不同間隙比和邊界層厚度下的近壁面彈性圓柱渦激振動開展了三維數值模擬研究。分析了各工況下的流場特性和圓柱振動響應特性。研究結果表明：

(1) 圓柱發生渦激振動時，跨中部分泄渦強度較高，尾渦較寬，與壁面邊界層的相互作用較強；兩端的情況則相反。

(2) 近壁圓柱的振動幅值與圓柱和邊界層之間的相對位置密切相關。當圓柱處于邊界層之外時，在間隙比不變條件下，隨著邊界層厚度的增大，順流向振幅增大，而橫流向振幅減小；當圓柱部分沒入邊界層時，在間隙比不變條件下，隨著邊界層厚度的增大，順流向振幅增大，而橫流向振幅小幅增大；當圓柱全部淹沒在邊界層中時，在間隙比不變條件下，隨著邊界層厚度的增大，圓柱的順流向和橫流向振幅均大幅減小。

(3) 圓柱順流向和橫流向振動同頻，均呈現一階振型。當圓柱位于邊界層外或部分處于邊界層中時，圓柱振動頻率幾乎不受邊界層厚度和間隙比影響，頻率均集中于St=0.222附近。當圓柱完全位于邊界層中時，頻率明顯減小，集中于St=0.199附近。

(4) 在近壁面條件下，圓柱振動軌跡呈現水滴形或橢圓形，隨著近壁圓柱逐漸沒入邊界層中，圓柱順流向和橫流向振動相位差逐漸增大，在無邊界層和薄邊界層工況中，圓柱振動表現為逆時針軌跡，在厚邊界層工況中，圓柱振動表現出順時針軌跡。