振蕩流作用下圓柱體結構的水動力特性研究

楊沖霄，袁昱超，薛鴻祥，唐文勇

(上海交通大學 海洋工程國家重點實驗室，上海 200240)

隨著我國海洋油氣資源開發逐漸向深海邁進，海洋工程結構物的安全性能越來越得到重視。海洋立管系統作為連接水面浮式裝置和海底設備的導管，在復雜的載荷作用下，極易出現碰撞、振動及疲勞等破壞形式，是深海工程裝備中最薄弱的環節[1]。在來流作用下，立管的交替泄渦現象產生振蕩的流體力，引起結構的渦激振動[2]。早期國外學者對該問題開展了大量的試驗研究[3-5]。其中Gopalkrishnan[3]對Re為10 000下的圓柱受迫振動進行了模型試驗，并首次構建了受迫振動流體力系數庫。在數值方法中，隨著計算機技術的發展，計算流體力學方法(computational fluid dynamics， 簡稱CFD)應運而生，為了提高計算效率，將立管簡化為二維圓柱是較為常見的研究手段。Muhamad和Krish[6]使用RANS k-ω湍流模型對Re=10 000來流下的圓柱進行單自由度受迫振動數值模擬，結果表明漩渦尾跡對圓柱振動頻率較為敏感，當振動頻率接近Strouhal頻率時，會發生漩渦的相位變換。王亞非[7]數值模擬了雙自由度彈性支撐圓柱體的渦激振動問題，結果表明振蕩流擴大了自激振動的鎖定范圍，并且使得順流向振動幅值大大增加，甚至超過了橫向所能激發的最大幅值。付博文等[8]基于切片理論，通過使用徑向基函數法作為OpenFOAM中的動網格策略，模擬了長細比為1 000的柔性立管在橫流向和順流向的振動，數值模擬重現了高階主控模態及主控模態的頻繁變換等大長細比柔性立管的渦激振動特性。王凱鵬[9]通過對均勻來流和剪切來流兩種來流形式下的靜止圓柱、橫向受迫振動圓柱繞流及圓柱渦激振動問題展開數值模擬，首次系統地分析了來流形式對圓柱繞流和圓柱渦激振動問題的影響，并指出來流形式是很重要的一個影響因素。總體來說目前的研究主要集中在定常流。

實際生產過程中，海洋浮式結構物在風、浪、流的聯合作用下會帶動立管在水中往復運動，此時立管的遭遇流場可等效為振蕩流場[10]。相比定常流，振蕩流場引起的圓柱繞流問題更為復雜，需要更加深入的研究，目前已有學者獲得了相關成果。Zhao等[11]模擬了振蕩流和均勻流共同作用下圓柱的渦激振動，探討了流動比例a對圓柱的響應，并發現了在a=0.8，Vr=7時，漩渦在一個振蕩流周期中經歷2S、2P和2T三種泄渦模式。鄧躍[12]對低雷諾數時單自由度彈性支撐圓柱在均勻流和振蕩流共同作用下的受迫和自激振動進行了數值模擬，結果發現當有振蕩流參與時，結構的鎖定范圍和振動幅值等都有明顯的變化。鄧迪等[13]采用OpenFOAM開源軟件對在靜水中做周期性振蕩運動的二維剛性圓柱渦激振動進行數值模擬，發現圓柱的橫向振動加劇了升力系數的變化，使得泄渦方向和圓柱表面漩渦分離點的位置發生了明顯變化。

目前對于振蕩流的研究主要集中在低雷諾數或小KC數工況，而已有認知表明，不同雷諾數及KC數對圓柱繞流特性影響顯著[14-16]。文中重點研究高雷諾數條件并兼顧大KC數工況。首先，對均勻流下雷諾數10 000時的圓柱受迫振動進行了數值模擬，通過與試驗結果進行對比，驗證了基于CFD方法復現受迫振動試驗并研究流體力系數的可行性。之后開展振蕩流下圓柱受迫振動數值仿真，分析不同KC數下圓柱的水動力特性及漩渦形態，并歸納出大KC數和小KC數下升力系數和漩渦發放特點。

1 數值分析方法

1.1 控制方程

考慮不可壓縮流場，采用雷諾時均方法(Reynolds Averaged Navier-Stokes， 簡稱RANS)，結合剪切應力輸運模型(k-ω SST)求解N-S方程，控制方程：

(1)

(2)

k-ω SST湍流模型的輸運方程：

(3)

(4)

式中：μt為渦黏性；Sij為平均速度應變率張量；τtij表示雷諾應力的渦黏性模型；σk、β*、σω、σω2均為經驗參數；Pω代表交錯擴散項。F1為混合函數，在近壁處采用Wilcox k-ω模型，邊界層邊緣和自由剪切層采用k-ε模型，中間通過F1函數實現過渡。

1.2 水動力系數識別方法

在Gopalkrishnan[3]的試驗中，圓柱進行受迫振動，位移函數：

y(t)=Asin(2πf0t)

(5)

式中：A為振動幅值；f0為振動頻率；t為運動時間。

由此得到圓柱運動的速度函數為：

v(t)=2πf0Acos(2πf0t)

(6)

當圓柱以式(5)進行振動時，垂直來流方向的升力可表示為:

L=L0sin(2πf0t+φ0)+Lssin(2πfst+φs)

(7)

式中：L0、Ls為升力幅值；φ0、φs為相位角；fs代表Strouhal頻率，其中下標為0表示與圓柱振動頻率相關，下標為s表示與Strouhal頻率相關。

當圓柱處于鎖定狀態時，Strouhal頻率消失，并且圓柱并未以頻率fs在振動，即該成分不參與到流體與結構的能量傳遞中，故在實際分析時，升力得到簡化：

L=L0sin(2πf0t+φ0)

(8)

升力幅值及相位角可通過對時歷曲線進行傅里葉擬合得到：

(9)

計算可得參數：

(10)

(11)

(12)

升力幅值和相位角可表示為：

(13)

(14)

升力系數由下式計算得到：

(15)

式中：l為圓柱高度，在二維計算中取單位高度；D為圓柱直徑；U為遭遇流速。

升力系數還可進一步分解得到激勵力系數和慣性力系數。升力系數中和圓柱運動速度同相位的部分定義為激勵力系數：

CL_V0=CL0sinφ0

(16)

激勵力系數為正代表能量從流體輸入到圓柱結構中，此時易發生圓柱的渦激振動，激勵力系數為負代表能量從結構輸出到流體中。

升力系數中和圓柱運動加速度同相位的部分定義為慣性力系數：

CL_A0=CL0(-cosφ0)

(17)

無因次振幅定義為：

A*=A/D

(18)

無因次頻率定義為：

f*=f0U/D

(19)

振蕩流定義如下：

u(t)=2πfBsin(2πft)

(20)

式中：B為振蕩流幅值；f為振蕩流頻率。

在均勻流中，采用雷諾數定義流體特征：

Re=UD/ν

(21)

振蕩流中，引入柯萊根—卡彭特數(Keulegan-Carpenter number，KC)表征來流：

KC=Umax/fD

(22)

Umax=2πfB

(23)

式中：Umax為振蕩流最大速度。

1.3 計算網格和邊界條件

計算域如圖1所示。X軸平行來流方向，Y軸垂直來流方向。圓柱直徑D=0.025 4 m，上下邊距圓柱中心8D，水平方向50D。由于圓柱周圍流場變化較為劇烈，在圓柱中心3.5D范圍內進行加密，如圖2所示。左端為速度入口，均勻流流速0.4 m/s；上下兩端為對稱邊界；右端為自由出流；圓柱采用無滑移固壁條件；壓力—速度耦合采用SIMPLE算法；時間項采用二階隱式積分方法；對流項采用二階迎風離散格式。

圖1 計算域Fig. 1 Computational domain

圖2 局部加密Fig. 2 Partial encryption

2 均勻流工況試驗驗證

2.1 網格收斂性驗證

首先對定常流速0.4 m/s、Re=10 000下的靜止圓柱進行了數值模擬，給出了普通網格、加密網格及其他學者計算結果的對比，如表1所示。

表1 兩套網格參數Tab. 1 Two sets of grid parameters

根據對比可知，兩套網格的計算結果均與試驗數據對應較好，在綜合考慮計算精度和計算成本的前提下，選擇普通網格完成振蕩流流動計算。

2.2 流體力系數對比

計算了均勻流振幅比為0.3時各無因次頻率下的流體力系數，并與Gopalkrishnan[3]的試驗數據進行了對比，升力系數、激勵力系數、慣性力系數及相位角結果如圖3所示。

圖3 流體力系數對比Fig. 3 Comparison of fluid force coefficient

由圖3可知，隨著無因次頻率的增加，升力系數幅值呈現逐漸增加的趨勢，文中模擬結果在趨勢和數值上和試驗結果吻合較好。激勵力系數變化較為復雜，在無因次頻率0.17處出現正峰值，在0.21處出現負峰值，模擬結果基本可復現這一現象，但數值上略有差距。慣性力系數變化與無因次頻率呈現負相關關系，文中CFD模擬結果較好的吻合了試驗數據。相位角在低頻率處為負值，在頻率0.15附近發生相位突變現象，之后隨著無因次頻率增加逐漸趨于0，CFD模擬結果也基本符合這些特征。但在低頻率處，文中模擬結果與試驗存在一定誤差，這可能是因為Gopalkrishnan所進行的試驗結果包含了復雜的三維效應，故二維模型暫未全面反映三維試驗結果。文獻[18]針對二維數值模擬與試驗誤差也提出了相似的分析。

綜上所述，文中所采用的網格劃分和CFD設置方法復現了圓柱受迫振動的試驗結果，在水動力系數模擬上具有較高的準確性，可運用到振蕩流中開展研究。

2.3 均勻流漩渦形態分析

均勻流工況漩渦發放形態如圖4所示，選取一個脫落周期的漩渦發放結果。圖4(a)時，圓柱向上方運動，尾渦在圓柱下流向產生，并逐漸向后延伸。當圓柱向下運動到圖4(b)位置時，漩渦從圓柱的下尾渦末端脫落，此時上尾渦的末端也分離出了即將脫落的漩渦。圖4(c)時，漩渦從上尾渦末端分離。脫落周期持續0.37 s，與圓柱此時的振動周期是吻合的。可以看到，漩渦在一個周期內呈現上下交替脫落現象，為典型的2S脫落模式。

圖4 均勻流下圓柱漩渦脫落(A*=0.5，f*=0.17)Fig. 4 Vortex shedding under uniform flow(A*=0.5，f*=0.17)

3 振蕩流結果與分析

3.1 工況設置

為了與前文均勻流形成對比，振蕩流流速幅值定為0.4 m/s，振蕩流可表示為：

v=0.4×sin(2πft)

(24)

基于第2節的研究基礎，通過改變振蕩流振蕩周期，研究不同KC數下的二維圓柱受迫振動特性。其中，振蕩流工況設置情況如表2所示。為研究圓柱振動幅度和振動頻率對計算結果的影響，同一振幅比下選擇8組無因次頻率進行計算，無因次振幅比設置6組，共計算240個振蕩流工況。

表2 振蕩流工況設置Tab. 2 Setting of oscillation flow condition

3.2 流體力系數分析

提取0.50D和1.00D兩組振幅比的計算結果進行分析，升力系數幅值變化情況見圖5(a)。在低振動頻率時，升力系數呈現較小的值，隨著振動頻率增加，升力系數幅值也逐漸增加。在該工況下，不同KC數對升力系數幅值影響不大。激勵力系數變化情況見圖5(b)。激勵力系數變化較為復雜，且對KC數的敏感性較高。此時不同KC數下的激勵力系數隨無因次頻率變化基本呈現先增加后減小的趨勢，峰值區出現在0.18附近，相比均勻流的0.17略有變化，符合文獻[12]的研究結論。在振動頻率0.15～0.20之間，KC數對激勵力系數影響較大，該區間也屬于均勻流試驗測得的鎖定發生的關鍵區域。由圖5(b)可知，小KC數下，激勵力系數變化較大。除在峰值區，KC=31.5時激勵力系數基本保持在較低水平。這可能是因為此工況下流速及流向變化較快，漩渦尚未脫落便遭遇反向流速，無法形成穩定的漩渦脫落周期。KC數增大后，激勵力系數也逐漸增加，但隨著KC數變大，流態逐漸趨于定常流，此時激勵力系數值將保持相近。

圖5 低振幅比(0.50D)流體力系數對比Fig. 5 Contrast of fluid force coefficients at low amplitude ratio (0.50D)

圖6給出了振動幅值1.00D時，升力系數幅值和激勵力系數隨KC數的變化情況。隨著振動幅值的增加，KC數對流體力系數的影響效應并不顯著，這與小振幅工況時是不同的。圖6(b)表明，激勵力系數隨著無因次頻率的增加逐漸減小。大KC數工況相對較為穩定，小KC數工況KC=31.5時振蕩流周期變化較快，漩渦脫落相對不穩定，規律性較弱，故在低無因次頻率下呈現出了下降趨勢。但整體來看，各工況下激勵力系數均小于0，此時流體對結構振動起阻尼作用，能量由結構傳向流體。激勵力系數由升力系數幅值和相位角計算得到，由Gopalkrishnan[3]均勻流試驗可知，低振幅比與高振幅比時相位角變化具有明顯不同。振蕩流工況下高低振幅比時相位角變化趨勢也不同，1.00D時相位角均保持為負值，故激勵力系數變化趨勢與0.50D存在較大差別。

圖6 高振幅比(1.00D)流體力系數對比Fig. 6 Contrast of fluid force coefficients at high amplitude ratio (1.00D)

綜上可知，在小振動幅值時，當圓柱的無因次頻率處于渦激振動鎖定區間內，KC數對激勵力系數影響較大，且大、小KC數下呈現不同的規律性。在大振動幅值時，KC數對該雷諾數下流體力系數的影響將逐漸減小，激勵力系數保持負值。

3.3 漩渦形態分析

基于振幅比0.50D計算結果，分析KC數對漩渦發放形態的影響。前文分析可得，振動頻率在0.15～0.20之間流體力系數受KC數變化較為敏感，因此選擇頻率0.17作為分析重點。將5組KC數分為兩類，其中，小KC數工況為組1和組2(KC=31.5和KC=63.0)，大KC數工況為組3、組4和組5(KC=126.0、KC=252.0和KC=503.9)。由于大小兩類工況下不同KC數的漩渦形態具有相似性，故大KC數和小KC數各選擇一組典型形態進行分析。

3.3.1 大KC數工況

圖7給出KC=503.9、無因次頻率0.17、一個振蕩流周期內漩渦脫落發展狀態。隨著流速的逐漸增大，尾渦在圓柱壁面上產生，并逐漸向圓柱右側延伸。當到達圖7(a)時，有漩渦從尾渦后方脫落。圖7(b)對應流速達到峰值，在由(a)到(b)的過程中，速度進一步增加，圓柱后方的尾渦也逐漸伸長。此階段內伸長的尾渦與之前脫落的漩渦連在一起，并未出現顯著漩渦脫離現象。經過圖7(b)后，流速逐漸下降，此時圓柱后方的漩渦開始上下交替脫落，呈現出典型的2S泄渦模式，如圖7(c)所示。圓柱運動到圖7(d)時，此時圓柱后方依舊有漩渦脫落，但已不再是2S模式，并且漩渦的大小和強度均開始下降。

圖7 大KC數下漩渦脫落(A*=0.5，f*=0.17，KC=503.9)Fig. 7 Vortex shedding at large KC number (A*=0.5，f*=0.17，KC=503.9)

當流速繼續下降時，脫落的漩渦也逐漸在圓柱后方消散。圖7(e)顯示此時流速很小，圓柱壁面已無法產生尾渦。之后流速反向并逐漸增加。漩渦開始在圓柱的左側產生，但此時尾渦較短，漩渦在離圓柱較近的位置脫落，如圖7(f)所示。流速進一步增加，尾渦逐漸延伸，但無法觀察到脫落的漩渦。這一階段類似圖7(a)至7(b)。圖7(g)時流速剛過最大值，漩渦開始從尾渦上脫落，并呈現出2S泄渦模式，與流速正向不同的是，這一階段2S模式的持續時間要更小。圖7(h)后，流速逐漸下降，圓柱的尾渦長度慢慢變短，強度逐漸減弱。

3.3.2 小KC數工況

圖8為KC數為31.5時的尾渦演化。圖8(a)時，流速剛經過最大值，尾渦在圓柱后側產生延長。隨著流速下降，從圖8(b)可以發現有單個漩渦從圓柱壁面上脫落，但由于沒有長時間的單側流向，難以觀察到漩渦一個接一個在圓柱的下流向消散的現象。圖8(c)時，流速反向，脫落的漩渦被反向流速帶回到圓柱附近，導致圓柱周圍漩渦分布復雜化，這也是小KC數時的典型狀態。之后流速增加，尾渦在圓柱的左側產生。

圖8 小KC數下漩渦脫落(A*=0.5，f*=0.17，KC=31.5)Fig. 8 Vortex shedding at low KC number (A*=0.5，f*=0.17，KC=31.5)

在f*=0.17下，振蕩流隨著KC數的增加，漩渦脫落的模式各有不同。小KC數下流體速度變化較快，漩渦從圓柱壁面脫落后難以保持穩定泄渦，當流速反向后，脫落的漩渦又被沖向圓柱，使得圓柱周圍壓力變化更加復雜。在大KC數時，漩渦以典型的2S模式從圓柱壁面脫落，且隨著速度下降，漩渦強度逐漸下降，脫落的漩渦也可在流速反向前消散完成。

3.4 升力系數時歷分析

為更好地分析漩渦產生發展過程和升力系數之間的關系，提取了A*=0.5，f*=0.17對應的升力系數時歷曲線，并對照漩渦發放過程，探討二者之間的內在聯系。

圖9展示了大KC數時的升力系數時歷曲線。3種大KC數工況均觀察到了振幅調制現象。圖9(c)給出了KC=503.9時升力系數變化情況。該曲線可以看到兩個升力系數增大區，即振蕩流速度幅值附近(72 s和88 s)。當運動到64 s時，振蕩流速度由0開始逐漸增大，升力系數也逐漸增加。當圓柱運動至72 s時，進入到2S發放模式，升力系數進入第一個峰值區。80 s時，流速降低為0，可以觀察到此時升力系數要小于峰值區。流速反向后，速度開始增加，再次進入到2S發放模式，此時升力系數同樣達到峰值區。

圖9 大KC數升力系數曲線Fig. 9 Amplitude curve of lift coefficient at large KC number

KC=126.0時，峰值區出現在18 s和22 s附近，并且可以看到第二個峰值區的范圍要略大于第一個，且峰值更高。KC=252.0時，峰值區出現在36 s和44 s，該工況下振幅調制現象較為對稱。

圖10給出了小KC數及均勻流工況下升力系數時歷曲線。當KC數較小時，流體速度變化較快，在一個振蕩流周期內無法觀察到明顯的兩次振幅調制現象，但振幅仍隨時間呈現波動規律。圖10(c)給出了均勻流下10個漩渦脫落周期內升力系數變化情況，升力系數基本不隨時間發生變化。

圖10 小KC數及均勻流升力系數曲線Fig. 10 Amplitude curve of lift coefficient at low KC number

4 結 語

基于CFD方法對振蕩流下的振動圓柱進行了數值模擬，分析了不同KC數工況圓柱的水動力特性及泄渦形態，得到以下結論：

1) 采用CFD方法復現均勻流圓柱受迫振動試驗的研究手段是可行的，對比流體力系數可知，CFD模擬結果在趨勢和數值上均與試驗結果吻合較好。

2) 振蕩流下，在低振幅比、低振動頻率時，振動圓柱的激勵力系數會隨著KC數改變發生顯著變化。小KC數時，流速變化較快，激勵力系數受KC數影響較大，隨著KC數逐漸增加，流態趨向均勻流，激勵力系數數值差異減小。隨著圓柱振幅比增加，不同KC數下的流體力系數基本保持一致。

3) 大KC數工況下由于流體在單向上可保持較長時間，漩渦有足夠時間發展和脫落，可觀察到明顯的2S泄渦模式；小KC數工況時，流體變化加快，脫落的漩渦還未完全消散便被反向流速沖回圓柱，漩渦發放規律性較弱。

4) 大KC數工況下的升力系數時歷曲線可在同一振蕩流周期中觀察到兩次明顯的振幅調制現象，分別對應正向、反向兩個流速階段；小KC數工況時，流速頻繁改變大小和方向，升力系數時歷不像均勻流下穩定，但變化不具備明顯規律性。