串列不等直徑雙圓柱海流能發電振子渦激振動數值模擬

白 旭， 陳 云

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

0 引言

渦激振動是一種常見的流固耦合現象，即流體經過圓柱體時會產生周期性脈動力，進而引發柱體振動，當圓柱體的自振頻率與旋渦脫落頻率相近時，兩者會發生共振，從而產生較大振幅。 研究表明，利用圓柱體渦激振動產生的振幅進行發電，能夠減少不可再生能源消耗量，具有重要的現實意義[1]。

渦激振動的研究多采用數值模擬和試驗分析的方法，在試驗分析中，雷諾數Re≥104，在數值模擬中，雷諾數Re≤102。 有研究表明，渦脫落過程產生的大尺寸尾流在一定程度上呈高、低雷諾數（102量級）相似的特征[2]。 高雷諾數的主要現象和特征可以從低雷諾數問題的研究中得到。 在渦激振動發電裝置的研究方面，美國密歇根大學發 明 的VIVACE （Vortex Induced Vibration for Aquatic Clean Energy） 能量轉換裝置的發電原理是將圓柱渦激振動產生的位移通過變速裝置帶動發電機產生電能[3]。 研究表明，在獲取能量方面，多柱體較單柱體有較大優勢[4]。 針對渦激振動多柱體的研究主要從雙圓柱入手，進而揭示多柱體的振動規律。

目前， 針對雙圓柱渦激振動的研究主要集中在等直徑串列、并列雙圓柱旋渦脫落形式、振動特性等方面。 對于直徑不等的多圓柱繞流及渦激振動現象的研究相對較少。 Wang H[5]對上、下游圓柱直徑比為2 的串列圓柱體進行了數值模擬（50≤Re≤200），其中上游圓柱靜止，下游圓柱做垂直流向的單自由度運動或同時做順流向、 垂直流向的雙自由度運動，研究結果表明，下游圓柱做垂直流向的單自由度運動時產生的振幅比單圓柱振動和兩個等直徑圓柱體同條件振動產生的振幅均要大。 Lam K M[6]對并列布置、串列布置和交錯布置方式下的上、 下游圓柱直徑比為2 的雙圓柱進行了試驗研究，其中上游大圓柱固定，下游小圓柱彈性支撐， 研究結果表明， 在串列及并列布置方式下， 彈性支撐的下游小圓柱的流體誘發振動會受到抑制，而在交錯布置方式下，流體誘發振動會被放大。Wang Y T[7]對串列不等直徑雙圓柱進行了數值模擬（100≤Re≤150），模擬結果表明，直徑比（d/D）和間距比（G/D）對升、阻力系數、圓柱周圍的壓力分布和渦旋脫落頻率均有重要影響。

綜上可知，針對不等直徑圓柱體渦激振動的研究主要集中于雙圓柱靜止繞流或上游圓柱靜止，下游圓柱做單自由度運動的情況，而有關不等直徑雙圓柱均彈性支撐同時做單自由度運動的研究較少。因此，本文計算分析了上、下游圓柱彈性支撐同時沿垂直流向做單自由度振動的情況，旨在觀察該工況下間距比和直徑比對雙圓柱體渦激振動響應特性及尾渦特征的影響，從而為海流能發電振子系統的直徑選取及布置提供理論依據。

1 計算模型

本文計算條件為雷諾數Re=300，位于亞臨界雷諾數區間，該范圍內的尾渦以固定頻率周期進行脫落，便于觀察分析。 涉及的計算參數：來流速度U=0.14 m/s，空氣密度ρ=1.225 kg/m3，流體的動力黏度系數υ 為1.48×10-5m2/s，上、下游圓柱的間距為G，間距比G/D=2～6，上、下游圓柱的質量之比均為2.4， 阻尼比ζ=0.001， 上游圓柱的直徑D=0.032 m，下游圓柱的直徑d=0.2D～0.8D。

圖1 為計算域及邊界示意圖。 圓柱圓心距離進口邊界為15D，距離上、下邊界均為15D。 進口采用速度入口，出口采用壓力出口，上、下邊界采用對稱邊界，圓柱為壁面邊界。 限定上、下游圓柱均做垂直流向單自由度振動。

2 數值計算方法

2.1 控制方程

不可壓縮粘性牛頓流體的控制方程為Navier-Stokes 方程。 笛卡爾坐標系中連續性方程和動量方程分別為

式中：u=（Ui，υ）T為流體的速度矢量；Ui為某時刻的速度分量，m/s；P 為壓力，Pa。

柱體振動方程可以表示為

2.2 重疊網格技術

在重疊網格方法中， 可以使用兩種類型的網格 （可以針對不同區域獨立生成） 來離散流場。第一種類型的網格是主要網格（背景網格），它是一個全局笛卡爾網格（整個計算域的H 型網格），該網格固定在空間中。 第二種類型的網格是與主要網格重疊， 并可以靜止或移動的次要網格。對于每個網格分別求解控制方程，并且使用插值算法在主網格和次網格之間有效地傳輸和交換數據[8]。

在本文中，次要網格（O 型網格） 覆蓋了圓柱周圍的區域。 網格之間的數據轉換是通過被稱為邊緣點的特定點來完成的。主要網格（背景網格） 的邊緣點是放置在圓柱體外部和其表面附近的節點， 而次要網格的邊緣點是放置在每個次要網格外部邊界線上的節點。 另一組點稱為孔點，在網格中切割并與固體表面相交。為了將信息從網格中的邊緣點傳輸到其他網格中的節點， 本文采用搜索算法指定對象網格中該邊緣點的周圍節點[9]，[10]。

圖2 離散化的計算網格Fig.2 Discretized computational grid

2.3 可靠性驗證

為了驗證本文模型的可靠性， 將利用本文模型計算得到的單圓柱及串列雙圓柱靜止繞流的結果與文獻結果進行對比。

2.3.1 單圓柱驗證

為了驗證程序的正確性， 對Re=200 情況下的單圓柱勻速繞流進行了數值模擬。 計算得到的升力系數Cl和斯托羅哈爾數St 如表1 所示。 由表1 可以看出， 本文的計算結果與前人的計算結果吻合較好，驗證了程序的正確性。

表1 單圓柱繞流結果對比Table 1 Comparison of single cylinder flow results

2.3.2 等直徑串列雙圓柱驗證

為進一步驗證計算程序的可靠性， 本文對串列等直徑雙圓柱繞流進行了數值模擬（G/D 分別為1.5，2，3 和4）。 計算得到的阻力系數Cd1、阻力系數Cd2和斯托羅哈爾數St 如表2 所示。 由表2可以看出， 本文的計算結果與前人的計算結果非常接近。 因此，驗證了本文網格劃分、時間步長設置和數值求解格式等程序的有效性和準確性。

表2 雙圓柱繞流結果對比Table 2 Comparison of the results of flow around a double cylinders

3 結果與分析

3.1 振幅響應

選取每個周期內振幅的均方根（Yrms）作為振幅幅值， 利用振幅比A*=Yrms/D 來分析振幅響應。 在不同間距比下，上、下游圓柱的振幅響應隨直徑比的變化曲線如圖3 所示。

圖3 不等直徑比下振幅響應示意圖Fig.3 Schematic diagram of amplitude response under unequal diameter ratio

從圖3 可以看出： 上游大圓柱和下游小圓柱的振幅比最大值均在G/D=2 時取得，其中上游大圓柱在G/D=2，d/D=0.6 時取得振幅比A*的最大值（約為0.39D）；下游小圓柱在G/D=2，d/D=0.2 時取得振幅比A*的最大值 （約為3.2 d）。 在小間距比下，圓柱間的干涉作用較強，上游圓柱的自由剪切層包圍整個下游圓柱， 增強了下游小圓柱的振幅響應； 當d/D 大于0.6，G/D=2 時，隨著直徑比的逐漸增大，上游大圓柱的振幅比A*逐漸減小；在其他間距比下，上游大圓柱的振幅比A*均隨著直徑比（d/D＞0.6）的增大而緩慢增加；隨著直徑比的逐漸增大，下游小圓柱的渦激振動阻力逐漸增大， 不同間距比下的振幅比A*均逐漸減小。在串列不等直徑雙圓柱系統中，在小間距比下，直徑比的變化對系統振幅響應的影響較大。

3.2 受力特性

3.2.1 升力特性

在不同間距比下，上、下游圓柱的升力系數均方根Clrms 隨直徑比的變化曲線如圖4 所示。

圖4 不同直徑比下升力系數均方根的變化曲線Fig.4 Root mean square variation curves of lift coefficient under unequal diameter ratio

從圖4（a）可以看出：對于上游大圓柱，當G/D=2，直徑比較小時，Clrms 較小，隨著直徑比的逐漸增加，Clrms 變化明顯， 并在d/D=0.6 后出現陡增趨勢，說明耦合不等直徑的雙圓柱振動時，圓柱間有較高流速的射流，隨著直徑比的逐漸增加，下游圓柱對上游圓柱的影響逐漸增強； 當G/D=3，4，6 時，隨著直徑比的增加，Clrms 變化緩慢，且Clrms 隨著G/D 的增大而增大； 當G/D=6 時，Clrms 的變化趨近于單圓柱的振動情況。 從圖4（b）可以看出：對于下游小圓柱，當d/D＜0.35 時，Clrms 在G/D=3 時取得最大值； 當d/D＞0.35 時，G/D=4 下的Clrms 始終大于其他間距比下的Clrms；隨著直徑比的逐漸增大，下游小圓柱承受的來自上游圓柱形成的脫落渦的沖擊力逐漸增大， 上游大圓柱對下游小圓柱的干擾作用逐漸增強，下游小圓柱的Clrms 逐漸增大。 在不同間距比（G/D=2 除外）下，下游小圓柱的Clrms 明顯高于上游大圓柱的Clrms。

3.2.2 阻力特性

在不同間距比下，上、下游圓柱的阻力系數均方根Cdrms 隨直徑比的變化曲線如圖5 所示。

圖5 不同直徑比下阻力系數均方根的變化曲線Fig.5 Root mean square variation curve of resistance coefficient under unequal diameter ratio

從圖5（a）可以看出：對于上游大圓柱，當G/D 一定時，隨著直徑比的逐漸增大，Cdrms 出現較強的波動性，這說明當間距比一定時，下游小圓柱的直徑變化對上游大圓柱的影響較大；當G/D=3，4，6 時，隨著直徑比的逐漸增大，Cdrms 的變化趨勢相似，但與G/D=2 時的變化趨勢相反；當直徑比較小時， 上游圓柱的自由剪切層重附在下游小圓柱上，下游圓柱處在上游圓柱的尾渦范圍內，對上游圓柱的影響不明顯。 從圖5（b）可以看出：對于下游小圓柱，在不同的間距比下，Cdrms 均隨著直徑比的增大而呈現出增大趨勢； 當直徑比較小時，下游圓柱處于上游圓柱的繞流回流區段，受到與原水流方向相反的作用力， 使得Cdrms 較小；隨著直徑比的逐漸增大，圓柱間的干擾作用減小， 下游小圓柱的振動狀態逐漸趨近于單圓柱振動狀態，使得Cdrms 逐漸增大。

3.3 頻域特征

工程中一般用斯特羅哈爾數（Strouhal，St）來確定旋渦脫落的頻率。 St 的計算式為

式中：fS為旋渦脫落的頻率，Hz。

計算串列不等直徑雙圓柱繞流時發現， 大小圓柱的渦脫落頻率相同，但均比單圓柱繞流小。圖6 給出了d/D=0.6 時不同間距比對應的升力信號的功率譜密度。

圖6 升力的功率譜密度分布示意圖Fig.6 Schematic diagram of power spectral density distribution of lift

從圖6 中可以看出：當G/D=2，3，4，6 時，上游大圓柱對應的旋渦脫落頻率分別為0.908 7，0.931 4，0.933 2，0.934 9，將其帶入式（7），計算得到St 分別為0.207 7，0.212 9，0.213 3，0.213 7；隨著間距比的增大，旋渦脫落頻率逐漸增大，同時St也有略微增大的趨勢， 但St 總體保持在0.21 左右，這說明St 有一定的魯棒性且尾流中旋渦泄放規則有序。 本文的仿真結果與文獻[15]，[16]中的仿真結果相接近，證明了本文仿真結果的有效性。

3.4 尾流特征

當G/D=2 時，不同直徑比下的串列雙圓柱的瞬時渦量如圖7 所示。

從圖7 中可以看出：當間距比相同時，隨著直徑比的增大，尾流形態不斷發生改變。 當d/D=0.2時，大圓柱分離的剪切層附著于下游小圓柱表面，大、小圓柱后的渦脫落形態較為穩定，呈現雙旋渦脫落形態，尾流雙旋渦呈現近似平行狀態，旋渦脫落不明顯，呈現相互連接的狀態；當d/D=0.4 時，大圓柱分離的剪切層附著于下游小圓柱上表面，剪切層接觸小圓柱表面部分減少，大、小圓柱后均出現穩定的渦脫落現象，旋渦脫落明顯；當d/D=0.6 時，大圓柱分離的剪切層不再附著于下游小圓柱下表面， 大、 小圓柱后均出現穩定的渦脫落現象，旋渦脫落明顯，穩定單旋渦數量增加；當d/D=0.8 時，大圓柱分離的剪切層不再附著于下游小圓柱表面，大、小圓柱后均出現穩定的渦脫落現象，旋渦脫落明顯，尾流渦之間的連接明顯減少，圓柱之間的干涉效應逐漸消失，基本呈現單渦形態。

4 結論

①隨著直徑比的逐漸增大， 大圓柱的振幅比變化較大，小圓柱的振幅比逐漸減小；存在直徑比范圍（d/D＜0.7）使得大、小圓柱的振幅比取得最大值；當間距比較小時，直徑比的變化對串列不等直徑圓柱系統的振幅響應影響較大。

②在間距比相同的情況下， 直徑比的增大會使下游圓柱的渦激振動響應增強； 隨著直徑比的逐漸增大， 升力系數均方根和阻力系數均方根總體呈現出逐漸增加的變化趨勢， 下游圓柱對上游圓柱的影響也隨著直徑比d/D 的增大而增強。

③在不同直徑比和間距比下，大、小圓柱的渦脫落頻率相同，但都較單圓柱繞流小；旋渦脫落頻率和St 均隨著間距比的增大而增大。