尾流干涉下串聯圓柱渦激振動數值模擬研究*

王曉凱 婁 敏

(中國石油大學(華東)石油工程學院)

0 引 言

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

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

1 數值方法

1.1 Overset Mesh

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

圖1 重疊網格構成圖Fig.1 Composition of overset mesh

利用重疊網格技術可以處理具有小縫隙部件的相對網格運動，與傳統網格相比，重疊網格對網格拓撲要求低，且網格生成相對容易，同時可以避免傳統動網格技術所面臨的網格動態更新時易出現的負體積現象，在網格運動期間始終能夠保持很高的網格質量，實現了局部結構網格在非結構網格中的使用。Overset Mesh在計算流體力學工程實際應用中的流程如圖2所示。

圖2 重疊網格在數值計算中的流程Fig.2 Numerical calculation process of overset mesh

1.2 控制方程

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

U=0

(1)

(2)

漩渦泄放對圓柱體作用力如圖3所示。由圖3可見，流體流經圓柱時會在其兩側產生漩渦泄放現象，在橫流向方向產生升力FL造成橫向運動，同時在順流向方向產生脈動阻力FD造成順流向的運動。升力與脈動阻力的表達式分別為：

(3)

(4)

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

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

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

圖4 圓柱振動模型Fig.4 Cylindrical vibration model

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

(5)

(6)

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

2 數值模型

2.1 模型參數

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

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

圖5 流場區域計算模型Fig.5 Computational model of flow field area

圖6 網格劃分Fig.6 Mesh generation

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

2.2 模型驗證

為檢驗上述計算模型與網格劃分的可靠性，當雷諾數為200時，分別進行單圓柱繞流和間距比為4的串聯雙圓柱繞流數值模擬，得到相應的升力和阻力系數，模擬結果如圖7所示。

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

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

(7)

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

表2 單圓柱繞流結果對比Table 2 Comparison between simulation results of flow around a single cylinder

表3 雙圓柱繞流結果對比Table 3 Comparison between simulation results of flow around two cylinders

3 計算結果

3.1 圓柱位移響應

對于圓柱結構的渦激振動問題，當圓柱的渦旋脫落頻率fs與固有頻率fn接近時，一般認為發生了“鎖定”現象。圖8給出了串聯兩圓柱的頻率比fs/fn隨約化速度Ur(表示來流速度和結構振動頻率之間的無量綱參數)的變化關系。從圖8可判斷兩圓柱的“鎖定”區間為6≤Ur≤11。

圖8 頻率比隨Ur的變化關系Fig.8 Variation of frequency ratio with Ur

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

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

3.2 圓柱受力特性

圖10為串聯雙圓柱在渦激振動過程中升力系數和阻力系數隨Ur的變化關系。

圖10 升力系數和阻力系數隨Ur的變化關系Fig.10 Variation of lift resistance coefficient with Ur

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

3.3 圓柱運動軌跡

圖11為上游圓柱和下游圓柱的運動軌跡隨Ur的變化關系。

圖11 圓柱運動軌跡隨Ur的變化關系Fig.11 Variation of cylinder’s moving trajectories with Ur

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

3.4 圓柱尾渦脫落模式

圖12為圓柱尾渦脫落模式隨Ur的變化關系。

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

從圖12可以看出，當Ur<5時，上游圓柱無明顯渦脫現象，因為此時來流速度較小，在圓柱兩側未能產生周期性漩渦脫落，所以串聯圓柱的尾流在下游圓柱后方匯集形成一條渦帶，使得圓柱系統尾流區內的漩渦以圖12a的模式進行整體脫落。隨著Ur不斷增大，從圖12c可見，Ur=8時尾渦脫落呈現“2S”模式，即每個周期內泄放兩個漩渦。當進入“鎖定”區間內，上游圓柱后方開始出現漩渦脫落，下游圓柱的漩渦脫落長度同時增加，上游圓柱分離的剪切層與下游圓柱側面發生碰撞，兩圓柱的漩渦在下游圓柱表面發生融合，使得漩渦強度提升，尾渦脫落模式因此發生變化。從圖12d可見，Ur≥12時脫落模式由“2S”轉變為“2P”。從上游圓柱一側脫落的漩渦與下游圓柱表面發生碰撞，撞擊后的漩渦一部分與下游圓柱同側漩渦融合，另一部分跟隨下游圓柱反向脫落的漩渦一起脫落，即每個周期內從下游圓柱兩側各釋放一對逆向旋轉的漩渦，最終導致尾渦脫落模式再次發生轉變。

4 結 論

(1)Overset-Mesh技術操作方便，在網格運動期間能夠保持很高的網格質量，可以較好地用于計算流體力學工程實際中。

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

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

(4)在“鎖定”區間內，上游圓柱脫落的漩渦與下游圓柱發生碰撞，同性質漩渦融合使得漩渦強度增大，尾渦脫落模式因此轉變，下游圓柱振幅相應增大。串聯圓柱的尾渦脫落模式隨Ur增大會發生顯著變化，進而深刻影響圓柱的運動響應和受力情況。