雙自由度渦激振動數值模擬方法研究

孫麗萍， 張 旭， 倪問池

(哈爾濱工程大學 深海工程技術研究中心，哈爾濱 150001)

雙自由度渦激振動數值模擬方法研究

孫麗萍， 張 旭， 倪問池

(哈爾濱工程大學 深海工程技術研究中心，哈爾濱 150001)

當前雙自由度渦激振動數值模擬的研究模擬結果精度并不十分理想，尤其是對于最大振幅的響應的估計。在Jauvtis與Williamson的實驗中，觀測到低質量比圓柱的最大振幅可達1.5D，并將其命名為“超上端分支”。然而，當前研究很少有模擬結果能達到1.3D，甚至無法捕捉到超上端分支，其原因很可能是湍流模型本身的缺陷以及數值模擬參數設置的不合理。針對上述問題，基于開源軟件OpenFOAM，采用了改進的k-ε湍流模型，研究了加速度對渦激振動響應的影響，優化了數值模擬方法。通過振幅、相位、軌跡、尾渦等響應的分析，證明在適當的加速度下，運用改進k-ε模型可對雙自由度渦激振動實現較精確的數值模擬。

渦激振動；k-ε湍流模型；加速度；數值模擬

隨著海洋工業的發展，渦激振動問題日益受到關注。數值模擬是研究渦激振動的重要方法，關于雙自由度渦激振動數值模擬的研究成果已有很多[1]，然而，已有的模擬結果精度并不理想。特別是對于最大振幅的響應的估計，在Jauvtis等[2]的實驗中，觀測到低質量比圓柱的最大振幅可達1.5D，并將其命名為“超上端分支”[3]，而從當前的文獻中看，很少有模擬結果能達到1.3D，甚至無法捕捉到超上端分支。

Pan等[4]運用SST模型對低質量比剛性圓柱進行了雙自由度渦激振動數值模擬，未能捕捉到上端分支。同樣，Sanchis等[5]也未能在雙自由度圓柱數值模擬中捕捉到上端分支。谷家揚等[6]對低質量比圓柱進行了單自由度以及雙自由度的渦激振動數值模擬，其橫向運動最大幅值均未超過1，也沒能捕捉到超上端分支。林琳等[7]運用RNG以及SST湍流模型分別對質量比為2.4的圓柱進行單自由度渦激振動數值模擬，其結果顯示，SST湍流模型的計算精度高于RNG模型，但其橫向運動最大幅值的計算結果仍較實驗值明顯偏小。Srinil等[8]對細長桿件進行了不同頻率比的雙自由度渦激振動數值模擬，并與實驗值進行了對比。根據其模擬結果可以發現，數值模擬結果的橫向最大振幅均小于實驗結果。

究其原因，可分為兩方面。首先，標準湍流模型在模擬分離流時會產生誤差。Medic[9]認為，標準雷諾平均模型由于自身不足，在計算不穩定分離流時會產生誤差，并會抑制漩渦的分離。Younis[10]經過對大量實驗數據分析之后提出，在模擬渦激振動的過程中，應當考慮從周期性的時均流場傳輸至湍流脈動場的能量傳輸過程。另外，對于數值模擬的方法，Guilmineau等[11]應用了“掃頻法”，并捕獲了上端分支。

針對上述問題，本文借鑒了Younis與Guilmineau的思路，基于開源軟件OpenFOAM，采用改進的k-ε湍流模型，研究了加速度對渦激振動響應的影響，優化了數值模擬方法，并對低質量比圓柱的振幅、相位、軌跡、尾渦等響應進行了分析。

1 湍流模型

1.1 改進k-ε湍流模型

標準兩方程k-ω模型是工程中常用的湍流模型之一，它運用兩個額外的輸運方程來表示的湍流的特性。這使得兩方程k-ε模型能考慮到歷史影響，如對流和擴散的湍動能。其具體形式如下：

輸運方程

湍動能：

(1)

耗散率：

(2)

渦黏性

(3)

k的產生項

(4)

式中，各系數的取值如表1所示。

表1 湍流模型系數Tab.1 Coefficients of turbulence closures

根據已有文獻，標準 RANS 法對于模擬渦激振動這類強分離流動的結果難以令人滿意[12-13]。對此，Younis等經過對大量實驗數據分析，提出以下觀點：在渦激振動的湍流流場中，伴隨著漩渦脫落，能量將從周期性的時均流場傳輸至湍流脈動場中。并指出，渦激振動流場中，能譜函數應該表示成如下形式：

E(κ,t)=(A0+A(t))κs

(5)

式中：κ為波數；s為匹配指數；A(t)為附加項，且在平穩流動中應當取0。

根據湍動能k與能譜函數的關系[14]

(6)

(7)

可以推導出適用于渦激振動的耗散率ε的表達式為

(8)

根據式(8)，應當對標準k-ε模型耗散率方程中的系數C1ε進行如下修正。

(9)

式中：Q為時均流場中每個單位體積中的湍動能；Ct為經驗常數，取為0.38。

2 數值模擬

2.1 計算域及網格

本次數值模擬所采用的網格如圖1所示。其中，坐標系定義：以圓柱中心點為原點，來流方向為x軸方向，垂直來流方向為y軸方向，z軸方向為圓柱展向。根據Willamson的實驗數據，設定D=0.038 1 m，計算域劃分為：-12D

圖1 計算域整體網格Fig.1 Global mesh for circular cylinder

網格劃分：圓柱周圍3D的范圍內采用O 型網格，在其他區域使用六面體網格。網格數量根據雷諾數不同進行相應的調整，改變貼近圓柱表面的第一層網格厚度，保證y+<1。邊界條件為：圓柱表面采用無滑移固壁，出口采用壓力梯度為零，圓柱繞流模型入口采用均勻來流，渦激振動模型運用自編程序，入口采用先勻加速，然后勻速的條件。其余均為對稱邊界。本次數值模擬采用了兩套網格，其具體信息如表2所示。

時間步長根據庫朗數確定。庫朗數的表達形式如下：

Courant number=u·Δt/Δx

(10)

式中：u為流速；Δt為時間步長；Δx為網格尺寸。在數值模擬過程中，根據不同的流速調整時間步長，使得庫朗數最大值為0.2左右。

表2 網格參數信息Tab.2 Parameters of numerical grids usedin numerical simulation

3 計算結果及分析

3.1 加速度對渦激振動響應的影響

G2號網格被應用于探究加速度對渦激振動響應的影響，具體分為2個階段：第一階段是勻加速階段。在這個階段中，通過自編程序，令流場速度從零增加到目標值。第二階段是勻速階段。當流速達到目標值后，保持恒定，令圓柱繼續振動，直到運動達到穩定。隨后，提取穩定段的橫向振幅和頻率的響應進行分析，并與實驗結果進行對比。

本文考察的重點為渦激振動的超上端分支響應。因此，在本節中，目標速度取為約化速度Ur=6，對應于超上端分支的穩定階段。加速時間設定為從10～200 s，對應于加速度的值為0.012 5～0.25每約化時間(Ut/D)。此外，勻速工況(流場初始速度直接設為Ur=6)也被包括其中。

不同加速度情況下，穩定階段圓柱的振幅以及頻率響應與圖2與圖3所示。其中，振幅取的是橫向最大無量綱振幅，即橫向最大振幅與直徑的比值，可表示為Aymax/D。頻率響應取的是橫向頻率比fy/fn，即橫向振動頻率與固有頻率的比值。此外，fst為固定圓柱的斯特勞哈爾頻率。

從圖2中可以看出，當加速時間小于100 s，對應于加速度大于0.025每無量綱時間時，穩定段的最大振幅響應幾乎與勻速工況相同，約為0.8D左右，明顯小于實驗值。

當加速時間大于100 s時，最大振幅響應隨著加速時間的增加而增大，直到加速時間為150 s左右，對應于加速度大小為0.017每無量綱時間以后，振幅達到最大值。此后，最大無量綱振幅保持不變，約為1.1左右，與實驗值接近。

圖3為不同加速時間下的橫向振動頻率響應。從圖中可以看出，當加速時間低于150 s，對應于加速度大于0.017每無量綱時間，橫向振動頻率接近斯特勞哈爾fst，表明未發生鎖定現象。當加速時間大于150 s時，頻率比變為1左右，接近圓柱在靜水中的固有頻率，表明此時發生了“鎖定”現象,這也是振動進入上端分支的特征之一。

綜合振幅與頻率的分析可以得出，數值模擬過程中，流場的加速度大小對響應結果有重要影響，只有當加速時間大于150 s，即加速度小于0.017每無量綱時間時，雙自由度渦激振動數值模擬才能捕捉到上端分支。

圖2 橫向最大振幅隨加速時間響應曲線

Fig.2 The maximum non dimensional transverse amplitude responses with different acceleration magnitudes

圖3 橫向振動頻率隨加速時間響應曲線

Fig.3 The transverse vibration frequency ratio responses, under different acceleration conditions

3.2 渦激振動結果分析

根據渦激振動響應隨加速度的變化規律，在雙自由度圓柱渦激振動數值模擬過程中，取加速度大小為0.015每無量綱時間。所得的最大振幅響應曲線如圖4所示。

從圖4中可以看出，采用改進k-ε湍流模型，在加速度0.015每無量綱時間的條件下，捕捉到了明顯的上端分支。最大振幅達到了1.4D左右，與實驗值較為接近。此外，U*=6.8時的渦泄頻率譜也在圖4中給出。可以看出，當振幅達到最大值時，渦泄頻率約為0.4 Hz，與固有頻率一致，此時振動頻率鎖定于固有頻率，發生鎖定現象。相較于Sanchis等未能捕捉到上端分支，且最大值約等于1D的數值模擬結果，本文采用的數值模擬方法大大提升了渦激振動模擬精度。

圖4 不同約化速度下的最大振幅響應Fig.4 Response of amplitude versus reduced velocity U*

此外，如上文所述，根據已有的文獻，當雷諾數較大時，運用標準RANS湍流模型模擬渦激振動，會產生誤差，并抑制漩渦的分離，對應于約化速度大于10以后，振幅衰減過快。從圖4的結果看，改進k-ε湍流模型的計算結果與實驗值較為接近，并未出現衰減過快的現象。可見改進k-ε湍流模型對模擬不穩定分離流動具有良好的效果。

總而言之，數值模擬的振幅響應與實驗結果吻合度較高，只是在從上端分支跳轉到下段分支對應的約化速度上略有差異，這也許是實際實驗中的三維效應以及其他干擾因素造成的。

圖5為典型約化速度下升力系數、阻力系數與位移歷時曲線。其中，Ur=3、6、7分別對應著初始分支、超上端分支以及下端分支。需要注意的是，本次數值模擬包括勻加速以及勻速兩個階段，對于Ur=3、6的情況，圖中取的為勻速穩定階段的數據，對于Ur=7的情況，圖中取的為Ur=7附近的從加速到勻速直到穩定的過程的數據。

從圖5中可以看出，當Ur從3變為6，對應于運動從初始分支跳躍到上端分支后，圓柱的升力系數、阻力系數與位移顯著增大，結合3.1節中的頻率分析，可以得出，此時振動發生了鎖定現象，圓柱的橫向振動被鎖定于固有頻率附近，產生共振，使得振幅迅速增加。

從圖5(c)中可以看出，當約化速度到達7附近時，圓柱升力系數、阻力系數與振幅突然顯著降低，預示著振動響應由超上端分支跳躍到下端分支。

振動過程中位移與升力的相位角的變化可以通過分析位移與升力系數曲線得到。分析圖5(a)與圖5(b)可得，當振動處于初始分支以及上端分支時，位移與升力系數的變化趨勢基本一致，相位角約為0°。此時升力對振動起促進作用，使得振幅不斷增加。

當Ur達到7附近時，如圖5(c)，對應著從超上端分支跳轉到下端分支，隨著振幅與升力系數的突然降低，位移與升力的相位角也突然從0°變化到180°左右。此后，升力與位移保持同相位，對振動起抑制作用，因此振幅大大減小，這與渦激振動基本理論與實驗現象吻合，可見本文采用的數值模擬方法能很好地捕捉渦激振動升力、阻力、位移以及相位的變化規律。

(a) Ur=3

(b) Ur=6

(c) Ur=7圖5 典型約化速度下升力系數、阻力系數與位移歷時曲線

Fig.5 Time histories of lift coefficientCt, drag coefficientCdand displacementy/Dat critical reduced velocities

本次數值模擬的振幅與尾渦響應如圖6(a)所示。此外，Jauvtis與Williamson實驗中觀測到的軌跡以及尾渦的DPIV結果被展示在圖6(b)中。

從圖6中可見，起初，橫向振幅很小，軌跡呈橫條形。隨著速度的增加，在初始分支的穩定階段，軌跡為傳統的“8字形”。當進入超上端分支以后，橫向與縱向振動的相位角發生變化，軌跡形狀逐步從“8字形”過渡到“新月形”。當運動跳轉到下端分支時，縱向振幅遠小于橫向振幅，軌跡呈“瘦8字形”。數值模擬的軌跡變化規律與實驗結果完全吻合。

(a) 數值模擬結果

(b) 實驗結果圖6 不同約化速度下的軌跡與尾渦形態

Fig.6 The trajectories and vorticity contours compared with Jauvtis and Williamson’s results (2004)

從尾渦形態上看，在初始分支中，尾渦的形態為“2S模式”，即在每個周期內，兩個獨立的尾渦均勻形成并從物體后方交替脫落。當振動處于上端分支時，特有的“2T模式”在數值模擬中被清楚地捕捉到。此時，物體后方每個周期發放出兩個渦對，每個渦對包含三個渦，并且三個渦的旋轉方向不完全相同。當振動跳轉到下端分支以后，尾渦形態變為“2P模式”，在每個周期內，物體后方發放出兩個渦對，并且渦對中的兩個漩渦旋轉方向相反。數值模擬的結果與DPIV實驗結果一致。

4 結 論

本文針對雙自由度渦激振動數值模擬中，難以捕捉到上端分支，以及下端分支振幅衰減過快的現象，提出了兩點改進措施，包括改進傳統的湍流模型以及采用合理的加速度值。

首先，針對標準雷諾平均模型，在計算不穩定分離流時會產生誤差，并會抑制漩渦的分離的缺陷，借鑒Younis在模擬渦激振動的過程中，應當考慮從周期性的時均流場傳輸至湍流脈動場的能量傳輸過程的思想，并基于OpenFOAM開源軟件，改進了標準k-ε湍流模型。

其次，探究了流場加速度與振幅響應之間的關系。發現當約化速度從0～6，加速時間小于150 s，對應于加速度大于0.017每無量綱時間時，橫向振幅較小且振動頻率接近fst,無法進入超上端分支。只有當加速度小于0.017每無量綱時間時，才能捕捉到超上端分支。

最后，運用改進的k-ε湍流模型模擬圓柱雙自由度渦激振動。數值模擬最大振幅響應與實驗結果基本吻合，最大振幅達到了1.4D。此外，還捕捉到了升力系數、阻力系數、位移曲線的變化趨勢，相位角的突變，軌跡的變化規律以及特有的“2T模式”尾渦。

運用改進k-ε模型，在適當的加速度下，可對雙自由度渦激振動實現較精確的數值模擬。

[1] 范杰利, 黃維平. 細長立管兩向自由度渦激振動數值研究[J]. 振動與沖擊, 2012, 31(24):65-68.

FAN Jieli, HUANG Weiping. Numerical simulation of 2-DOF vortex-induced vibration of a long riser[J]. Journal of Vibration & Shock, 2012, 31(24):65-47.

[2] JAUVTIS N, WILLIAMSON C H K. The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[J]. Journal of Fluid Mechanics, 2004, 509:23-62.

[3] WILLIAMSON C. A high-amplitude 2T mode of vortex-induced vibration for a light body in XY motion[J]. European Journal of Mechanics-B/Fluids, 2004, 23(1):107-114.

[4] PAN Z Y, CUI W C, MIAO Q M. Numerical simulation of vortex-induced vibration of a circular cylinder at low mass-damping using RANS code[J]. Journal of Fluids & Structures, 2007, 23(1):23-37.

[5] SANCHIS A, S?EVIK G, GRUE J. Two-degree-of-freedom vortex-induced vibrations of a spring-mounted rigid cylinder with low mass ratio[J]. Journal of Fluids & Structures, 2008, 24(6):907-919.

[6] 谷家揚, 楊建民, 肖龍飛. 兩種典型立柱截面渦激運動的分析研究[J]. 船舶力學, 2014(10):1184-1194.

GU Jiayang, YANG Jianmin, XIAO Longfei. Study on vortex induced motion of two typical different cross-section columns[J]. Journal of Ship Mechanics, 2014(10):1184-1194.

[7] 林琳, 王言英. 不同湍流模型下圓柱渦激振動的計算比較[J]. 船舶力學, 2013，17(10):1115-1125.

LIN Ling, WANG Yanying. Comparison between different turbulence models on vortex induced vibration of circular cylinder[J]. Chuan Bo LI Xue/journal of Ship Mechanics, 2013, 17(10):1115-1125.

[8] SRINIL N, ZANGANEH H, DAY A. Two-degree-of-freedom VIV of circular cylinder with variable natural frequency ratio: Experimental and numerical investigations[J]. Ocean Engineering, 2013, 73(8):179-194.

[9] MEDIC G. Etude mathématique des modèles aux tensions de Reynolds et simulation numérique d’écoulements turbulents sur parois fixes et mobiles[D]. Bibliogr, 1999.

[10] YOUNIS B A, PRZULJ V P. Computation of turbulent vortex shedding[J]. Computational Mechanics, 2006, 37(5):408-425.

[11] GUILMINEAU E, QUEUTEY P. Numerical simulation of vortex-induced vibration of a circular cylinder with low mass-damping in a turbulent flow[J]. Journal of Fluids & Structures, 2004, 19(4):449-466.

[13] JUNG Y W, PARK S O. Vortex-shedding characteristics in the wake of an oscillating airfoil at low Reynolds number[J]. Journal of Fluids & Structures, 2005, 20(3):451-464.

[14] POPE S B. Turbulent flows[M]. cambridge university press, cambridge, U.K. 2000, 771.

Numericalsimulationmethodfor2-DOFvortex-inducedvibration

SUN Liping， ZHANG Xu， NI Wenchi

(Deepwater Engineering and Technology Research Center，Harbin Engineering University，Harbin 150001， China)

The existing simulation results for 2-DOF vortex-induced vibration are not satisfied, especially, the maximum amplitude estimation of its responses. In tests of Jauvtis & Williamson, the maximum amplitude of a cylinder with a low mass ratio could reach 1.5D, this phenomenon was named “super upper branch”. But, few simulation results of current studies can reach 1.3Deven unable to capture super upper branch. The reasons were likely to be defects of the turbulence model itself and the unreasonable setting of numerical simulation parameters. Aiming at the problems mentioned above, an improvedk-εturbulence model was used to study the effect of acceleration on responses of vortex-induced vibration here based on the software OpenFOAM and optimize the numerical simulation method. By analyzing responses of amplitude, phase, trajectory and tail vortex, it was shown that 2-DOF vortex-induced vibration can be numerically simulated more accurately with the improvedk-εturbulence model under appropriate accelerations.

vortex-induced vibration (VIV);k-εturbulence model; acceleration; numerical simulation

國家工信部海洋工程裝備科研項目(G014614002)

2016-06-14 修改稿收到日期：2016-10-19

孫麗萍 女，教授,博士生導師

張旭 女，博士生，1991年5月。E-mail：zhangxu_heu@163.com

U661.1

A

10.13465/j.cnki.jvs.2017.23.004