亞臨界雷諾數下單圓柱渦激振動研究

周力，邱中秋，袁亞帥，宗智

大連理工大學 船舶工程學院，遼寧 大連 116024

0 引 言

海洋立管結構大多以柔性管形式存在，在洋流的作用下自由振動。然而由于海底地形復雜多變，鋪設在海底的某段立管兩端容易受到約束而形成一定的跨度，在洋流作用下振動并產生破壞。研究中通常把這段立管作為彈性安裝的剛性管考慮，利用自激振動或強迫振動的方法對其相關水動力特性進行探討[1-2]。

Williamson 等通過強迫振動實驗研究對振動立管尾渦脫落的結構模式進行了分類，并依據渦脫模式劃分了不同的區域[3]。Peppa 等研究分析了低雷諾數下來流和立管結構之間的能量傳遞關系，發現在較低無量綱振幅比時尾渦脫落模式為“2S”（S 表示單個旋渦脫落），較高無量綱振幅比時尾渦脫落模式逐漸變得復雜[4]。Meneghini 等將立管的振動頻率和幅值作為變量，對Re＝200 時的橫向振動立管尾渦脫落頻率與“鎖定”現象進行了研究，確定了該雷諾數下發生“鎖定”時的最大振幅值及振動頻率的變化范圍[5]。Morse 等[6-7]通過圓柱的強迫振動和自激振動組合實驗發現了新的尾渦脫落模式，位于“2S”和“2P”模式的過渡區（P 表示一對旋轉方向相反的旋渦脫落），并將其定義為“2P0”模式。

王凱鵬等[8]基于緊致插值法對均勻流(Re=200)中圓柱橫向強迫振動問題進行了二維數值模擬研究。朱永健等[9]通過對定常流中橫向振動圓柱的升力突變現象研究，發現隨著強迫振動振幅的增加，圓柱脈動升力系數出現了突變衰減，圓柱靜止情況下的渦脫頻率與受迫振動頻率控制的波動出現了相位逆變，使得圓柱與流體的能量傳遞出現了逆向改變。

從海洋工程應用角度來看，彈性支撐立管強迫振動研究的主要目的是揭示自激振動更深層次的規律，為海洋立管渦激振動預測模型提供可靠的實驗數據[10]。事實上，大多數渦激振動（vortexinduced vibration，VIV）預測程序是基于強迫振動實驗數據形成數據庫后建立相應的預報模型，以滿足工程應用需要，如VIVANA 和SHEAR7[11-12]。

自激振動能更直觀地觀察圓柱受到流體力作用發生VIV 時的現象，但強迫振動能更深層次地揭示流體力與振動相互作用下的水動力特性[13]。強迫振動和自激振動的流場變化及振動響應存在很大的差異，目前鮮有文獻對2 類實驗的內部聯系進行深入研究?；诖?，本文將結合強迫振動數值計算與自激振動實驗2 種方法，將強迫振動的動力響應、渦脫模式的轉變與該圓柱自激振動出現“鎖定”時的最大振幅聯系起來，期望為深入討論二者之間的聯系提供參考。

1 數值模型

1.1 控制方程

以二維不可壓縮黏性流體模擬流場運動，其控制方程Navier-Stokes 方程組在坐標系下表示為：

本文采用的求解器為 STAR-CCM+ 軟件中的Realizablek-ε 湍流模型。標準k-ε 湍流模型的主要優點是計算快速、穩定，計算結果合理，適用于高雷諾數的流動，但不建議用于分離流動。Realizablek-ε 湍流模型是標準k-ε 湍流模型的一個改進模型，可以較好地模擬帶有強壓力梯度的邊界層流動和分離流問題，由Shih 等[14]提出，其流湍動能及其耗散率輸運方程為：

1.2 計算域和邊界條件

STAR-CCM+ 軟件可以有效地實現圓柱在均勻來流下發生橫向強迫振動的模擬。首先對圓柱作減運算，對包含減運算在內的重疊網格施加不同的場函數，可實現模擬圓柱在不同頻率和振幅下的運動。流經圓柱的均勻流速度為U0。

計算域如圖1 所示，來流速度方向與x軸平行，圓柱橫向振動方向與y軸平行，圓柱直徑D=0.05 m，圓柱中心上下兩側20D處為對稱邊界；左端距離圓柱中心20D處為速度進口，右端距離圓柱中心50D處為壓力出口；圓柱表面為無滑移壁面條件。

圖1 圓柱受迫振動計算域示意圖Fig. 1 Schematics of the forced-oscillation computational domain

圓柱垂直于來流方向沿y軸做橫向強迫振動，瞬時位移：

式中：A為圓柱橫向振動振幅；fe為圓柱強迫振動的激振頻率。

圓柱強迫振動過程中沿y軸方向的瞬時速度：

圓柱開始振動時的瞬時速度vy(0)=2πfeA。

2 計算可行性驗證

2.1 固定圓柱驗證計算

網格劃分及其數量會直接影響到計算結果的準確性。為此，采用多邊形和棱柱層的網格劃分方法，對圓柱橫向振動的延展區域和尾部區域進行了局部加密，在圓柱周圍區域設置包裹圓柱運動的運動重疊網格，重疊區域大小為2D×5D。網格劃分如圖2 所示，計算模型是具有較高精度的結構化網格。采用壁面函數法進行圓柱體近壁面處理，圓柱體表面第1 層網格滿足y+≈1 條件。由Re=20 000，計算得到來流速度U0=0.356 m/s。流體的密度取為997 kg/m3，其動力黏度取為0.895×10?3kg/(s·m)，則可計算得到緊鄰壁面第1 層網格的厚度為?s=0.04 mm。設置棱柱層數n=10，增長率q=1.3，由δ=(?s(1?qn))/(1?q)，則可計算得到邊界層總厚度為1.67 mm。

圖2 固定網格劃分示意圖Fig. 2 Schematics of fixed mesh division

對網格、時間步長進行了無關性驗證，在保證計算精度的同時，應節約計算成本，最終選取網格數為39 萬，時間步長為0.01 s。為了驗證計算模型的可靠性，對雷諾數Re=20 000 均勻流中固定圓柱和振動圓柱的繞流問題進行了數值計算，如圖3 所示，通過觀察升力系數的能量譜密度曲線，確定圓柱渦街脫落的特征頻率為1.46 Hz。

圖3 Re=20 000 固定圓柱升力系數能量譜密度Fig. 3 Energy spectral density of lift coefficient of fixed cylinderwith Re=20 000

根據公式：

其中，fov為圓柱尾渦脫落頻率，可得到該固定圓柱在Re=20 000 時的斯特勞哈爾數St=0.204 8，與文獻[15]中St=0.2 基本保持一致。

2.2 圓柱強迫振動計算

在相同來流條件下對圓柱進行強迫振動數值計算，將振幅比A/D=0.5，強迫振動頻率比范圍0.5≤fe/fs≤1.5（fs為斯特勞哈爾頻率）時的數值計算結果與喻晨欣等[16]的數值結果和Sarpkaya[17]的實驗結果進行對比驗證。從圖4 可以看出，本文計算結果與文獻[16-17]的最大升力系數幅值變化趨勢保持一致，隨著頻率比的增大，均呈現不斷增大的趨勢，并在fe/fs=1 附近出現了跳躍[18-19]。由于本文與參考工況對應的雷諾數存在差異，對應的頻率比也有所不同，導致結果存在偏差，與實際情況相符。以上算例驗證了本文模型可用于求解亞臨界雷諾數下圓柱的強迫振動問題。

圖4 A/D=0.5 時最大升力系數Clmax 隨fe/fs 變化曲線Fig. 4 Clmax variation versus fe/fs with A/D=0.5

3 計算結果分析

3.1 強迫振動數值結果分析

本文采用外加激勵的方法，使圓柱在均勻來流下做類似自激振動的運動。通過控制振幅A和激振頻率fe，可以實現使圓柱位移時域曲線呈簡諧運動的規律，無量綱振幅比A/D和激振頻率比fe/fn如表1 所示（fn為結構固有頻率，考慮附加水質量）。渦激振動達到“鎖定”時對應的振幅存在閾值[3]，最小振幅比取A/D=0.3。通過自由衰減振動實驗測得fn=0.714 2 Hz。

表1 無量綱頻率比及振幅比范圍Table 1 Non-dimensional frequency ratio and amplitude ratio

考慮到強迫振動過程中圓柱受到流體力和激振力的雙重作用，假定在數值計算中得到的升力為Fm，完全由流體力提供的升力為Fl，強迫振動對圓柱升力的影響用慣性力Fi表示，三者關系可用式（9）來表示：

當流體力Fl與慣性力Fi方向相同時，Fm>Fl；當Fl與Fi方向相反時，Fm

圖5 為頻率比fe/fn=0.83 時不同振幅比下，單周期內圓柱振動速度vy曲線及相應升力系數Cl曲線。從圖中可以看到，振動速度曲線始終遵循簡諧運動規律，但由于受到強迫振動帶來的慣性力影響，單周期內升力系數曲線不再規則變化。當振幅比A/D足夠大時，可觀察到升力系數曲線單周期內包含多個峰值。

圖5 單周期內振動速度及升力系數變化曲線Fig. 5 Oscillating speed and lift coefficient varies within a single period

同時可以發現，隨著振幅比A/D的增大，在圓柱上升階段（0～0.5T），最大振動速度（0.25T）對應的升力系數值會由正值轉變為負值；而在圓柱下降階段（0.5T～T），最大振動速度（0.75T）對應的升力系數值由負值轉變為正值。強迫振動研究中通常取圓柱固定（A/D=0）時的升力系數為正值，所以在接近A/D=0 一端，升力系數應取正值。在本文圓柱強迫振動數值計算中，將在上升階段圓柱振動速度達到最大時（0.25T）的位置選為升力系數監測位置。

圖6 為VIV 預測程序SHEAR7 中的保守模型，各無量綱頻率比下的Cl-A/D曲線通常由強迫振動實驗所得數據點（0，Cl0），（（A/D）*，Clmax）、（（A/D）max，0）通過二次函數擬合得到，其中（A/D）*和（A/D）max分別表示升力系數為零和最大值時對應的無量綱振幅比[19]。該模型呈現升力系數Cl隨著振幅比A/D的增大先上升后下降、由正升力系數逐漸轉變為負升力系數的趨勢。

圖6 升力系數隨無量綱振幅比變化曲線Fig. 6 Variation of lift coefficient with dimensionless A/D

流體力具有增大圓柱振幅趨勢的作用，而振動則會限制圓柱的振幅繼續增大[20]。在兩者的相互作用下，當升力系數為正時，說明振動對流體力的限制能力有限，此時若繼續增大來流速度，立管振幅仍有增大趨勢，因此渦激振動未達到“鎖定”狀態時圓柱振幅會隨著來流速度的增大而增大。隨著振幅比的增大，振動對流體力的限制能力也進一步增強，直至平衡位置處的升力系數為零，二者達到平衡，振幅達到最大且不再有增大的趨勢。值得一提的是，在渦激振動中，振動的限制作用不會超過流體力的激勵作用，即升力系數不可能取負值，所以達到“鎖定”狀態時，在一定流速范圍內，即使增加流速，圓柱的振幅仍能保持最大值不變。

圓柱在外加激勵作用下發生強迫振動，通過改變振幅比A/D，觀察監測位置處的升力系數，可擬合出各激振頻率比fe/fn下的升力系數隨振幅比A/D的變化曲線。對各工況下數值計算所得結果進行三階多項式擬合，并控制擬合條件Cl0=0.168（Re=20 000 時固定圓柱的升力系數），可得擬合曲線。在VIV 預測程序SHEAR7 中輸入相關參數條件，同樣可得預報的擬合曲線。將同一頻率比下的2 種擬合曲線進行對比分析，結果如圖7所示。

由圖7 可見，各激振頻率下數值計算結果擬合曲線與SHEAR7 擬合曲線的變化趨勢高度吻合，均呈現出升力系數隨著振幅比A/D增大先上升后下降的變化規律。傳統強迫振動計算中升力系數的取值方法往往未呈現該規律，說明選取上升階段最大振動速度時的升力系數，對于建立VIV 預測模型是適合的。

為了進一步了解數值計算結果擬合曲線與渦激振動振幅響應之間的聯系，從頻率比fe/fn=1 附近的零升力系數點著手分析。如圖8 所示，可觀察到圓柱平衡位置處的零升力系數點均保持在振幅比A/D=0.8 附近。通過類比SHEAR7 預測模型進行分析，可知在亞臨界雷諾數下該圓柱渦激振動出現“鎖定”狀態時最大振幅應保持在0.8D左右。由此可對該圓柱在自激振動頻率比fov/fn≈1時的最大響應振幅做出預測，其中fov為自激振動時圓柱的尾渦脫落頻率。

圖8 各工況下Cl 隨振幅比A/D 變化曲線Fig. 8 Variation of lift coefficient with A/D under all conditions

同時由圖8 可見，在同一振幅比下，隨著頻率比的增大，升力系數逐漸增大。當頻率比fe/fn=1.42 時，相對于其他組，升力系數大幅升高，發生了突增現象[12]，相關機理有待進一步探討。

3.2 自激振動實驗對比分析

強迫振動數值計算結果表明，監測位置處升力系數為零時對應的振幅比均保持在A/D=0.8 附近。在自激振動工況下，預測彈性安裝時該圓柱在亞臨界雷諾數下發生渦激振動的最大振幅為0.8D左右，且此時處于“鎖定”狀態（自激振動頻率比fov/fn≈1）。為了驗證預測結果，設計并進行了圓柱自激振動實驗。

圖9 為自激振動實驗裝置圖，振動圓柱嵌合在振動框架中，并彈性安裝于固定裝置。彈簧具有一定的的彈簧系數，振動框架可帶動圓柱沿兩側的滑軌做垂向往復運動。為了盡量減小三維效應，固定圓柱兩端的U 型架端口盡可能采用厚度較薄、質量較輕的端板。實驗所采用圓柱的直徑與數值計算中保持一致，彈簧系數K=230.48 N/m，均勻來流速度的取值范圍為0.13～0.34 m/s，對應的約化速度Ur取值范圍為3.7～9.4 m/s。

圖9 圓柱渦激振動實驗裝置Fig. 9 Experimental device setup of cylindrical vortex-induced vibration

通過改變來流速度，監測得到圓柱振動幅度及相應的自激振動頻率比fov/fn，可觀察到振動圓柱發生“鎖定”時的最大振幅，以及對應的頻率比范圍。圖10 為自激振動實驗中振幅比A/D和頻率比fov/fn隨約化速度Ur的變化情況。由圖可見，在fov/fn接近1 時，振幅比大幅升高；且在一定的fov/fn范圍內，振幅比A/D不發生大的改變，說明發生了“鎖定”現象。同時可觀察到“鎖定”范圍對應的約化速度使圓柱的最大振幅達到0.8D左右。

圖10 圓柱振幅比及自激振動頻率比隨約化速度變化圖Fig. 10 Changes in the amplitude ratio and self-induced vibration frequency ratio of the cylinder at different reduced speeds

3.3 渦脫模式轉變分析

Williamson 等[3]對振動圓柱尾流中漩渦脫落模式進行了實驗研究，并按尾渦脫落時的結構形式劃分了不同的區域，發現“鎖定”區域附近的主要渦脫落模式為2S，2P 和P+S。隨著無量綱振幅比的增加，在渦脫模式轉變臨界線附近，通常會從2S 模式轉變為P+S 或2P 模式。實驗研究表明，在低雷諾數下（Re<300），P+S 模式會取代2P模式，以至于整個區域只出現P+S 模式，而未出現2P 模式。而在較大雷諾數時，P+S 和2P 模式之間的邊界與臨界線基本一致?？梢姡S著雷諾數的改變，同一區域內渦脫模式發生轉變后不確定呈現2P 模式或P+S 模式，但轉變臨界線基本保持不變。

在同一強迫振動頻率比下，選取振幅比范圍為A/D=0.6～1.0 內同一時刻的渦量圖。如圖11 所示，各頻率比工況中的振幅比由A/D=0.6 至A/D=1.0 依次增大。可以發現，在振幅比小于0.8 時，如A/D=0.7 時紅色線圈所標記渦對，圓柱尾渦脫落模式為P+S 模式；當振幅比大于0.8 時，如A/D=0.9 時紅色線圈所標記渦對，圓柱尾渦脫落模式為2P 模式。說明在圓柱強迫振動中A/D=0.8 附近是尾渦脫落模式發生轉變的臨界線，監測位置零升力系數點與尾渦脫落模式轉變同時發生。

圖11 各頻率比下渦量隨振幅比變化圖Fig. 11 Variation of vorticity with amplitude ratio at various frequency ratios

4 結 論

本文研究了圓柱強迫振動時升力系數響應值的選取位置，并進行相關機理分析，對圓柱發生渦激振動“鎖定”時的最大振幅進行了預報和對比驗證，對平衡位置升力系數為零時的渦脫模式轉變進行了探討。通過分析，得到如下結論：

1） 亞臨界雷諾數下圓柱強迫振動時的升力系數響應與自激振動相比存在較大差異，當強迫振動中圓柱速度達到最大時，監測到的升力僅由流體力提供，將圓柱在上升階段的振動速度達到最大時（0.25T）對應的升力系數作為計算結果是合理的。

2） 數值模擬計算結果擬合曲線與VIV 預測程序SHEAR7 的相關擬合曲線吻合較好，且變化規律保持高度一致，升力系數均隨著振幅比的增大，呈現出先增大后減小的趨勢，提供了一種以強迫振動數值計算結果建立Cl-A/D模型用于渦激振動響應振幅預測的方法。

3） 數值計算結果表明，圓柱在強迫振動過程中平衡位置處的零升力系數點均保持在振幅比A/D=0.8 附近，采用自激振動實驗方法對預測結果進行了對比驗證。實驗結果表明，當自激振動頻率比fov/fn接近于1 時，振幅比大幅升高，且在一定的自激振動頻率比fov/fn范圍內，振幅比A/D沒有較大改變，說明發生了“鎖定”現象。同時確定了圓柱在“鎖定”時對應的最大振幅為0.8D左右，證明了該預報方法具有可行性。

4） 對平衡位置處的零升力系數與相應的圓柱尾渦脫落模式進行了探討。計算結果表明：在振幅比小于0.8 時，圓柱尾渦脫落模式為P+S 模式；當振幅比大于0.8 時，圓柱尾渦脫落模式為2P 模式。強迫振動圓柱Cl-A/D擬合曲線零升力系數對應的振幅比與圓柱在渦激振動中“鎖定”時的最大響應振幅比基本一致，且圓柱尾渦脫落模式在此振幅比下發生了轉變。

未來可通過對不同雷諾數工況下的圓柱振動進行計算，對圓柱在均勻來流下發生“鎖定”的頻率比區間進行預報。