李　駿，李　威

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

基于SSTk-ω湍流模型的二維圓柱渦激振動(dòng)數(shù)值仿真計(jì)算

李駿，李威

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

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

關(guān)鍵詞：渦激振動(dòng)；SST k-ω湍流模型；立管

0引言

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

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

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

本文參照Khalak 和 Williamson 的實(shí)驗(yàn)?zāi)Ｐ停x圓柱的質(zhì)量比m*=2.4 ，質(zhì)量阻尼比 m*ξ=0.013；折合速度U*=Uinlet/(fnD)從2.0增加13.9；相應(yīng)的雷諾數(shù)從1 700增加到11 600。采用SSTk-ω湍流模型研究彈性支撐剛性圓柱在自由來流下的橫向渦激振動(dòng)問題，為評(píng)判該模型在這一研究的正確性，所有數(shù)值仿真結(jié)果都將與Khalak 和 Williamson 的相關(guān)實(shí)驗(yàn)結(jié)果作對(duì)比。

1計(jì)算模型

1.1　控制方程

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

(1)

(2)

(3)

式中：渦流粘度μt由湍流模型給出；δij為克羅內(nèi)克函數(shù)；k為湍流動(dòng)能。

(4)

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

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

因此，圓柱在其橫向方向的無量綱運(yùn)動(dòng)方程為：

(5)

式中：U*為折合速度；ζ為系統(tǒng)結(jié)構(gòu)的阻尼比；m*為質(zhì)量比；y為無量綱的橫向位移；CL為橫向的升力系數(shù)。

比較兩組治療前后呼吸困難程度評(píng)分（MMRC）、6 min步行距離（6MWD）、肺功能狀態(tài)（FEV1預(yù)計(jì)值）等。

當(dāng)圓柱渦激振動(dòng)進(jìn)入“鎖定區(qū)域”，尾渦的渦脫頻率與圓柱振動(dòng)頻率完全相同，此時(shí)圓柱的位移為：

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分別為升力系數(shù)的振幅及它在速度方向與加速度方向上的分量；ωex為圓柱振動(dòng)圓頻率；φ為升力系數(shù)和響應(yīng)位移之間的相位角。

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

(11)

(12)

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

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

Cvortex=CL0-Cpotential，

(13)

其中，

(14)

類似地，相應(yīng)的無量綱運(yùn)動(dòng)方程還可以表示為：

(15)

式中：ma為圓柱的附加質(zhì)量；ζ′為包含了附加質(zhì)量系統(tǒng)結(jié)構(gòu)的阻尼比；φvortex為旋渦升力系數(shù)與響應(yīng)位移之間的相位角，可通過之前提到的相似方法求解得到。

采用Newmar k-β方法求解圓柱的橫向振動(dòng)響應(yīng)位移，同時(shí)通過求解二維的RANS方程獲得相應(yīng)的流體力系數(shù)。

1.2　網(wǎng)格劃分及邊界條件

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

圖2　計(jì)算域和網(wǎng)格模型Fig.2　Computational field and grid model

2結(jié)果與分析

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

2.1　圓柱振動(dòng)響應(yīng)幅值比A*及頻率比f*

在對(duì)圓柱渦激振動(dòng)的研究中，圓柱的響應(yīng)振幅比和頻率比始終是最關(guān)注的2個(gè)因素。圖3給出了用SSTk-ω湍流模型計(jì)算得到的振幅比和頻率比分別與實(shí)驗(yàn)結(jié)果的對(duì)比(本文實(shí)驗(yàn)結(jié)果數(shù)值全部以空心圓點(diǎn)表示)。

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

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

圖3(b)給出了不同折速下對(duì)應(yīng)的響應(yīng)頻率比，在鎖定區(qū)內(nèi)部，圓柱體的實(shí)際振動(dòng)頻率fex與固定圓柱體的泄渦頻率fst分離，同時(shí)頻率比穩(wěn)定在1.15附近，而在解鎖區(qū)，圓柱體的實(shí)際振動(dòng)頻率fex與固定圓柱體的渦脫頻率fst相同，這與前人的實(shí)驗(yàn)結(jié)果大致相同。

2.2　圓柱尾渦渦脫模式

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

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

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

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

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

2.3　流體力系數(shù)和相位角

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

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

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

3結(jié)語

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

1)圓柱橫向振動(dòng)上端分支最大響應(yīng)幅值為0.747，略小于實(shí)驗(yàn)值，可能原因是實(shí)際實(shí)驗(yàn)?zāi)Ｐ拖蚨S的簡化,而在頻率響應(yīng)計(jì)算中，圓柱鎖定區(qū)的頻率比穩(wěn)定在1.15附近，與實(shí)驗(yàn)值較為接近。

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

3)在初始分支向上端分支的過渡中，渦相角φvortex發(fā)生突變，且尾渦的泄渦模式也發(fā)生改變，由2S模式變?yōu)?P模式；而在上端分支過渡到下端分支時(shí)φvortex并沒有發(fā)生突變，渦脫模式也未變化，說明φvortex是否發(fā)生突變是評(píng)判泄渦模式是否發(fā)生變化的標(biāo)準(zhǔn)。

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

借鑒，而且還為進(jìn)一步的實(shí)驗(yàn)和理論研究提供了參考。

參考文獻(xiàn)：

[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 - ) ，男，碩士研究生，主要從事海洋結(jié)構(gòu)物設(shè)計(jì)與制造工作。

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

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

doi：10.3404/j.issn.1672-7649.2015.02.006

中圖分類號(hào)：U661.44

文獻(xiàn)標(biāo)識(shí)碼：A