質量比對近壁面兩向自由度圓柱渦激振動的影響

劉旭菲， 陳威霖， 及春寧

(1. 浙江水利水電學院 水利與環境工程學院，杭州 310018; 2. 天津大學 水利工程仿真與安全國家重點實驗室，天津 300072)

過去幾十年對圓柱繞流和渦激振動(vortex-induced vibrations，VIV)的研究多限定在無限區域內的自由來流條件[1-7]。然而，實際工程多涉及有壁面存在的情況，如海床上的輸油管道、電纜以及近壁的熱交換管等。壁面的存在使得近壁面一側剪切層的發展受到抑制，會改變圓柱的受力和振動響應。因此，對近壁面圓柱渦激振動的研究具有重要的工程價值和科學意義。

對近壁面圓柱繞流而言，無量綱壁面邊界層厚度(δ/D，D為圓柱直徑)、間隙比(G/D，G為圓柱下表面與壁面之間的距離)以及雷諾數(Re)等對圓柱周圍流場以及脫落旋渦與壁面邊界層的相互作用有顯著影響，并因此影響圓柱受到的流體力[8-9]。根據G/D的不同，Price等[10]將近壁面圓柱繞流分成了4個不同的區域：當G/D=0.125時，旋渦脫落被抑制，僅上側剪切層在較遠下游出現小幅擺動；當G/D=0.250～0.375時，圓柱下側剪切層開始與上側剪切層發生相互作用，使得圓柱上側發生旋渦脫落；當0.501.0時，壁面的影響幾乎可以忽略，近壁圓柱旋渦脫落與孤立圓柱的情況相同。類似地，Lei等[11]發現當G/D=0.2～0.3時，圓柱旋渦脫落被抑制，而Buresti等[12]試驗研究發現當G/D<0.4時，圓柱旋渦脫落被完全抑制。此外，隨δ/D的增加，旋渦脫落被抑制的臨界G/D增加。Zdravkovich[13]發現當圓柱位于壁面邊界層之外時，G/D對圓柱阻力的影響較弱，而當圓柱位于壁面邊界層之內時，阻力隨G/δ的增加而減小。

當圓柱可以自由振動時，圓柱與壁面剪切層發生復雜的相互作用，近壁面的影響變得更加顯著。Wang等[14]研究發現，壁面效應會顯著改變圓柱的振幅和頻率，使得圓柱與壁面邊界層之間發生顯著的非線性相互作用。Zhao等[15]數值模擬研究發現，即使間隙比很小(G/D=0.002)，圓柱的渦激振動仍可以被激發。根據圓柱與壁面作用的不同，他們將脫渦模式分成了單渦模式、碰撞反彈前脫渦模式以及碰撞反彈后脫渦模式。Yang等[16]試驗研究發現，當G/D>0.66時，圓柱的振動響應明顯不同于G/D<0.3的情況。Li等[17]發現壁面抑制了圓柱下側的剪切層，使得流向和橫向振動具有相同的振動頻率。Tham等[18]數值模擬研究發現，當Re=100和m*=10時，隨G/D的減小，鎖定區域變窄。當G/D<1.0時，交替脫落旋渦模式被明顯抑制；當G/D≤0.6時，在初始和下端分支之間出現了第三分支。Li等[19]通過二維和三維的數值模擬研究了Re=100,300的近壁面圓柱渦激振動，其中間隙比固定為G/D=0.9。研究發現，隨δ/D的增加，圓柱的振動開始于更大的折合流速上，共出現了4種不同的尾流模式：W2S(A) ——旋渦從圓柱兩側脫落，且下側旋渦可以向下游移動；W2S(B)——與W2S(A)相似，但下側旋渦更弱，且很快消失；1S ——僅圓柱上側脫落旋渦；NS——無旋渦從圓柱脫落。Barbosa等[20]研究發現，當G/D>2.0時，壁面對鎖定區域內的振動影響可以忽略；當G/D=0.75～2.00時，壁面邊界層使得圓柱振幅減小，但流場基本保持對稱；當G/D<0.75時，振動的圓柱在某些情況下與壁面發生碰撞。對近壁面圓柱渦激振動軌跡的研究發現，近壁面圓柱渦激振動的軌跡呈橢圓形，與孤立圓柱的8字型軌跡形成了鮮明的對比[21-23]。陳鎣等[24]通過試驗研究了均勻流下強迫振動的近壁面圓柱的水動力特性，其中雷諾數為Re=2×105。他們發現圓柱受到的平均阻力隨G/D的減小而下降；近壁面的存在對圓柱的能量傳遞有重要的影響。楊兵等[25]通過水槽試驗研究了固壁條件下單圓柱的流向振動，發現隨著來流流速的增加圓柱的流向振動經歷發生、發展及最后消失的過程。劉俊等[26]通過試驗研究了近壁面渦激振動觸發和停振的臨界流速，發現觸發的臨界流速隨G/D的減小而增加。

以上的研究表明，近壁面圓柱渦激振動受雷諾數和間隙比等因素的影響顯著。然而，孤立圓柱渦激振動的研究表明，質量比也是影響渦激振動的重要因素[27-31]。但已有的研究較少關注質量比對近壁面圓柱渦激振動的影響。為此，本文選取3個質量比(m*=2，10和20)，對質量比對近壁面圓柱渦激振動的影響進行深入研究，以期對近壁面圓柱渦激振動機理進行更深入的探討。

1 數值方法及其驗證

1.1 控制方程

流固耦合問題的數值模擬采用浸入邊界法[32]，控制方程如下

(1)

?·u=0

(2)

式中：u為速度；t為時間;p為壓強;ρ為流體密度;υ為運動黏滯系數；?為梯度算子；f為附加體積力矢量，代表流固耦合邊界條件。

對以上控制方程采用二階的Admas-Bashforth時間格式進行離散，可得守恒形式如下

(3)

?·un+1=0

(4)

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

(5)

式中：I(φ,Xi)和D(Φ,x)為引入的插值函數，φ為定義在流體網格節點x上的物理量，如流體的速度u、壓強p等,Φ為定義在物面邊界點Xi上的物理量，如固體速度Vn+1等;上標n+1，n+1/2，n，n-1為時間步;Δt為時間步長。

針對傳統浸入邊界法施加邊界條件精度不高的情況，Ji等提出了基于嵌入式迭代的浸入邊界法，將浸入邊界法嵌入到壓強泊松方程的迭代求解中，利用壓強的中間解比初始值更接近真實值的特點，迭代修正附加體積力，在不顯著增加額外計算耗時的前提下，提高整個算法的求解精度。有關浸入邊界法的細節，請參考Ji等的研究，此處不再贅述。

對做兩自由度運動的剛性圓柱，其運動方程可以用下述方程來描述

(6)

1.2 模型設置

在近壁面圓柱渦激振動的數值模擬中，計算域的流向和橫向寬度分別為Lx和H。如圖1所示，Lx=Lu+Ld，其中Lu=9D和Ld=80D；H=H1+H2，其中H1=1.1D和H2=97D。此外，為保證模擬結果的準確性，在圓柱周圍設置一個大小為8D×8D的加密區域，加密區域內流向和橫向網格大小均為1/64D。加密區域外的網格大小以一定的比例增加[33-36]。計算域邊界條件設置如下。入口為Dirichlet邊界條件，出口為Neumann邊界條件，上邊界為自由可滑移邊界，下邊界為不可滑移邊界。無量綱時間步長為ΔtU∞/D=0.005。需要說明的是本文所采用網格精度和時間步長均與Chen等研究中的一致，參數收斂性得到了充分驗證，讀者可自行查閱，此處不再贅述。

圖1 近壁面圓柱渦激振動的計算域和邊界條件設置Fig.1 Computational domain and boundary conditions for VIV of a near-wall cylinder

1.3 數值模擬驗證

以近壁面圓柱渦激振動為例，驗證本文數值方法的準確性。為保證結果的可比性，這里采用與驗證算例相等的計算域：圓柱中心距離入口為29D，距離出口為45D，圓柱中心距離上邊界為10D。間隙比G/D=0.6，折合流速Ur=6.0，質量比m*=10，阻尼比為零。圓柱可兩向自由度振動。如表1所示，本文的模擬結果與Tham等吻合很好，最大的誤差僅為2.4%，說明了本文數值方法的準確性。表1中，橫向無量綱最大振幅定義為Ymax/D=(ymax-ymin)/2D, 其中:ymax為最大位移和ymin為最小位移；流向無量綱均方根振幅Xrms/D，其中Xrms為流向位移的均方根值；Clrms和Cdrms分別為升力系數和阻力系數均方根值。

表1 近壁面單圓柱渦激振動的數值模擬結果與驗證算例[18]的對比Tab.1 Comparison of VIV of a near-wall cylinder with Ref [18]

2 結果和討論

2.1 質量比對振動響應的影響

圖2給出了不同質量比下圓柱橫向和流向振幅隨Ur的變化情況。為與孤立圓柱(無壁面)進行對比，圖2給出了m*=2和Re=100條件下兩向自由振動孤立圓柱的結果。本節旨在研究質量比對近壁面圓柱渦激振動的影響，故未給出m*=10,20條件下孤立圓柱渦激振動的結果?？傮w上來看，隨著m*的增加，振幅呈現下降的趨勢，且圓柱振動也開始于更大的折合流速。當m*=2時，振動從Ur=3.0開始，而當m*=10和20時，振動開始于Ur=3.6，這與孤立圓柱的情況類似[38-39]。m*=10和20工況下橫向和流向振幅差別不大。兩種工況下，圓柱的橫向最大振幅分別為Ymax/D=0.59和Ymax/D=0.56。圓柱的流向振動表現為兩個峰值，且第一個峰值明顯低于第二個峰值，這與低雷諾數條件下孤立圓柱的單峰響應不同，而與高雷諾數條件下孤立圓柱流向振動的雙峰響應[40]類似。近壁圓柱的順流向和橫流向振動均出現遲滯，隨著m*的增加，遲滯區域寬度先減小后增大，并在m*=10時達到最小。遲滯區域后，m*=2工況下圓柱的振幅緩慢下降，最終趨于零；而m*=10和20工況下，圓柱的振幅維持在零上，圓柱不再振動。需要說明的是：盡管個別工況下圓柱的橫流向振幅大于圓柱與壁面之間的間隙，但由于圓柱振動的平衡位置發生了上移，圓柱并未與壁面發生碰撞。

圖3給出了不同質量比下圓柱運動軌跡圖。如圖3(a)所示，當m*=2時，軌跡均為雨滴形。此時，圓柱的橫向和流向振動同頻。隨著Ur的增加，軌跡的上部逐漸趨向于下游，由于旋渦脫落時向下游振動的圓柱與旋渦之間的作用時間增加，旋渦產生的更大激勵作用，可激發較大振幅的振動[41]。而遲滯發生之后，振動軌跡上部趨向上游，旋渦脫落時圓柱與旋渦之間的作用時間減少，旋渦產生的激勵作用也降低，圓柱振幅下降。當m*=10和20時，軌跡基本上保持與m*=2工況相同的趨勢，但是也出現了兩個明顯的不同：一是由于流向振動2倍頻成分占比增加，在Ur較低時，軌跡表現出類8字型；二是當振幅較大時，軌跡更加趨向于橢圓形。需要說明的時：本文所有工況的圓柱振動軌跡在上部均為順時針方向，即圓柱由上游向下游通過上頂點。出于簡潔，圖3僅標記了m*=4,Ur=4工況的軌跡方向，其余工況的軌跡方向與之相同。

圖2 不同質量比條件下圓柱橫向和流向振幅隨折合流速的變化(VIV為m*=2的孤立圓柱)Fig.2 Variation of the transverse and in-line vibration amplitudes with the reduced velocity at different mass ratios (VIV represents the results of an isolated cylinderwith m*=2)

圖3 不同質量比條件下近壁圓柱運動軌跡隨Ur的變化Fig.3 Variation of the vibration trajectoryat different mass ratios

2.2 質量比對脫渦頻率的影響

圖4給出了不同質量比下無量綱脫渦頻率(St)隨Ur的變化情況。當圓柱的振動開始后，St隨著Ur的增加逐漸下降，直到遲滯時發生跳躍。m*=2工況下，尾流中旋渦脫落一直存在，但由于脫渦模式的變化，遲滯后的St比遲滯前的明顯較大，且穩定在St=0.15附近。m*=10和20工況下，遲滯后尾流不脫落旋渦。此外，對比不同質量比下的St值可以發現，在相同的折合流速下，m*較大時，圓柱的脫渦頻率較高。當脫渦頻率(圓柱振動頻率)與圓柱自然頻率相近時，鎖定現象出現。通過與自然頻率對比，在m*=2工況下，鎖定出現在Ur=3.0～6.5，而m*=10和20工況下，鎖定分別出現在Ur= 3.5～7.3和Ur=3.5～7.6上。可見，隨著m*的增加，圓柱渦激振動的鎖定區間寬度增加。

圖4 不同質量比條件下圓柱脫渦頻率隨Ur的變化(VIV為m*=2的孤立圓柱)Fig.4 Variation of the vortex shedding frequency with the reduced velocity at different mass ratios (VIV represents the results of an isolated cylinder with m*=2)

2.3 質量比對尾流模式的影響

圖5給出了不同質量比下圓柱尾流隨Ur的變化情況。由于m*=10和20工況下圓柱尾流類似，這里僅給出m*=2和10的情況。當Ur較低時，振幅較小，圓柱下側的剪切層出現明顯擺動，但并沒有旋渦脫落，此時的尾流為1S模式，如圖5(ai)和圖5(bi)所示。隨著Ur增加，振幅增加，圓柱下側剪切層與壁面邊界層發生相互作用，下側開始發生旋渦脫落，此時圓柱下側脫落的旋渦迅速分裂為兩個，形成為 C+S(C指兩個同相旋轉的旋渦)模式[42]，如圖5(aii)～圖5(aiii)和圖5(bii)～圖5(biv)所示。但是由于下側剪切層受到壁面的抑制，脫落的旋渦僅存在于較短距離的下游。隨著Ur繼續增大，振幅進一步增加，剪切層的擺動進一步增強，圓柱下側有一個旋渦脫落，此時的尾流為2S模式[43]，如圖5(aiv)～圖5(aviii)和圖5(bv)～圖5(bvii)所示。當m*=2且Ur=5.8(Ur減小)和6.4時，由于邊界層完全由間隙間通過，從下側剪切層脫落的漩渦被顯著拉長。遲滯發生后，m*=2工況下，振幅隨Ur的增加而下降，剪切層的擺動逐漸減弱，圓柱下側剪切層也不再脫落旋渦，尾流呈1S模式，如圖5(aix)所示；m*=10工況下，圓柱振幅迅速下降到零，圓柱的上下剪切層變不擺動，形成穩定尾流，如圖5(bix)所示。

圖5 不同質量比和折合流速條件下圓柱尾渦的變化Fig.5 Wake patterns at different mass ratios and reduced velocities

2.4 遲滯現象的機理解釋

在研究的質量比范圍內，近壁面圓柱渦激振動均出現了遲滯現象。本節嘗試對遲滯現象給出機理解釋。如圖6所示，增折合流速時，順流向平衡位置偏移逐漸增大，在遲滯發生時出現小幅跳躍。另一方面，遲滯發生前，橫流向平衡位置偏移隨著折合流速的增加逐漸增加，然而遲滯發生時，平衡位置跳躍到距離壁面更近的位置上，圓柱與壁面之間的間隙明顯變小。而后，隨著折合流速的增大，平衡位置偏移逐漸減小。而減折合流速時，圓柱順流向和橫流向平衡位置偏移呈現出與上述相反的趨勢。

圖6 當m*=10時圓柱橫向和流向平衡位置偏移隨Ur的變化Fig.6 Variation of the transverse and in-line shifts of the balanced position with the reduced velocity at m*=10

橫流向振動平衡位置隨折合流速偏移的這種變化，直接導致了壁面邊界層與圓柱及其剪切層相互作用模式的不同。圖7給出了m*=10和Ur=7.2工況下增、減折合流速時尾流在一個振動周期內的變化情況。如圖7(a)所示，增折合流速時，圓柱可以達到距離壁面更遠的地方，當圓柱從最高處向壁面移動時，圓柱下側剪切層可以充分發展。由于壁面之間的距離較大，間隙流也相對較強，壁面邊界層通過圓柱下方。在壁面邊界層的推動下，圓柱下側剪切層包裹圓柱的下側以及底部，使得圓柱下側的壓強更低，進而促進了圓柱向壁面運動。當圓柱距離壁面較近時，壁面邊界層重附著于圓柱的上側，與圓柱上側剪切層相融合，使得圓柱上側脫落的旋渦更強。由于壁面邊界層給圓柱上側剪切層提供渦量，使得圓柱上側的壓強更低。因此，當圓柱遠離壁面時，圓柱受到了更強的激勵作用。由于以上兩種激勵作用的影響，圓柱發生較大幅度的振動。如圖7(b)所示，減折合流速時，壁面邊界層不再與圓柱發生直接的相互作用，僅從圓柱下側通過，抑制了圓柱下側剪切層的發展，圓柱下側剪切層和壁面邊界層穩定不擺動，而上側剪切層僅在較遠下游出現擺動，這就解釋了為什么圓柱的振動較小。

綜合以上分析可知，增折合流速時，壁面邊界層會在圓柱向下運動時通過圓柱與壁面之間的間隙，而在圓柱向上運動時重復著于圓柱上側表面上，這種不穩定的相互作用模式對圓柱振動起到促進作用。而減折合流速時，壁面邊界層與振動的圓柱之間沒有直接作用，并且限制了圓柱下側剪切層的發展，對圓柱振動起到抑制作用?？偨Y來說，圓柱振動的遲滯現象與壁面邊界層在增、減折合流速條件下的雙穩態有關。需要說明的是，該遲滯現象僅出現在圓柱距離壁面較近的工況下。Chen等指出，當圓柱距離固壁較遠時，壁面邊界層與運動圓柱之間的雙穩態現象消失，相應地，該遲滯現象也不再出現。進一步增大圓柱與壁面的距離(如壁面無限遠的孤立圓柱)，在特定條件下圓柱振動的遲滯現象又重新出現，但孤立圓柱振動遲滯的成因與尾渦模式切換等有關，這與近壁圓柱振動遲滯的物理機制不同。

圖7 m*=10和Ur=7.2條件下一個振動周期內圓柱尾渦的變化Fig.7 Vortex shedding process in a vibration cycle at m*=10 and Ur=7.2

3 結 論

本文研究了質量比對近壁面兩向自由度圓柱渦激振動響應、脫渦頻率以及尾流模式的影響，解釋了近壁面渦激振動遲滯現象的機理。研究發現，質量比對近壁面圓柱渦激振動有著顯著的影響。隨著質量比的增加，圓柱的振動開始于更高的折合流速上，而且圓柱的振幅也更小一些。受壁面邊界層影響，圓柱流向和橫向的振動頻率相等。因此，在絕大多數情況下，圓柱振動軌跡為雨滴形。但當m*=10和20時，流向振動的兩倍頻成分占比增加，在折合流速較小時，呈現出類8字形軌跡。隨著質量比的增加，運動軌跡更加趨向于橢圓型。此外，還發現當振動軌跡上部偏向下游時，振動圓柱與脫落旋渦的作用時間更長，旋渦對圓柱振動的促進作用更大，而軌跡上部偏向上游時則恰好相反。

當振幅較小時，各質量比條件下，圓柱下側剪切層出現明顯擺動，但并無旋渦形成。對應地，尾流出現了1S模式。隨著振幅的增加，圓柱下側剪切層脫落兩個旋渦，形成了C+S模式。振幅更大時，圓柱下層剪切層與壁面產生強烈的相互作用，圓柱下側剪切層僅脫落一個旋渦，對應的尾流為2S模式。在遲滯區域之后，m*=2工況下，尾流為1S模式，而m*=10和20工況下，尾流為穩定模式。對增、減折合流速條件下的尾流模式分析發現，壁面邊界層重復著的雙穩態性是導致近壁面圓柱渦激振動出現遲滯的原因。當增折合流速時，壁面邊界層周期性地重復著于圓柱上側，對圓柱上側剪切層的發展起到了顯著的促進作用，并進而激勵圓柱產生更大的振幅；當減折合流速時，壁面邊界層從圓柱的下側通過，抑制了下側剪切層的發展和圓柱振動。

Vol．41 No．12 2022