基于一種固體區(qū)域迭代算法的圓柱渦激振動數(shù)值計算

章大海，石凡奇，王 君，郝木明

（中國石油大學（華東）化學工程學院，山東青島266580）

基于一種固體區(qū)域迭代算法的圓柱渦激振動數(shù)值計算

章大海，石凡奇，王 君，郝木明

（中國石油大學（華東）化學工程學院，山東青島266580）

利用Fluent平臺的用戶自定義程序（UDF）以及動網(wǎng)格模型，實現(xiàn)了圓柱運動方程的一種迭代求解算法，分別對層流、湍流狀態(tài)下，彈性支承圓柱體在一定約化速度下的渦激響應進行了數(shù)值模擬，探討了不同阻尼比對渦激響應的影響。結(jié)果表明：采用該迭代求解算法對彈性支承圓柱渦激振動的預測結(jié)果較為合理；隨著阻尼比的逐漸增加，初始支振幅、升阻力系數(shù)時程曲線將由多頻率拍振，最終變?yōu)閱我活l率主導的振動，且渦激振幅逐漸減?。怀速|(zhì)量-阻尼比聯(lián)合參數(shù)m*ζ外，阻尼比ζ本身也應作為一個重要的渦激影響參數(shù)單獨進行考量。

渦激振動；流固耦合；阻尼比；數(shù)值模擬

0 引 言

流體流經(jīng)管道時在管體兩側(cè)形成的交替渦脫會誘發(fā)管道的周期性振動，而管道的振動又會擾亂流場渦的脫落，這種流體-固體耦合作用的現(xiàn)象稱為“渦激振動”（Vortex-induced Vibration，VIV）。長輸管道，特別是海洋立管的渦激振動是當今工業(yè)界研究的熱點和難點。立管渦激振動預測水平的提高，能直接減少海洋平臺建造成本。加強渦激振動的研究對石油工業(yè)的意義不言而喻。

由于立管現(xiàn)場全比例實驗極其耗費資源，研究者主要從數(shù)值模擬和縮比尺實驗進行渦激振動研究。相關(guān)綜述可參考Sarpkaya[1]，Williamson[2-3]，Gabbai[4]等。最早且最著名的實驗研究是Feng[5]在英屬哥倫比亞大學風洞所作的自激振動實驗。Feng對質(zhì)量比m*=248的圓柱體在橫向運動的響應進行了研究。其中橫向振幅表現(xiàn)出兩個響應分支，初始分支和下端分支，振幅幅值為0.6D。以Williamson為首的實驗小組[6-9]對低質(zhì)量比m*=O（10）量級的振動系統(tǒng)進行了一系列的實驗研究，其目的在于探究渦激振動的內(nèi)在機理，并為渦激振動的鎖定、遲滯現(xiàn)象找到了合理的解釋。Willamson小組的實驗結(jié)果表明：（1）振動幅值由質(zhì)量-阻尼聯(lián)合參數(shù)控制。當質(zhì)量-阻尼組合參數(shù)m*ζ較大時，呈現(xiàn)初始支和下支兩支響應，而該參數(shù)較小時，呈現(xiàn)三支響應（出現(xiàn)了上支）；（2）旋渦脫落在初始支為2S模態(tài)，在上支和下支為2P模態(tài)，初始支與上支之間存在“遲滯性過渡”，而上支與下支之間為“間隙性轉(zhuǎn)換”；（3）當m*ζ為定值時，鎖定區(qū)域的速度范圍隨m*的增大而減小。Zhao[10]將彈性支撐圓柱體自激振動得到的振動曲線作為受迫振動的輸入數(shù)據(jù)，結(jié)果表明，二者可得到幾乎完全相同的流體力信息和尾渦模式，更為直觀地顯現(xiàn)了兩種振動實驗的聯(lián)系。

受限于人們對湍流的有限認知以及計算機技術(shù)發(fā)展水平，至今仍不能對渦激振動進行精確模擬。目前，彈性支撐圓柱體的渦激振動的CFD模擬主要有以下途徑：基于有限體積法的直接數(shù)值模擬（DNS）、大渦模擬（LES）、雷諾時均方法（RANS）、以及最近提出的離散渦方法（DVM）和格子Boltzmann方法（LBM）。林琳等[11-12]利用RANS方法比較了固定圓柱繞流外流場與振動圓柱外流場的區(qū)別，認為SST k-ω模型的模擬結(jié)果優(yōu)于RNG k-ε模型；何長江[13]采用動網(wǎng)格的動態(tài)鋪層模型和滑移網(wǎng)格模型以及用戶自定義程序，將計算結(jié)構(gòu)響應的Newmark代碼嵌入FLUENT軟件，建立了二維圓柱橫流向渦激振動的計算模型。潘志遠[14]利用RANS方法結(jié)合一定的網(wǎng)格策略模擬了水流當中圓柱體橫向自激振動與受迫振動。黃智勇[15]在其研究基礎(chǔ)上模擬了圓柱的二自由度自激振動。研究表明，在質(zhì)量比低于3.5時，兩自由度圓柱的橫向振幅比單自由度情況大。董婧等[16]采用離散渦方法和四階Runge-Kutta法計算了低雷諾數(shù)條件下圓柱體的渦激振動，分析了決定圓柱振動的影響參數(shù)及尾渦結(jié)構(gòu)與氣動力的關(guān)系，觀察到“拍”和“相位開關(guān)”等現(xiàn)象。陶實[17]應用格子Boltzmann方法對彈性支承圓柱渦激振動進行了數(shù)值模擬，研究Sg數(shù)、質(zhì)量比、頻率比對渦激響應的影響，結(jié)果表明質(zhì)量比對振幅的抑制效果最顯著；且頻率比足夠大時，會呈現(xiàn)“8”字型運動軌跡。

以上模擬方法各有優(yōu)勢，能較為準確地捕捉渦激振幅和頻率響應。但在固體區(qū)域控制方程的處理均無一例外采用Runge-Kutta法或Newmark方法。本文利用Fluent平臺內(nèi)嵌的用戶自定義程序（User-defined Functions,UDF）以及動網(wǎng)格模型，實現(xiàn)了圓柱運動方程的一種迭代求解算法，研究了彈性支承圓柱體在一定約化速度范圍內(nèi)的渦激響應，以及不同阻尼比對渦激響應的影響。計算結(jié)果表明，與現(xiàn)有模擬方法相比，本文采用的迭代求解算法也能對彈性支承圓柱渦激振動做出相同精度的預測；隨著阻尼比的逐漸增加，初始支振幅、升阻力系數(shù)時程將由多頻率拍振，最終演變?yōu)閱我活l率主導的振動，且渦激振幅逐漸減小。因而本文認為，除了Williamson小組[6-9]所建議的質(zhì)量-阻尼比聯(lián)合參數(shù)m*ζ外，阻尼比ζ本身也應作為一個重要的渦激影響參數(shù)單獨進行考量。

1 計算模型

1.1 幾何模型

彈性支承圓柱振動系統(tǒng)簡化為圖1所示的物理模型。整個計算域取為如圖1所示的長方形區(qū)域。圓柱直徑D，圓柱中心距離上，下及左邊界的距離為10D，距離右邊為20D，以保證充分瀉渦；總網(wǎng)格數(shù)約2.0×104個，且為保證圓柱附近流場計算的精度，緊貼圓柱布置隨圓柱運動的結(jié)構(gòu)化網(wǎng)格，厚度約1D，如圖2所示。其中，圓管表面第一層網(wǎng)格的高度為0.01D。流體密度ρ=1 000 kg/m3，粘度ν=1.01× 10-6m2/s。湍流情況采用SST k-ω模型，網(wǎng)格左端為速度入口，右端為壓力出口；上下邊界為對稱邊界條件。由于渦激振動涉及到計算網(wǎng)格的運動，利用動網(wǎng)格的重構(gòu)和彈簧光順方法，實現(xiàn)網(wǎng)格的運動和更新。

圖1 幾何模型示意圖Fig.1 Model sheme of elastically mounted cylinder

圖2 圓柱附近網(wǎng)格Fig.2 Grid near cylinder

1.2 結(jié)構(gòu)動力學方程

對于圖1中的單自由度（single degree of freedom,sdof）彈簧-質(zhì)量-阻尼系統(tǒng)，其控制方程為：

其中：y，m，c，k分別為圓柱位移，圓柱質(zhì)量，結(jié)構(gòu)阻尼系數(shù)和系統(tǒng)剛度系數(shù)。FL（t）為系統(tǒng)所受橫向流體力

進行無量綱化后可得

1.3 計算流程

基于Fluent平臺，利用其用戶自定義函數(shù)（UDF）實現(xiàn)固體和流體控制方程的耦合。首先將（3）式簡化為如下形式：

根據(jù)（7）、（8）式利用迭代思想就可進行結(jié)構(gòu)方程求解程序的編寫。編寫代碼時，給定初始速度為0 m/s，在每時間步只需提取流體力FL，就可對結(jié)構(gòu)運動方程予以求解?；贔luent的UDF內(nèi)嵌機制，將求解代碼嵌入Fluent，實現(xiàn)流體與結(jié)構(gòu)控制方程的耦合；同時利用動網(wǎng)格的彈簧光順（spring-based smoothing）及局部重構(gòu)（local remeshing）模型對網(wǎng)格進行更新，即可對圓柱渦激振動問題進行時域范圍求解。

2 彈性支承圓柱的渦激響應

考慮層流和湍流兩種情況下彈性支承圓柱的渦激響應。為與Williamson小組的文獻保持一致，圖3～8采用約化速度表達式U*=U∞/fnwD進行處理。層流情況模型直徑D=0.01 m，質(zhì)量比m*=2.8，水中的固有頻率fnw=0.35 Hz，約化速度U*=3～16，對應雷諾數(shù)范圍為100～600。湍流情況D=0.038 1 m，質(zhì)量比m*=8.6，fnw=0.4 Hz，約化速度U*=3～12，對應雷諾數(shù)范圍為1 800～8 000，計算結(jié)果如圖3～6所示。

圖3給出了本文的計算結(jié)果與實驗數(shù)據(jù)的對比。湍流情況下，本文得出的最大振幅出現(xiàn)在U*=5.6處，最大振幅值約為A/D=0.52，振幅變化趨勢與Govardhan[7]的實驗結(jié)果相符，但本文模擬出初始支的最大振幅比實驗結(jié)果偏小。產(chǎn)生本文結(jié)果的重要原因可能是本文采用的弱耦合算法。弱耦合算法在一個時間步內(nèi)分開來求解流場和固體區(qū)域的控制方程，流體、固體區(qū)域物理量（流體力、固體位移）通過流固交界面插值計算來進行交換。因此，在耦合界面上的動量是不完全守恒的[19]，這可能是造成本文算例振幅偏大的最重要的因素。值得注意的是，在相同約化速度下，層流狀態(tài)和湍流狀態(tài)下的渦激振幅響應有明顯差別。圖3清楚地表明，湍流狀態(tài)下的振幅響應比層流狀態(tài)大，且鎖定區(qū)域?qū)臒o量綱速度范圍也更大。因此，單獨用約化速度不能作為渦激響應的衡量標準，還應加上標明流場本身性質(zhì)的雷諾數(shù)。圖4給出了約化速度U*=5.6時，對應鎖定狀態(tài)下流體力系數(shù)與圓柱振幅時程。可以看出：在鎖定現(xiàn)象發(fā)生時，圓柱的升力系數(shù)曲線和圓柱振幅曲線的波峰、波谷對應，二者相位一致，表明流體向圓柱輸送能量，出現(xiàn)類似共振的“鎖定現(xiàn)象”。阻力系數(shù)頻率約為升力系數(shù)兩倍。以上結(jié)論均與實驗結(jié)論一致，因而本文采用的迭代算法是有效的。

圖3 振幅曲線[7，18]Fig.3 Amplitude response as a function of reduced velocity

圖4 U*=5.6，振幅及流體力系數(shù)曲線Fig.4 Time traces of cylinder amplitude and fluid forces,U*=5.6

圖5 層流狀態(tài)下圓柱體單個振動周期內(nèi)流場渦量圖，U*=5.6Fig.5 Vorticity magnitude contours within a vortex shedding period,laminar flow,U*=5.6

圖6 湍流狀態(tài)下圓柱體單個振動周期內(nèi)流場渦量圖，U*=7.5Fig.6 Vorticity magnitude contours within a vortex shedding period,turbulent flow,U*=7.5

本文層流狀態(tài)下的模擬捕捉到2S的瀉渦模式，如圖5所示；圖6給出的湍流狀態(tài)下（對應U*=7.5）圓柱一個振動周期內(nèi)渦量等值線圖清晰地捕捉到了2P瀉渦模式。Williamson小組[5-9]實驗結(jié)果表明：低質(zhì)量比彈性支撐圓柱渦激振動中與振動響應初始分支對應的尾渦為2S模式，而與上端分支和下端分支對應的尾渦則為2P模式。本文結(jié)果與之相符，這也驗證了計算方法的有效性。

3 阻尼比對渦激響應的影響

研究層流和湍流情況下阻尼比對渦激響應的影響。層流、湍流工況均取渦激響應初始支，U*=5.6。采用阻尼比為ζ=0，4.6×10-3，9.2×10-3，9.2×10-2，其余實驗參數(shù)與上節(jié)相同。模擬結(jié)果如圖7～8所示。

圖7 層流情況下不同阻尼比圓柱的渦激響應時程及對應頻譜。頻譜圖中兩條虛線分別是系統(tǒng)在水中的自然頻率fnw=0.35和空氣中的自然頻率fn=0.41Fig.7 Amplitude and fluid force response and corresponding spectrum under different damping ratios,laminar case. Dash lines indicate the natural frequencies of the system in water and air respectively

圖8 湍流情況下不同阻尼比圓柱的渦激響應時程及對應頻譜。頻譜圖中兩條虛線分別是系統(tǒng)在水中的自然頻率fnw=0.40和空氣中的自然頻率fn=0.42Fig.8 Amplitude and fluid force response and corresponding spectrum under different damping ratios,turbulent case. Dash lines indicate the natural frequencies of the system in water and air respectively

圖7、圖8表明，阻尼比對渦激響應有顯著影響。在振動響應的初始支，層流情況下，阻尼比ζ=0時，升力系數(shù)、渦激振幅時程表現(xiàn)出多頻率拍振現(xiàn)象，最大振幅在0.50左右；阻尼比ζ=4.6×10-3，9.2× 10-3時，振幅、升力時程表現(xiàn)出兩個頻率主導的類似差拍振動現(xiàn)象。升力系數(shù)和渦激振幅時程的頻譜分析顯示，ζ=4.6×10-3時，二者主要集中在fex=0.439、0.409兩個頻率振動，最大振幅約0.51；阻尼比ζ增大為9.2×10-3時，升力和振幅主要集中在頻率fex=0.439、0.427振動，兩頻率相互靠近，最大振幅0.49；而當阻尼比繼續(xù)增大到9.2×10-2時，升力系數(shù)和振幅集中在單一頻率振動，最大振幅為0.30。湍流情況下的情況與此類似，隨著阻尼比的增大，振動響應越來越明顯地表現(xiàn)為靠近fnw的單頻振動，頻率遠離fnw的振動成分逐漸消失。可以看出，隨著阻尼比的逐漸增加，初始支的振幅、升阻力系數(shù)時程由多頻率拍振，逐漸演變?yōu)閮深l率主導的類似差拍振動，最終變?yōu)閱我活l率主導的振動，且振幅逐漸減小。因而本文認為，除了Williamson小組[6-9]所建議的質(zhì)量-阻尼比聯(lián)合參數(shù)m*ζ外，阻尼比ζ本身也應作為一個重要的渦激影響參數(shù)單獨進行考量。

4 結(jié) 語

本文跟蹤了圓柱繞流現(xiàn)象的最近研究成果，利用Fluent平臺內(nèi)嵌的用戶自定義程序（UDF）以及動網(wǎng)格模型，實現(xiàn)了圓柱運動方程的另一種迭代求解算法，研究了彈性支承圓柱體在一定約化速度范圍內(nèi)的渦激響應，以及不同阻尼比對渦激響應的影響。計算結(jié)果表明：

（1）與現(xiàn)有固體區(qū)域求解方法相比，本文采用的迭代求解算法也能對彈性支承圓柱渦激振動做出相同精度的預測；

（2）隨著阻尼比的逐漸增加，渦激響應初始支振幅、升阻力系數(shù)時程將由多頻率拍振，最終演變?yōu)閱我活l率主導的振動，且渦激振幅逐漸減小。因而本文認為，除了Williamson小組[6-9]所建議的質(zhì)量-阻尼比聯(lián)合參數(shù)m*ζ外，阻尼比ζ本身也應作為一個重要的渦激影響參數(shù)單獨進行考量。

[1]Sarpkaya T.A critical review of the intrinsic nature of vortex-induced vibrations[J].Journal of Fluids and Structures,2004, 19(4):389-447.

[2]Williamson C H K,Govardhan R.Vortex-Induced Vibrations[J].Annual Review of Fluid Mechanics,2004,36(1):413-455.

[3]Williamson C H K,Govardhan R.Govardhan.A brief review of recent results in vortex-induced vibrations[J].Journal of Wind Engineering and Industrial Aerodynamics,2008,96(6-7):713-735.

[4]Gabbai R D,Benaroya H.An overview of modeling and experiments of vortex-induced vibration of circular cylinders[J]. Journal of Sound and Vibration,2005,282(3-5):575-616.

[5]Feng C C.The measurement of vortex induced effects in flow past stationary and oscillating circular and dissection cylinders [D].Canada:University of British Columbia,1968.

[6]Morse T L,Williamson C H K.The effect of Reynolds number on the critical mass phenomenon in vortex-induced vibration[J].Physics of Fluids,2009,21(4):045105.

[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,420:85-130.

[8]Morse T L,Williamson C H K.The effect of Reynolds number on the critical mass phenomenon in vortex-induced vibration[J].Phisics of Fluids,2009,21(045105):1-8

[9]Khalak A,Williamson C H K.Dynamics of a hydroelastic cylinder with very low mass and damping[J].Journal of Fluids and Structures,1996,10:455-472.

[10]Zhao J,Nemes A,Lo Jacono D,et al.Comparison of fluid forces and wake modes between free vibration and tracking motions of a circular cylinder[C]//18th Australian Fluid Mechanics Conference.Launceton,Australia,2012.

[11]林 琳,王言英.自激振動圓柱體湍流場及發(fā)展變化的研究[J].振動與沖擊,2013,32(14):164-173. Lin Lin,Wang Yanying.Computation of a turbulence flow and its development around a self-induced vibrating circular cylinder[J].Journal of Vibration and Shock,2013,32(14):164-173.

[12]林 琳,王言英.不同湍流模型下圓柱渦激振動的計算比較[J].船舶力學,2013,17(10):1115-1125. Lin Lin,Wang Yanying.Comparison between different turbulence models on vortex induced vibration of circular cylinder [J].Journal of Ship Mechanics,2013,17(10):1115-1125.

[13]何長江,段忠東.二維圓柱渦激振動的數(shù)值模擬[J].海洋工程,2008,26(1):57-63. He Changjiang,Duan Zhongdong.Numerical simulation of vortex induced vibration on 2D circular cylinders[J].Ocean Engineering,2008,26(1):57-63.

[14]潘志遠.海洋立管渦激振動機理與預報方法研究[D].上海:上海交通大學,2005. Pan Zhiyuan.Vortex induced vibration of marine riser and response prediction[D].Shanghai:Shanghai Jiao Tong University,2005.

[15]黃智勇,潘志遠,崔維成.兩向自由度低質(zhì)量比圓柱體渦激振動的數(shù)值計算[J].船舶力學,2007,11(1):1-9. Huang Zhiyong,Pan Zhiyuan,Cui Weicheng.Numerical simulation of VIV of a circular cylinder with two degrees of freedom and low mass-ratio[J].Journal of Ship Mechanics,2007,11(1):1-9.

[16]董 婧,宗 智,李章銳,等.兩自由度運動圓柱繞流的離散渦方法模擬[J].船舶力學,2012,16(1-2):9-20. Dong Jing,Zong Zhi,Li Zhangrui,et al.Numerical simulation of flow around a cylinder of two degrees of freedom motion using the discrete vortex method[J].Journal of Ship Mechanics,2012,16(1-2):9-20.

[17]陶 實,周 力,宗 智.基于格子Boltzmann方法的圓柱兩自由度渦激振動的研究[J].水動力學研究與進展,2013, 28(2):167-175. Tao Shi,Zhou Li,Zong Zhi.Vortex-induced vibration of circular cylinder with two degrees of freedom based on Lattice Boltzmann simulation[J].Chinese Journal of Hydrodynamics,2013,28(2):167-175.

[18]Anagnostoplutos P,Bearman P W.Response characteristics of a vortex-excited cylinder at low reynolds numbers[J].Journal of Fluids and Structures,1992,6:39-50.

[19]陳 鋒,王春江,周 岱.流固耦合理論與算法評述[J].空間結(jié)構(gòu),2012,18(4):55-63. Chen Feng,Wang Chunjiang,Zhou Dai.Review of theory and numerical methods of fluid-structure interaction[J].Spatial Structures,2012,18(4):55-63.

Numerical simulation of VIV of a cylinder using an iteration method in calculating cylinder motion

ZHANG Da-hai,SHI Fan-qi,WANG Jun,HAO Mu-ming
(College of Chemical Engineering,China University of Petroleum,Qingdao 266580,China)

An iteration method is adopted in numerical simulation,to solve the motion of an elastically mounted cylinder of VIV based on Fluent UDF codes.Cases of laminar flow and turbulent flow are conducted respectively,within a range of reduced velocity,and cases with different damping ratios are analyzed. Comparison with experimental data verifies the validity of the method.The resuts show that the lift coefficient and the amplitude of the cylinder will transform from multi-frequency oscillation to single frequency oscillation,and the maximum amplitude will decrease as damping ratio increases.Therefore damping ratio ζ should be regarded as an independent parameter in addition to the well known mass-damping parameter m*ζ.

Vortex-Induced Vibration;solid-fluid interaction;damping ratio;numerical simulation

O237

：Adoi:10.3969/j.issn.1007-7294.2016.04.006

1007-7294（2016）04-0430-09

2015-07-13

山東省自然科學基金項目（ZR2012EEQ011）

章大海（1978-），男，博士，副教授，E-mail:dhzhang@upc.edu.cn；石凡奇（1989-），男，碩士研究生，E-mail:andy20110920@hotmail.com。