振蕩流中圓柱體橫流渦激振動特性數值分析

王坤鵬， 遲慶海， 張一兆

(江蘇科技大學 船舶與海洋工程學院，江蘇 鎮江 212003)

海洋立管是連接海底和海洋平臺的細長柔性結構。在來流作用下，立管兩側的渦脫落將導致立管在橫流和順流2個方向產生渦激振動，當渦脫落頻率和立管振動頻率鎖定時，渦激振動會被放大，從而加劇立管的疲勞破壞。因此，有必要對渦激振動問題開展深入的研究。

回顧過去近半個世紀的時間里，對于渦激振動的研究已經取得了諸多成果，具體可參閱文獻[1-4]。彈性支撐圓柱體是研究渦激振動的典型結構，許多學者采用該結構開展了圓柱體單自由度[5-6]和雙自由度[7-8]的渦激振動實驗研究。隨著計算機和計算流體力學(CFD)方法的不斷發展，越來越多的學者開始關注基于CFD的渦激振動研究，該方法直接求解N-S方程，可以有效捕捉流體力和渦脫落模態等[9-12]。

在實際海洋環境中，波浪和頂端平臺運動帶動立管與流體之間產生的相對振蕩流會導致渦激振動的產生。Sumer等[13]、Kozakiewicz等[14]和王俊高等[15]采用實驗的方法對振蕩流中圓柱體的渦激振動特性開展研究，發現了間歇性振動和模態躍遷等現象，推動了振蕩流中渦激振動數值研究的發展。鄧迪等[16]采用切片法對王俊高的實驗模型進行了數值模擬，有效捕捉到了振蕩流中渦激振動的現象。鄧躍[17]選取KC數為0.56和2.51的振蕩流，研究了均勻流和振蕩流疊加下圓柱體的渦激振動，分析了圓柱體的振動幅值和鎖定區間。Zhao等[18]研究了振蕩流KC數為10和20時圓柱體的渦激振動，發現當最大折合速度小于8時，單一頻率振動，而當最大折合速度大于8時，多頻率振動，且振動周期會受到渦脫數量和方向的影響。整體來看，針對振蕩流中圓柱體渦激振動的研究采用的KC數較小，且缺少渦脫落模態和流體力之間關系的分析。

本文基于Fluent軟件建立了振蕩流中圓柱體橫流渦激振動的二維數值模型。該模型通過兩端邊界的簡諧運動實現振蕩流，并采用了RANS方程和SSTk-ω湍流模型模擬流場。采用該模型對已公布的實驗模型進行分析，驗證了數值模型的有效性。基于該模型開展了振蕩流KC數為25和50這2種情況下圓柱體的橫流渦激振動模擬，重點研究了圓柱體的激發模態、流體力和渦脫落特征。

1 振蕩流數值方法

盡管在湍流狀態下，渦激振動的渦脫落具有一定的三維特征，但是已公布的研究表明二維數值模型可以較好的模擬流體力和渦脫落模態[9]，因此本文采用二維模型進行模擬。圖1為二維圓柱體渦激振動的計算流域，圓柱體直徑為D，流向長度為80D，橫向寬度為25D，圓柱體位于計算流域的中心位置。通過左右兩側邊界的來回運動產生振蕩流，其速度表達式為：

U(t)=Umcos(2πfost)

(1)

式中：Um為邊界往返運動的最大速度；fos為往返的頻率；t為時間。振蕩流的KC數定義為KC=Um/(fosD)。

圖1 計算流域

對于振蕩流中流場，本文采用非定常不可壓縮的RANS控制方程模擬。在考慮圓柱體運動帶來的網格變化時，RANS方程的表達式為：

(2)

(3)

(4)

式中：vt為湍流粘性系數；k為湍流動能；δij為克羅內克爾符號。由于雷諾應力項的存在，本文結合SSTk-ω湍流模型求解方程(3)，該湍流模型廣泛應用于渦激振動的研究，且預報結果理想[10]。

圓柱體兩側的渦脫落產生的升力可以激勵圓柱體橫流方向的振動。圓柱體振動的控制方程可表示為:

(5)

式中：y是圓柱體垂直于來流方向上的位移；m、c、k分別為單位長度圓柱體的質量、結構阻尼與彈簧剛度；Fy表示作用在圓柱體上的升力，通過圓柱體表面的壓力積分得到。本文采用Newmark-β方法求解方程(5)。

圓柱體和左右邊界的運動通過UDF傳入到Fluent，兩者的運動會引起網格的變化，本文采用擴散光順的網格光順方法，通過求解擴散方程得到網格節點的位移，從而實現網格運動。網格節點的計算位移公式為：

兩大流派的研究各有其優勢、貢獻和不足。工作意義流派的研究見長于分析個人行為動機和過程。通過深入觀察多種職業員工的工作重塑現象，此類研究發現了各種不同的行為策略，歸納出了覆蓋范圍較廣的工作重塑行為方式，洞察了認知重塑對個人工作意義和工作身份的重要影響。這些研究深化了人們對工作重塑理論的認識，也發展了理論模型。同時，由于此類視角的研究往往采取更加耗時的質性方法，對研究者的時間、精力和分析理解能力都要求較高，所以相應的研究數量也很有限。

·(γUj)=0

(6)

式中γ為擴散系數，本文取值1.5。

2 驗證分析

本節采用振蕩流中的渦激振動實驗模型[13]對數值模型進行驗證，并與Zhao等[18]的數值結果進行比較。實驗模型的相關參數為：KC=10，D=0.02 m，fN=0.77 HZ。計算流域內設置了3個速度監測點A+、B+、C+，來分析兩端邊界的簡諧振動產生的流場，如圖1所示。圖2為部分區域計算網格，圓柱體表面有200個節點，邊界層厚度為0.001D，保證無量綱壁面距離y+小于1.0，無量綱步長τ為0.002，這里定義τ=(tUm/D)。

圖2 圓柱體周圍網格

圖3(a)是Ur=10，Um=0.154 m/s時，3個速度監測點沿流向的速度隨時間變化曲線，Ur為最大折合速度，定義為Ur=UmfN/D。可以發現A+、B+、C+的速度變化相同且最大值為0.154 m/s，與邊界Um相同；另外，從圖3(b)可以看出，除了圓柱體的附近區域，流場的速度均勻分布。因此，邊界往返運動可以有效地產生振蕩流場。

圖3 流場速度

圖4是本文預報的平均幅值比(Ay/D)與Zhao等[18]模擬的數值結果以及Sumer[13]的實驗結果的曲線圖。可以看出，在Ur小于5時，本文驗證結果與Zhao等[18]的結果相近，均高于實驗結果；但在Ur大于5時，本文結果介于文獻[18]和實驗結果之間。與實驗結果相比，數值預報產生的鎖定區間偏大，這不排除實驗誤差的原因。因此，需要后續的實驗研究和更有效的三維數值模型開展進一步驗證。

圖4 KC=10振蕩流中不同Ur下的振幅比

3 算例分析與討論

本節對KC=25和KC=50，最大折合速度Ur從2～15的情況進行了數值模擬，并對位移、激發頻率、流體力和渦脫落特征進行了分析。圓柱體直徑D為0.04 m，質量比為2.2，靜水中的固有頻率fN為0.98 Hz，雷諾數Re由最大流速確定，取值范圍為3 185～23 208。

3.1 位移和頻率

圖5是KC=25、最大折合速度Ur為5、8和15時,圓柱體在y方向的時間位移歷程曲線。可以看出在Ur=5和8時，振動較為規律；而在Ur=15時，振幅曲線上部出現2次峰值，下部幅值相對較大，同樣的現象在Zhao等[18]模擬的KC=10的情況也有出現。

圖5 KC=25圓柱體y方向位移歷程曲線

圖6是KC=50、最大折合速度Ur為5、8和15時，圓柱體在y方向的時間位移歷程曲線。觀察到Ur=5時曲線變化與KC=25、最大折合速度Ur為8的情況相近；而在Ur=15時，2次峰值不再只出現于正y方向，而是上下交替出現。綜合KC=25和50的情況來看，在Ur較小時，圓柱體在振蕩流中的渦激振動比較規律。

圖6 KC=50圓柱體y方向位移歷程曲線

圖7 KC=25時的頻譜

圖8 KC=50時的頻譜

下面主要討論隨著最大折合速度的增加，平均振幅比Ay和各激發模態的振幅及對應頻率比的變化情況。每一個Ur取3個主要激發模態，振幅定義為Ay1，Ay2和Ay3，對應的頻率依次為f1、f2和f3。當KC=25時，結合圖9(a)中最大激發模態的振幅和頻率比變化可以發現，鎖定從Ur=4開始， 在Ur=10時結束；當KC=50時，從圖10(a)中可以看到，鎖定區間為Ur從4～12。可見，隨著KC數的增大，鎖定區間有變大趨勢。從圖9(a) 和10(a)的幅值比包絡線可以看出，KC數對渦激振動中的最大幅值影響較小，均在0.8倍直徑附近。從圖9(b) 和10(b)中頻率比的變化可以發現，當處于鎖定區間時，幅值最大模態的頻率(主導頻率)鎖定在固有頻率附近；鎖定結束后，當KC=25時，主導頻率隨著最大折合速度的增大而增大，但是當KC=50時，主導頻率仍然在固有頻率附近，其原因可能是大KC數對應大振蕩周期，所以低流速持續的時間較長，容易產生低頻渦脫落，從而產生低頻流體力，如圖11(b)中虛線圈出部分，而離固有頻率較近的低頻流體力更容易激發圓柱的振動；另外，KC數較大時渦脫落產生的高頻流體力的頻率間隔要小于KC數較小的情況，如圖11所示，而相近的高頻流體力會相互干擾，不利于鎖定。

圖9 KC=25不同最大折合速度下的振幅和頻率比

圖10 KC=50不同最大折合速度下的振幅和頻率比

圖11 Ur=15時的流體力幅值譜

3.2 渦脫落和流體力

圓柱體的受力和振動可以通過觀察圓柱體的渦脫落和流場變化來解釋，由于云圖較多，現結合圓柱體在y方向的受力和位移曲線選取幾組典型渦量變化等值線圖進行分析，圖中X=x/D，Y=y/D，橫向箭頭表示流向，豎直箭頭表示圓柱體運動方向。

圖12和圖13是KC=25、Ur=5時，不同時刻的渦量等值線圖、圓柱體的流體力和位移曲線圖。從圖12(a)和12(b)中可以發現，圖中渦1和渦2是上一周期脫落的漩渦，在流向改變后，分別從圓柱體的上方和下方經過，此時對應圖13中在t=81.6 s和82.08 s左右出現流體力峰值的時候。從圖12(c)和12(d)中可以發現，流速經過峰值(或谷值)后，圓柱體的渦脫模式發生明顯改變，前者為流速谷值之前，后者為谷值之后，此時對應圖13中A區域受力曲線出現小幅波動，相同的情況如圖中的虛線箭頭。

圖12 不同時刻的渦量變化等值線

圖13 振蕩流速和圓柱體y方向不同時刻的受力位移曲線

當KC=50，Ur=5時，圖14中出現了“2P”形式的渦脫模式，但這種渦脫模式并不是對稱的，如圖14(b)和14(c)所示，這種改變引起了圖15中流體力曲線在B區域出現持續的小幅波動。變化較穩定的A區域，有較大的流體力，原因同圖13一樣，也與脫落的渦在流向改變后經過圓柱體有關。另外還可以看出，圖15中位移曲線的峰值略遲于力曲線的峰值，但與流速的最大值差不多同時出現，在圖13中也出現同樣的現象。

在Ur=15時，KC=25振蕩流的流速及對應的圓柱體受力和位移見圖16。在圖17(b)中，振蕩流方向開始變化，即將脫落的渦1沿著圓柱體上表面滑動到另一側，導致圖16中受力曲線A區域的小幅波動。KC=50時的渦量云圖及流體力特征與KC=25時相似，圖18 (a)的渦量云圖與圖19中受力曲線A區域的小幅波動對應。由圖16和圖19中可以發現，流體力的最大值和流速的峰值基本同時出現，與上述Ur=5的結果不同。

圖14 KC=50,Ur=5時不同時刻的渦量變化等值線圖

圖15 振蕩流速和圓柱體y方向不同時刻的受力位移曲線

圖16 振蕩流速和圓柱體y方向不同時刻的受力位移曲線

圖18 KC=50,Ur=15時不同時刻的渦量變化等值線圖

在Ur=15時，盡管KC=25和50對應的渦脫落方式相似，但KC=50時渦數量偏多，這是因為大KC數振蕩流產生渦脫落的流速時間跨度較大。在流向反轉之后，這些脫落的渦反向運動經過圓柱體，影響圓柱體的受力和振動。這里主要討論KC=50的情況，如圖18(b)，上半周期脫落的渦2主要從圓柱體上側“相交”而過，并被圓柱體分為2個大小不等的渦，渦3主要從圓柱體下方“相切”而過，渦4從圓柱體上方“相切”而過。這些渦經過圓柱體的時刻對應流體力的幅值1、2和3，均高于后續出現的幅值，綜合之前的情況說明流向反轉后經過圓柱體的渦可以增加圓柱體的流體力。由于圓柱體的受力出現1、2和3及后續的高頻幅值，導致渦激振動位移曲線出現B和C 2個區域的高頻振動，振動幅度較小的原因是高頻力的頻率遠離靜水中的固有頻率。

4 結論

1) 振蕩流能激發圓柱體的多模態振動，且KC數越大多模態特性越明顯；另外，處于鎖定區間時，主導模態的幅值明顯高于其他模態。

2) 隨著KC數的增大，振蕩流中的鎖定區間有變寬趨勢，一個振蕩流周期內從圓柱體脫落的渦數也會增加，而KC數的增大似乎對渦激振動幅值比包絡線的最大值影響較小。

3) 在流速方向發生改變時，圓柱體已脫落的渦會隨著流體從圓柱體的上下兩側交替經過，增加圓柱體的流體力。

4) 在大KC數和大Ur的情況下，最大流速附近出現多個周期的高頻流體力，并引起圓柱體的小幅高頻振動。