潮流能是一种蕴藏在海洋中的可再生能源，广泛分布于世界各地。中国拥有广阔的海域，潮流能资源丰富^[1]。目前，国内主要是利用水轮机将潮流能转化成机械能来开发潮流能，而水轮机式潮流能发电装置能够适用的最小水深为2~5 m^[2]。然而，通过使用涡激振动发电装置，可以在水深不超过2 m的低流速海域收集能量。因此，涡激振动发电可以进一步利用水深度较浅的潮流能区域，增加可开发的潮流能量。

国内、外针对涡激振动的研究主要是采用实验研究与数值模拟的手段。在实验研究方面，早在1996年，Khalak等^[3]和Williamson等^[4]设计了一套刚性圆柱振子的涡激振动实验装置，研究圆柱振子的振动特性，并发现了“2T”模态。2007年，Stappenbelt等^[5]使用简易涡激振动实验装置，研究了单自由度和双自由度下不同质量比的圆柱振子振幅响应运动规律。2008年，密歇根大学的Bernitsas等^[6]利用涡激振动原理提出并建立了涡激振动发电装置(Vortex-induced vibration for aquatic clean energy converter，VIVACE)的物理模型，并在美国海洋可再生能源实验室(Marine Renewable Energy Laboratory, MRELab)对其进行了一系列试验。2011年，Lee等^[7]设计了第一代虚拟阻尼/弹簧(Virtual damper/spring，Vck)系统，用虚拟元件替代涡激振动发电装置的阻尼和弹簧，在短时间内进行了大量试验，并与具有物理弹簧和阻尼器的VIVACE系统的振幅响应进行比较，验证了Vck涡激振动发电系统的性能。2016年，Sun等^[8]在第一代Vck的基础上设计开发了第二代Vck系统，不同于第一代的是，第二代Vck系统有一套独立的虚拟嵌入式程序，能够使用虚拟信号对涡激振动装置进行研究，提高研究效率。在数值模拟方面，2009年，徐枫等^[9]使用数值模拟的方法，在单、双自由度的基础上，研究了不同截面形状的圆柱振子的流致振动现象。2016年，白旭等^[10]使用尾流振子模型，对不同质量比下的双自由度圆柱振子涡激振动的振幅进行了数值模拟，并得到了不同质量比下圆柱振子顺流向和横流向的振幅响应。2018年，中国海洋大学的蔡云雯等^[11]开发设计了双圆柱振子涡激振动潮流能转换装置，并使用CFD软件探索了圆柱振子间距和雷诺数对下游圆柱振子涡激振动的影响。

当前研究主要基于固定的弹簧刚度，且大部分涡激振动能量转化装置只能在特定流速范围的水流中进行能量利用，无法满足在不同海况下圆柱振子的振幅最大化的要求。为了使涡激振动发电系统能够对能量进行高效利用，本文将通过数值模拟来探究不同刚度圆柱振子在时变流速条件下不同刚度圆柱振子的振幅响应、频率响应和尾涡特征变化。

1 数值模拟 1.1 单自由度涡激振动仿真模型

当流体流经圆柱振子时，圆柱振子会出现横流向振动响应和顺流向振动响应。由于质量比较大的圆柱振子在流场中顺流向的振幅很小，可忽略不计^[12]，因此本文中仅专注于研究横流向振动响应。为简化分析，文中采用弹簧-质量-阻尼系统模型，如图 1所示。该模型中：U代表流速；m代表圆柱振子的质量；D代表圆柱振子的直径；C代表系统结构阻尼系数；K代表弹簧的刚度系数。

1.2 控制方程

在圆柱绕流问题中，假设流体为黏性不可压缩流体，服从雷诺平均纳维-斯托克斯(Reynolds average Navier-Stockes，RANS)方程：

$ \frac{{\partial {{\overline u }_i}}}{{\partial {x_i}}} = 0，$

(1)

$ \frac{{\partial {\rho _{\rm{f}}}{{\overline u }_i}}}{{\partial t}} + \frac{{\partial {\rho _{\rm{f}}}\overline {{u_i}{u_j}} }}{{\partial {x_j}}} = - \frac{{\partial \overline P }}{{\partial {x_i}}} + \mu {\nabla ^2}{\overline u _i} + \frac{{{A_{\rm{R}}}}}{{\partial {x_j}}}。$

(2)

式中：ρ_f表示黏性不可压缩流体密度；u_i和u_j分别表示i和j方向上速度的时间平均值；x_i和x_j分别表示i和j方向上的分量；P、t和μ分别表示平均压力、时间和运动黏度；A_R表示雷诺应力张量。

基于涡激振动发电装置的1-DOF二维模型，系统的控制方程可写为

$ m\ddot y + C\dot y + Ky = {F_y}。$

(3)

式中F_y是圆柱振子在y方向上所受到的流体力。

$ {F_y} = \frac{1}{2}{C_y}\rho {U^2}DL。$

(4)

式中：C_y是升力系数；ρ是流体的密度；U是水流速度；D是圆柱直径；L是圆柱长度。

结合式(3)和(4)，单自由度涡激振动系统的运动控制方程可写作

$ \ddot y + 2\xi {\omega _0}\dot y + \omega _0^2y = \frac{{{F_y}}}{m}。$

(5)

式中：ξ为阻尼比，若令阻尼为C且弹簧刚度为K，则$\xi = \frac{C}{{2\sqrt {K/m} }}$；ω₀为圆柱振子的固有频率，若令弹簧刚度为K，则$ {\omega _0} = \sqrt {K/m} $；m为圆柱质量。

1.3 计算域和边界条件

在单自由度涡激振动数值模拟中，为确保圆柱振子的运动不受边界的影响，并能利用清晰的尾迹进行流场分析，一般计算域的高度需不小于20D，需要预留出的尾迹区域应不小于22.5D^[12]。

经综合考量后，圆柱振子仿真计算域的设计尺寸如图 2所示。圆柱振子的中心距入口边界(Inlet)和出口边界(Outlet)分别为10D和30D，距上、下边界(Symmetry)均为10D(上、下边界间距为20D)，圆柱振子周围的嵌套边界(Overset)直径为3D。

1.4 网格划分

本文使用能够解决边界运动问题的重叠网格技术对单自由度圆柱振子涡激振动进行仿真。首先划分包裹圆柱振子的前景网格和表示流体域的背景网格，如图 3(a)所示。然后对两个网格重叠区域的嵌套边界(Overset)处进行插值，插值网格如图 3(b)所示。最后，由于用户自定义程序(UDF)能够使信息在两个网格之间互相传递，因此实现了圆柱振子的双向流固耦合仿真。

重叠网格技术容易实现动网格，但前提是网格必须均匀且高质量，以确保信息传输的准确性。为保证良好的网格质量，对整个流体域采用结构网格进行划分，并设置圆柱振子表面的边界层网格如图 4所示。

1.5 网格无关性验证

为了研究网格划分和时间步长对单自由度圆柱振子涡激振动模型仿真结果的影响，采用三种网格密度对圆柱振子模型进行划分，并比较了多种情况下的最大无量纲振幅、升力系数均方根和阻力系数平均值。

6种工况下所得到的结果如表 1所示，其中约化速度U_r为无量纲常数。当时间步长减小时，数据趋于稳定。考虑到计算精度和计算量，本文选取第5种划分方式，即时间步长为0.002 s，圆周节点数为200的网格划分方式。

1.6 模型方法验证

为了验证单自由度涡激振动模型，本文采用了Stappenbelt等^[5]的实验数据：圆柱直径D=0.055 4 m；圆柱质量(包括封闭流体)m=3.502 kg；系统阻尼比ξ=0.006；系统刚度为453 N/m；质量比m^*=2.36。将这些数据输入到单自由度仿真模型中进行数值模拟，从而仿真得到圆柱振子的最大无量纲振幅，然后同实验结果进行对比验证，如图 5所示。

通过对比结果可知，质量比在2.36的情况下，随着约化速度的提高，无量纲振幅先增大后减小，共振区域大致出现在约化速度U_r的范围为6~8，与文献[5]的实验结果相比，振幅的变化曲线同文献[5]的实验数据基本吻合，验证了模型的可靠性。

2 仿真结果及分析

本文使用CFD和重叠网格技术，基于Modir等^[13]的实验模型进行研究，具体数据如表 2所示。研究包括定常流速和时变流速两种情况，其中定常流速下的研究用于确定系统的最佳刚度，并为时变流速的研究提供数据参考。时变流速的研究主要分析不同刚度下圆柱振子的振动特性和尾涡特征。

首先，进行定常流速下的研究，在不同流速下，分别对弹簧刚度为125、365和495 N/m的涡激振动系统进行数值模拟，得到圆柱振子无量纲振幅A/D随流速的变化规律，具体结果如图 6所示。对于以上3种弹簧刚度的圆柱振子，其最大无量纲振幅对应的流速随弹簧刚度的减小而逐渐降低，并且每种流速下都存在最佳的弹簧刚度。

其次，依据定常流速的研究结论，调整弹簧刚度系数以控制圆柱振子。具体的设置方案分为2种工况：流速随时间下降(工况1)；流速随时间上升(工况2)。每种工况包括5个阶段，每个阶段对应着不同的流速和弹簧刚度，如表 3和4所示。

2.1 振幅响应

本文对单自由度涡激振动系统的圆柱振子进行了双向流固耦合仿真，得到了两种工况下圆柱振子的位移时程曲线，如图 7和8所示。

工况1下的圆柱振子的位移时程曲线如图 7(a)和(b)所示：第一阶段(0~40 s)的流速为0.8 m/s，且在接近30 s时圆柱振子从静止状态进入共振状态，振幅达到0.044 m；40 s时进入第二阶段，流速降低至0.6 m/s，圆柱振子振幅随流速的降低而减小；80 s时进入第三阶段，弹簧刚度系数调整为365 N/m，共振的流速区间发生变化，圆柱振子振幅由0.04 m提高到约0.042 m；120 s时进入第四阶段，流速变为0.4 m/s，圆柱振子振幅明显降低；160 s时进入第五阶段，弹簧刚度系数调整为125 N/m，此时涡激振动系统在0.4 m/s的流速下，圆柱振子振幅重新增加，约为0.04 m。

工况2下圆柱振子的位移时程曲线如图 8(a)和(b)所示：第一阶段(0~40 s)，圆柱振子经历起振并最终达到共振状态，振幅达到0.04 m；40 s时进入第二阶段，流速变为0.6 m/s，此时流速增加，圆柱振子进入下分支，振幅减小；80 s时进入第三阶段，弹簧刚度系数调整为365 N/m，涡激振动系统的共振区间变大，在0.6 m/s的流速下实现共振，圆柱振子振幅显著增加，且增加到0.043 m；120 s时进入第四阶段，流速提高到0.9 m/s，圆柱振子再次进入下分支，振幅减小；160 s时进入第五阶段，弹簧刚度系数变为125 N/m，涡激振动系统再次达到共振状态，最大振幅约0.04 m。

2.2 频率响应

Bernitsas等^[6]设计开发了一套VIVACE涡激振动发电装置，并用此装置研究了发电功率的影响因素。研究结果表明，圆柱振子振幅和振动频率的乘积越高，越有利于发电装置对于能量的高效利用。因此，对圆柱振子位移时程曲线进行傅里叶分析，得到5个阶段下圆柱振子的频率-振幅图，如图 9和10所示。

在工况1(见图 9)下：当流速由0.8 m/s升至0.6 m/s或由0.6 m/s升至0.4 m/s时，圆柱振子振动频率逐渐降低，频率比也随之减小。此外，随着弹簧刚度的减小，圆柱振子振动频率也会减小。然而，随着涡激振动系统中圆柱振子固有频率的降低，频率比反而增大，其数值接近于1，此时振幅和振动频率的乘积较大。

在工况2(见图 10)下：当流速由0.4 m/s升至0.6 m/s时，圆柱振子的振动频率明显提高，而振幅减小；当流速由0.6 m/s升至0.9 m/s时，圆柱振子的振幅降低，振动频率区间减小，振动频率大小基本保持不变。弹簧刚度的变化会导致圆柱振子振动频率的变化，进而使振幅增大。此时，随着涡激振动系统中圆柱振子固有频率的改变，当频率比接近1时达到较高的振幅和频率乘积较大。

2.3 尾涡脱落模式

在时变流场中，弹簧刚度的变化对圆柱振子尾涡脱落模式产生影响，如图 11和12所示。

工况1(见图 11)下：第一阶段，圆柱振动稳定后，其尾涡脱落模式为“2P”模式，如图 11(a)所示；当圆柱振子由第二阶段进入第三阶段时，弹簧刚度由495 N/m降至365 N/m，尾涡之间的距离增加，但尾涡脱落模式并未发生变化，如图 11(b)和(c)所示；第四阶段中，圆柱振子的尾涡脱落模式为“2S”模式，即在一个振动周期内，圆柱振子尾流方向形成2个尾涡，如图 11(d)所示；第五阶段，弹簧刚度由365 N/m降至125 N/m，圆柱振子的尾涡产生分离，形成了两对方向相反的尾涡，尾涡脱落模式变为“2P”模式，即在一个振动周期内，圆柱振子的尾流方向形成有2对尾涡, 如图 11(e)所示。对于单自由度涡激振动系统而言，不同的尾涡脱落模式对应着不同的振动响应分支，其中“2P”模式对应着更大振幅的响应。

工况2(见图 12)下：第一阶段，圆柱振子流速处于刚度为125 N/m的共振区间，其尾涡脱落模式为“2P”模式，如图 12(a)所示。第二阶段，流场流速发生改变，导致圆柱振子脱离了共振区间，进入了下分支，圆柱振子的尾涡进行了融合，尾涡脱落模式转变为“2S”模式，如图 12(b)所示。第四阶段，流速再次发生改变，导致圆柱振子每一对尾涡的间距明显减小，处于“2P”模式的临界区域，见图 12(d)所示；第三、五阶段，随着弹簧刚度增大，圆柱振子的共振区间得以增大，使得其在0.6和0.9 m/s流速时的尾涡间距增大，并重新进入共振状态，尾涡脱落模式转变为“2P”模式，如图 12(c)和(e)所示。

3 结论

本文基于CFD模型和重叠网格技术，对时变流场中不同刚度的圆柱振子涡激振动进行了流固耦合分析，并探究了单自由度涡激振动系统的振动响应规律，得出以下结论：

(1) 弹簧刚度范围为125~495 N/m时，单自由度圆柱振子可以通过调节弹簧刚度系数来维持不同流速下的共振状态，并保持其横流向振幅的最大值。

(2) 改变弹簧刚度可以使圆柱振子的振动频率靠近系统的固有频率1 s^-1，以此维持共振状态下的频率比获得最大的发电率。

(3) 通过调整单自由度涡激振动系统的弹簧刚度，可使圆柱振子的尾涡脱落模式由“2S”模式转变为“2P”模式，增大尾涡间距，从而使圆柱振子的振动维持在最大振幅处。