林杉杉， 李 莉， 安然然

(沈陽化工大學(xué) 計算機科學(xué)與技術(shù)學(xué)院， 遼寧 沈陽 110142)

海洋面積占據(jù)地球表面積的75 %，其中蘊藏著豐富的資源及能量.在海洋工程研究中，當(dāng)水流經(jīng)一些非流線型結(jié)構(gòu)時，會在其尾部發(fā)生邊界層分離，形成周期性交替脫落的漩渦，當(dāng)脫渦頻率與結(jié)構(gòu)體的固有頻率相等時，則會產(chǎn)生渦致振動(vortex-induced vibration，VIV).渦致振動過程中交替脫落的漩渦會在垂直于流體來流方向產(chǎn)生周期變化的作用力，從而激勵彈性支撐的壓電俘能器產(chǎn)生振動并發(fā)電.

目前，對于單、雙圓柱的繞流及渦致振動的數(shù)值及實驗問題已取得了大量的研究成果.鄒琳等[1]研究了串列雙圓柱渦致振動的數(shù)值模擬情況，研究發(fā)現(xiàn)隨著兩圓柱間距比的增大，下游圓柱進入鎖定區(qū)間的折合速度逐漸減小，下游圓柱對上游圓柱的影響逐漸減小.然而三圓柱做為圓柱群中的基本組成單元，對其渦致振動的研究則較少.當(dāng)流體流過圓柱時，由于各個圓柱的相互作用及相互干擾使得圓柱的受力情況更加復(fù)雜，徐有恒等[2]進行了三圓柱在風(fēng)洞中的實驗，對壓力系數(shù)的分布情況進行了研究，徐楓等[3]對正三角形排列的圓柱繞流情況進行了數(shù)值模擬.雖然多年來國內(nèi)外對于渦致振動的研究非常多，并且取得了大量研究成果，但是將渦致振動作為可利用的能源進行俘能技術(shù)的研究卻還剛剛起步，在數(shù)值計算和實驗研究方面的成果非常少.Molino-Minero-Re等[4]對圓柱外加單懸臂梁的壓電裝置進行了水槽內(nèi)渦激振動的能量收集實驗，實驗證明在不同尺寸的圓柱時實驗裝置的產(chǎn)能情況，其最大可以產(chǎn)生0.3 μW的電能.Abdelkefi[5]對剛性圓柱橫自由度方向安裝壓電能量換能器的渦致振動能量獲取問題進行了研究，確定了負載電阻對鎖頻區(qū)域的影響.

本文利用ADINA對串列三圓柱渦致振動情況進行數(shù)值模擬，并對結(jié)果進行分析.將由于渦致振動而使柔性圓管產(chǎn)生的位移作為機械能的輸入對機電耦合情況進行數(shù)值模擬，探索渦致振動的影響因素以及產(chǎn)生電壓的情況.

1 數(shù)值計算模型及方法

1.1 壓電發(fā)電結(jié)構(gòu)及計算模型

針對李莉等課題組提出的渦致振動壓電發(fā)電結(jié)構(gòu)(見圖1)進行研究[6]，該結(jié)構(gòu)將壓電懸臂梁沿柔性圓管軸向放置于圓管內(nèi)部，利用流體流經(jīng)圓管時產(chǎn)生的漩渦脫落帶動圓管垂直來流方向做周期性的振動，從而使壓電結(jié)構(gòu)由于振動而產(chǎn)生電能.

圖1 流致振動壓電發(fā)電結(jié)構(gòu)

將2個參數(shù)相同的壓電能量收集結(jié)構(gòu)串列，置于剛性圓柱后方，計算模型如圖2所示.原點位于上游圓柱的中心處，x軸為流體來流方向，z軸沿圓管垂直向上，y軸垂直于xz面.其中L為2個圓管圓心之間的距離，D為圓管的外直徑.利用計算機數(shù)值模擬計算剛性圓柱后的2個柔性圓管振動情況.

圖2 三圓管計算模型

1.2 數(shù)值分析方法

渦致振動為非線性流固耦合運動，由于柔性圓管內(nèi)置了壓電懸臂梁，且懸臂梁平面與來流方向一致，因此懸臂梁在來流方向起到了桁梁作用，這使得柔性圓管沿來流方向的運動很小，從而可以忽略不計.對于流體的運動，假設(shè)流體為不可壓縮的無黏性流體，則滿足納維-斯托克斯(N-S)方程

(1)

式中：dv為x、y、z軸上的速度；f為質(zhì)量力；ρ為流體密度;pl為流體的壓力.

模態(tài)分析法[7]是研究圓柱結(jié)構(gòu)渦致振動特性的一種常用方法.假設(shè)柔性圓管的振動位移是由一系列固有模態(tài)的線性疊加而成，則根據(jù)圓管上不同離散點的應(yīng)變即可求得圓管的位移情況.高度為H的圓管的振動位移y(z，t)可按模態(tài)分解為

(2)

式中：ωn(t)為各個模態(tài)權(quán)重函數(shù);t為時間;φn為振型函數(shù);z為圓管高度.

對于一端固定一端自由的圓管，其振型函數(shù)可以表示成

(3)

式中n為n階模態(tài).由式(2)和式(3)可得圓管在振動時各點的振動位移

z∈(0,H).

(4)

則位移的二階導(dǎo)數(shù)為

z∈(0,H).

(5)

根據(jù)材料力學(xué)曲率與應(yīng)變之間的關(guān)系可得

z∈(0,H).

(6)

式中：ε為測量獲得應(yīng)變;R為圓管的半徑;s為模態(tài)疊加個數(shù).通過計算式(5)和式(6)即可求得某個時刻的權(quán)重函數(shù)，將權(quán)重系數(shù)結(jié)果代入式(4)即可求得任一點的位移.將圓柱上所有點的運動位移進行合并計算，即可得到圓柱整體的振動情況，這個計算過程需要通過數(shù)值模擬軟件實現(xiàn).此處利用ADINA對串列三圓柱進行數(shù)值模擬，分析柔性圓管在渦致振動情況下的位移形變情況.

2 流固耦合數(shù)值模擬與結(jié)果分析

利用ADINA的流固耦合模塊(FSI)對串列三圓柱進行數(shù)值模擬.采用長方形流體域中第1個圓柱為剛性圓柱，剛性固定在流體域中，起到流致振動中阻流體的作用，第2和第3個圓柱為柔性圓管的串列方式，柔性圓管底端固定，上端在x、y、z方向自由.流固耦合數(shù)值模擬的具體參數(shù)如表1所示.

表1 流固耦合數(shù)值模擬相關(guān)參數(shù)

根據(jù)實際的海底洋流速度，選取入口速度為1.1 m/s.在三圓柱繞流過程中，上游圓柱對下游圓柱的干涉一直存在，但下游圓柱對上游圓柱的干涉隨著間距比的增大逐漸減小，為了減小圓柱間的相互干涉對實驗結(jié)果的影響，實驗研究了圓柱之間間距變化對振動狀態(tài)的影響，結(jié)果發(fā)現(xiàn)當(dāng)圓柱之間的間距在(3～4)D之間時振動效果較好，因此選取間距比L/D=4進行數(shù)值模擬[1].為了提高計算精度，圓柱周圍的網(wǎng)格劃分較密集.流體域和固體域的網(wǎng)格劃分示意圖如圖3和圖4所示.

圖3 流體域網(wǎng)格劃分

圖4 圓柱網(wǎng)格劃分

經(jīng)過計算和后處理流體域求解文件，得到流體域xy截面的壓力云圖如圖5所示.

圖5 流體的壓力云圖

水流經(jīng)串列圓柱，在流動過程中沿著圓柱面產(chǎn)生法向力，即流體流經(jīng)剛性圓柱阻流體后產(chǎn)生穩(wěn)定的交替脫落的漩渦，漩渦與柔性圓管相互耦合，使柔性圓管發(fā)生形變.

經(jīng)過計算和后處理固體域的求解文件，得到柔性圓管的形變情況如圖6所示.剛性圓柱在流體流動過程中起到阻流體作用，固定不動，因此在數(shù)據(jù)分析時不進行分析.選取柔性圓管在不同時間段的振動情況進行分析，若無特殊說明下文中論述的第1、2個圓管都是阻流體下游順流方向的第1、2個柔性圓管.

圖6 固體域模型的應(yīng)力云圖

串列三圓柱的流固耦合仿真結(jié)果顯示：流速由零逐漸增加的初始階段，兩個圓管未發(fā)生明顯的形變，而后第2個柔性圓管較先發(fā)生較大的擾動，第一個圓管擾動較??；隨著時間的推移，兩個圓管在垂直于來流(即y軸)方向的振動位移逐漸變大，達到最大值后又逐漸減小.圖6只能反映圓管的大致振動狀態(tài)，并不能得到具體的形變情況，因此選取圓柱頂端的一點進行具體的形變大小即位移分析，其t-y軸位移曲線如圖7、圖8所示.

圖7 第1個圓管時間-位移曲線

由于圓管之間的相互作用，第2個柔性圓管在實驗進行2 s左右開始起振，而第1個圓管的起振時間則在2.5 s左右，滯后于第2個圓管.然后兩個圓管都在5 s時達到最大振幅，兩個圓管的振幅分別為0.969 mm、1.34 mm，即兩個圓管均已受到漩渦的影響并產(chǎn)生了振動，兩個圓管的振動均較穩(wěn)定，其位移時程曲線為較規(guī)則的正弦曲線，為理想的壓電懸臂梁載荷輸入狀態(tài).

圖8 第2個圓管時間-位移曲線

為了進一步分析流速對串列三圓柱結(jié)構(gòu)渦致振動情況的影響，采用控制變量法只改變流體速度(0.4～1.4 m/s)，其他參數(shù)不變進行多組仿真實驗.考慮到海底洋流的實際流速，數(shù)值模擬中最高流速設(shè)置為1.4 m/s.由于剛性圓柱阻流體固定不動，所以在研究形變位移時不予以考慮.經(jīng)過多組仿真實驗，得到每組仿真結(jié)果中兩個柔性圓管上部任一點的最大振動位移(即柔性圓管的最大形變位移)與流速之間的關(guān)系，如圖9所示.由圖9可以看出:兩柔性圓管的最大振幅先隨著流速的增大不斷增大，流速為1.1 m/s時圓管的響應(yīng)振幅達到最大值，而后隨著流速的增大圓管的振幅開始逐漸減小.這說明在流速為1.1 m/s時脫渦頻率與圓管自身固有頻率較為接近，產(chǎn)生共振現(xiàn)象，使柔性圓管的響應(yīng)振幅達到最大值，此時柔性圓管處于鎖頻的穩(wěn)定振動狀態(tài).本文中圓管流速-響應(yīng)振幅的變化趨勢為先逐漸增加，達到最大值后逐漸減小，與Khalak等[8]研究結(jié)果趨勢相同，進一步證實了流體在流經(jīng)鈍體后產(chǎn)生規(guī)律的漩渦，與柔性圓管發(fā)生共振現(xiàn)象，此時響應(yīng)振幅達到最大值.在同一時刻，第2個柔性圓管的響應(yīng)振幅始終大于第1個柔性圓管的響應(yīng)振幅，這表明第2個柔性圓管內(nèi)的壓電懸臂梁能產(chǎn)生更大的電能.因此在進行壓電仿真時采用第2個圓管的響應(yīng)振幅作為輸入負載進行仿真.

圖9 流速-圓管最大振幅關(guān)系

3 壓電耦合仿真分析

采用壓電雙晶片懸臂梁結(jié)構(gòu)做為機電轉(zhuǎn)換結(jié)構(gòu)，銅片的兩側(cè)為壓電層并聯(lián)結(jié)構(gòu)，壓電懸臂梁的一端固定在柔性圓管的底部.當(dāng)柔性圓管在渦致振動的作用下發(fā)生形變時，壓電結(jié)構(gòu)也隨之振動，由于壓電耦合作用而產(chǎn)生電荷.雙晶片懸臂梁的相關(guān)參數(shù)如表2所示.

由于柔性圓管在渦致振動時垂直于水流方向的橫向振動位移要遠大于順流方向位移，而且能量收集結(jié)構(gòu)中置于柔性圓筒內(nèi)的壓電梁平面與水流方向平行且相對位置固定，進一步減小了圓筒順流方向位移.因此，壓電耦合仿真時除將壓電懸臂梁軸向底端邊界條件設(shè)置為固定外，還將順流方向壓電梁的位移設(shè)置成對稱的邊界條件.

表2 雙晶片懸臂梁的相關(guān)參數(shù)

為提高數(shù)值模擬效率，將壓電結(jié)構(gòu)自由端由于渦致振動而產(chǎn)生振幅響應(yīng)作為機電耦合的初始輸入值進行機電耦合數(shù)值模擬，得到的壓電懸臂梁電壓分布云圖如圖10所示.

圖10 機電耦合電壓分布云

由圖10可以看出：壓電懸臂梁上在距離固定端很近的位置產(chǎn)生的電壓最大，隨著與固定端距離的變大，電壓逐漸減小.另外，為了清楚地得到壓電懸臂梁中的電壓分布情況，將壓電懸臂梁的中線平均分為50個點，然后將壓電懸臂梁上各點產(chǎn)生的電壓都耦合在壓電懸臂梁中線的各個點上，得到懸臂梁中線上50個點的電壓分布圖如圖11所示.由圖11可知:當(dāng)流速為1.1 m/s時第2個柔性圓管內(nèi)壓電懸臂梁產(chǎn)生的電壓從自由端到固定端呈逐漸減小趨勢，在距離固定端很近的位置產(chǎn)生了最大電壓值(68.54 V).這是由于柔性圓管的彎曲振動帶動壓電懸臂梁自由端彎曲振動，而壓電懸臂梁的另一端固定，使得距離固定端很近的位置產(chǎn)生的形變最大，因此產(chǎn)生的電壓也最大；相反隨著與固定端距離的增大，懸臂梁產(chǎn)生的彎曲形變逐漸減小，從而產(chǎn)生的電壓也逐漸變小.

圖11 懸臂梁中線電壓分布

將流固耦合分組仿真得到的第2個柔性圓管最大形變值作為壓電懸臂梁機電耦合仿真負載進行分組仿真，得到了懸臂梁最大開路輸出電壓與流速關(guān)系曲線，如圖12所示.

圖12 流速-最大開路電壓關(guān)系

由圖12可以看出：流速在0.4～0.6 m/s時，電壓峰值的上升比較平緩；在0.7～1.1 m/s時，電壓的峰值上升速度較快；到1.1 m/s時達到峰值68.54 V，后又減小，與圖9第2個柔性圓管的位移變化規(guī)律一致.

4 結(jié) 論

利用ADINA對前置剛性圓柱阻流體的串列三圓柱壓電能量收集結(jié)構(gòu)進行了流固耦合和機電耦合的仿真研究，分析了流速在0.4～1.4 m/s范圍時，柔性圓管振動形變幅值和內(nèi)置懸臂梁產(chǎn)生的電壓變化情況.仿真結(jié)果表明串列三圓柱壓電能量收集結(jié)構(gòu)中第2個柔性圓管產(chǎn)生的形變位移始終大于第1個圓管的形變位移，當(dāng)流速為1.1 m/s時柔性圓管的形變幅值達到峰值，此時圓管處于鎖頻的穩(wěn)定振動狀態(tài)，其振動規(guī)律呈明顯的周期性.根據(jù)流固耦合的仿真結(jié)果，將第2個柔性圓管頂端的最大形變值作為輸入負載進行機電耦合仿真，得到了懸臂梁的電壓分布情況，當(dāng)流速為1.1 m/s時，懸臂梁的開路輸出電壓達到最大值為68.54 V，其與流速的關(guān)系曲線與流固耦合仿真得到的最大形變值與流速的關(guān)系一致.而流固耦合仿真得到的柔性圓管最大形變值與Khalak的渦致振動實驗結(jié)果一致，證明了仿真數(shù)據(jù)的可靠性.從而說明前置阻流體的串列三圓柱壓電能量收集結(jié)構(gòu)能在渦致振動條件下產(chǎn)生持續(xù)的周期性變化的電壓，該結(jié)論為實際研制渦致振動能量收集結(jié)構(gòu)提供了理論基礎(chǔ).