林杉杉， 李 莉， 安然然

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

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

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

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

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

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

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

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

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

圖2 三圓管計(jì)算模型

1.2 數(shù)值分析方法

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

(1)

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

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

(2)

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

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

(3)

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

z∈(0,H).

(4)

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

z∈(0,H).

(5)

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

z∈(0,H).

(6)

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

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

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

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

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

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

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

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

圖5 流體的壓力云圖

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

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

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

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

圖7 第1個(gè)圓管時(shí)間-位移曲線

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

圖8 第2個(gè)圓管時(shí)間-位移曲線

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

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

3 壓電耦合仿真分析

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

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

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

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

圖10 機(jī)電耦合電壓分布云

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

圖11 懸臂梁中線電壓分布

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

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

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

4 結(jié) 論

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