串列不等直徑雙圓柱海流能發(fā)電振子渦激振動(dòng)數(shù)值模擬

白 旭， 陳 云

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

0 引言

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

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

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

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

1 計(jì)算模型

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

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

2 數(shù)值計(jì)算方法

2.1 控制方程

不可壓縮粘性牛頓流體的控制方程為Navier-Stokes 方程。 笛卡爾坐標(biāo)系中連續(xù)性方程和動(dòng)量方程分別為

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

柱體振動(dòng)方程可以表示為

2.2 重疊網(wǎng)格技術(shù)

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

在本文中，次要網(wǎng)格（O 型網(wǎng)格） 覆蓋了圓柱周圍的區(qū)域。 網(wǎng)格之間的數(shù)據(jù)轉(zhuǎn)換是通過(guò)被稱為邊緣點(diǎn)的特定點(diǎn)來(lái)完成的。主要網(wǎng)格（背景網(wǎng)格） 的邊緣點(diǎn)是放置在圓柱體外部和其表面附近的節(jié)點(diǎn)， 而次要網(wǎng)格的邊緣點(diǎn)是放置在每個(gè)次要網(wǎng)格外部邊界線上的節(jié)點(diǎn)。 另一組點(diǎn)稱為孔點(diǎn)，在網(wǎng)格中切割并與固體表面相交。為了將信息從網(wǎng)格中的邊緣點(diǎn)傳輸?shù)狡渌W(wǎng)格中的節(jié)點(diǎn)， 本文采用搜索算法指定對(duì)象網(wǎng)格中該邊緣點(diǎn)的周圍節(jié)點(diǎn)[9]，[10]。

圖2 離散化的計(jì)算網(wǎng)格Fig.2 Discretized computational grid

2.3 可靠性驗(yàn)證

為了驗(yàn)證本文模型的可靠性， 將利用本文模型計(jì)算得到的單圓柱及串列雙圓柱靜止繞流的結(jié)果與文獻(xiàn)結(jié)果進(jìn)行對(duì)比。

2.3.1 單圓柱驗(yàn)證

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

表1 單圓柱繞流結(jié)果對(duì)比Table 1 Comparison of single cylinder flow results

2.3.2 等直徑串列雙圓柱驗(yàn)證

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

表2 雙圓柱繞流結(jié)果對(duì)比Table 2 Comparison of the results of flow around a double cylinders

3 結(jié)果與分析

3.1 振幅響應(yīng)

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

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

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

3.2 受力特性

3.2.1 升力特性

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

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

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

3.2.2 阻力特性

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

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

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

3.3 頻域特征

工程中一般用斯特羅哈爾數(shù)（Strouhal，St）來(lái)確定旋渦脫落的頻率。 St 的計(jì)算式為

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

計(jì)算串列不等直徑雙圓柱繞流時(shí)發(fā)現(xiàn)， 大小圓柱的渦脫落頻率相同，但均比單圓柱繞流小。圖6 給出了d/D=0.6 時(shí)不同間距比對(duì)應(yīng)的升力信號(hào)的功率譜密度。

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

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

3.4 尾流特征

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

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

4 結(jié)論

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

②在間距比相同的情況下， 直徑比的增大會(huì)使下游圓柱的渦激振動(dòng)響應(yīng)增強(qiáng)； 隨著直徑比的逐漸增大， 升力系數(shù)均方根和阻力系數(shù)均方根總體呈現(xiàn)出逐漸增加的變化趨勢(shì)， 下游圓柱對(duì)上游圓柱的影響也隨著直徑比d/D 的增大而增強(qiáng)。

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