振蕩流中圓柱體橫流渦激振動(dòng)特性數(shù)值分析

王坤鵬， 遲慶海， 張一兆

(江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003)

海洋立管是連接海底和海洋平臺的細(xì)長柔性結(jié)構(gòu)。在來流作用下，立管兩側(cè)的渦脫落將導(dǎo)致立管在橫流和順流2個(gè)方向產(chǎn)生渦激振動(dòng)，當(dāng)渦脫落頻率和立管振動(dòng)頻率鎖定時(shí)，渦激振動(dòng)會(huì)被放大，從而加劇立管的疲勞破壞。因此，有必要對渦激振動(dòng)問題開展深入的研究。

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

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

本文基于Fluent軟件建立了振蕩流中圓柱體橫流渦激振動(dòng)的二維數(shù)值模型。該模型通過兩端邊界的簡諧運(yùn)動(dòng)實(shí)現(xiàn)振蕩流，并采用了RANS方程和SSTk-ω湍流模型模擬流場。采用該模型對已公布的實(shí)驗(yàn)?zāi)Ｐ瓦M(jìn)行分析，驗(yàn)證了數(shù)值模型的有效性?；谠撃Ｐ烷_展了振蕩流KC數(shù)為25和50這2種情況下圓柱體的橫流渦激振動(dòng)模擬，重點(diǎn)研究了圓柱體的激發(fā)模態(tài)、流體力和渦脫落特征。

1 振蕩流數(shù)值方法

盡管在湍流狀態(tài)下，渦激振動(dòng)的渦脫落具有一定的三維特征，但是已公布的研究表明二維數(shù)值模型可以較好的模擬流體力和渦脫落模態(tài)[9]，因此本文采用二維模型進(jìn)行模擬。圖1為二維圓柱體渦激振動(dòng)的計(jì)算流域，圓柱體直徑為D，流向長度為80D，橫向?qū)挾葹?5D，圓柱體位于計(jì)算流域的中心位置。通過左右兩側(cè)邊界的來回運(yùn)動(dòng)產(chǎn)生振蕩流，其速度表達(dá)式為：

U(t)=Umcos(2πfost)

(1)

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

圖1 計(jì)算流域

對于振蕩流中流場，本文采用非定常不可壓縮的RANS控制方程模擬。在考慮圓柱體運(yùn)動(dòng)帶來的網(wǎng)格變化時(shí)，RANS方程的表達(dá)式為：

(2)

(3)

(4)

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

圓柱體兩側(cè)的渦脫落產(chǎn)生的升力可以激勵(lì)圓柱體橫流方向的振動(dòng)。圓柱體振動(dòng)的控制方程可表示為:

(5)

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

圓柱體和左右邊界的運(yùn)動(dòng)通過UDF傳入到Fluent，兩者的運(yùn)動(dòng)會(huì)引起網(wǎng)格的變化，本文采用擴(kuò)散光順的網(wǎng)格光順方法，通過求解擴(kuò)散方程得到網(wǎng)格節(jié)點(diǎn)的位移，從而實(shí)現(xiàn)網(wǎng)格運(yùn)動(dòng)。網(wǎng)格節(jié)點(diǎn)的計(jì)算位移公式為：

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

·(γUj)=0

(6)

式中γ為擴(kuò)散系數(shù)，本文取值1.5。

2 驗(yàn)證分析

本節(jié)采用振蕩流中的渦激振動(dòng)實(shí)驗(yàn)?zāi)Ｐ蚚13]對數(shù)值模型進(jìn)行驗(yàn)證，并與Zhao等[18]的數(shù)值結(jié)果進(jìn)行比較。實(shí)驗(yàn)?zāi)Ｐ偷南嚓P(guān)參數(shù)為：KC=10，D=0.02 m，fN=0.77 HZ。計(jì)算流域內(nèi)設(shè)置了3個(gè)速度監(jiān)測點(diǎn)A+、B+、C+，來分析兩端邊界的簡諧振動(dòng)產(chǎn)生的流場，如圖1所示。圖2為部分區(qū)域計(jì)算網(wǎng)格，圓柱體表面有200個(gè)節(jié)點(diǎn)，邊界層厚度為0.001D，保證無量綱壁面距離y+小于1.0，無量綱步長τ為0.002，這里定義τ=(tUm/D)。

圖2 圓柱體周圍網(wǎng)格

圖3(a)是Ur=10，Um=0.154 m/s時(shí)，3個(gè)速度監(jiān)測點(diǎn)沿流向的速度隨時(shí)間變化曲線，Ur為最大折合速度，定義為Ur=UmfN/D?？梢园l(fā)現(xiàn)A+、B+、C+的速度變化相同且最大值為0.154 m/s，與邊界Um相同；另外，從圖3(b)可以看出，除了圓柱體的附近區(qū)域，流場的速度均勻分布。因此，邊界往返運(yùn)動(dòng)可以有效地產(chǎn)生振蕩流場。

圖3 流場速度

圖4是本文預(yù)報(bào)的平均幅值比(Ay/D)與Zhao等[18]模擬的數(shù)值結(jié)果以及Sumer[13]的實(shí)驗(yàn)結(jié)果的曲線圖?？梢钥闯?，在Ur小于5時(shí)，本文驗(yàn)證結(jié)果與Zhao等[18]的結(jié)果相近，均高于實(shí)驗(yàn)結(jié)果；但在Ur大于5時(shí)，本文結(jié)果介于文獻(xiàn)[18]和實(shí)驗(yàn)結(jié)果之間。與實(shí)驗(yàn)結(jié)果相比，數(shù)值預(yù)報(bào)產(chǎn)生的鎖定區(qū)間偏大，這不排除實(shí)驗(yàn)誤差的原因。因此，需要后續(xù)的實(shí)驗(yàn)研究和更有效的三維數(shù)值模型開展進(jìn)一步驗(yàn)證。

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

3 算例分析與討論

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

3.1 位移和頻率

圖5是KC=25、最大折合速度Ur為5、8和15時(shí),圓柱體在y方向的時(shí)間位移歷程曲線?？梢钥闯鲈赨r=5和8時(shí)，振動(dòng)較為規(guī)律；而在Ur=15時(shí)，振幅曲線上部出現(xiàn)2次峰值，下部幅值相對較大，同樣的現(xiàn)象在Zhao等[18]模擬的KC=10的情況也有出現(xiàn)。

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

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

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

圖7 KC=25時(shí)的頻譜

圖8 KC=50時(shí)的頻譜

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

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

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

圖11 Ur=15時(shí)的流體力幅值譜

3.2 渦脫落和流體力

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

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

圖12 不同時(shí)刻的渦量變化等值線

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

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

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

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

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

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

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

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

4 結(jié)論

1) 振蕩流能激發(fā)圓柱體的多模態(tài)振動(dòng)，且KC數(shù)越大多模態(tài)特性越明顯；另外，處于鎖定區(qū)間時(shí)，主導(dǎo)模態(tài)的幅值明顯高于其他模態(tài)。

2) 隨著KC數(shù)的增大，振蕩流中的鎖定區(qū)間有變寬趨勢，一個(gè)振蕩流周期內(nèi)從圓柱體脫落的渦數(shù)也會(huì)增加，而KC數(shù)的增大似乎對渦激振動(dòng)幅值比包絡(luò)線的最大值影響較小。

3) 在流速方向發(fā)生改變時(shí)，圓柱體已脫落的渦會(huì)隨著流體從圓柱體的上下兩側(cè)交替經(jīng)過，增加圓柱體的流體力。

4) 在大KC數(shù)和大Ur的情況下，最大流速附近出現(xiàn)多個(gè)周期的高頻流體力，并引起圓柱體的小幅高頻振動(dòng)。