圆柱形涡致振动压电能量收集阵列的仿真分析

林杉杉， 李 莉， 安然然

(沈阳化工大学 计算机科学与技术学院， 辽宁 沈阳 110142)

海洋面积占据地球表面积的75 %，其中蕴藏着丰富的资源及能量.在海洋工程研究中，当水流经一些非流线型结构时，会在其尾部发生边界层分离，形成周期性交替脱落的漩涡，当脱涡频率与结构体的固有频率相等时，则会产生涡致振动(vortex-induced vibration，VIV).涡致振动过程中交替脱落的漩涡会在垂直于流体来流方向产生周期变化的作用力，从而激励弹性支撑的压电俘能器产生振动并发电.

目前，对于单、双圆柱的绕流及涡致振动的数值及实验问题已取得了大量的研究成果.邹琳等[1]研究了串列双圆柱涡致振动的数值模拟情况，研究发现随着两圆柱间距比的增大，下游圆柱进入锁定区间的折合速度逐渐减小，下游圆柱对上游圆柱的影响逐渐减小.然而三圆柱做为圆柱群中的基本组成单元，对其涡致振动的研究则较少.当流体流过圆柱时，由于各个圆柱的相互作用及相互干扰使得圆柱的受力情况更加复杂，徐有恒等[2]进行了三圆柱在风洞中的实验，对压力系数的分布情况进行了研究，徐枫等[3]对正三角形排列的圆柱绕流情况进行了数值模拟.虽然多年来国内外对于涡致振动的研究非常多，并且取得了大量研究成果，但是将涡致振动作为可利用的能源进行俘能技术的研究却还刚刚起步，在数值计算和实验研究方面的成果非常少.Molino-Minero-Re等[4]对圆柱外加单悬臂梁的压电装置进行了水槽内涡激振动的能量收集实验，实验证明在不同尺寸的圆柱时实验装置的产能情况，其最大可以产生0.3 μW的电能.Abdelkefi[5]对刚性圆柱横自由度方向安装压电能量换能器的涡致振动能量获取问题进行了研究，确定了负载电阻对锁频区域的影响.

本文利用ADINA对串列三圆柱涡致振动情况进行数值模拟，并对结果进行分析.将由于涡致振动而使柔性圆管产生的位移作为机械能的输入对机电耦合情况进行数值模拟，探索涡致振动的影响因素以及产生电压的情况.

1 数值计算模型及方法

1.1 压电发电结构及计算模型

针对李莉等课题组提出的涡致振动压电发电结构(见图1)进行研究[6]，该结构将压电悬臂梁沿柔性圆管轴向放置于圆管内部，利用流体流经圆管时产生的漩涡脱落带动圆管垂直来流方向做周期性的振动，从而使压电结构由于振动而产生电能.

图1 流致振动压电发电结构

将2个参数相同的压电能量收集结构串列，置于刚性圆柱后方，计算模型如图2所示.原点位于上游圆柱的中心处，x轴为流体来流方向，z轴沿圆管垂直向上，y轴垂直于xz面.其中L为2个圆管圆心之间的距离，D为圆管的外直径.利用计算机数值模拟计算刚性圆柱后的2个柔性圆管振动情况.

图2 三圆管计算模型

1.2 数值分析方法

涡致振动为非线性流固耦合运动，由于柔性圆管内置了压电悬臂梁，且悬臂梁平面与来流方向一致，因此悬臂梁在来流方向起到了桁梁作用，这使得柔性圆管沿来流方向的运动很小，从而可以忽略不计.对于流体的运动，假设流体为不可压缩的无黏性流体，则满足纳维-斯托克斯(N-S)方程

(1)

式中：dv为x、y、z轴上的速度；f为质量力；ρ为流体密度;pl为流体的压力.

模态分析法[7]是研究圆柱结构涡致振动特性的一种常用方法.假设柔性圆管的振动位移是由一系列固有模态的线性叠加而成，则根据圆管上不同离散点的应变即可求得圆管的位移情况.高度为H的圆管的振动位移y(z，t)可按模态分解为

(2)

式中：ωn(t)为各个模态权重函数;t为时间;φn为振型函数;z为圆管高度.

对于一端固定一端自由的圆管，其振型函数可以表示成

(3)

式中n为n阶模态.由式(2)和式(3)可得圆管在振动时各点的振动位移

z∈(0,H).

(4)

则位移的二阶导数为

z∈(0,H).

(5)

根据材料力学曲率与应变之间的关系可得

z∈(0,H).

(6)

式中：ε为测量获得应变;R为圆管的半径;s为模态叠加个数.通过计算式(5)和式(6)即可求得某个时刻的权重函数，将权重系数结果代入式(4)即可求得任一点的位移.将圆柱上所有点的运动位移进行合并计算，即可得到圆柱整体的振动情况，这个计算过程需要通过数值模拟软件实现.此处利用ADINA对串列三圆柱进行数值模拟，分析柔性圆管在涡致振动情况下的位移形变情况.

2 流固耦合数值模拟与结果分析

利用ADINA的流固耦合模块(FSI)对串列三圆柱进行数值模拟.采用长方形流体域中第1个圆柱为刚性圆柱，刚性固定在流体域中，起到流致振动中阻流体的作用，第2和第3个圆柱为柔性圆管的串列方式，柔性圆管底端固定，上端在x、y、z方向自由.流固耦合数值模拟的具体参数如表1所示.

表1 流固耦合数值模拟相关参数

根据实际的海底洋流速度，选取入口速度为1.1 m/s.在三圆柱绕流过程中，上游圆柱对下游圆柱的干涉一直存在，但下游圆柱对上游圆柱的干涉随着间距比的增大逐渐减小，为了减小圆柱间的相互干涉对实验结果的影响，实验研究了圆柱之间间距变化对振动状态的影响，结果发现当圆柱之间的间距在(3～4)D之间时振动效果较好，因此选取间距比L/D=4进行数值模拟[1].为了提高计算精度，圆柱周围的网格划分较密集.流体域和固体域的网格划分示意图如图3和图4所示.

图3 流体域网格划分

图4 圆柱网格划分

经过计算和后处理流体域求解文件，得到流体域xy截面的压力云图如图5所示.

图5 流体的压力云图

水流经串列圆柱，在流动过程中沿着圆柱面产生法向力，即流体流经刚性圆柱阻流体后产生稳定的交替脱落的漩涡，漩涡与柔性圆管相互耦合，使柔性圆管发生形变.

经过计算和后处理固体域的求解文件，得到柔性圆管的形变情况如图6所示.刚性圆柱在流体流动过程中起到阻流体作用，固定不动，因此在数据分析时不进行分析.选取柔性圆管在不同时间段的振动情况进行分析，若无特殊说明下文中论述的第1、2个圆管都是阻流体下游顺流方向的第1、2个柔性圆管.

图6 固体域模型的应力云图

串列三圆柱的流固耦合仿真结果显示：流速由零逐渐增加的初始阶段，两个圆管未发生明显的形变，而后第2个柔性圆管较先发生较大的扰动，第一个圆管扰动较小；随着时间的推移，两个圆管在垂直于来流(即y轴)方向的振动位移逐渐变大，达到最大值后又逐渐减小.图6只能反映圆管的大致振动状态，并不能得到具体的形变情况，因此选取圆柱顶端的一点进行具体的形变大小即位移分析，其t-y轴位移曲线如图7、图8所示.

图7 第1个圆管时间-位移曲线

由于圆管之间的相互作用，第2个柔性圆管在实验进行2 s左右开始起振，而第1个圆管的起振时间则在2.5 s左右，滞后于第2个圆管.然后两个圆管都在5 s时达到最大振幅，两个圆管的振幅分别为0.969 mm、1.34 mm，即两个圆管均已受到漩涡的影响并产生了振动，两个圆管的振动均较稳定，其位移时程曲线为较规则的正弦曲线，为理想的压电悬臂梁载荷输入状态.

图8 第2个圆管时间-位移曲线

为了进一步分析流速对串列三圆柱结构涡致振动情况的影响，采用控制变量法只改变流体速度(0.4～1.4 m/s)，其他参数不变进行多组仿真实验.考虑到海底洋流的实际流速，数值模拟中最高流速设置为1.4 m/s.由于刚性圆柱阻流体固定不动，所以在研究形变位移时不予以考虑.经过多组仿真实验，得到每组仿真结果中两个柔性圆管上部任一点的最大振动位移(即柔性圆管的最大形变位移)与流速之间的关系，如图9所示.由图9可以看出:两柔性圆管的最大振幅先随着流速的增大不断增大，流速为1.1 m/s时圆管的响应振幅达到最大值，而后随着流速的增大圆管的振幅开始逐渐减小.这说明在流速为1.1 m/s时脱涡频率与圆管自身固有频率较为接近，产生共振现象，使柔性圆管的响应振幅达到最大值，此时柔性圆管处于锁频的稳定振动状态.本文中圆管流速-响应振幅的变化趋势为先逐渐增加，达到最大值后逐渐减小，与Khalak等[8]研究结果趋势相同，进一步证实了流体在流经钝体后产生规律的漩涡，与柔性圆管发生共振现象，此时响应振幅达到最大值.在同一时刻，第2个柔性圆管的响应振幅始终大于第1个柔性圆管的响应振幅，这表明第2个柔性圆管内的压电悬臂梁能产生更大的电能.因此在进行压电仿真时采用第2个圆管的响应振幅作为输入负载进行仿真.

图9 流速-圆管最大振幅关系

3 压电耦合仿真分析

采用压电双晶片悬臂梁结构做为机电转换结构，铜片的两侧为压电层并联结构，压电悬臂梁的一端固定在柔性圆管的底部.当柔性圆管在涡致振动的作用下发生形变时，压电结构也随之振动，由于压电耦合作用而产生电荷.双晶片悬臂梁的相关参数如表2所示.

由于柔性圆管在涡致振动时垂直于水流方向的横向振动位移要远大于顺流方向位移，而且能量收集结构中置于柔性圆筒内的压电梁平面与水流方向平行且相对位置固定，进一步减小了圆筒顺流方向位移.因此，压电耦合仿真时除将压电悬臂梁轴向底端边界条件设置为固定外，还将顺流方向压电梁的位移设置成对称的边界条件.

表2 双晶片悬臂梁的相关参数

为提高数值模拟效率，将压电结构自由端由于涡致振动而产生振幅响应作为机电耦合的初始输入值进行机电耦合数值模拟，得到的压电悬臂梁电压分布云图如图10所示.

图10 机电耦合电压分布云

由图10可以看出：压电悬臂梁上在距离固定端很近的位置产生的电压最大，随着与固定端距离的变大，电压逐渐减小.另外，为了清楚地得到压电悬臂梁中的电压分布情况，将压电悬臂梁的中线平均分为50个点，然后将压电悬臂梁上各点产生的电压都耦合在压电悬臂梁中线的各个点上，得到悬臂梁中线上50个点的电压分布图如图11所示.由图11可知:当流速为1.1 m/s时第2个柔性圆管内压电悬臂梁产生的电压从自由端到固定端呈逐渐减小趋势，在距离固定端很近的位置产生了最大电压值(68.54 V).这是由于柔性圆管的弯曲振动带动压电悬臂梁自由端弯曲振动，而压电悬臂梁的另一端固定，使得距离固定端很近的位置产生的形变最大，因此产生的电压也最大；相反随着与固定端距离的增大，悬臂梁产生的弯曲形变逐渐减小，从而产生的电压也逐渐变小.

图11 悬臂梁中线电压分布

将流固耦合分组仿真得到的第2个柔性圆管最大形变值作为压电悬臂梁机电耦合仿真负载进行分组仿真，得到了悬臂梁最大开路输出电压与流速关系曲线，如图12所示.

图12 流速-最大开路电压关系

由图12可以看出：流速在0.4～0.6 m/s时，电压峰值的上升比较平缓；在0.7～1.1 m/s时，电压的峰值上升速度较快；到1.1 m/s时达到峰值68.54 V，后又减小，与图9第2个柔性圆管的位移变化规律一致.

4 结 论

利用ADINA对前置刚性圆柱阻流体的串列三圆柱压电能量收集结构进行了流固耦合和机电耦合的仿真研究，分析了流速在0.4～1.4 m/s范围时，柔性圆管振动形变幅值和内置悬臂梁产生的电压变化情况.仿真结果表明串列三圆柱压电能量收集结构中第2个柔性圆管产生的形变位移始终大于第1个圆管的形变位移，当流速为1.1 m/s时柔性圆管的形变幅值达到峰值，此时圆管处于锁频的稳定振动状态，其振动规律呈明显的周期性.根据流固耦合的仿真结果，将第2个柔性圆管顶端的最大形变值作为输入负载进行机电耦合仿真，得到了悬臂梁的电压分布情况，当流速为1.1 m/s时，悬臂梁的开路输出电压达到最大值为68.54 V，其与流速的关系曲线与流固耦合仿真得到的最大形变值与流速的关系一致.而流固耦合仿真得到的柔性圆管最大形变值与Khalak的涡致振动实验结果一致，证明了仿真数据的可靠性.从而说明前置阻流体的串列三圆柱压电能量收集结构能在涡致振动条件下产生持续的周期性变化的电压，该结论为实际研制涡致振动能量收集结构提供了理论基础.