考慮渦泄模式的低質量比圓柱渦激運動數值模擬

孫洪源，何棟梁，林海花，常 爽，于福臨

(1. 山東交通學院 船舶與輪機工程學院，山東 威海 264209；2. 中國海洋大學 工程學院，山東 青島 266100)

0 引 言

隨著海洋浮式結構的廣泛應用，Spar 海洋平臺、Spar 型浮式風電基礎等所包含的立柱式結構，在一定流速下發生渦激振動的現象。為了區別，將這種長徑比比較低、漂浮于水面的剛性結構物的渦激振動現象稱為渦激運動。渦激運動涉及諸多科目，如流體力學、計算流體力學、振動學、結構力學和統計學等[1]。

海洋結構物渦激運動（Vortex induced motion，VIM）這一熱點問題越來越受到深海工程技術人員和研究人員的關注，它是典型鈍體繞流中升力和拖曳力所導致的直接后果[2]。目前國際上已有一些學者對渦激運動進行了研究，其中包含對海洋平臺等結構物進行渦激運動集體預報等[3-5]。深水Spar 平 臺、浮筒等均具有特殊的深吃水立柱式結構，當水流經過時，圓柱體后方發生周期性渦旋分離，致使結構受到沿流向的拖曳力以及垂直于流向的升力[6]。

本文重點研究的是低質量比圓柱與浮式圓柱渦激運動特性之間的關系。將運動系統簡化為質量（m）-彈簧（k）-阻尼（ζ）系統，分析浮式圓柱運動的控制方程并通過4 階Runge-Kutta 法離散，得到一段時間內步長的速度和位移。通過軟件的自定義功能編程并嵌入到Fluent 軟件中，可以實現圓柱體的運動以及流場的更新，流場更新后再次反作用于圓柱，實現圓柱與流體的流固耦合。這已成為安全性及運動性能評估方面必須考慮的重要原因[7]。

1 數值計算模型

1.1 渦激運動的主要研究方法

隨著計算機技術的發展，計算流體力學（CFD）方法逐漸應用到研究浮式海洋平臺渦激振動的問題中[8]，采用數值模擬方法彌補了模型試驗對于一些多參數的敏感性分析和流場捕捉方面成本太高的不足。為了保證有限元仿真的可靠性，在建模時需要考慮模型的力學性質，采用合理的單元類型和應變方程去模擬實際的結構[9]。

渦激運動的研究在于規定流體的參數，因為它與許多參數有關，這些參數都會直接或者間接影響這個流體力。

1.2 渦激運動的相關參數

1）雷諾數：Re=Uc·(D/v)；

2）斯托哈爾數：S t=fv·(D/Uc)；

3）約化速度：Ur=Uc·(TSWAY/D)；

4）無量綱振幅：A/D(=[Amax?)Amin]/2D；

5）質量比：m?=m/πρD2L/4。

式中：Uc為水來流速度；D為圓柱動力直徑；V為水粘性系數；L為圓柱高度；fv為圓柱渦泄模式頻率；TS WAY為圓柱橫搖周期；A為 圓柱運動幅值；ρ為流體密度。

1.3 渦激運動求解過程

本文將彈性支撐圓柱的渦激運動系統簡化為質量-彈簧-阻尼系統，如圖1 所示。

圖1 彈性支撐圓柱渦激運動模型Fig. 1 Vortex-induced motions model of elastically supported cylinder

均勻流中，圓柱渦激運動控制方程為：

式中：t為 時間；m為質量；Cx=4πmζx fnx和Cy=4πmζy fny分別為順流向和橫流向運動的阻尼系數，其中 ζ為阻尼比；Kx和Ky分別為順流向和橫流向的系統剛度；x和y分別為順流向和橫流向渦激運動位移；Fd(t)為阻力；Fl(t)為升力。

本文計算區域與流場網格劃分情況采用圓柱直徑D=0.1m，計算區域大小為 20D×4 0D，速度入口距離圓心位置為10D，壓力出口距離圓心位置為 30D，上下邊界為滑移型距離圓心1 0D，柱體表面邊界為無滑移型。為了考慮渦激運動數值模擬的動網格應用，圓柱周圍設置了直徑 7D的隨體網格，隨體網格采用結構網格，其他區域采用非結構網格。流場計算應用Fluent 求解器，選擇SSTK?ω湍流模型，近壁面中設置了10 層邊界層，邊界層高度為0.14，網格高度符合y+≈1（y+是無量綱化的壁面距離）。動量方程中的壓力-速度耦合應用SIMPLEC 算法。為了減少數值的耗散，動量、湍流動能、耗散率應用2 階迎風格式。

網格劃分既要保證足夠的計算精確程度，又要兼顧運算的時間成本。在正式計算之前需要測試網格的密度，以獲得最佳效果。由于在數值模擬中采用了動網格技術，所以單純的測試圓柱繞流并不能獲得計算精度與時間的最優方案。本文在質量比m?=2.6，兩向固有頻率相等fnx=fny=0.38 Hz 的條件下，測試了在約化速度Ur=5 的圓柱渦激運動，結果見表1。考慮計算精度與時間成本，本文選擇了網格b 作為最優方案。

表1 網格測試結果Tab. 1 Grid test results

2 實驗結果與分析

2.1 計算工況

為了驗證數值模擬方法的可靠性，選擇質量比為2.6 的工況與Jauvtis 和Williamson[10]實驗結果進行對比。后期對質量比2.6 的圓柱與質量比為1 的浮式圓柱進行渦激運動特性研究比較，具體工況信息如表2 所示。

表2 計算工況表Tab. 2 Calculation condition table

2.2 渦泄模式的比較

圖2 為不同響應分支所對應的渦旋泄放模式：左側方為數值模擬結果，右側方為Jauvtis 和Williamson的實驗結果。

在初始渦旋泄放模式有2 種，分別為S S模式和2S模式。在初始分支的開始階段，由于順流向幅值大于橫流向幅值處于主導地位，渦泄模式表現為特有的S S（streamwise symmetric vortices shed per cycle）模式，即每個周期泄放一對對稱的渦旋，如圖2（a）所示。隨著流速增加，橫流向振幅一旦大于順流向振幅時，渦旋泄放即表現為 2S（two single vortices shed percycle）模式，如圖2（b）所示。

在上端分支，當橫向振幅達到最大值時，旋渦泄放出現了 2T（two triplets of vortices per cycle）模式，如圖2（c）所示。

在下端分支，渦旋泄放表示為 2P（two pair vortices shed per cycle）模式，即1 個周期泄放2 個渦，如圖2（d）所示。

圖3 給出了 2T模式在1 個周期的旋渦圖，當圓柱處于橫流向最大位移時，速度最小但加速度最大，在這種較大加速度作用下產生了3 個渦。其中一個渦的旋轉方向與另外2 個相反，且這個渦的強度小于其余2 個，數值計算較好模擬了這一特有現象。

圖2 不同響應分支對應的渦泄模式Fig. 2 Vortex release modes corresponding to different response branches

圖3 2T 渦泄模式1 個周期內的渦旋圖Fig. 3 2T Vortex diagram within one cycle of vortex release mode

2.3 幅值比較

圖4給出了不同約化速度下無量綱振幅與實驗結果的對比。由圖4（a）可見，數值模擬與實驗所得的無量綱橫流向振幅Ay?曲線變化趨勢是一致的，同樣表現出：初始分支I（Initial branch）、上端分支U（Upper branch）、下端分支L（Lower baranch）以及去同步化階段D（Desynchronization）。特別在初始分支、下端分支和去同步化階段，數值模擬與實驗結果吻合良好，能夠真實反映實驗結論。然而在上端分支，數值模擬沒有達到實驗中的最大橫流向振幅。同樣對于圖4（b），數值模擬與實驗所得的無量綱順流向振幅Ax?曲線趨勢大致吻合，只是在上端分支的順流向振幅有所低估。出現這種現象是因為在上端分支時，圓柱的運動幅值相對較大，大大減弱了圓柱的軸向相關性，而在二維圓柱數值模擬中，內部假設了其軸向是完全相關的[11]。另外，Fluent 數值模擬中采用了雷諾平均法，未考慮流場中的隨機擾動[12]。通過比較可以得出結論：應用數值計算可以較好的模擬圓柱渦激運動，但在最大振幅峰值方面有所低估。

圖4 不同約化速度下無量綱振幅與實驗結果的對比Fig. 4 Comparison of dimensionless amplitude and experimental results at different reduction speeds

2.4 軌跡比較

圖5給出了不同約化速度下，數值模擬與Jauvtis和Williamson 實驗所得的運動軌跡比較。數值模擬所得到的初始分支、上端分支和下端分支的圓柱運動軌跡與實驗結果基本相似，再次驗證了數值模擬方法的可靠。

2.5 考慮渦泄模式的浮式圓柱渦激運動特性

隨著約化速度的不斷增加，順流向渦激運動頻率均為橫流向渦激運動頻率的2 倍關系。Ur=4.47 處于渦激運動初始分支的后部，對比圖2（b），同樣在低阻尼比條件下，質量比2.6 的圓柱表現為 2S模式，而質量比為1 的浮式圓柱表現為 2P模式。圖6 給出了此約化速度下一個周期的渦旋泄放圖。

Ur=7.11 處于渦激運動的上端分支，與圖2（c）一樣表現出 2T模式。順流向渦激運動頻率為橫流向渦激運動頻率的2 倍關系。

圖5 不同約化速度下圓柱運動軌跡Fig. 5 Cylindrical trajectory under different reduced speeds

圖6 浮式圓柱渦激運動在 Ur =4.47 時一個周期渦泄圖（m*=1，）Fig. 6 A periodic vortex release diagram of floating cylindrical vortex induced motion at Ur=4.47（m*=1，）

Ur=11.84 處于渦激運動的下端分支，渦泄模式與圖2（d）一樣表現為 2P模式。Ur=20.26 處于下端分支與去同步化分支過渡區域，橫流向渦激運動頻率有2 個，分別是0.66 Hz 和1.61 Hz，0.66 Hz 與下端分支中的鎖定頻率相近，1.61 Hz 與該流速下的固定圓柱繞流頻率相等。雖然橫流向運動表現為2 個頻率，但由于在去同步化分支中渦激運動幅值非常小，渦泄模式表現為2S模式。

表3 給出了質量比為1 的浮式圓柱橫流向、順流向渦激運動響應譜，圖7 給出了各自對應的渦泄圖。順流向與橫流向頻率比為2 倍關系，初始分支和去同步化分支渦泄模式為 2S模式，其他響應分支為 2P模式。

可以看出它的上端分支的范圍很小，與初始分支難以區分，所以上端分支中的 2P渦泄模式表現并不明顯。與圖2 對比可見，同樣在低質量比條件下，浮式圓柱的渦旋泄放表現出了不同的形式。

表3 浮式圓柱X/D、Y/D 渦激運動響應頻譜圖（m*=1）Tab. 3 Spectral diagram of X/D and Y/D vortex induced motion of floating cylinder（m*=1）

圖7 浮式圓柱X/D，Y/D 渦激運動渦泄圖Fig. 7 Vortex-exhaust diagram of floating cylinder X/D and Y/D vortex-induced motion

3 結 語

本文應用CFD 數值模擬方法對二維圓柱在不同流速下的渦激運動問題進行研究。

1）通過與Jauvtis&Williamson 在2004 年的經典實驗結果的對比，數值計算可以較好模擬圓柱渦激運動，但在最大振幅方面有所低估。初始分支、上端分支和下端分支的運動軌跡與實驗結果基本相似。數值模擬再現了實驗中的S S， 2S， 2T和 2P渦泄模式。

2）浮式圓柱渦激運動的初始分支，表現為2 種不同的渦泄模式，分別為P+S模式和 2P模式。而質量比為2.6 的圓柱，在初始分支中渦旋泄放表現為S S模式和 2S模式。