基于IBM 法的低雷諾數下渦激振動高質量比效應的研究

郭 濤，張紋惠,2，王文全，羅竹梅

(1. 昆明理工大學建筑工程學院工程力學系，昆明 650500；2. 中國能源建設集團云南省電力設計院有限公司，昆明 650000；3. 四川大學水力學與山區河流開發保護國家重點實驗室，成都 610065；4. 昆明理工大學能源與動力工程系，昆明 650093)

流體-結構互動問題是自然界中十分普遍的物理現象，也是一個經典的流體力學問題。當流體流過不具有流線型的鈍體結構時，會在固體壁面形成黏性邊界層，在一定雷諾數范圍內該邊界層容易與壁面發生分離，在結構尾部形成交替式的周期性脫落漩渦，產生周期性的上下拖曳力，也就是升力、阻力，從而誘發結構振動。渦激振動(Vortex-induces vibration, VIV)是促使結構產生疲勞破壞的重要原因，其危害和影響十分廣泛。例如：空氣動力學中，機翼在流動空氣中的振動問題[1-3]，導彈無側滑大攻角分行時，背風面分離渦的不對稱性導致的附加偏航力問題；水力機械中，導葉、葉片在水流作用下的擺動和振動問題[4-5]；橋梁工程中，震驚世界的塔科馬大橋顫振風毀事件(1940)[6]。近期的虎門大橋渦激振動事件(2020 年5 月5 日)[7]，拉索風雨激振問題[8-11]等；高層建筑物中，如深圳市賽格廣場大廈風致渦激共振造成的有感振動事件(2021 年5 月18 日)[12]；以及海洋立管系統(如生產立管、輸油、輸氣立管等)和系泊系統[13-14]領域十分常見的渦激振動現象等。因此，結構在流場中的渦激振動抑制以及從渦激振動中獲取能量的研究[15-22]問題成為工程界和學術界廣泛關注的問題。

1968 年Feng[23]以空氣為介質的高質量比渦激振動實驗，奠定了后續渦激振動的研究基礎。自此，各國學者對圓柱繞流渦激振動的質量比效應開展了廣泛的研究[24-39]。Mittal 等[24]通過有限元數值計算得出，低質量比(m?=4.25)單圓柱會出現“弱鎖定”現象，而較高質量比(m?=25)未發現“弱鎖定”現象。Tu 等[25]發現低質量比單圓柱不僅會出現“弱鎖定”現象，同時還會伴有大幅橫流向振動。而且Tu 等[26]和Prasanth 等[27]在低質量比串聯雙圓柱中也發現了“弱鎖定”現象。Jiang等[28]、楊驍等[29]和杜曉慶等[30]的研究也表明：質量比對雙圓柱渦激振動有顯著影響。谷家揚等[31]和陳正壽等[32]基于CFD 法研究發現低質量比單圓柱渦激振動時，鎖定區間和橫流向振幅均隨質量比的增大而減小。卞正寧等[33]基于SST 湍流模型研究發現高雷諾數下，高質量比的圓柱渦激振動會出現高幅分支。李小超等[34]研究發現質量比控制和影響了鎖定頻率及最大能量轉換率的無量綱流速范圍。實驗方面，Williamson 等[35-36]進行了一系列不同質量比m?(0.5～20)下的剛性圓柱體渦激振動和尾渦脫落實驗，為后人提供了豐富的參考數據，實驗表明：在鎖定區間和振動幅值方面，低質量比情況表現出來的結果與高質量比情況有著明顯差異，而且質量比-阻尼組合m?ζ對振動響應也有重要影響。

截至目前，對較高質量比m?情況下的渦激振動研究較少，1992 年Anagnostopoulos 等[37]做的比較經典的彈性支撐剛性圓柱體渦激振動實驗，質量比非常高，達到m?=148。趙劉群等[38]和鄧躍等[39]采用ALE 和CFD 法對該實驗參數下的低雷諾數圓柱體VIV 進行了數值計算，得到了非鎖定區間、鎖定區間和“拍頻”現象。除此以外，很多數值研究工作涉及到的質量比m?并不高，從0.5～25 不等。2021 年5 月18 日，深圳市賽格廣場大廈在低速風場(風速發生時風速約為10 m/s)作用下產生了強烈的有感振動問題。經多方論證認為是由于樓頂桅桿的第四階非對稱風致渦激振動(2.12 Hz 頻率)和大廈及桅桿動力特性改變的耦合，誘發結構的高階彎扭組合振動[40]。以此問題為出發點，為進一步澄清渦激振動的高質量比效應，本文以彈性支撐形式下的“質量-彈簧-阻尼系統”單圓柱體結構(便于和前人研究作對比)為研究對象，基于銳利界面浸入邊界法，采用C++編程計算了低雷諾數(Re=80～150)流場中，不同高質量比(m?=14.8 ～280 )和阻尼比( ζ=0.0012 ～0.036)作用下的結構振動響應、結構受力、“鎖定區間”、相位突變以及尾渦脫落問題。研究了質量比對結構振動響應特性、升阻力系數、渦脫頻率等水力不穩定現象的影響規律，分析了尾渦脫落模態的演變過程，探討了導致結構振動的流固耦合機制。

1 計算方法

為避免復雜貼體網格的更新和畸形對動邊界流場計算效率、精度的影響，以充分掌握結構場的受力特性，本文采用浸入邊界法進行流固耦合計算。浸入邊界法(Immersed boundary method，IBM)最初由Peskin[41]提出，并成功應用到模擬人類心臟血液流動中的瓣膜運動中。早期的IBM 法，是將固體邊界離散為一組由彈簧連接的單元，適用于柔性體[42]，但變形不易通過流場物理量求得，不具有一般性。Goldstein 等[43]發展了虛擬邊界法來求解剛性邊界問題，主要使用反饋力法，配合譜方法來漸進地施加所需的邊界條件，也稱反饋力法。后來Mohd-Yusof[44]提出了直接構造受力區的速度場來確保滿足邊界條件的直接力法。Fadlum等[45]、Le 等[46]、Wu 等[47]、Uhlmann[48]、王文全等[49]和郭濤等[50]將其應用于流固耦合研究，獲得了與貼體網格下計算結果一致性較好的結果。以上浸入邊界法為擴散界面法(Diffused-interface method)，均通過計算力源項來施加浸入邊界條件。與擴散界面法相對的另一類方法為銳利界面法(Sharpinterface method)，該方法將力源項的影響通過界面單元處的速度插值直接施加到計算域，即用插值方程替代浸入界面附近處的流體控制方程，實施較簡單。本文采用一種改進的銳利界面浸入邊界法[51]模擬復雜邊界的流激振動問題。該方法將整個物理區域劃分成純流體區域以及包含固體的次流體區域。在流體次區域，構造“虛擬點—受力點—垂足點”的計算結構。利用雙線性插值法對虛擬點最近的4 個歐拉網格點進行插值計算，得到虛擬點的速度，并通過強制滿足固體邊界的無滑移條件計算出受力點的速度，將其作為流動邊界條件，實現復雜動邊界的流動數值模擬。

1.1 控制方程

將整個物理區域(包括流體、固體)看作不可壓縮牛頓流體的粘性流動。在整個計算域內，流體采用Euler 描述，固體采用Lagrange 描述。其連續性方程和動量方程可表示為：

1.2 時間推進

采用二階時間分布投影法對時間進行離散，來求解控制式(1)和式(2)，具體步驟如下：

首先，在次區域內通過算法判斷和標識流體節點、固體節點，確定受力點。在任意的流體網格點周圍可以找到相鄰的4 個網格點，若相鄰的這4 個網格點中只要有一個為固體點，那么認為該點為考慮力源項影響的流體網格點即受力點，如圖1(a)所示；然后，構造如圖1(b)所示的“虛擬點—受力點—垂足點”計算結構：過受力點F 向固體邊界做垂線交固體邊界于垂足點B，反向延長垂線到虛擬點I，使點I～點F的距離等于受力點F到垂足點B的距離。利用線性插值來確定受力點F的速度uf。其中，垂足點B在固體邊界上，當圓柱為固定邊界時，其速度恒等于零。本文算例中，固體為具有單自由度彈性支承的圓柱體，發生渦激振動后，固體邊界存在被動振動的橫流向位移。所以，B點的初始速度為零，之后每一時間步下的速度由式(17)求得。

圖1 固體邊界處理Fig. 1 Solid boundary

虛擬點I的速度利用其周圍4 個Euler 網格點速度進行雙線性插值得到：

式中：un為距離虛擬點最近的4 個Euler 網格上的點(三點不可在同一直線上)； φi為常系數：

式中：

式中，下標1～4 分別指距離虛擬點I最近的4 個Euler 網格點。

3)考慮新的壓力項，更新流場，求得下一步的速度u，即：

式中，pn+1為通過解耦壓力Poisson 方程式(11)得到的下一步壓力值。

1.3 空間離散

對于對流項，設 φ為流場內的任意變量，對于流體區域內部點，采用如下格式進行離散：

在邊界位置，對流項采用一階迎風格式進行離散。

對于擴散項，中間節點處采用中心差分格式，即：

在邊界上，采用一階迎風格式。

在對壓力的離散過程中，采用五點差分離散壓力泊松方程：

對于Neumann 邊界條件，在邊界點處采用增設虛點的方法[52]。對式(14)等號右邊項，內部節點采用中心差分格式：

邊界節點采用一階迎風格式，如在左邊界上(i=0)，有：

2 渦激振動算例

2.1 計算模型

如圖2 所示，計算域為一矩形區域，長×寬=15D×10D，流體次區域為邊長3D的正方形，次區域左側距流場入口的距離為3.5D。次區域中央浸沒的剛性圓柱設置為“質量-彈簧-阻尼結構”的彈性支撐形式，如圖3 所示。圓柱直徑為D，圓心坐標為(0, 0)。計算域左側為均勻來流入口邊界，采用Dirichlet 邊界條件，即u=U，v=0；上下兩側均為無穿透邊界；右側為自由出流邊界，采用Neumenn 邊界，即 ?u/?x=0 ， ?v/?x=0。整個流場采用一套間距為 Δx=Δy=0.025的均勻四邊形網格，固體邊界采取與流體網格尺度相等的等間距離散，時間步長為 Δt=0.002 ，時間歷程t=160 s，采用本文算法計算機時為29 h (32 G 內存的普通PC 機)。

圖2 整體計算域及邊界條件Fig. 2 Computational domain and boundary conditions

圖3 渦激振動計算模型Fig. 3 Computational model of vortex-induced vibration

在流體作用下，圓柱發生被動形式的剛性縱向振動，采用一般的運動方程進行描述：

式中：y為圓柱縱向位移；m為結構的質量；c為結構的阻尼；k為結構的剛度；Fy為縱向的流體力(升力)。 計算過程中所涉及到的參數均做了無量綱化處理，見表1。

表1 各參數的無量綱化處理Table 1 The dimensionless parameter

將表1 中的無量綱化參數代入式(17)中得到無量綱化的運動方程，即：

2.2 不同來流情況對流場的影響

圖4 和圖5 為固定參數質量比m?=148，阻尼比 ζ=0.0012，約化速度U?=1 時，渦脫頻率fs、振幅比Y隨雷諾數(Re=80～150)的變化規律。

圖4 渦脫頻率 fs/fn 隨Re 數的變化規律Fig. 4 Variation of vortex shedding frequency ( fs/fn)with Reynolds number

圖5 振幅比Y 隨Re 數的變化規律Fig. 5 Variation of amplitude (Y) ratio with Reynolds number

從圖4 可以看出，渦脫頻率fs隨著雷諾數Re的變化趨勢可以分成三個階段。初期：fs隨著Re數的增大而增大，fs/fn<1這一區域是“弱鎖定”區間。第二階段：當Re數進入某一區域后，fs保持不變開始穩定下來，“渦脫頻率fs與被圓柱的固有頻率fn所捕獲，fs與fn基本相同”，這一區域便是渦激振動中的“鎖定區間”，振動頻率單一，而振幅迅速增加。后期：當Re數超過這一區間后，渦脫頻率將“擺脫”圓柱的固有頻率，fs再次隨著Re數的增大而增大，“鎖定現象”消失。實驗[37]所得鎖定區間為Re=105～125；ALE 方法[38]約為Re=98～105，偏離較遠；CFD+動網格方法[39]得到的鎖定區間為Re=90～125，區間較廣，上限值與實驗相同；而本文方法所得區間為Re=100～130，區間更加接近實驗值。從圖5 中也可看出，在鎖定區間內，渦激振動強度將達到最大，最大振動幅值均在0.5D左右。由于計算方法的差異，達到最大振幅所對應的Re有所不同，但都在鎖定區間內靠近Re數較小的一側。本文方法和實驗所得最大振幅對應的Re數都為110。

從圖6 中也可看出升力隨著雷諾數Re的增大而增大。初期Re數較小時，流場對圓柱所產生的升力系數也比較小，而影響結構渦激振動縱向位移的主要作用力就是升力，因此振幅Y也很小，基本趨近于0。進入鎖定區間后，升力增加幅度有限，而振幅卻急劇增加(最大達十幾倍)，此時結構振動不在依賴于升力，結構發生了共振，Re=110時，橫流向振幅達到最大；當Re=130 時，結構開始擺脫鎖定狀態，振動開始變弱，而且此時橫流向振幅與升力系數發生了明顯的“相位突變”現象；當Re=150 時，結構已經完全擺脫鎖定狀態，此時不僅渦激振動變弱、消失，振幅與升力系數仍然處于反相位狀態，而且升力系數隨時間的推進逐漸出現了“拍頻”現象。

圖6 不同雷諾數下升力系數CL 及振幅比Y 時程曲線Fig. 6 The time-history curve of lift coefficients (CL) and amplitude ratio (Y) with different Reynolds number

圖7 為鎖定區間內，當振幅達最大時，一個周期內尾跡渦的演化過程。從圖7 中可知。圓柱后面正、負漩渦交替脫落，向下游發展，呈現出典型的卡門渦街現象，其渦脫落模態為2S 型，圓柱表面的邊界層分布、分離剪切層的卷起和圓柱后方的分離區域清晰可見。與固定繞流不同的是，隨著圓柱的上下振動，漩渦中心與圓柱中心不再總是處于一條水平線上。

圖7 一個振動周期T 內尾跡渦隨時間的演化(Re=110，范圍： -8≤w≤8，間隔1，圖中實線為正值，虛線為負值)Fig. 7 The evolution of vorticity Isolines with time in a vibration period (Re=110, Range:-8≤w≤8, Spacing 1, Solid lines are positive isolines,Dashed lines are negative isolines)

2.3 質量比 m? 、阻尼比 ζ 和質量-阻尼比組合m?ζ對渦激振動的影響

圖8 給出了不同質量比結構下，振動振幅隨雷諾數Re的變化規律。

由圖8 可知，不同高質量比m?下，結構振幅曲線趨勢基本相同。當m?參數發生成倍變化時，曲線的變化并不明顯。結構的鎖定區間都位于Re=100～130，而且發生共振時對應的雷諾數均為110，位于鎖定區間靠近Re較小的一側，最大振幅也都在0.5D左右，無太大變化，即m?對振幅的影響并不大。現如今，新能源已被廣泛利用，從渦激振動中獲取海洋能的多少與結構受流場力的作用而產生的振幅大小息息相關。因此，若只考慮減小結構質量比m?的方法，對提高結構振動幅度的效果有限。

圖8 不同質量比 m?對比振幅比Y 的影響Fig. 8 Variation of amplitude ratio (Y) with mass ratio ( m?)and Reynolds number

圖9 為相同質量比m?=14.8 、不同阻尼比ζ(考慮了振動系統阻尼、彈簧剛度的綜合系數)對渦激振動幅值的影響。

圖9 不同阻尼比 ζ對振幅比Y 的影響Fig. 9 Variation of amplitude ratio (Y) with damping ratio ( ζ)and Reynolds number

由圖9 可知，在質量比相同的情況下，不同阻尼比結構系統的振動響應隨雷諾數的變化趨勢大致相同。 ζ的變化并沒有對渦激振動的鎖定區間造成影響；雖然鎖定區間和變化趨勢相同，但當ζ發生變化時，在鎖定區間內，圓柱發生渦激振動的最大響應幅值是有差別的，隨著 ζ的增加，振幅比Y逐漸減小，振動受到了抑制。所以，質量比m?相同的情況下，減小阻尼比ζ，無法使圓柱的鎖定區間增大，但是對提高結構振動幅度的效果明顯。

圖10 討論了高質量比m?=14.8和較高質量比m?=148結構系統，在質量阻尼比相同情況下(m?ζ=0.1776)，渦激振動振幅隨雷諾數變化的關系曲線。

圖10 質量阻尼比對縱向振幅的影響Fig. 10 Variation of amplitude ratio (Y) with mass damping ratio ( m?ζ) and Reynolds number

從圖10 可以得到，當m?ζ保持不變時，圓柱的渦激振動幅值基本相同，最大值還是在0.5D左右；但是，鎖定區間卻發生了較大變化。相對低的質量比結構系統，由非鎖定狀態進入鎖定區間的雷諾數更低，鎖定狀態持續時間長，擺脫鎖定狀態也晚。因此，質量-阻尼比相同時，質量比低的結構系統發生渦激振動時的鎖定區間更廣，雷諾數大致為Re=90～140，比m?=148 的較高質量比結構系統提高了1.67 倍，能更早的進入鎖定區間。而且，質量比低時，結構系統發生最大振動響應時對應的雷諾數也小。

3 結論

本文利用一種銳利界面浸入邊界法對二維彈性支承圓柱流固耦合問題進行了數值計算，并討論了不同參數(雷諾數Re、阻尼比 ζ 、質量比m?)對圓柱縱向渦激振動的影響，通過對比分析，得到以下結論：

(1)通過與其它文獻和實驗結果的對比，驗證了本方法的準確性和有效性。

(2)對于圓柱繞流問題，當圓柱發生渦激振動時，在鎖定區間內，振動強度將達到最大。本文計算結果與實驗結果[37]均表明最大振動幅值均在0.5D左右，此時所對應的雷諾數都為Re=110，都在鎖定區間內靠近雷諾數較小的一側。實驗所得鎖定區間為Re=105～125，本文方法所得區間為Re=100～130。

(3)質量比m?對鎖定區間、最大振動幅值以及發生共振時對應的雷諾數都沒有太大影響；同樣，阻尼比 ζ對鎖定區間并無太大影響。但是，對鎖定區間內的最大振動幅值是有影響的，隨著 ζ的增加，振動受到了抑制。因此，相反若是對于海洋能發電等新能源領域，應該考慮從減小振動系統阻尼、彈簧剛度等方面著手，以便獲取更多海洋能。

(4)質量阻尼比m?ζ相同情況下，雖然振動振幅不會有明顯變化，但是質量比低的結構系統鎖定區間更廣。發生最大振動響應時對應的雷諾數也小。