考慮渦泄模式的低質(zhì)量比圓柱渦激運(yùn)動(dòng)數(shù)值模擬

孫洪源，何棟梁，林?；?，常 爽，于福臨

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

0 引 言

隨著海洋浮式結(jié)構(gòu)的廣泛應(yīng)用，Spar 海洋平臺(tái)、Spar 型浮式風(fēng)電基礎(chǔ)等所包含的立柱式結(jié)構(gòu)，在一定流速下發(fā)生渦激振動(dòng)的現(xiàn)象。為了區(qū)別，將這種長(zhǎng)徑比比較低、漂浮于水面的剛性結(jié)構(gòu)物的渦激振動(dòng)現(xiàn)象稱為渦激運(yùn)動(dòng)。渦激運(yùn)動(dòng)涉及諸多科目，如流體力學(xué)、計(jì)算流體力學(xué)、振動(dòng)學(xué)、結(jié)構(gòu)力學(xué)和統(tǒng)計(jì)學(xué)等[1]。

海洋結(jié)構(gòu)物渦激運(yùn)動(dòng)（Vortex induced motion，VIM）這一熱點(diǎn)問(wèn)題越來(lái)越受到深海工程技術(shù)人員和研究人員的關(guān)注，它是典型鈍體繞流中升力和拖曳力所導(dǎo)致的直接后果[2]。目前國(guó)際上已有一些學(xué)者對(duì)渦激運(yùn)動(dòng)進(jìn)行了研究，其中包含對(duì)海洋平臺(tái)等結(jié)構(gòu)物進(jìn)行渦激運(yùn)動(dòng)集體預(yù)報(bào)等[3-5]。深水Spar 平 臺(tái)、浮筒等均具有特殊的深吃水立柱式結(jié)構(gòu)，當(dāng)水流經(jīng)過(guò)時(shí)，圓柱體后方發(fā)生周期性渦旋分離，致使結(jié)構(gòu)受到沿流向的拖曳力以及垂直于流向的升力[6]。

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

1 數(shù)值計(jì)算模型

1.1 渦激運(yùn)動(dòng)的主要研究方法

隨著計(jì)算機(jī)技術(shù)的發(fā)展，計(jì)算流體力學(xué)（CFD）方法逐漸應(yīng)用到研究浮式海洋平臺(tái)渦激振動(dòng)的問(wèn)題中[8]，采用數(shù)值模擬方法彌補(bǔ)了模型試驗(yàn)對(duì)于一些多參數(shù)的敏感性分析和流場(chǎng)捕捉方面成本太高的不足。為了保證有限元仿真的可靠性，在建模時(shí)需要考慮模型的力學(xué)性質(zhì)，采用合理的單元類型和應(yīng)變方程去模擬實(shí)際的結(jié)構(gòu)[9]。

渦激運(yùn)動(dòng)的研究在于規(guī)定流體的參數(shù)，因?yàn)樗c許多參數(shù)有關(guān)，這些參數(shù)都會(huì)直接或者間接影響這個(gè)流體力。

1.2 渦激運(yùn)動(dòng)的相關(guān)參數(shù)

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

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

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

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

5）質(zhì)量比：m?=m/πρD2L/4。

式中：Uc為水來(lái)流速度；D為圓柱動(dòng)力直徑；V為水粘性系數(shù)；L為圓柱高度；fv為圓柱渦泄模式頻率；TS WAY為圓柱橫搖周期；A為 圓柱運(yùn)動(dòng)幅值；ρ為流體密度。

1.3 渦激運(yùn)動(dòng)求解過(guò)程

本文將彈性支撐圓柱的渦激運(yùn)動(dòng)系統(tǒng)簡(jiǎn)化為質(zhì)量-彈簧-阻尼系統(tǒng)，如圖1 所示。

圖1 彈性支撐圓柱渦激運(yùn)動(dòng)模型Fig. 1 Vortex-induced motions model of elastically supported cylinder

均勻流中，圓柱渦激運(yùn)動(dòng)控制方程為：

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

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

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

表1 網(wǎng)格測(cè)試結(jié)果Tab. 1 Grid test results

2 實(shí)驗(yàn)結(jié)果與分析

2.1 計(jì)算工況

為了驗(yàn)證數(shù)值模擬方法的可靠性，選擇質(zhì)量比為2.6 的工況與Jauvtis 和Williamson[10]實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。后期對(duì)質(zhì)量比2.6 的圓柱與質(zhì)量比為1 的浮式圓柱進(jìn)行渦激運(yùn)動(dòng)特性研究比較，具體工況信息如表2 所示。

表2 計(jì)算工況表Tab. 2 Calculation condition table

2.2 渦泄模式的比較

圖2 為不同響應(yīng)分支所對(duì)應(yīng)的渦旋泄放模式：左側(cè)方為數(shù)值模擬結(jié)果，右側(cè)方為Jauvtis 和Williamson的實(shí)驗(yàn)結(jié)果。

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

在上端分支，當(dāng)橫向振幅達(dá)到最大值時(shí)，旋渦泄放出現(xiàn)了 2T（two triplets of vortices per cycle）模式，如圖2（c）所示。

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

圖3 給出了 2T模式在1 個(gè)周期的旋渦圖，當(dāng)圓柱處于橫流向最大位移時(shí)，速度最小但加速度最大，在這種較大加速度作用下產(chǎn)生了3 個(gè)渦。其中一個(gè)渦的旋轉(zhuǎn)方向與另外2 個(gè)相反，且這個(gè)渦的強(qiáng)度小于其余2 個(gè)，數(shù)值計(jì)算較好模擬了這一特有現(xiàn)象。

圖2 不同響應(yīng)分支對(duì)應(yīng)的渦泄模式Fig. 2 Vortex release modes corresponding to different response branches

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

2.3 幅值比較

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

圖4 不同約化速度下無(wú)量綱振幅與實(shí)驗(yàn)結(jié)果的對(duì)比Fig. 4 Comparison of dimensionless amplitude and experimental results at different reduction speeds

2.4 軌跡比較

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

2.5 考慮渦泄模式的浮式圓柱渦激運(yùn)動(dòng)特性

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

Ur=7.11 處于渦激運(yùn)動(dòng)的上端分支，與圖2（c）一樣表現(xiàn)出 2T模式。順流向渦激運(yùn)動(dòng)頻率為橫流向渦激運(yùn)動(dòng)頻率的2 倍關(guān)系。

圖5 不同約化速度下圓柱運(yùn)動(dòng)軌跡Fig. 5 Cylindrical trajectory under different reduced speeds

圖6 浮式圓柱渦激運(yùn)動(dòng)在 Ur =4.47 時(shí)一個(gè)周期渦泄圖（m*=1，）Fig. 6 A periodic vortex release diagram of floating cylindrical vortex induced motion at Ur=4.47（m*=1，）

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

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

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

表3 浮式圓柱X/D、Y/D 渦激運(yùn)動(dòng)響應(yīng)頻譜圖（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 渦激運(yùn)動(dòng)渦泄圖Fig. 7 Vortex-exhaust diagram of floating cylinder X/D and Y/D vortex-induced motion

3 結(jié) 語(yǔ)

本文應(yīng)用CFD 數(shù)值模擬方法對(duì)二維圓柱在不同流速下的渦激運(yùn)動(dòng)問(wèn)題進(jìn)行研究。

1）通過(guò)與Jauvtis&Williamson 在2004 年的經(jīng)典實(shí)驗(yàn)結(jié)果的對(duì)比，數(shù)值計(jì)算可以較好模擬圓柱渦激運(yùn)動(dòng)，但在最大振幅方面有所低估。初始分支、上端分支和下端分支的運(yùn)動(dòng)軌跡與實(shí)驗(yàn)結(jié)果基本相似。數(shù)值模擬再現(xiàn)了實(shí)驗(yàn)中的S S， 2S， 2T和 2P渦泄模式。

2）浮式圓柱渦激運(yùn)動(dòng)的初始分支，表現(xiàn)為2 種不同的渦泄模式，分別為P+S模式和 2P模式。而質(zhì)量比為2.6 的圓柱，在初始分支中渦旋泄放表現(xiàn)為S S模式和 2S模式。