海洋能是目前世界上为数不多的可再生能源种类，世界各国对其开发利用给予了足够的重视。潮流能作为海洋能的重要组成部分之一，有能量密度高、可预测性强及相对稳定等优点，据统计全球中可利用开发的潮流能总量达3×10¹¹ kW^[1-4]。根据潮流能捕获形式的不同，目前主要有叶轮式和非叶轮式潮流能装置^[5]，其中叶轮式潮流发电装置是目前较为流行的不能方式，而根据其叶轮旋转方向，可以将潮流能叶轮分为水平式和垂直式叶轮^[6-8]。近几十年来，垂直轴潮流能叶轮在“理想潮流环境”中的水动力研究理论及方法已经比较成熟。Coiro^[9]针对Koblod垂直轴叶轮，建立了专门的双盘面多流管模型，进而计算了其该水轮机的效率。Camporeale等^[10]和张亮等^[11]计算了不同型式水轮机的推力、效率等相关性能参数，分别使用流管法和模型实验对其进行了研究，发现利用流管法计算水轮机性能时准确度较高。再此基础上，Li^[12]提出了离散涡模型，可以用来计算多个叶轮之间的相互干扰。计算流体力学(computational fluid dynamics, CFD)方法在垂直轴叶轮方面的应用已经取得了很多研究成果。李志川等^[13]将CFD方法运用到垂直轴叶轮水动力计算中，并分别讨论了不同偏角形式水轮机，重点研究了叶片数、安装角等对叶轮的水动力性能影响，并通过模型试验证明了CFD方法在垂直轴叶轮研究上的可行性。王凯等^[14-15]采用ANSYS软件对两叶片立轴叶轮的艏摇和纵摇运动进行了模拟分析，研究了艏摇和纵摇对浮式立轴水轮机的水动力性能影响。张万超等^[16]运用CFD方法，模拟了垂直轴叶轮在均匀流下的受迫摇荡运动，研究了不同摇荡频率下的水动力载荷及相应水动力系数的变化特性。由于旋转机械下的波浪载荷模拟异常困难，波流作用下水轮机水动力的研究非常有限。波流联合作用下的，垂直轴叶轮旋转的同时会伴随波浪摇荡，其运动响应非常复杂，使用常规的CFD方法及软件很难对此问题进行直接模拟。流载荷作用下，叶轮发生旋转运动，而与其固结的浮式载体在波浪作用下会发生摇荡运动。

为了便于研究波浪中垂直轴潮流能叶轮的运动响应及水动力问题，本文将其复杂模型进行简化，假定叶轮在均匀流中作受迫微幅简谐运动，研究不同横摇幅值对浮式立轴潮流能叶轮的水动力影响。

1 计算流体力学数值模拟 1.1 理论基础

真实浪流环境下，浮式立轴叶轮的运动如图 1所示，为便于研究，定义密度为ρ的流域中来流速度V的方向为x轴正向。

针对浮式垂直轴叶轮在浪流作用下的运动工况，为了方便分析，定义无量纲参数为：

叶轮推力系数：

$ {C_{Fx}} = {F_x}/(0.5\rho {V^2}DH) $

(1)

叶轮侧向力系数：

$ {C_{Fy}} = {F_y}/(0.5\rho {V^2}DH) $

(2)

叶轮能量利用率系数：

$ {C_P} = Q{\omega _0}/(0.5\rho {V^3}DH) $

(3)

叶尖速比：

$ \lambda = R\omega /V $

(4)

无因次横摇角速度：

$ \bar \omega = uL/V = LA\omega {\rm{cos}}\left( {\omega t} \right)/V $

(5)

无因次横摇角加速度：

$ \bar a = \dot uL/V = - A{\omega ^2}{L^2}{\rm{sin}}\left( {\omega t} \right)/{V^2} $

(6)

式中：D、R和ω₀分别为叶轮直径、半径和旋转角速度；ω、L、u和A分别为横摇频率、臂长、速度与幅值；Z、θ、C和H分别为叶片数量、位置角、弦长、展长；F_x、F_y和Q分别为叶轮受到的推力、侧向力及转矩。

式(1)、(2)中水动力系数可以分为3项：

$ \left\{ \begin{array}{l} {C_{Fx}} = {C_x} + {n_x}\bar \omega + {m_x}\bar a\\ {C_{Fy}} = {C_y} + {n_y}\bar \omega + {m_y}\bar a \end{array} \right. $

(7)

式(7)中n_x、n_y、m_x、m_y用三角级数形式展开:

$ \left\{ \begin{array}{l} {n_x} = \frac{{\delta {C_{Fx}}}}{{\delta \bar \omega }} = {\rm{ }}{n_{x, 0}} + \sum\limits_{k = 1}^\infty {{n_{x, k}}} {\rm{sin}}\left( {kz\theta \left( t \right) + {\psi _{\omega x, k}}} \right)\\ {n_y} = \frac{{\delta {C_{Fy}}}}{{\delta \bar \omega }} = {\rm{ }}{n_{y, 0}} + \sum\limits_{k = 1}^\infty {{n_{y, k}}} {\rm{sin}}(kz\theta \left( t \right) + {\psi _{\omega y, k}})\\ {m_x} = \frac{{\delta {C_{Fx}}}}{{\delta \bar a}} = {\rm{ }}{m_{x, 0}} + \sum\limits_{k = 1}^\infty {{m_{x, k}}} {\rm{sin}}(kz\theta \left( t \right) + {\psi _{ax, k}})\\ {m_y} = \frac{{\delta {C_{Fy}}}}{{\delta \bar a}} = {\rm{ }}{m_{y, 0}} + \sum\limits_{k = 1}^\infty {{m_{y, k}}} {\rm{sin}}(kz\theta \left( t \right) + {\psi _{ay, k}}) \end{array} \right. $

(8)

最后推力系数可以写成下列级数形式:

$ \begin{array}{l} {C_{Fx}} = {C_x} + {n_x}\bar \omega + {m_x}\bar a = \\ \;\;\;\;\;\;\;\;{C_{x, 0}} + \sum\limits_{k = 1}^\infty {{C_{x, k}}} {\rm{sin}}(kz\theta \left( t \right) + {\psi _{x, k}}) + \\ \;\;\;\;\;\;\;\;\left\{ {{n_{x, 0}} + \sum\limits_{k = 1}^\infty {[{n_{x, k}}{\rm{sin}}(kz\theta \left( t \right) + {\psi _{\omega x, k}})]} } \right\}\bar \omega + \\ \;\;\;\;\;\;\;\;\left\{ {{m_{x, 0}} + \sum\limits_{k = 1}^\infty {[{m_{x, k}}{\rm{sin}}(kz\theta \left( t \right) + {\psi _{ax, k}})]} } \right\}\bar a \end{array} $

(9)

侧向力系数可以写成下列级数形式：

$ \begin{array}{l} {C_{FY}} = {C_y} + {n_y}\bar \omega + {m_y}\bar a = \\ \;\;\;\;\;\;\;\;\;{C_{y, 0}} + \sum\limits_{k = 1}^\infty {{C_{y, k}}{\rm{sin}}(kz\theta \left( t \right) + {\psi _{y, k}}) + } \\ \;\;\;\;\;\;\;\;\;\left\{ {{n_{y, 0}} + \sum\limits_{k = 1}^\infty {[{n_{y, k}}{\rm{sin}}(kz\theta \left( t \right) + {\psi _{\omega y, k}})]} {\rm{ }}} \right\}\bar \omega + \\ \;\;\;\;\;\;\;\;\;\left\{ {{m_{y, 0}} + \sum\limits_{k = 1}^\infty {[{m_{y, k}}{\rm{sin}}(kz\theta \left( t \right) + {\psi _{ay, k}})]} {\rm{ }}} \right\}\bar a \end{array} $

(10)

式(10)和(11)中的C_Fx和C_Fy可以通过CFD软件计算得到，再通过最小二乘法拟合得到公式中的各项。

1.2 立轴水轮机模型

本文使用翼型为NACA0018的垂直轴三维叶轮模型，叶轮直径D为0.8 m，叶片数Z=2，展长H=0.6 m，弦长C=0.12 m，叶片属于固定偏角形式，安装角为0°。

图 2所示为立轴叶轮横摇运动的计算网格模型。总网格计算域为长35D, 宽和高都为25D的长方形域。叶片周围第1层网格高度为0.000 1，满足第1层网格Y⁺ < 5(Y+是湍流模型中靠近壁面距离的无量纲量), 旋转域、横摇域和总域之间用滑移网格连接，在计算中网格不存在变形，确保了计算过程中的网格质量，提高了计算精度。3个域之间的相对运动可以模拟立轴叶轮在浪流作用下的横摇运动。

图 3为立轴叶轮横摇运动计算模型的边界条件设置，设置方法可以参考文献[16]。

2 数值结果分析 2.1 CFD方法的可行性

图 4为均匀来流下立轴叶轮效率的试验值和CFD计算值得对比图，从图中可以看出实验值和计算值有较高的吻合度，误差较小，在速比2.5时效率达到最大值。通过对比试验和计算工况，计算值和实验值之间的误差可能是由实验中壁面效应和机械摩擦造成的，但是其不影响CFD计算方法在立轴叶轮数值模拟中的准确性。

2.2 横摇幅值对叶轮推力和侧向力的影响

本文研究的立轴叶轮横摇模型主要对比在来流速度1.0 m/s，横摇频率1.2 rad/s, 速比2.5工况下的3个不同横摇幅值(5°、10°和15°)的叶轮推力和侧向力时历曲线的变化规律。

图 5(a)和(c)为立轴浮式叶轮在均匀流无横摇运动下的推力和侧向力的时历曲线，其包络线表现为直线形式。图 5(b)和(d)为浪流作用下的立轴叶轮发生横摇运动时的时历曲线，从图中可以明显看出时历曲线存在周期性波动。特别是侧向力的上下包络线波动更加明显，叶轮的横摇运动方向和叶轮的侧向力方向重合因此侧向力更加明显。

图 6和图 7分别为3个横摇幅值下叶轮推力和侧向力的上下包络线对比图，其中0°代表的是立轴叶轮在均匀来流中无摇荡运动的叶轮推力和侧向力的上下包络线，呈直线形式。3个横摇幅值下的上下包络线呈周期性波动，波动频率和横摇频率一致，这是因为叶轮在浪流联合作用下发生横摇运动，从而引起叶轮推力和侧向力的周期性波动，因此力的波动频率和横摇频率一致。从4幅图中可以看出，推力和侧向力的上下包络线波动幅值都是随着横摇幅值增大而增大。

为了进一步对比横摇运动下立轴叶轮推力和侧向力的变化规律，分别定义推力和侧向力的平均值和波动幅值：

$ \left\{ \begin{array}{l} {{\bar C}_{Fx}} = \int_{{t_1}}^{{t_2}} {{C_{Fx}}{\rm{d}}t/({t_2} - {t_1})} \\ {{\bar C}_{Fy}} = \int_{{t_1}}^{{t_2}} {{C_{Fy}}{\rm{d}}t/({t_2} - {t_1})} \\ {{\hat C}_{Fx}} = {C_{Fx{\rm{max}}}} - {C_{Fx{\rm{min}}}}\\ {{\hat C}_{Fy}} = {C_{Fy{\rm{max}}}} - {C_{Fy{\rm{min}}}} \end{array} \right. $

(11)

图 8中3个横摇幅值下的推力和侧向力的平均值和波动值对比图。图 8(a)可以看出随着横摇幅值增大，横摇运动下叶轮推力和侧向力平均值基本不变，说明横摇运动对叶轮受力的平均值没有影响。图 8(b)中推力和侧向力的波动幅值随着横摇幅值增大而增大，这是因为叶轮横摇角速度会由于横摇幅值增大而增大，角速度增大从而导致叶轮受力增大。

2.3 不同横摇幅值的叶轮涡量图及叶片压力云图

图 9为叶轮横摇运动下3个不同横摇幅值的同一时刻涡量图。随着横摇幅值增大，叶轮尾流中涡量也越大，这是因为横摇幅值增大引起叶轮处合速度增大，从而引起尾涡增大。图 9(c)中叶轮下游尾涡明显多于图 9(a)中的尾涡。

在同一水平高度下分别定义横摇运动中叶轮叶片，漏出叶片长度短的为1号叶片，长的为2号叶片。横摇幅值15°的叶轮2号叶片明显高于5°横摇幅值2号叶片高度，这是叶轮横摇幅值不同造成的。随着横摇幅值增加，叶轮下叶片漏出的展长也越大。

图 10对应图 9中的2号叶片表面压力。从图 10中可知，叶片上压力分布形状相似但是压力值不同，随着横摇幅值增大，叶片表面压力逐渐增大。这是因为随着横摇幅值增大，叶轮同一时刻的横摇速度会增大，横摇速度和来流速度叠加导致叶轮盘面的合速度增大，叶轮横摇过程中受到的力也相应增大。

3 不同横摇幅值下的水动力系数

根据前面CFD计算的横摇幅值为5°、10°和15°工况下的叶轮推力和侧向力时历曲线，利用最小二乘法根据式(10)和式(11)对其进行拟合分析，得到所需的各项系数，并分析其变化规律。表 1为横摇叶轮推力C_Fx根据式(10)得到的各项系数，表 2为横摇叶轮侧向力C_Fy根据式(11)得到的各项系数。

从表 1可知，C_Fx展开系数中，C_{x, 0}、C_{x, 1}和ψ_{x, 1}表示均匀水动力项系数，n_{x, 0}、n_{x, 1}和ψ_{ωx, 1}表示阻尼项系数，m_{x, 0}、m_{x, 1}和ψ_{ax, 1}表示附加质量项系数。均匀水动力项各系数相差不大，说明不同横摇幅值下叶轮的均匀水动力值相当。ψ_{x, 1}、ψ_{ωx, 1}表明均匀水动力项和阻尼项发生周期性波动。从各项系数来看，随着横摇振幅的增加，阻尼项各项系数绝对值波动幅值不大，总体上微幅增加。附加质量项各系数绝对值在数值上远小于均匀水动力项和阻尼项。因此，横摇运动下叶轮推力的组成成分中均匀水动力和阻尼力占主要作用，附加质量力作用较小。

表 2中横摇侧向力展开各系数中，C_{y, 0}、C_{y, 1}和ψ_{y, 1}表示侧向力均匀水动力系数，n_{y, 0}、n_{y, 1}和ψ_{ωy, 1}表示侧向力阻尼力系数，m_{y, 0}、m_{y, 1}和ψ_{ay, 1}表示侧向力附加质量力系数。随着横摇幅值增加，均匀水动力各系数变化不大，说明横摇运动下叶轮受力中均匀水动力项不受横摇幅值影响。ψ_{y, 1}、ψ_{ωy, 1}和ψ_{ay, 1}侧向力中均匀水动力、阻尼力和附加质量力都发生周期性波动，且各项都占据一定比例。

4 结论

1) 立轴叶轮在浪流共同作用发生横摇运动时，叶轮的推力和侧向力发生周期性波动，波动频率和横摇频率一致。

2) 横摇运动下叶轮推力和侧向力波动幅值随着横摇幅值增大而增大。

3) 立轴叶轮在浪流作用下的横摇运动的受力平均值不受横摇幅值影响，但受力的波动幅值随横摇幅值增大而增大。

4) 横摇运动下叶轮推力和侧向力系数的各分解项发生周期性波动，推力项中附加质量力所占比例远小于均匀水动力和阻尼力，差一个数量级。侧向力项中3项所占比例在同一数量级。