世界经济持续高速发展，人类对于能源的需求也急剧增加，常规化石能源不断枯竭，同时环境恶化及节能减排的压力等问题日益突出，世界各国都在不断寻找可替代的可再生清洁能源。潮流能作为绿色清洁可再生能源，因其具有较好的可再生性、可预测性和可观的蕴藏量等特点，逐渐成为世界各国研究开发的重点。基于涡激振动的水流动能转换装置(VIVACE)^[1]因其可在低流速水流条件下有效提取水流动能，得到了学术界和工程界的广泛关注。后来Bernitsas团队^[2]采用在圆柱表面安装砂纸以改变表面粗糙度的方法来增强涡激振动，研究过程中发现通过合理的改变圆柱表面粗糙度的分布，不仅可以增强涡激振动，还可使柱体在“频率锁定”结束之后，使柱体发生驰振，从而拓宽了柱体高幅值振动的流速范围。驰振^[3]是流动分离产生的流体动力负阻尼引起的结构失稳式振动，具有大振幅、低频率的特点，主要与截面形状、水流攻角、流速以及质量、刚度、阻尼等力学特性有关，涡激振动自限制特性，而发生驰振时结构振幅随流速的增加不断增大。因此，Bernitsas等^[4]将VIVACE装置的概念拓展为基于流激振动的能量转换装置。

柱体截面几何形状是影响流激振动的一个重要因素，柱体流激振动特性会因截面形状的不同而存在很大不同^[5]。目前绝大部分关于柱体流激振动的研究主要集中在圆柱上，而对其他形状柱体流激振动的研究还比较少。徐枫等^[6]针对不同建筑物在风作用下的动力特性，对圆柱、方柱、正三角形柱和正六边形柱在低雷诺数(Re=200)下的流激振动进行了数值模拟研究。邓见等^[7]对方柱绕流及流激振动开展了数值模拟研究。Nemes等^[8]对不同攻角下方柱的流激振动进行了研究。针对VIVACE装置，Ding等^[9]、李恒^[10]、Haque^[11]等开展了不同截面形状柱体的流激振动数值模拟研究。燕翔等^[12]、Lian等^[13]围绕不同截面形状柱体的流激振动能量转换特性开展了实验研究。

虽然针对不同截面形状柱体流激振动问题，已有一些学者开展了研究工作，但不同的研究者的结论存在较大差异。例如，Lian等^[13]的实验结果表明圆柱与三角柱的能量转换效果高于方柱及菱形柱，但李恒^[10]研究所得的柱体能量转换效率为菱形柱最高，其次为方柱、三角柱、梯形柱。影响流激振动的参数众多，由于其中很多参数是相关的，在研究中很难保持所有的参数一致，因而上述学者的研究在改变柱体截面形状的同时，一些重要影响参数同时也发生了变化，这使得通过比较不同截面形状柱体结果而得到的结论具有较大的不确定性。例如，Ding等^[9]和Haque^[11]的研究中，不同截面形状柱体，其质量比不尽相同，但根据以往的研究^[14-15]，质量比对频率锁定的发生、频率锁定的约化流速范围以及振幅都有重要影响。

综上所述，虽然已有一些学者对不同截面形状柱体流激振动现象及能量获取效果进行了研究，但是仍然存在很多问题需要进一步深入研究，特别是高雷诺数、高阻尼下柱体流激振动特性及能量转换特性的研究还非常匮乏。本文选取圆形、方形、类梯形、梯形和T字形5种不同截面形式的柱体进行研究，开展高雷诺数、高阻尼下柱体流激振动数值模拟。为了更好地比较不同截面形状的柱体的流激振动特性和能量转换特性，对不同形状柱体的参数进行精心设计，使不同截面形状柱体的质量比、阻尼比、同一流速下的雷诺数和约化速度等重要无量纲参数均保持一致。

1 流体与结构模拟参数

圆柱、方柱、类梯形柱、梯形柱和T字形柱5种不同截面形状的柱体见图 1所示，图中D为截面特征长度，U为来流速度。为了使各柱体的主要无量纲参数保持一致，对各截面形状特征进行分析，在保持截面特征长度D不变的情况下，同一流速对应的雷诺数$\mathit{Re}\left( {Re = \frac{{UD}}{v}, v} \right.$为流体的运动粘度)也相同。在采用同一材料的情况下，即柱体的材料密度不变，所有柱体单位长度的结构质量m和附加质量m_a的比值将不变，即不同截面形状柱体的质量比m^*(m^*=m/m_a)保持不变；在采用同一材料及选取相同弹簧刚度的情况下，通过使各柱体截面面积相等可使不同截面形状柱体的自振频率保持一致，进而使同一流速下约化速度保持不变。对于方柱，由于难以在特征长度D不变的情况下，使柱体截面面积与圆柱的相同，本文通过调整弹簧刚度使方柱的自振频率与其他截面形状的柱体保持一致。

基于分析，各柱体截面几何尺寸如图 1所示，其中取D=0.1 m，其他尺寸根据使各柱体截面面积相等确定。各柱体材料密度均相同ρ_s=1 880 kg/m³，水的密度取ρ=1 000 kg/m³，从而质量比m^*为1.88(取静水附加质量系数为1.0)。除方柱的弹簧刚度有所调整之外，其他形状柱体的弹簧刚度均保持相同，为$\sqrt[2]{{{{21}_1}}} = 600{\rm{N}}/{\rm{m}}$，方柱的弹簧刚度为$\sqrt[2]{{{{21}_1}}} = 764{\rm{N}}/{\rm{m}}$，使所有柱体的水下自振频率f_n均为0.820 3 Hz，进而使同一流速对应的约化速度U_r(U_r=U/f_nD)均相同。为了更贴近于实际高阻尼情况，取阻尼比ζ为0.12。考虑到我国河流和海洋中的水流流速绝大多数较低，计算流速范围取为0.2~1.3 m/s，相应的Re范围为20 000~130 000(取υ=1.0×10^-6)。

2 流激振动数值模型

采用连续性方程与雷诺平均N-S方程对流场进行模拟求解，其方程组表示为：)

$ \partial {u_i}/\partial {x_i} = 0 $

(1)

$ \frac{{\partial {u_i}}}{{\partial t}} + {u_j}\frac{{\partial {u_i}}}{{\partial {x_j}}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}} + v\frac{{{\partial ^2}{u_i}}}{{\partial x_j^2}} - \frac{{\partial \overline {u_i^\prime u_j^\prime } }}{{\partial {x_j}}} $

(2)

式(1)、(2)为张量表达式，下标等于x和y，分别表示来流方向和垂直于来流方向；ρ为流体密度; p为压力；u_x、u_y分别为x、y方向的流体平均速度；u′ _x、u′ _y分别为沿x、y方向的流体脉动速度。考虑到数值模拟涉及到复杂的截面几何形状，对柱体近壁边界、边界层及自由剪切层的处理要求较高，故采用SSTk-ω湍流模型^[16](剪切应力输运模型)对方程组进行封闭。SST k-ω湍流模型是Menter^[16]结合k-ω湍流模型和k-ε湍流模型开发而成，能够较好地处理近壁小尺度问题。

柱体振动采用二阶振动微分方程描述为：

$ m\ddot y + c\dot y + \lambda y = {F_y} $

(3)

式中：m为振动结构质量；c为系统阻尼；λ为弹簧刚度；$y、\dot y、\ddot y$分别为结构振动的位移、速度和加速度；F_y为作用在柱体上的流体升力。

选定50D×20D的区域为计算域，柱体中心距离左侧入流边界10D，距离右侧出流边界的距离为40D，到上下边界的距离均为10D。考虑到是湍流环境，入流边界在给定来流速度的同时，还给定了湍流动能、流速及耗散率；出流边界给定压力和零压力梯度；上下边界为滑移固壁边界，柱体表面为无滑移的流固耦合边界，即满足：位移协调条件d_f-d_s=0、力平衡条件τ_f-τ_s=0和无滑移的速度平衡条件u-u_s=0，其中，d_f、d_s分别为流体和固体位移量，τ_f、τ_s分别为流体和固体的应力，u、u_s分别为流体和固体的速度。

采用有限体积法对计算域进行离散。考虑到雷诺数较高，在进行计算域离散时，对柱体周围及尾流区等参数梯度变化较大的区域进行了网格加密，网格过密势必使计算量和所需的存储量增大，使时间步长减小，造成计算耗时不能接受。各柱体计算域的网格划分原则基本一致，但考虑到截面形状存在较大差异，对柱体边界附近长宽为20D的正方形区域内采用了不同离散的分块方法，因而网格分布规则和数量也都不同。通过试算，圆柱、方柱、类梯形柱、梯形柱和T字形柱的流体域网格数分别确定为52 000、27 500、48 000、37 800、33 800。图 2为圆柱和T字形柱边界附近的网格划分情况。流固耦合的求解采用双向迭代耦合求解方法，为保证柱体发生较大位移时流体计算域的网格质量，采用任意拉格朗日欧拉(arbitrary Lagrange-Euler，ALE)法根据流固耦合求解过程中边界或物体的运动、变形，对每个时间步迭代之前的流体计算域网格进行更新和重构。

为了验证数值模型的正确性，将数值计算结果与Lee等^[17]在高阻尼、高雷诺数下的圆柱涡激振动实验结果进行比较，实验圆柱直径D=0.088 9 m，弹簧刚度λ=400 N/m，质量比m^*=1.927，阻尼比ζ=0.12。数值计算结果与实验结果的比较见图 3，可以看出数值计算得到振幅和频率随流速的变化规律均与实验较为一致，且最大相对误差均小于15%，两者吻合较好。

3 能量转换

假设柱体的振动为简谐振动，柱体振动位移可表示为：

$ y(t) = A\sin (2{\rm{ \mathit{ π} }}ft) $

(4)

式中：A为振幅；f为振动频率。单个振动周期内流体对柱体所做的功，即柱体获取的功率为：

$ P_{c}=\frac{1}{T} \int_{0}^{T} F_{y} \dot{y} \mathrm{d} t $

(5)

根据式(3)等号两边积分应相等，有：

$ {P_c} = \frac{1}{T}\int_0^T {(m\ddot y + c\dot y + \lambda y)} \dot y{\rm{d}}t = 8{{\rm{ \mathit{ π} }}^3}m\zeta {(Af)^2}{f_n} $

(6)

式中：ζ为系统阻尼比；$\zeta = \frac{c}{{2\sqrt {m\lambda } }};{f_n}$为静水中的自振频率。流体流动的功率可表示为：

$ {P_w} = \rho {U^3}DL/2 $

(7)

式中L为柱体长度。于是，能量转换效率为：

$ \eta = \frac{{{P_c}}}{{{P_w}}} = \frac{{16{{\rm{ \mathit{ π} }}^3}m\zeta {{(Af)}^2}{f_n}}}{{\rho {U^3}DL}} $

(8)

4 模拟结果与分析

不同截面形状柱体的模拟计算参数已在第1节中给出，图 4为不同截面形状柱体的振幅比和频率比随约化速度的变化。从图中可以看出，圆柱、类梯形柱、梯形柱和T字形柱均在U_r=2.4附近起振，方柱在U_r=3.7附近起振。圆柱的幅值曲线在U_r=6.1附近存在一个峰值，具有明显的自限制特性，表现出典型的涡激振动特征。对于其他柱体，振幅比均随U_r的增大而增大，不存在峰值，表现为驰振的振动特征，方柱、类梯形柱、梯形柱及T字形柱的振幅比在计算流速范围内总体上表现为不断增大的趋势，似乎不存在极限值，不过T字形柱、梯形柱和方柱分别在U_r=12.2(Re=100 000)、9.8(Re=80 000)和11.0(Re=90 000)附近出现了一个减小的趋势，经分析，尾涡的脱落形态在此时发生了变化，如从图 10可以看出，雷诺数由70 000增加至90 000时梯形柱尾涡脱落模式由2P模式变为了2S+2P模式。在U_r < 12.2范围内同一流速下T字形柱振幅比较其他柱体的为大，且随U_r的增大速率最大，类梯形柱和梯形柱次之，方柱较为平缓。在U_r < 6.1时，同一流速下各柱体的振幅比较为接近，U_r>7.5后方柱、类梯形柱、梯形柱和T字形柱的振幅比明显较圆柱的为大。由图 4(b)可以看出，在同一流速下，圆柱的频率比较其他柱体的为大。圆柱的频率比在发生涡激振动时小于1.0，即振动频率较自振频率小，随着约化速度的增大，频率比逐渐增大，并大于1.0，即振动频率大于自振频率，在U_r=6.1时达到1.4，之后在6.1 < U_r < 11.0的范围内基本保持不变，有趣的是在U_r=11.0附近频率比突然迅速增大至2.3左右，在U_r=13.4时达到了2.5，圆柱振动频率远大于自振频率，这一现象在Khalak等^{[15, 18]}和Govardhan等^[19]的低质量比、低阻尼涡激振动实验中亦有所观察到。其他4种柱体的频率比曲线与圆柱较为相似，但均小于1.0，即振动频率小于自振频率，均未发生共振，频率比在U_r>6.1后基本保持不变。

柱体的振动是由于作用在柱体上的升力引起的，升力系数C_y定义为：

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

(9)

式中F_y为垂向方向上的升力，随时间脉动，如图 5。图 6给出了不同截面形状柱体的升力系数幅值随U_r的变化曲线，图中C_{y, max}为升力系数幅值。各柱体的升力系数幅值曲线形状较为相似，均存在一个峰值，即随着约化速度的增加升力系数幅值先增大后减小，不同的是，圆柱、方柱和T形柱升力系数幅值曲线较类梯形柱和梯形柱的尖，即增大或减小的速率均较快。同一流速下T字形柱升力系数幅值均较其他柱体的为大，圆柱、方柱、类梯形柱、梯形柱和T字形柱的最大升力系数幅值分别为4.07、2.58、2.88、3.04和4.81。

图 7~11给出了5种柱体在Re=30 000(U_r=3.7)、50 000(U_r=6.1)、70 000(U_r=8.5)和90 000(U_r=11.0)4个雷诺数下的尾涡形态。Re=30 000时，除了方柱没有出现明显的漩涡脱落之外，其他柱体尾涡脱落均呈现2S模式，即一个振动周期内截面上下各有一个单涡脱落；Re=50 000时，圆柱和方柱的尾涡脱落模式为2S，类梯形柱尾涡以2P模式脱落，即一个振动周期内截面上下各有一对涡脱落，梯形柱尾涡以S+P模式脱落，即一个振动周期内含有一个单涡和一对涡，T字形柱尾涡脱落模式为2P；Re=70 000时，圆柱尾涡脱落模式为2P，方柱尾涡呈2S模式脱落，类梯形柱尾涡以S+P模式脱落，梯形柱尾涡以2P模式脱落，T字形柱尾涡脱落模式为2P；Re=90 000时，圆柱尾涡模式为2S模式，方柱尾涡以S+P模式脱落，类梯形柱、梯形柱及T字形柱尾涡呈明显的2S+2P模式脱落形态，即一个振动周期内有6个涡脱落。

由式(7)、(8)可知，柱体由于振动所获取的机械能功率与柱体质量、阻尼、自振频率、振幅和频率有关，而能量转换效率除了与上述因素有关外，还与流体密度、来流速度以及柱体在水流方向的投影面积有关。图 12分别给出了不同截面形状柱体的获能功率及转换效率与U_r的关系曲线。由于各柱体频率在U_r>6.1后基本保持不变，功率曲线的形状与振幅曲线较为一致，圆柱功率在U_r=6.1附近出现峰值，其他柱体的功率均随流速(或Re、U_r)的增大而增大。在U_r < 12.2的范围内同一流速下T字形柱获得的功率均较其他柱体获得的为大。虽然方柱的质量比其他柱体的为大(方柱截面面积较大)，但由于振幅和频率较其他柱体的小，使得其在大部分流速下的功率较其他柱体的小。虽然圆柱最大功率较T字形柱、类梯形柱及梯形柱的小，但其产生最大功率的流速较其他柱体的小，即在保持质量、阻尼和自振频率相等情况下，圆柱产生最大功率所需的流速较其他柱体小，不过需要注意的是，在圆柱产生最大功率的流速下T字形柱的功率较圆柱的大。由图 12(b)可以看出，各柱体能量转换效率曲线均存在一个峰值，5种柱体中T字形柱的最大能量转换效率为最大，达到了42.5%，其他依次为圆柱、类梯形柱、梯形柱、方柱，分别为27%、15.4%、11.45%、3.8%。虽然除圆柱之外的4种柱体的功率均随流速(或Re、U_r)的增大而增大，但由于流速的增大使流体功率迅速增大，从而导致能量转换效率在达到峰值后逐渐减小。圆柱和方柱的最大能量转换效率均出现在U_r=6.1附近，类梯形柱和梯形柱的均出现在U_r=9.8附近，T子形柱的出现在U_r=7.3附近，表明为了达到最大能量转换效率，圆柱和方柱所需要的流速最小，其次为T字形柱。由式(8)可知，能量转换效率与振幅和频率的乘积的二次方成正比，与流速的三次方成反比。圆柱振动为涡激振动，由于其频率较大，且振幅在U_r < 6.1时与其他柱体的较为接近，因而在低流速下较容易达到最大能量转换效率。另一方面，T字形柱振动为驰振，振幅较大，频率随较圆柱的小，但较其他柱体的为大，因而也较容易在低流速下达到最大能量转换效率。

5 结论

1) 圆柱的振动表现为典型涡激振动特征，具有明显的自限制特性。其他4种柱体的振动均存在驰振，振幅随流速增大而增大。在同一流速下T字形柱的振幅比较其他柱体大，最大达到了2.0，其随U_r的增大速率也最大，方柱最为平缓。在同一流速下，圆柱的频率比较其他柱体的为大，各柱体频率比在U_r < 6.1时均随约化速度的增大而增大，圆柱最大频率比达到了1.5左右，其他柱体的频率比均小于1.0。

2) 各柱体的升力系数幅值曲线形状较为相似，均存在一个峰值，圆柱、方柱和T形柱升力系数幅值曲线较类梯形柱和梯形柱的尖。在U_r < 12.2范围内同一流速下T字形柱升力系数幅值均较其他柱体的为大。数值模拟观察到了2S、2P、S+P和2S+2P等尾涡脱落模式。

3) 类梯形柱和梯形柱的截面形状较为接近，其振动特征也较为相似，表明截面形状较小的变化并不会引起振动特征的剧烈变化。

4) 功率曲线的形状与振幅曲线较为一致，圆柱功率在U_r=6.1附近出现峰值，其他柱体的功率均随流速(或Re、U_r)的增大而增大。在U_r < 12.2的范围内同一流速下T字形柱获得的功率均较其他柱体获得的为大。圆柱产生最大功率所需的流速最小。各柱体能量转换效率曲线均存在一个峰值，5种柱体中T字形柱的最大能量转换效率为最大，达到了42.5%，其他依次为圆柱、类梯形柱、梯形柱、方柱，分别为27%、15.4%、11.45%、3.8%。为了达到最大能量转换效率，圆柱和方柱所需要的流速最小，其次为T字形柱。从分析结果来看，T字形柱和圆柱对于低流速下的能量转换较为有利。