流致振动广泛地存在于许多工程领域中^[1]，作为储量巨大的清洁能源，海流能与风或海浪相比，具有季节性、周期性，其流动状况高度可预测。但近海岸洋流流速往往较低，导致难以利用传统方式获取海洋能量。流致振动捕能技术正在成为可再生能源发电的可行解决方案，需要开展重大的研究、开发和部署计划，才能在该领域取得真正的产业成功^[2]。近年来，许多学者开始研究利用流致振动获取低流速下海洋动能的发电方式。

流致振动是指流体流经固体时会对固体表面施加交替相间的流体力，使得固体发生往复运动，而固体的往复运动又改变流体流态，进而改变作用于固体表面的流体力^[3]。尽管有关流致振动发电的研究并未形成系统的研究体系^[4]，但仍取得了成果。Belloli等^[5]比较了临界流动状态下涡激振动与亚临界振动的特征。Sun Hai等^[5-7]研究具有端部弹性约束的刚性振子的流致振动，对幅度响应，频率响应，能量收集和效率的实验结果进行了介绍和讨论，并使用2个汽缸串联进行实验，给出了能量收集和效率的实验结果、推导出最优功率的包络。罗竹梅^[8]提出了一种结合并列和串列构形的矩形布置四圆柱振动系统，并提出机械耦合的概念描述这种系统的动态特性，将由4个圆柱组成的振动系统整体视为一个质量-弹簧-阻尼系统。Kim等^[9]证明了交替提升技术(ALT)转换器是一种能量吸收装置，可有效地在速度低至1.0~1.5 m/s的水流中工作，这是目前水车和涡轮机无法达到的速度范围。燕翔^[10]分析了电机参数、机械参数及振子截面参数对三棱柱振动模式、能量利用的影响，明确了三棱柱在流致振动能量利用上的优势。练继建等^[4]剖析了研究当中存在的各类问题，并对流致振动能量利用的研究方向提出了建议。这些研究集中在物理结构、流速等参数对振动的影响上，尚未存在对振动捕能效果的组合优化的研究。本文工作是研究在现有非理想条件下，所能达到的能量转化率最大值，即最优捕能效率值。

1 振动能量获取的数学物理模型 1.1 能量的获取过程

对于本文所建立的流致振动模型，其运动中产生的可利用能量经历了2次削弱的过程:1)振动系统获取流体动能时，有一部分流体动能耗散在振子尾流粒子的碰撞、摩擦中，还有一部分流体未参与做功；2)振动系统获取的能量会有一部分损失在振动系统机械部分的摩擦中，剩下的部分才是可以加以利用的机械能，也即输入发电机的能量。能量转化过程如图 1所示。

1.2 振子物理模型及受力分析

计算采用的振子直径D=0.1 m，展向长度L=1 m，振子表面光滑，由2支拉压弹簧连接在支架上，弹簧刚度为K=500 N/m。模拟质量比m^*取1.1、1.4、1.7、2.0这4个数值，阻尼比ξ从0.03~0.21每隔0.03取值一次，总计取7个值。对以上28种情况，约化速度从2~12，间隔0.1取值计算。

在流致振动过程中，振子在竖直方向受3个力：即弹簧拉压力F_spr、结构及水流的阻力F_damp、流体动压力F_Fluid。

如图 2所示，当振子处于平衡位置以下且运动方向向下时，其加速度方向向上。故可以得到等式：

$ {m_{{\rm{osc}}}}a + {F_{{\rm{spr}}}} + {F_{{\rm{damp}}}} = {F_{{\rm{Fluid}}}} $

(1)

整理得：

$ \begin{array}{l} ({m_{{\rm{osc}}}} + {m_a}){Y^{\prime \prime }} + ({c_{{\rm{fluid}}}} + {c_{{\rm{structure}}}}){Y^\prime } + KY = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}\rho U\varGamma {\rm{cos}}(2\pi {f_{{\rm{st}}}}t) \end{array} $

(2)

式中:c_Fluid为流体附加阻尼；c_structure为振动结构的总阻尼；Γ为离散点涡强度；U为水流的来流速度；D、L为圆柱振子的直径、长度；m_a为振子的附加质量；流体力F_Fluid根据文献[13]讨论得到。

1.3 振幅比A^*、频率比f^*计算

引入涡强方程Γcos(2πf_stt)=βUDq、无量纲位移y=Y/D、无量纲时间τ=tω_st，化简式(2)得到振子运动方程；引入范德波尔方程描述无量纲漩涡强度q，得到升力荷载变化方程。二式相互耦合，联立得^[14]：

$ \left\{ {\begin{array}{*{20}{l}} {{y^{\prime \prime }} + \left( {2\xi \delta + \frac{\gamma }{\mu }} \right){y^\prime } + {\delta ^2}y = \frac{{\rho {D^2}}}{{{m_{{\rm{ osc }}}} + {m_a}}}\frac{1}{{8{\pi ^2}S_t^2}}\frac{{{C_{L0}}}}{2}q}\\ {{q^{\prime \prime }} + \varepsilon ({q^2} - 1){q^\prime } + q = h{y^{\prime \prime }}} \end{array}} \right. $

(3)

式中:δ=f_n/f_st是结构水中自振频率与斯特劳哈尔泻涡频率的比值；流体阻尼参数γ与平均拖曳力系数C_D成正比^[15-16]，与文献[15-16]中的模型相比本文流体阻尼较小，取0.5;μ=(m_osc+m_a)/ρD²；斯特劳哈尔数S_t取0.193；脉动参考升力系数C_L0取0.3^[17-18]；对m^*≤13的中低质量比振子，范德波尔参数ε=0.018e^0.21m*；通过与Stansby的实验数据比较，流体振子经验参数h取12。

振幅比A^*是无量纲位移y的极值，振幅比f^*是振子振动频率与振子在水中的自振频率的比值。使用四阶龙格库塔法求解式(3)，时间步长取0.1，得到不同约化速度下振子随周期振动的振幅比A^*，对震荡的时历曲线进行傅里叶展开，求得频率比f^*。

1.4 功率计算 1.4.1 流体功率计算

水流做功的功率可以根据水流对振子的压力乘以水流流速求得。根据伯努利方程，水流产生的动压强为(ρU²)/2。常见的水流功率的计算方式有2种，即水流压强作用面S取振子与流体直接接触面的投影D×L，或取振子振动过程扫过的面积(2A+D)×L。这两者并无本质区别，只是对流体动能利用的关注角度不同。本文采用后一种方案，且因连杆在水流方向的投影面积相比圆柱振子的投影面积可以忽略不计，忽略水流对连杆部分的作用力，取水流压强的作用面为S=(2A+D)×L。水流功率表示为：

$ {P_{{\rm{ Fluid }}}} = \frac{1}{2}\rho {U^2} \times S \times U = \frac{1}{2}\rho {U^3}(2A + D)L $

(4)

1.4.2 振子捕能功率计算

将式(2)流体附加阻尼项移到等式右侧，得到：

$ \begin{array}{*{20}{c}} {({m_{{\rm{osc}}}} + {m_a}){y^{\prime \prime }} + {c_{{\rm{structure}}}}{y^\prime } + Ky = }\\ {\frac{1}{2}\rho U\varGamma {\rm{cos}}(2\pi {f_{st}}t) - {c_{{\rm{Fluid}}}}{y^\prime }} \end{array} $

(5)

式(5)右侧即水流力与流体摩擦的差值，左右相等，故左侧部分积分即可得到振动系统捕获的能量，而后再除以周期时长，即得到流体的平均有效功率。振子捕能包含振动系统的机械能、振动系统机械摩擦损耗2部分。其数值上与与水流对振子做的功相等。对能量及功率的计算为：

$ \begin{array}{*{20}{c}} {{P_{{\rm{ Mesh }}}} = \frac{1}{{{T_{{\rm{ osc }}}}}}{W_{{\rm{ Mesh }}}} = \frac{1}{{{T_{{\rm{ osc }}}}}}\int\limits_0^{{T_{{\rm{ ose }}}}} {(({m_{{\rm{ osc }}}} + {m_a}){y^{\prime \prime }} + } }\\ {{c_{{\rm{ structure }}}}{y^\prime } + Ky){y^\prime }{\rm{d}}t} \end{array} $

(6)

根据简谐理论，振子运动位移y随时间呈简谐变化，即y=Asin(2πf_osct)，其中振子的振动主频为f_osc，其值为T_osc=1/f_osc。在稳定振动中，液体与振子发生共振，根据结构动力学共振的规律，此时的振动主频与升力主频相等，即f_osc=f_Fluid。

根据三角函数的奇偶性、周期性，含有y、y″的项积分后为0，代入式y、y′、y″，化简为：

$ \begin{array}{*{20}{c}} {{P_{{\rm{ Mesh }}}} = \frac{1}{{{T_{{\rm{ osc }}}}}}\int\limits_0^{{T_{{\rm{ osc }}}}} {{c_{{\rm{ structure }}}}} {y^{\prime 2}}{\rm{d}}t = }\\ {\frac{1}{{{T_{{\rm{ osc }}}}}}\int\limits_0^{{T_{{\rm{ osc }}}}} {{c_{{\rm{ structure }}}}} 4{\pi ^2}f_{{\rm{ osc }}}^2{A^2}{\rm{co}}{{\rm{s}}^2}(2\pi {f_{{\rm{ osc }}}}t){\rm{d}}t = 2{\pi ^2}f_{{\rm{ osc }}}^2{A^2}{c_{{\rm{ structure }}}}} \end{array} $

(7)

1.5 捕能效率计算

由于流量守恒，通过发电装置的流体速度不会发生突变，电机转换的能量仅由压强差提供，因此流体能量损失天然地存在一个上限，即贝茨极限。在桨叶式发电装置中贝茨极限取59.3%。此外，因为尾流粒子间、流体固体间的摩擦，振子不能充分利用水流本身的能量损失，因此振动系统对水流能量的利用率还要小于贝茨极限。定义振子捕能功率与流体损失功率的比值为捕能效率，表达式为：

$ \eta = \frac{{{P_{{\rm{ Mesh }}}}}}{{{P_{{\rm{ Fluid }}}}}} = \frac{{2{\pi ^2}f_{{\rm{ osc }}}^2{A^2}{c_{{\rm{ structure }}}}}}{{\frac{1}{2}\rho {U^3}(2A + D)L}} = \frac{{4{\pi ^4}\xi ({m^*} + 1){{({f^*}{A^*})}^2}}}{{U_r^3(2{A^*} + 1)}} $

(8)

式中：质量比m^*、阻尼比ξ为结构参数；振幅比A^*、频率比f^*为流固耦合参数；约化速度U_r为流体参数。流固耦合参数与流体、固体参数均相关，无法直接获得，计算过程见2.3节。

2 捕能效率的影响分析 2.1 模型可靠性验证及计算方案

Khalak等^[19]进行了m^*=2.4、m^*ξ=0.013的低质量比、低阻尼比物理实验；罗竹梅^[8]采用SST k-ω高雷诺数湍流模型在相同的质量比、阻尼比下进行数值模拟。本文在相同条件下，使用尾流振子模型计算约化速度2 < U_r < 12情况下振子的振幅比与频率比，并把计算结果与上述文献结果代入式(8)，将得到的3组捕能效率进行对比，从而验证计算模型的可靠性。其中Khalak and Williamson实验与罗竹梅模拟的各流速下频率比振幅比数值由已知数据点线性插值求得。

振子在不同约化速度下的3种捕能效率曲线如图 4所示。从图中可以看出，本文模拟的曲线走势与文献结果较好吻合，最大捕能效率介于文献[19]与文献[8]的实验、模拟结果，且极为接近。如表 1所示，模拟的最大捕能效率与文献结果的误差分别是7.8%和23.3%，这说明了本文模拟方法的可靠性。模拟振动在U_r=4.5时获得最大捕能效率2.58%，而参考文献分别在U_r=4.5和U_r=4.6时对应取得最大捕能效率。相比下，模拟计算的整体捕能趋势均略微滞后，这可能是范德波尔经验参数取值偏小的缘故。

根据模拟结果，振子共振稳定后的横向振动响应近似正弦规律的谐波振动^[8]，取模拟计算的最后19个周期，以分析质量比、阻尼比对振子振幅比、频率比、捕能效率的影响，并总结其相互关系。

2.2 质量比、阻尼比对振子振幅比、频率比的影响

随约化速度的增加，不同质量比、阻尼比条件下的振子振幅比先增加后降低，质量比、阻尼比越低，振幅越大。记取最大振幅比的约化速度为最佳U_r，则振子的最佳U_r集中在6.3~6.9。最佳U_r值同阻尼比几乎没有关联性，但会随着质量比的增大稍稍提前。此外振幅比的变化率与质量比无关，但随着阻尼比的减小，振幅比增大得越来越快。

频率比同约化速度成正相关，其增大变化率经历了陡峭-平缓-陡峭的过程。平缓阶段4 < U_r < 8，此时振子处于锁定阶段，泻涡频率受振子振动频率的控制，两者发生共振，因而此阶段振子的振幅比最大。阻尼比越小，振子受锁定的影响越大。取最大捕能效率的U_r也位于锁定区域内。

由图 5(a)知，平缓区间内的频率比与质量比无关。此外，不同阻尼比条件下频率比随约化速度的变化可以分成2个区间。区间一为3.3 < U_r < 5.1，此区域内阻尼比越小频率比越大；此后各阻尼比小的频率比值迅速接近，当U_r>5.5，振子阻尼比越大频率比越大。

2.3 不同质量比、阻尼比下捕能效率的分析

不同质量比、阻尼比下振子捕能效率随约化速度的变化情况如图 6所示。

1) 总体而言，当约化速度处于4.5~6.4时，多数振子捕能效率η>10%，能量利用率较高。最大捕能效率见图 6(a)，m^*=1.1、ξ=0.15时，捕能效率η=14.27%，此时U_r=5.1。

2) 在趋势上，捕能效率η随约化速度的增加先增加再降低，在约化速度U_r=5.3附近取得最大值，定义η取最大值对应的U_r为最优约化速度，此刻捕能效率称为最优效率η_b。效率值η在远离取最优U_r处变化平缓，在最优U_r邻域也有较小变化率形成的缓和段，但在邻域外侧随U_r变化的程度尤为激烈。同图 5的锁定区间比较，最优邻域恰好位于稳定共振区域，而位于进入、退出共振段区域的邻域外侧，捕能效率变化较大。

最优U_r的位置随阻尼比的增大微微增大，与质量比的取值无关。

3) 最优效率η_b随质量比的变化见图 7，阻尼比ξ较小时，η_b与m^*正相关；ξ较大时，η_b与m^*反相关。这是由于随着阻尼比的增大，振幅比迅速变小，同时频率比在U_r=5.1附近几乎没有变化。

4) 振子阻尼比较大时，最优效率η_b随m^*的变化越剧烈，从图 7可看出，对于不同的m^*，均存在ξ使η_b取极大值。如式(8)所示，这是因为阻尼比较小时，阻尼比本身对效率的影响较大；随着阻尼比增大，较大的阻尼比导致振幅比迅速减小，振幅比成为影响效率η的主导因子。不同质量比下取得最优效率的阻尼比数值可以由式(10)确定。

2.4 最优阻尼比分析

使捕能效率最大的阻尼比随质量比的变化而改变，为获取最优捕能效率，需要在模拟前估算最优阻尼比数值。用最小二乘法对η_b作关于m^*、ξ的多项式曲面拟合，从而求解最优阻尼比, 如图 8所示。

曲面方程：

$ \begin{array}{*{20}{l}} {{\eta _b} = 0.707{\kern 1pt} {\kern 1pt} 6{\xi ^2}{m^*} + 13.54{\xi ^3} - 0.419{\kern 1pt} {\kern 1pt} 2\xi {m^*} - }\\ {9.734{\xi ^2} + 0.021{\kern 1pt} {\kern 1pt} {\kern 1pt} 38{m^*} + 2.207\xi + 0.015{\kern 1pt} {\kern 1pt} {\kern 1pt} 35}\\ {(2 \le {U_r} \le 12;1.1 \le {m^*} \le 2;0.03 \le \xi \le 0.3)} \end{array} $

(9)

拟合的复相关系数R-square=0.975 3，均方根误差RMSE=0.003 007，误差平方和SSE=0.000 298 4，故式(9)可以准确地描述质量比、阻尼比发生改变时，最优捕能效率η_b的变化趋势。

对式(9)右侧取关于ξ的偏导，并使其为零，得到关于m^*、ξ的等式(10)，即某一m^*下取得最佳效率η_b的阻尼比ξ求解公式：

$ \begin{array}{*{20}{c}} {\xi = \frac{{(19.468 - 1.415{m^*})}}{{81.24}} - }\\ {\frac{{\sqrt {2.003{m^{*2}} + 11.011{m^*} + 20.41} }}{{81.24}}} \end{array} $

(10)

分别代入m^*=1.1到2，得ξ=0.148、0.138、0.128、0.117，符合对图 6的预测。这也间接验证了式(9)的可靠性。

将m^*=1.1、ξ=0.148代入式(9)，得η_b=14.50%，即计算区域内的最大捕能效率。

3 结论

1) 振子的捕能效率峰值早于振幅比的峰值出现，这是由于振幅比、频率比的增长速率较约化速度的增长速率变化较小的缘故。

2) 根据模拟结果，振子在U_r=5.1、m^*=1.1、ξ=0.15时，捕能效率取最大值η_b=14.27%；根据经验公式，振子在m^*=1.1、ξ=0.148时，捕能效率取最大值η_b=14.50%。

3) 捕能效率随约化速度的变化曲线存在最大值，最佳捕能效率对应的约化速度U_r位于区间[5.0，5.3]。

4) 在式(9)所处的低质量比区间上，最优捕能效率与质量比成反相关，质量比越小捕能效率越高；但阻尼比存在使捕能效率达到最大的极值，其值随着质量比的增大，由0.093降低至0.068。