陳東陽

仿真秀科普作者

(資料圖片)

1、柱體結構渦激振動定義

對于海洋工程、風工程上普遍采用的圓柱形斷面結構物，流體繞過柱體時會產生交替發放的瀉渦，這種交替發放的瀉渦又會在柱體上生成順流向及橫流向周期性變化的脈動壓力。

如果此時柱體是彈性支撐的，或者柔性管體允許發生彈性變形，那么脈動流體力將引發柱體(管體)的周期性振動，這種規律性的柱狀體振動反過來又會改變其尾流的瀉渦發放形態。這種流體-結構物相互作用的問題被稱作”渦激振動”(Vortex-Induced Vibration :VIV)，它是典型的一種流固耦合振動現象。

2、柱體結構渦激振動動力學模型及仿真方法簡介

進行二維彈性支撐柱體的數值模擬是研究海洋工程和風工程中柱體結構渦激振動現象和機理的重要手段。

根據牛頓第二定律，2-DOF彈性支撐的柱體運動的控制方程可以寫為：

(1)

式中， m為圓柱體的質量， c為結構阻尼系數， k為結構剛度系數。

式（1）又可以寫為：

(2)式（2）中，柱體固有頻率

，

阻尼比

對于渦激振動的數值模擬大多以柱體的橫向振動研究為主，同時考慮橫向和來流向的耦合振動研究相對較少。其主要原因在于，對柱體結構渦激振動數值模擬，必須保證柱體周圍網格質量非常好，才能有效預測渦激振動響應。

一般使用網格質量高的結構化網格，但是當柱體結構發生較大振動位移時，周圍流場網格會發生畸變，甚至產生負網格，導致計算失敗。如果再考慮流向耦合振動，計算難度將非常大，計算成功率也將明顯降低。

如果采用非結構化網格，且采用網格重構技術，可以吸收柱體較大的振動位移，但是非結構化網格質量相比結構化網格質量較差，且采用非結構化網格必定需要大大增加網格量，從而大大增加了計算時間。

此外，如果對于渦激振動抑制裝置設計研究，也就是在柱體結構表面形狀變復雜的情況下，網格劃分難度和計算難度也將大大增加。因此，尋求一種既可以保證網格質量，又能不大幅度增加網格數量，且可以避免網格畸變或者負網格問題的方法十分重要。

基于CFD商業軟件FLUENT和結構動力學原理，通過用戶自定義函數（UDF）及嵌套網格技術，建立了2-DOF彈性支撐柱體結構VIV數值模型。非定常不可壓縮流體RANS方程為：

(4)

(5)

式（5）中，

式中，為不可壓縮流體的密度；表示方向上的瞬時速度分量，為方向上速度脈動量，為速度的時間平均值；分別表示笛卡爾坐標系、時間、壓力、運動粘度；為湍流黏度，下標“t”表示湍流；為湍動能；是“Kronecker delta”符號，就是當時，，當時，。湍流模型選用SST k–ω湍流模型。通過計算流場，可以得到二維柱體表面的壓力分布，進而可以得到作用在二維柱體上的升力和阻力系數：

(6)

(7)

結合方程（2），2-DOF彈性支撐的柱體運動的控制方程可以寫為：

(8)（a） 兩自由度彈性支撐剛性柱體（b） 二維兩自由度彈性支撐剛性柱體VIV模型

圖1 2-DOF彈性支撐圓柱體VIV模型示意圖

(a) 背景網格與嵌套網格

(b) 挖洞和插值(c) 整個流體域的網格

圖2 2-DOF彈性支撐圓柱體流場計算網格

兩自由度彈性支撐剛性柱體在流體作用下的結構示意圖如圖1(a)所示，二維2-DOF振動柱體VIV模型示意圖如圖1(b)所示。一般柱體流場的尾跡區域需要大于等于22.5D（D為柱體直徑），整體局域高度一般需要大于等于20D，柱體振動才不受流體區域邊界的影響。

因此，綜合考慮計算條件的情況下，流場域的尺寸大小如圖1(b)中標注所示，尾跡區域30D，柱體前端和上下距離柱體都是10D。包圍這柱體的組分網格外邊界直徑大小為3D。流場入口邊界條件為速度入口，出口為壓力出口，上下壁面為滑移壁面，柱體表面即動邊界為無滑移壁面。

流場隨著柱體邊界的改變而改變，通過動網格技術來實現流場中柱體邊界的運動。嵌套網格技術是最新的動網格技術，主要適用于剛性邊界運動問題。如圖2所示，流場域網格劃分采用的是嵌套網格。如圖2（a）所示，背景網格和嵌套網格都使用結構化網格，靠近柱體表面部分為邊界層網格（），較好的保證了網格質量。

采用嵌套網格技術，可以無需擔心網格畸變以及負網格導致求解失敗等問題。同時，不會較多的增加計算量。嵌套網格即多重網格相互重疊組合成的一組網格。有可能存在兩套或者兩套以上的網格相互重疊。

嵌套網格求解的大致思路為：首先劃分包裹柱體的組分網格（組分網格數量為5262），和外流場的背景網格（背景網格數量為15731），求解器識別嵌套網格邊界，對被組分網格遮蔽的背景網格部分進行“挖洞”，然后對嵌套區域邊界單元進行插值，將背景區域的邊界單元變量信息插值到嵌套區域的邊界單元（如圖2（b）所示），最后進行流場計算。整個流場的計算網格如圖2（c）所示。

對于流場的數值計算，時間項采用全隱式積分方法，對流項則采用二階迎風離散格式。控制方程中速度分量與壓力的耦合則采用COUPLED算法進行處理。初始條件為。

圖3 2-DOF彈性支撐圓柱體VIV計算流程圖

流場域求解基于CFD商業軟件FLUENT，根據邊界條件獲得流場和二維柱體表面的壓力、速度等信息。提取作用在柱體表面的力，然后代入柱體的結構運動方程，通過求解二維柱體的運動方程，得到當前時間步長下的柱體運動的位移和速度。同時利用得到的柱體位移和瞬時速度，更新流場網格，然后進行下一個時間步的計算。這個雙向流固耦合仿真過程是通過FLUENT軟件的用戶自定義函數（UDF）來實現的。

UDF中可以使用標準C語言的庫函數，也可使用FLUENT中預定義的宏。通過預定義宏可以獲得FLUENT計算過程中的流場數據。FLUENT中用戶自定義函數是通過DEFINE宏來實現的。基于CFD的2-DOF彈性支撐柱體VIV數值求解的計算流程圖如圖3所示。圖中的虛線框內為通過C語言編制的UDF程序來實現。

3、部分仿真結果展示

（a）不同約化速度下的振幅分布（b）頻率比隨約化速度變化圖

圖4 基于CFD模型的2-DOF柱體VIV計算結果

柱體渦激振動的最大幅值及頻率比隨不同約化速度的變化如圖4所示，與實驗數據對比，誤差較小，驗證了本文計算方法的正確性。

從圖4（a）可以看出，數值仿真出3種響應分支，在=3~4之間的時候，原始分支向上端分支轉變；當=5~6之間的時候，出現下端分支。在上端分支中振幅達到最大值0.98，而在下端分支中振幅最大值0.642；從=11的時候開始，圓柱體的響應位移又回落到一個很小的數值。

從圖4（b）可以看出，在頻率“鎖定”區間=4~10內，柱體的實際振動頻率與固定柱體的泄渦頻率分離，不再符合與Re數關系圖。同時，柱體的實際渦瀉頻率與柱體固有頻率比值穩定在1.15附近，而在解鎖區域，柱體的實際振動頻率與固定柱體的渦脫頻率相同，這與前人的實驗結果大致相同。

圖5 時，不同時刻的渦量云圖（周期T=3.34s）

圖5給出了=5時彈性支撐柱體的75s~78.5s的渦量云圖，包含了一個周期的運動，從圖中可以看出，=5時的渦脫模式為P+S模式（即一個渦脫周期內有一個單個渦+一對渦形成）。

Govardhan 和 Williamson 的實驗研究表明一般在柱體振幅較大時候渦脫模式為P+S或者2P（即在一個渦脫周期內有2對尾渦形成），在振幅較小的時候渦脫模式為2S（即在一個渦脫周期內有2個單獨的尾渦形成）。圖5中的黑虛線為柱體的原始位置，紅點為柱體當前時刻的中心位置，從圖中可以看出，柱體振動游走的軌跡是一個“8”字形。

以上是ANSYSFLUENT流固耦合仿真教程中的柱體結構渦激振動仿真案例。