基于SST k-ω湍流模型的二維圓柱渦激振動數值仿真計算

李　駿，李　威

(華中科技大學 船舶與海洋工程學院，湖北 武漢 430074)

基于SSTk-ω湍流模型的二維圓柱渦激振動數值仿真計算

李駿，李威

(華中科技大學 船舶與海洋工程學院，湖北 武漢 430074)

摘要：隨著海洋工程逐漸向深海發展，廣泛使用的柔性立管因為高頻振動很容易受到嚴重的損傷。因此對于這類細長柔性結構的渦激振動(VIV)研究是當今的一個熱點，同時隨著計算機的快速發展，CFD (計算流體力學) 技術成為研究渦激振動問題不可或缺的一種方法。本文采用雷諾平均納維爾-斯托克斯(RANS)方程，并結合SST k-ω 湍流模型，研究低質量比彈性支撐剛性圓柱體的渦激振動問題。從振幅響應、頻率響應、3個響應分支的水動力性能、尾渦模式等方面和Williamson相關實驗作對比。結果表明，SST k-ω 湍流模型能夠有效準確地模擬圓柱繞流的渦激振動。本文豐富了海洋工程的理論研究，為柔性立管的實際應用提供了一定的理論指導。

關鍵詞：渦激振動；SST k-ω湍流模型；立管

0引言

很多工程領域中都會涉及到渦激振動 (VIV) 問題，在海洋管道、立管或者其他水下結構中尤其明顯。近年來隨著柔性材料的廣泛應用，產生了大量關于低質量比的圓柱渦激振動研究成果[1-3]。

對于彈性支撐的剛性圓柱在自由來流下的橫向渦激振動研究，Williamson 的團隊通過一系列實驗[4-7]為這項研究做出了巨大貢獻。他們的實驗結果可以概括成以下3點：一是對于低質量比的彈性支撐剛性圓柱在橫向振動方向會產生3種明顯的振幅響應分支，分別是初始分支、上端分支和下端分支；二是在不同的響應分支，尾渦的脫落模式也各有不同；三是在相鄰響應分支的過渡點會伴隨著相應的相位角突變。

借助CFD (計算流體力學) 模擬渦激振動時，大致有3個方法處理N-S(納維爾-斯托克斯) 方程，分別是雷諾平均納維-斯托克斯(RANS)、直接數值模擬(DNS)和大渦模擬(LES)。相比其他2種方法，RANS方法已經足夠成熟，并且可以在較短的計算時間內得到相當精確的結果。Pan和Cui[8]在Khalak和Williamson實驗基礎上使用RANS方法來模擬二維數值模型。他們發現在渦脫模式和不同分支之間的轉變方面與實驗結果基本一致，但是模擬獲得的響應振幅遠低于實驗結果，尤其是最大響應振幅的時間點只有1個，這與實驗結果完全不一致。

本文參照Khalak 和 Williamson 的實驗模型，定義圓柱的質量比m*=2.4 ，質量阻尼比 m*ξ=0.013；折合速度U*=Uinlet/(fnD)從2.0增加13.9；相應的雷諾數從1 700增加到11 600。采用SSTk-ω湍流模型研究彈性支撐剛性圓柱在自由來流下的橫向渦激振動問題，為評判該模型在這一研究的正確性，所有數值仿真結果都將與Khalak 和 Williamson 的相關實驗結果作對比。

1計算模型

1.1　控制方程

對于非定常不可壓縮流體，其RANS方程可以寫作：

(1)

(2)

(3)

式中：渦流粘度μt由湍流模型給出；δij為克羅內克函數；k為湍流動能。

(4)

同時，結合“切片理論”，可將計算模型簡化至二維：首先將結構看成大量的沿立管軸向的切片，每個切片都可以當做1個彈性支撐的剛性圓柱，并且等效于只有1個自由度的線性彈簧質量阻尼系統。簡化過程如圖1所示。

圖1　三維立管簡化過程Fig.1　The simplification process of a 3D riser

因此，圓柱在其橫向方向的無量綱運動方程為：

(5)

式中：U*為折合速度；ζ為系統結構的阻尼比；m*為質量比；y為無量綱的橫向位移；CL為橫向的升力系數。

比較兩組治療前后呼吸困難程度評分（MMRC）、6 min步行距離（6MWD）、肺功能狀態（FEV1預計值）等。

當圓柱渦激振動進入“鎖定區域”，尾渦的渦脫頻率與圓柱振動頻率完全相同，此時圓柱的位移為：

y=A*sin(ωext)，

(6)

CL=CL0sin(ωext+φ)，

(7)

CL=CLvcos(ωext)+CLa[-sin(ωext)]，

(8)

CLv=CL0sinφ,CLa=-CL0cosφ，

(9)

φ=arctan(-CLv/CLa)。

(10)

式中：CL0，CLv及CLa分別為升力系數的振幅及它在速度方向與加速度方向上的分量；ωex為圓柱振動圓頻率；φ為升力系數和響應位移之間的相位角。

根據Parkinson[9]的研究，響應振幅比(A*=Ymax/D)和頻率比(f*=fex/fn)為

(11)

(12)

因此，一旦得到了響應振幅比和頻率比，就可通過上面2個方程反推得到CLv和CLa，以及相位角φ。

此外，根據Lighthill[10]的研究，圓柱體上的總流體力F可以分解為附加質量力Fpotential和旋渦力Fvortex。同樣，將上述力無量綱化后可得到：

Cvortex=CL0-Cpotential，

(13)

其中，

(14)

類似地，相應的無量綱運動方程還可以表示為：

(15)

式中：ma為圓柱的附加質量；ζ′為包含了附加質量系統結構的阻尼比；φvortex為旋渦升力系數與響應位移之間的相位角，可通過之前提到的相似方法求解得到。

采用Newmar k-β方法求解圓柱的橫向振動響應位移，同時通過求解二維的RANS方程獲得相應的流體力系數。

1.2　網格劃分及邊界條件

圖2(a)為仿真模型的計算域，圓柱直徑為D，計算域是一個寬為20D和長為30D的矩形。來流的邊界條件在圓柱體中心左邊10D的位置處，計算域右端則是流動出口邊界。對稱邊界則分別在圓柱體上下10D的位置處，無滑移邊界則是在圓柱體表面處。圖2(b)為計算模型的網格劃分，采用混合網格劃分方式，由小及大的網格分布不僅可以減少仿真計算時間，并且可以得到更精確的計算結果。圓柱周圍較密的網格可以精確模擬流動的尾渦情況。對于任意的折合速度，圓柱體外的第一圈網格始終劃分在y+≤0.5的范圍內。

圖2　計算域和網格模型Fig.2　Computational field and grid model

2結果與分析

采用SSTk-ω湍流模型模擬彈性支撐的剛性圓柱體的渦激振動問題。需要指出的是，在流體域的仿真計算過程中，無量綱時間步U∞Δt/D=0.005 。同時以二階隱式時間離散和二階迎風空間離散格式對控制方程進行離散，并采用Simplec算法耦合求解平均速度場合壓力場；而針對結構域的計算，則是通過對Fluent進行二次開發完成，即通過Fluent的用戶自定義函數 (UDF), 在C語言環境下編寫求解結構動力學運動方程的程序，來求解圓柱的振動響應；同時，結構變形引起的流場網格變化通過Fluent的動網格模型完成。

2.1　圓柱振動響應幅值比A*及頻率比f*

在對圓柱渦激振動的研究中，圓柱的響應振幅比和頻率比始終是最關注的2個因素。圖3給出了用SSTk-ω湍流模型計算得到的振幅比和頻率比分別與實驗結果的對比(本文實驗結果數值全部以空心圓點表示)。

從圖3(a)可以看出，SSTk-ω模型成功計算出3種響應分支。在折合速度U*=3.5的時候，原始分支向上端分支轉變；當折合速度U*=5.2時，出現下端分支。在上端分支中振幅達到最大值0.747，而在下端分支中振幅最大值0.642 ；當折合速度U*=11.0的時候，圓柱體的響應位移又回落到一個很小的數值，表明此時圓柱已進入了“渦脫區域”。

圖3　不同折合速度下的最大振幅及振動頻率(m*=2.4)Fig.3　Amplitude and frequency responses as a function of　reduced velocity (m*=2.4)

圖3(b)給出了不同折速下對應的響應頻率比，在鎖定區內部，圓柱體的實際振動頻率fex與固定圓柱體的泄渦頻率fst分離，同時頻率比穩定在1.15附近，而在解鎖區，圓柱體的實際振動頻率fex與固定圓柱體的渦脫頻率fst相同，這與前人的實驗結果大致相同。

2.2　圓柱尾渦渦脫模式

Govardhan和Williamson的研究結果[11]表明，圓柱渦激振動的不同響應分支對應的尾渦泄渦模式也各不相同：初始分支對應的泄渦是2S模式(即在一個渦脫周期內有2個單獨的尾渦形成)，而在上端分支和下端分支對應的渦脫都是2P模式(即在一個渦脫周期內有2對尾渦形成)，不同的是上端分支中2對渦的渦強大致相等，而下端分支中第2個渦強度較第1個渦強度更小。圖4所示為U*=3.0，4.8和7.3時一個泄渦周期內渦量變化的等值線圖。

折速U*=3.0時，位于初始分支。尾渦泄放并不同于經典的卡門渦街(2S模式)，而更像是P+S模式(每個渦脫周期中從圓柱一側泄放出2對渦，而在另一側泄放出1個單渦)。

折速U*=4.8時，位于上端分支。從圖4(b)可知,仿真模擬的泄渦模式是一個典型的2P渦脫模式，并與Williamson的實驗結果完全吻合。

圖4　一個渦脫周期內的渦量等值線圖Fig.4　Vorticity magnitude contours within a vortex　shedding period

但在下端分支，SSTk-ω湍流模型的模擬結果似乎也不是2P模式，這說明泄渦模式還與圓柱的振動形式有關。采用SSTk-ω湍流模型成功模擬出2P模式，這是由于SSTk-ω湍流模型的湍流粘度計算中考慮到了湍流切應力的傳播，而從可以有效的模擬流動分離。

2.3　流體力系數和相位角

如前所述，一旦圓柱的響應振幅比A*和頻率比f*的值已知，就可以通過相關方程求得總升力系數CL，漩渦力系數Cvortex，以及對應的總相位角φ和渦相角φvortex。圖5給出了相應的流體力系數和對應的相位角隨折合速度的變化關系。

圖5　不同折合速度下的升力系數和相位角(m*=2.4)Fig.5　Force and phase angle variations with U*(m*=2.4)

從圖中可看出，在初始分支中渦相位角φvortex接近0°，表明漩渦力系數Cvortex和附加質量力系數Cpotential同相，因此總升力系數CL值很大。在上端分支，盡管渦相位角φvortex接近180°導致Cvortex和Cpotential異相，但由于圓柱振幅很大，因此Cpotential很大，且遠遠超過了Cvortex的值，所以此時CL值也較大。但在下端分支，Cvortex和Cpotential大小相當且異相，故升力系數CL變得相當小。此外，在初始分支過渡到上端分支處，由于φvortex發生突變，因此泄渦模式發生變化，由2S模式變為2P模式；而在上端分支過渡到下端分支時，φvortex并沒有發生突變，因此泄渦模式也并未發生變化，依然保持在2P模式，這說明φvortex是否發生突變是評判泄渦模式是否發生變化的標準。

3結語

本文基于雷諾平均納維-斯托克斯(RANS)方程, 研究了低質量比m*=2.4時的圓柱橫向渦激振動問題，采用SSTk-ω湍流模型求解RANS方程，較好地再現了實驗中的很多現象。

1)圓柱橫向振動上端分支最大響應幅值為0.747，略小于實驗值，可能原因是實際實驗模型向二維的簡化,而在頻率響應計算中，圓柱鎖定區的頻率比穩定在1.15附近，與實驗值較為接近。

2)SSTk-ω湍流模型成功模擬出了處于上端分支的2P渦脫模式，這得益于該模型在計算湍流粘度中考慮到了湍流切應力的傳播，能有效模擬流動分離。

3)在初始分支向上端分支的過渡中，渦相角φvortex發生突變，且尾渦的泄渦模式也發生改變，由2S模式變為2P模式；而在上端分支過渡到下端分支時φvortex并沒有發生突變，渦脫模式也未變化，說明φvortex是否發生突變是評判泄渦模式是否發生變化的標準。

從上面的分析可看出，無論是振幅響應和頻率響應的計算，還是對升力系數和相位角的分析，SSTk-ω湍流模型的模擬結果較接近于真實的物理現象，并且采用SSTk-ω湍流模型成功地模擬出了2P模式，這些結論不僅為今后的數值模擬工作提供

借鑒，而且還為進一步的實驗和理論研究提供了參考。

參考文獻：

[1]BEARMAN P W.Circular cylinder wakes and vortex-induced vibrations[J].Journal of Fluids and Structures,2011(27):648-658.

[2]CATALANO P,WANG M,IACCARINO G.Numerical simulation of the flow around a circular cylinder at high Reynolds numbers[J].International Journal of Heat and Fluid Flow,2003(24)：463-469.

[3]GAO D M F,MINGHAM C G.Numerical and experimental investigation of turbulent flow around a vertical circular cylinder[C]//20th International Offshore and Offshore and Polar Engineering Conference Proceedings,2010：639-643.

[4]KHALAK A,WILLIAMSON C H K.Dynamics of a hydr-oelastic cylinder with very low mass and damping[J].Journal of Fluids and Structures,1996(10)：455-472.

[5]WILLIAMSON C H K.Advances in our understanding of vortex dynamics in bluff body wakes[J].Journal of Wind Engineering and Industrial Aerodynamics,1997:3-32.

[6]KHALAK A,WILLIAMSON C H K.Motions, forces and mode transitions in vortex dynamics in bluff body wakes[J].Journal of Wind Engineering and Industrial Aerodynamics,1999(13)：813-851.

[7]GOVARDHAN R,WILLIAMSON C H K.Modes of vortex formation and frequency response of a freely vibrating cylinder[J].Journal of Fluid Mechanics,2000：85-130.

[8]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 and Structures,2007(23):23-37.

[9]PARKINSON G V.Mathematical models of flow-induced vibrations of bluff bodies[J].Flow-Induced Structural Vibrations,1974:81-127.

[10]LIGHTHILL J.Fundamentals concerning wave loading on offshore structures[J].Journal of Fluid Mechanics,1986：667-681.

[11]GOVARDHAN R,WILLIAMSON C H K.Mean and fluctuating velocity fields in the wake of a freely-vibrating cylinder[J].Journal of Fluid and Structures,2000(15)：85-130.

Numerical simulation of vortex-induced vibration of a two-dimensional circularcylinder based on the SSTk-ωturbulent model

LI Jun1，LI Wei1,2

(Department of Naval Architecture and Ocean Engineering, Huazhong University of Science and

Technology,Wuhan 430074,China)

Abstract：Due to the great damage to widely utilized flexible structures in ocean engineering, vortex-induced vibration (VIV) of such long flexible marine structures is still a hot issue that needs more theoretical research, and CFD techniques become gradually indispensable to study the VIV problem. In this paper, two-dimensional Reynolds-averaged Navier-Stokes (RANS) equations are adopted to investigate transverse VIV of elastically mounted rigid cylinder with low mass-damping, and SST k-ω turbulent model is applied to solve the RANS equations. By comparing the cylinder displacement response, frequency response, vortex shedding modes of three different response branches, and other hydrodynamic coefficients with the previous research in detail, the numerical results indicate that SST k-ω model is invalid and accurate for VIV of the elastically mounted rigid cylinder. This investigation provides theoretical evidence for the numerical simulation of VIV of marine riser in engineering application.

Key words：vortex-induced vibration (VIV); SST k-ω turbulent model;riser

作者簡介：李駿( 1990 - ) ，男，碩士研究生，主要從事海洋結構物設計與制造工作。

收稿日期：2014-03-05； 修回日期： 2014-07-07

文章編號：1672-7649(2015)02-0030-05

doi：10.3404/j.issn.1672-7649.2015.02.006

中圖分類號：U661.44

文獻標識碼：A