尾流干涉下串聯(lián)圓柱渦激振動(dòng)數(shù)值模擬研究*

王曉凱 婁 敏

(中國(guó)石油大學(xué)(華東)石油工程學(xué)院)

0 引 言

T.KITAGAWA等[7]在亞臨界雷諾數(shù)下改變間距比，研究了上游圓柱的尾流對(duì)下游圓柱流體力的影響。B.S.CARMO等[8]對(duì)兩圓柱串聯(lián)流動(dòng)進(jìn)行仿真分析，上游圓柱保持固定，下游圓柱僅允許在橫流向自由振動(dòng)，與孤立圓柱對(duì)比結(jié)果表明，串聯(lián)圓柱的振動(dòng)最大位移顯著增大。郭曉玲等[9]對(duì)低雷諾數(shù)下(Re=150)等直徑串聯(lián)圓柱的流動(dòng)干涉進(jìn)行了數(shù)值模擬研究，發(fā)現(xiàn)圓柱間距比的改變會(huì)導(dǎo)致“鎖定”區(qū)間發(fā)生變化。ZHAO M.等[10]重點(diǎn)關(guān)注了大圓柱對(duì)小圓柱尾流與振動(dòng)的影響，在特定的間距范圍內(nèi)發(fā)現(xiàn)圓柱響應(yīng)會(huì)出現(xiàn)4個(gè)不同的分區(qū)。宋振華等[11-12]研究了順流向振動(dòng)對(duì)橫流向最大振幅和鎖定區(qū)間的影響，同時(shí)針對(duì)3根附屬桿抑制裝置對(duì)立管所受流體力的抑制效果、泄渦模式及抑制機(jī)理進(jìn)行了深入研究。F.J.HUERA-HUARTE等[13]改變上游圓柱與下游圓柱的直徑比，通過振幅和流體力系數(shù)的變化討論了尾流干涉的影響。婁敏等[14-15]在尾流干涉下考慮流固耦合作用對(duì)串聯(lián)圓柱進(jìn)行了仿真研究，發(fā)現(xiàn)隨著約化速度增大，渦激振動(dòng)遠(yuǎn)離“鎖定”狀態(tài)的多頻現(xiàn)象逐漸減弱；同時(shí)采用不同的抑振裝置對(duì)串聯(lián)立管渦激振動(dòng)抑制進(jìn)行了試驗(yàn)分析，綜合比較得出控制桿的抑制效果更佳。

從目前情況來看，對(duì)于圓柱渦激振動(dòng)的研究已經(jīng)趨于成熟，然而多數(shù)研究集中在單自由度情況下，同時(shí)涉及到小尺寸圓柱尾流干涉下的研究十分有限。為此，本文利用Fluent軟件中的SSTk-ω湍流模型，基于新穎的Overset-Mesh方法對(duì)流場(chǎng)建模處理，通過雷諾平均法求解黏性N-S方程，并運(yùn)用動(dòng)網(wǎng)格技術(shù)實(shí)現(xiàn)流固耦合，針對(duì)小尺寸串聯(lián)雙圓柱在雙自由度下的渦激振動(dòng)位移響應(yīng)、受力特性、運(yùn)動(dòng)軌跡以及渦脫模式進(jìn)行二維仿真研究。研究結(jié)果可以為深海油氣開發(fā)中立管的渦激振動(dòng)分析提供理論依據(jù)。

1 數(shù)值方法

1.1 Overset Mesh

Overset Mesh又被稱為重疊網(wǎng)格或者嵌套網(wǎng)格，重疊網(wǎng)格通常由一個(gè)背景網(wǎng)格和至少一個(gè)前景網(wǎng)格所組成，如圖1所示。重疊網(wǎng)格是將流體計(jì)算域劃分成多個(gè)簡(jiǎn)單的部件區(qū)域，然后各部件區(qū)域單獨(dú)劃分Parts網(wǎng)格，并將其嵌入到背景網(wǎng)格中，再通過Overset Interface連接重疊網(wǎng)格單元區(qū)域，在前景網(wǎng)格和背景網(wǎng)格的接連區(qū)域通過插值來實(shí)現(xiàn)流場(chǎng)信息的傳遞，迭代計(jì)算的過程中會(huì)不斷搜索重疊區(qū)域邊界，將真實(shí)的物理邊界識(shí)別出來并參與計(jì)算。

圖1 重疊網(wǎng)格構(gòu)成圖Fig.1 Composition of overset mesh

利用重疊網(wǎng)格技術(shù)可以處理具有小縫隙部件的相對(duì)網(wǎng)格運(yùn)動(dòng)，與傳統(tǒng)網(wǎng)格相比，重疊網(wǎng)格對(duì)網(wǎng)格拓?fù)湟蟮?，且網(wǎng)格生成相對(duì)容易，同時(shí)可以避免傳統(tǒng)動(dòng)網(wǎng)格技術(shù)所面臨的網(wǎng)格動(dòng)態(tài)更新時(shí)易出現(xiàn)的負(fù)體積現(xiàn)象，在網(wǎng)格運(yùn)動(dòng)期間始終能夠保持很高的網(wǎng)格質(zhì)量，實(shí)現(xiàn)了局部結(jié)構(gòu)網(wǎng)格在非結(jié)構(gòu)網(wǎng)格中的使用。Overset Mesh在計(jì)算流體力學(xué)工程實(shí)際應(yīng)用中的流程如圖2所示。

圖2 重疊網(wǎng)格在數(shù)值計(jì)算中的流程Fig.2 Numerical calculation process of overset mesh

1.2 控制方程

不可壓縮黏性牛頓流體的控制方程為N-S方程。在直角坐標(biāo)系下，基于雷諾平均的連續(xù)性方程和動(dòng)量方程分別為：

U=0

(1)

(2)

漩渦泄放對(duì)圓柱體作用力如圖3所示。由圖3可見，流體流經(jīng)圓柱時(shí)會(huì)在其兩側(cè)產(chǎn)生漩渦泄放現(xiàn)象，在橫流向方向產(chǎn)生升力FL造成橫向運(yùn)動(dòng)，同時(shí)在順流向方向產(chǎn)生脈動(dòng)阻力FD造成順流向的運(yùn)動(dòng)。升力與脈動(dòng)阻力的表達(dá)式分別為：

(3)

(4)

圖3 漩渦泄放對(duì)圓柱體作用力Fig.3 Acting force of vortex discharge on the cylinder

式中：CL為升力系數(shù)；CD為阻力系數(shù)，與雷諾數(shù)Re有關(guān)，無量綱；ρ為流體密度，kg/m3；D為圓柱體直徑，m。

早期對(duì)于渦激振動(dòng)的研究主要集中在較大質(zhì)量比的圓柱，順流向振幅極小，因此大多數(shù)研究選擇忽略順流向運(yùn)動(dòng)而采用主要觀察橫流向運(yùn)動(dòng)響應(yīng)特征的單自由度模型。雖然順流向的阻力FD通常要比橫流向的升力FL小得多，但由于阻力周期是升力的2倍，所以順流向阻力對(duì)圓柱結(jié)構(gòu)的影響不可忽略。對(duì)于本文采用的小質(zhì)量比圓柱，順流向運(yùn)動(dòng)會(huì)對(duì)橫流向造成影響，與橫流向的單自由運(yùn)動(dòng)相比，兩自由度運(yùn)動(dòng)能夠使橫流向獲得相對(duì)較大的幅值。圓柱振動(dòng)模型如圖4所示。

圖4 圓柱振動(dòng)模型Fig.4 Cylindrical vibration model

本文根據(jù)圓柱在流場(chǎng)中的受力形式和運(yùn)動(dòng)情況將振動(dòng)模型簡(jiǎn)化為雙自由度(2-DOF)彈簧-質(zhì)量-阻尼系統(tǒng)，便于研究順流向與橫流向之間的耦合運(yùn)動(dòng)。振動(dòng)模型控制方程為：

(5)

(6)

式中：x、y分別為圓柱順流向位移與橫流向位移，m；m為圓柱質(zhì)量，kg；c為結(jié)構(gòu)阻尼系數(shù)；k為系統(tǒng)剛度系數(shù)；FD(t)為圓柱所受順流向阻力，N；FL(t)為圓柱所受橫流向升力，N。

2 數(shù)值模型

2.1 模型參數(shù)

數(shù)值模擬中圓柱模型的基本參數(shù)如表1所示。根據(jù)T.K.PRASANTH等[16]對(duì)網(wǎng)格尺寸的研究，當(dāng)滿足尾流區(qū)長(zhǎng)度大于25D，整體高度區(qū)域大于20D時(shí)，圓柱體兩自由度渦激振動(dòng)的響應(yīng)將不再受流體區(qū)域邊界的影響。流場(chǎng)區(qū)域計(jì)算模型如圖5所示。選取D=0.018 m，兩圓柱為等徑串聯(lián)排布，圓柱間距比L/D=4(L表示兩圓柱圓心間的距離)，整個(gè)流場(chǎng)區(qū)域設(shè)為40D×20D，下游圓柱尾流區(qū)域長(zhǎng)度為26D。流場(chǎng)建模時(shí)基于Overset-Mesh技術(shù)進(jìn)行網(wǎng)格劃分，如圖6所示。從圖6可見，整個(gè)流場(chǎng)區(qū)域由矩形的藍(lán)色背景網(wǎng)格和2個(gè)圓環(huán)形的紅色前景網(wǎng)格組成，對(duì)圓柱體周圍和漩渦發(fā)散區(qū)域進(jìn)行加密處理，其他區(qū)域進(jìn)行線性加密，逐漸增大網(wǎng)格密度，合理控制網(wǎng)格總量。通過設(shè)定初始流速的方式來控制流場(chǎng)雷諾數(shù)的變化，在設(shè)置計(jì)算區(qū)域的邊界時(shí)，入口處采用速度邊界(Velocity-inlet)，出口處采用壓力邊界(Pressure-outlet)，上、下兩側(cè)面為對(duì)稱邊界條件(Symmetry)，流場(chǎng)與立管重疊部分為無滑移壁面邊界(Wall)。

表1 圓柱模型基本參數(shù)Table 1 Basic parameters of cylindrical model

圖5 流場(chǎng)區(qū)域計(jì)算模型Fig.5 Computational model of flow field area

圖6 網(wǎng)格劃分Fig.6 Mesh generation

本文數(shù)值模擬研究采用了SST(Shear Stress Transport)k-ω模型。SSTk-ω湍流模型屬于雷諾平均法，是一種混合兩方程的渦黏性分區(qū)模型，該湍流模型兼具了k-ω模型和k-ε模型的優(yōu)點(diǎn)。當(dāng)計(jì)算模型近壁面區(qū)域時(shí)，使用k-ω模型可以更好地解決邊界層問題；而在遠(yuǎn)離壁面的區(qū)域則通過k-ε模型能夠較好地模擬發(fā)展充分的湍流流動(dòng)，這使得SSTk-ω模型的應(yīng)用范圍變得更加廣泛。同時(shí)SSTk-ω模型考慮了湍流剪切應(yīng)力的傳播方式，并對(duì)ω方程的交叉擴(kuò)散進(jìn)行了合并，在求解負(fù)壓梯度問題時(shí)會(huì)有更好的表現(xiàn)，整體計(jì)算的精確度和可靠性也得到提升。

2.2 模型驗(yàn)證

為檢驗(yàn)上述計(jì)算模型與網(wǎng)格劃分的可靠性，當(dāng)雷諾數(shù)為200時(shí)，分別進(jìn)行單圓柱繞流和間距比為4的串聯(lián)雙圓柱繞流數(shù)值模擬，得到相應(yīng)的升力和阻力系數(shù)，模擬結(jié)果如圖7所示。

圖7 圓柱繞流模擬結(jié)果Fig.7 Simulation results of flow around the cylinder

將升力系數(shù)和阻力系數(shù)的相關(guān)數(shù)據(jù)通過MATLAB讀取，并進(jìn)行FFT(Fast Fourier Transform)傅里葉變換，得到漩渦脫落頻率fs，再通過fs計(jì)算斯特勞哈爾數(shù)Sr：

(7)

將所得結(jié)果與文獻(xiàn)[17-21]等數(shù)值模擬結(jié)果進(jìn)行對(duì)比，如表2和表3所示。從表2和表3可見，模擬得到的升力系數(shù)CL、阻力系數(shù)CD和斯特勞哈爾數(shù)Sr與相關(guān)文獻(xiàn)中的參考值均較為接近，確保了所選取的網(wǎng)格劃分方法、數(shù)值求解格式等參數(shù)設(shè)置的可靠性，可采用該模型繼續(xù)進(jìn)行流固耦合的計(jì)算。表3中的CD1表示上游阻力系數(shù)，CD2表示下游阻力系數(shù)。

表2 單圓柱繞流結(jié)果對(duì)比Table 2 Comparison between simulation results of flow around a single cylinder

表3 雙圓柱繞流結(jié)果對(duì)比Table 3 Comparison between simulation results of flow around two cylinders

3 計(jì)算結(jié)果

3.1 圓柱位移響應(yīng)

對(duì)于圓柱結(jié)構(gòu)的渦激振動(dòng)問題，當(dāng)圓柱的渦旋脫落頻率fs與固有頻率fn接近時(shí)，一般認(rèn)為發(fā)生了“鎖定”現(xiàn)象。圖8給出了串聯(lián)兩圓柱的頻率比fs/fn隨約化速度Ur(表示來流速度和結(jié)構(gòu)振動(dòng)頻率之間的無量綱參數(shù))的變化關(guān)系。從圖8可判斷兩圓柱的“鎖定”區(qū)間為6≤Ur≤11。

圖8 頻率比隨Ur的變化關(guān)系Fig.8 Variation of frequency ratio with Ur

圓柱渦激振動(dòng)響應(yīng)隨Ur的變化曲線如圖9a所示。順流向平均位移與圓柱直徑之比X/D隨約化速度的變化曲線如圖9b所示。從圖9a可見：采用振幅均方根與圓柱直徑之比Y/D對(duì)橫流向振動(dòng)響應(yīng)進(jìn)行描述，兩圓柱的橫流向無量綱振幅Y/D隨Ur發(fā)生振蕩變化，觀察發(fā)現(xiàn)橫流向振幅可明顯分為初始分支、上端分支和下端分支3個(gè)響應(yīng)區(qū)間；當(dāng)進(jìn)入“鎖定”區(qū)間內(nèi)，下游圓柱的振幅迅速增大，當(dāng)Ur=10時(shí)下游圓柱的Y/D取得最大值0.79，此時(shí)被“完全鎖定”。因?yàn)轭l率比不再隨約化速度的變化繼續(xù)增大，漩渦之間相互融合使得漩渦強(qiáng)度得到提升，所以在流體與圓柱之間產(chǎn)生了劇烈的非線性動(dòng)力相互作用，致使圓柱結(jié)構(gòu)產(chǎn)生大幅度振動(dòng)，形成上端分支的響應(yīng)區(qū)間。從圖9b可見：兩圓柱的順流向平均位移X/D呈現(xiàn)增大趨勢(shì)；由于隨著Ur的增大，流經(jīng)圓柱結(jié)構(gòu)的來流速度相應(yīng)增加，圓柱在順流向方向上受到的沖擊作用越來越強(qiáng)，故順流向平均位移持續(xù)增大。

圖9 圓柱渦激振動(dòng)響應(yīng)隨Ur的變化曲線Fig.9 Variation of cylinder's vortex-induced vibration response with Ur

3.2 圓柱受力特性

圖10為串聯(lián)雙圓柱在渦激振動(dòng)過程中升力系數(shù)和阻力系數(shù)隨Ur的變化關(guān)系。

圖10 升力系數(shù)和阻力系數(shù)隨Ur的變化關(guān)系Fig.10 Variation of lift resistance coefficient with Ur

由圖10a可見，兩圓柱的脈動(dòng)升力系數(shù)隨Ur的變化過程總體相似，表現(xiàn)為先增后減再增的三段式進(jìn)程，這說明約化速度Ur的改變對(duì)脈動(dòng)升力系數(shù)有著顯著的影響。由圖10b可見，當(dāng)Ur較小時(shí)，上游圓柱平均阻力系數(shù)始終大于下游圓柱，且下游圓柱的阻力峰值小于上游圓柱，而在Ur>6之后，兩圓柱的平均阻力系數(shù)逐漸減小，下游圓柱阻力系數(shù)短暫大于上游圓柱。下游圓柱受到上游圓柱尾流干涉的影響，圓柱周圍的實(shí)際來流速度小于上游圓柱，因此下游圓柱的平均阻力系數(shù)在低約化速度時(shí)均小于上游圓柱，同時(shí)峰值也會(huì)小于上游圓柱。在“鎖定”區(qū)間內(nèi)，下游圓柱橫流向振幅迅速增大，明顯大于上游圓柱，使得在來流方向上兩圓柱已不在同一平面，上游圓柱無法繼續(xù)對(duì)下游圓柱產(chǎn)生較大的影響，所以下游圓柱平均阻力系數(shù)在此期間短暫大于上游圓柱。隨著Ur繼續(xù)增大，最終兩圓柱的平均阻力系數(shù)趨于穩(wěn)定。

3.3 圓柱運(yùn)動(dòng)軌跡

圖11為上游圓柱和下游圓柱的運(yùn)動(dòng)軌跡隨Ur的變化關(guān)系。

圖11 圓柱運(yùn)動(dòng)軌跡隨Ur的變化關(guān)系Fig.11 Variation of cylinder’s moving trajectories with Ur

從圖11可發(fā)現(xiàn)：上游圓柱在初始分支和上端分支響應(yīng)區(qū)間的運(yùn)動(dòng)軌跡以“8”字形、“新月”形封閉曲線為主，周期性并不穩(wěn)定，這說明當(dāng)間距比L/D=4時(shí)下游圓柱對(duì)上游圓柱振動(dòng)模態(tài)的影響作用比較有限，僅引入些許的隨機(jī)性擾動(dòng)；當(dāng)Ur繼續(xù)增大，在離開“鎖定”區(qū)間后，橫流向振幅迅速減小，上游圓柱的運(yùn)動(dòng)軌跡逐漸向“橢圓”形轉(zhuǎn)變；下游圓柱的軌跡形狀隨Ur的增大在“橢圓”形、“8”字形和“新月”形等之間轉(zhuǎn)換，運(yùn)動(dòng)軌跡相比上游圓柱更加復(fù)雜。由于下游圓柱受到上游圓柱的遮蔽干擾以及尾流干涉作用的影響，故其呈現(xiàn)出不穩(wěn)定的非周期性運(yùn)動(dòng)軌跡，此外，約化速度也對(duì)圓柱的運(yùn)動(dòng)有明顯的影響。

3.4 圓柱尾渦脫落模式

圖12為圓柱尾渦脫落模式隨Ur的變化關(guān)系。

圖12 尾渦脫落模式隨Ur的變化關(guān)系Fig.12 Variation of wake shedding mode with Ur

從圖12可以看出，當(dāng)Ur<5時(shí)，上游圓柱無明顯渦脫現(xiàn)象，因?yàn)榇藭r(shí)來流速度較小，在圓柱兩側(cè)未能產(chǎn)生周期性漩渦脫落，所以串聯(lián)圓柱的尾流在下游圓柱后方匯集形成一條渦帶，使得圓柱系統(tǒng)尾流區(qū)內(nèi)的漩渦以圖12a的模式進(jìn)行整體脫落。隨著Ur不斷增大，從圖12c可見，Ur=8時(shí)尾渦脫落呈現(xiàn)“2S”模式，即每個(gè)周期內(nèi)泄放兩個(gè)漩渦。當(dāng)進(jìn)入“鎖定”區(qū)間內(nèi)，上游圓柱后方開始出現(xiàn)漩渦脫落，下游圓柱的漩渦脫落長(zhǎng)度同時(shí)增加，上游圓柱分離的剪切層與下游圓柱側(cè)面發(fā)生碰撞，兩圓柱的漩渦在下游圓柱表面發(fā)生融合，使得漩渦強(qiáng)度提升，尾渦脫落模式因此發(fā)生變化。從圖12d可見，Ur≥12時(shí)脫落模式由“2S”轉(zhuǎn)變?yōu)椤?P”。從上游圓柱一側(cè)脫落的漩渦與下游圓柱表面發(fā)生碰撞，撞擊后的漩渦一部分與下游圓柱同側(cè)漩渦融合，另一部分跟隨下游圓柱反向脫落的漩渦一起脫落，即每個(gè)周期內(nèi)從下游圓柱兩側(cè)各釋放一對(duì)逆向旋轉(zhuǎn)的漩渦，最終導(dǎo)致尾渦脫落模式再次發(fā)生轉(zhuǎn)變。

4 結(jié) 論

(1)Overset-Mesh技術(shù)操作方便，在網(wǎng)格運(yùn)動(dòng)期間能夠保持很高的網(wǎng)格質(zhì)量，可以較好地用于計(jì)算流體力學(xué)工程實(shí)際中。

(2)下游圓柱的“鎖定”區(qū)間與上游圓柱保持一致，并未受到流動(dòng)干涉的影響而產(chǎn)生后移現(xiàn)象，說明在間距比L/D=4條件下，上游圓柱對(duì)下游圓柱的尾流干涉作用相對(duì)較弱。

(3)在約化速度Ur=6附近存在一個(gè)臨界值Ur-crit，當(dāng)Ur>Ur-crit時(shí)，下游圓柱的橫流向振幅大于上游圓柱，而當(dāng)Ur

(4)在“鎖定”區(qū)間內(nèi)，上游圓柱脫落的漩渦與下游圓柱發(fā)生碰撞，同性質(zhì)漩渦融合使得漩渦強(qiáng)度增大，尾渦脫落模式因此轉(zhuǎn)變，下游圓柱振幅相應(yīng)增大。串聯(lián)圓柱的尾渦脫落模式隨Ur增大會(huì)發(fā)生顯著變化，進(jìn)而深刻影響圓柱的運(yùn)動(dòng)響應(yīng)和受力情況。