当今社会，化石能源在能源比例中仍然占据很大比例，但其可开采量越来越少.所以，世界上很多国家已经开始重视对绿色可再生能源的开发研究，同时潮流能也因为它的可持续性、高密度和可预测性的特点，正逐渐引起人们的重视^[1-4].

潮流能水轮机是主要的潮流能开发装置，根据主轴和来流方向的关系可以分为水平轴和垂直轴(立轴)两种^[5].在较深的海况中，浮式平台是水轮机最好的载体，因为浮式平台不仅可以作为水轮机、发电机及其他设备的载体，而且浮式平台的维护、拆卸和安装更加便利^[6].但是浮式平台在波浪的作用下会发生6个自由度的摇荡运动，这种工况会对水轮机水动力性能有很大影响.

浮式海洋结构物在波浪作用下的水动力系数变化一直是研究者关注的重点，但是大多数海洋浮式结构都有着复杂的外形以至于对它们水动力性能的研究存在很大阻力，所以早期学者最开始都是用最简单的结构来研究波浪下的水动力问题，如圆柱体和球体.Yeung^[7]选取有限水深下的圆柱体来研究其水动力性能，通过匹配其内部和外部问题的特征函数来研究圆柱的升沉、摇荡和翻滚3类辐射问题，然后给出了圆柱的阻尼系数和附加质量系数.结果表明, 升沉附加质量是对数奇异的，而且在低频时阻尼系数趋近一个常数.Lopes等^[8]研究了半潜式球体的脉动运动形式下的水动力性能，结果表明在接近共振频率时，振荡形式下的阻尼系数和附加质量系数比升沉、摇摆形式下的系数要高，而且水深对阻尼系数的影响仅限于低频范围内.另外很多研究人员利用多级方法、特征函数法或者其他方法对简单的结构物(圆柱、球)在波浪下的水动力性能做了详细的研究^[9-11].

但是真正的海洋结构物有着复杂的外形，按照上述方法对它们在波浪中的水动力性能问题的求解存在很大困难.因此CFD方法和实验方法成为研究波浪中复杂外形的海洋结构物水动力性能的主要方法.Pavon等^[12]利用实验和CFD方法对比分析半潜式浮式风机垂荡板的水动力系数，实验模型直径为1 m，是已知最大的实验模型.实验测量了阻尼系数和附加质量系数，并和CFD的计算结果对比，为浮式风机的概念设计提供参考标准.2010年，Galloway等^[13]研究了三叶片水平轴水轮机在波浪中的性能影响，发现推力波动最大增加了37%、力矩增加35%；2013年，Lust等^[14-15]做了相似实验发现波浪中叶轮性能波动现象明显，但这些研究也仅限于水平轴水轮机.

到目前为止，对浮式立轴水轮在波浪中水动力性能影响还处于探索阶段.本文利用CFD技术研究立轴水轮机发生强迫纵摇运动情况下的水动力性能变化情况，针对不同速比下的水轮机工况，并和均匀流中的水轮机水动力性能对比，找出纵摇对水轮机水动力的影响规律；并采用最小二乘法原理^[16]拟合纵摇水轮机的推力和侧向力系数时历曲线，得到波浪中的水轮机阻尼系数及附加质量系数.

1 数值模拟 1.1 理论基础

垂直轴叶轮来流方向为X轴正向，表 1为立轴潮流能水轮机的主要参数.

为了方便分析，定义量纲一的参数如下：

叶轮推力系数${C_{FX}} = \frac{{{F_X}}}{{0.5\rho {v^2}DH}}$，叶轮侧向力系数${C_{FY}} = \frac{{{F_Y}}}{{0.5\rho {v^2}DH}}$，叶轮的能量利用率系数${C_P} = \frac{{Q\omega }}{{0.5\rho {v^3}DH}}$，叶尖速比$\lambda = \frac{{R\omega }}{v}$，无因次纵摇角速度$\bar \omega ' = \frac{{\omega 'R'}}{v}$，无因次纵摇角加速度$\bar a' = \frac{{\dot {\omega '}{{R'}^2}}}{{{v^2}}}$，纵摇角度ξ=Asin(ω_Zt).其中：水轮机推力系数和侧向力系数可分为均匀水动力项、阻尼项和附加质量项.

$\begin{array}{*{20}{c}} {{C_{FX}} = C_x^0 + {n_{xx}}\bar \omega ' + {m_{xx}}\bar a',}\\ {{C_{FY}} = C_y^0 + {n_{yy}}\bar \omega ' + {m_{yy}}\bar a',}\\ {\bar \omega ' = \frac{{\bar \omega 'R'}}{v} = \frac{{\xi R'}}{v} = \frac{{A{\omega _z}{\rm{cos}}\left( {{\omega _z}t} \right)R'}}{v},}\\ {\bar a' = \frac{{\dot {\omega '}{{R'}^2}}}{{{v^2}}} = \frac{{\xi {{R'}^2}}}{{{v^2}}} = \frac{{ - A\omega _z^2{\rm{sin}}\left( {{\omega _z}t} \right){{R'}^2}}}{{{v^2}}}.} \end{array}$

阻尼项和附加质量项是叶轮位置角θ的函数，将它们可以展开三角级数形式为：

$\begin{array}{*{20}{c}} {{n_{xx}} = \frac{{\delta {C_{FX}}}}{{\delta \bar \omega '}} = n_{xx}^0 + \sum\limits_{k = 1}^\infty {n_{xx}^k{\rm{sin}}\left( {kz\theta \left( t \right) + \psi _{\omega 'x}^k} \right)} ,}\\ {{n_{yy}} = \frac{{\delta {C_{FX}}}}{{\delta \bar \omega '}} = n_{yy}^0 + \sum\limits_{k = 1}^\infty {n_{yy}^k{\rm{sin}}\left( {kz\theta \left( t \right) + \psi _{\omega 'y}^k} \right)} ,}\\ {{m_{xx}} = \frac{{\delta {C_{FX}}}}{{\delta \bar a'}} = m_{xx}^0 + \sum\limits_{k = 1}^\infty {m_{xx}^k{\rm{sin}}\left( {kz\theta \left( t \right) + \psi _{a'x}^k} \right)} ,}\\ {{m_{yy}} = \frac{{\delta {C_{FY}}}}{{\delta \bar a'}} = m_{yy}^0 + \sum\limits_{k = 1}^\infty {m_{yy}^k{\rm{sin}}\left( {kz\theta \left( t \right) + \psi _{a'y}^k} \right)} ,} \end{array}$

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

$\begin{array}{l} {C_{FX}} = C_x^0 + {n_{xx}}\bar \omega ' + {m_{xx}}\bar a' = C_{xx}^0 + \sum\limits_{k = 1}^\infty {C_{xx}^k{\rm{sin}}\left( {kz\theta \left( t \right) + } \right.} \\ \quad \quad \;\;\left. {\psi _{xx}^k} \right) + \left\{ {n_{xx}^0 + \sum\limits_{k = 1}^\infty {\left[ {n_{xx}^k{\rm{sin}}\left( {kz\theta \left( t \right) + \psi _{\omega 'x}^k} \right)} \right]} } \right\}\bar \omega ' + \\ \quad \quad \;\left\{ {m_{xx}^0 + \sum\limits_{k = 1}^\infty {\left[ {m_{xx}^k{\rm{sin}}\left( {kz\theta \left( t \right) + \psi _{a'x}^k} \right)} \right]} } \right\}\bar a'. \end{array}$

(1)

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

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

(2)

利用CFD方法模拟立轴水轮机在波浪中的纵摇运动，可以得到C_FX和C_FY的时历曲线，然后用最小二乘法拟合可以得到均匀水动力项、阻尼项和附加质量项.

1.2 立轴水轮机模型

立轴水轮机的试验研究一般是在循环水槽中进行的，水轮机被安在实验平台上，水轮机在水流冲击下转动，并通过设定循环水槽中水流流速来模拟不同的海流情况，实验总布置图如图 1(a)所示.试验平台具体布置如图 1(b)所示，传动轴、齿轮、传感器、发电机和刹车电机等依次安装在实验平台上.传动轴与一伞形齿轮相连，伞形齿轮通过蜗杆连接到转矩转速传感器，传动轴上方还装有一个小型刷环引电器；平台下方有一法兰盘，用于与水轮机模型相连，通过该实验平台可以测量出水轮机的效率等.哈尔滨工程大学循环水槽实验室对立轴水轮机做了一系列实验研究，实验段尺度：长×宽×深为(8.0 m×1.7 m×1.5 m)，工作流速：0.2~2.0 m/s, 如图 1(c)所示.实验所采用的水轮机模型参数为：v=1 m/s，D=0.8 m，H=0.6 m，C=0.12 m，Z=2，叶片翼型为NACA0018，实验模型如图 1(d)所示.

除了实验方法来研究立轴水轮机的水动力性能外，CFD方法研究水轮机的最大优势在于可以获得更多的流场瞬时信息，计算机的快速发展也推动着计算流体力学的发展.本文采用CFD方法模拟相同参数的水轮机来研究其水动力性能，并和实验结果对比分析，证明CFD方法对立轴水轮机研究的可行性.如图 2为CFD计算网格模型，叶轮旋转中心距离纵摇中心0.3 m(0.5H)，本文模拟立轴水轮机运动情况需要同时模拟叶轮的旋转运动和纵摇运动(ξ=Asin(ω_Zt))，所以将计算网格分为旋转域、纵摇域和外域3部分，如图 2(b)、(c)、(d)所示，各个域之间的相对运动都是旋转运动，且3个域之间采用滑移网格形式.

整个计算域网格都采用结构化网格，边界条件设置如下：模型的入口边界为速度入口；出口边界为压力出口，相对压力为0；外域的侧面边界和叶片表面分别为自由滑移壁面和无滑移壁面条件；旋转域和纵摇域分别给定旋转角速度和纵摇简谐运动；计算选取的湍流模型为SST模型，选择瞬态求解器，时间步长为叶轮自转1°所使用的时间，时间步长的选取文献[5, 7]的研究方法.

2 结果与分析 2.1 CFD方法的可行性证明

在来流速度v=1 m/s，叶轮不发生纵摇的工况下比较同一个水轮机在实验方法和CFD方法得出的效率—速比对比如图 3所示.从图 3中可以看出速比2.5时，两种方法的水轮机效率都达到峰值，不同速比下的效率都比较接近，微小的差异可能是由于实验设备限制、自由液面等造成的，但这一结果对比证明了CFD方法在立轴水轮机水动力性能研究方面的可行性.

2.2 纵摇对立轴水轮机效率的影响

来流速度v=1 m/s，分别模拟均匀流中不发生摇荡的工况和波浪中发生纵摇运动.假设水轮机的纵摇频率在常见的波浪频率范围内，规定纵摇频率为1.2 rad/s.测量两种运动工况下不同速比时的水动力性能.

如图 4为纵摇工况下和非纵摇工况下的水轮机效率-速比对比图.从图 4中看出纵摇对水轮机效率曲线趋势没有大的影响，都是在速比2.5左右效率达到峰值.纵摇工况下的效率值略大于非纵摇的值，特别在峰值处更加明显.

从图 5(a)可以看出，非纵摇水轮机的尾涡几乎是对称分布的，图 5(b)、(c)是纵摇时两个不同状态下的水轮机尾涡场.

从图中可以看出水轮机尾涡受到纵摇的影响，和非纵摇尾涡场有明显差别.在纵摇情况下，叶片周围水流场会受到纵摇的影响，纵摇工况诱导的速度和自身旋转的速度相互影响使叶片周围流场复杂化，从而引起叶片受力增大，切向力的增大会引起叶轮效率提高.

2.3 纵摇对水轮机推力侧向力的影响

水轮机的推力和侧向力是影响水轮机结构强度和水动力性能的重要参数，通过对比在均匀流中不发生纵摇的立轴水轮机和波浪中发生纵摇情况下的水轮机的推力和侧向力曲线(λ=1.5)，为水轮机的设计提供参考.

如图 6(a)、(b)所示，在均匀流中的立轴水轮机的推力和侧向力是按照其自身旋转频率规则波动，两条时历曲线波峰和波谷的极值几乎不变.在波浪中发生纵摇运动的水轮机推力和侧向力的波动既受到自身旋转频率的影响，又受叶轮纵摇频率的影响.如图 6(c)、(d)所示，纵摇工况下的叶轮推力和侧向力的上下包络线明显有很大波动.

为了对比分析纵摇运动不同速比下的水轮机推力和侧向力波动规律，分别选取速比λ=1.0、1.5、2.0、2.5、3.0，比较它们的推力侧向力包络曲线，如图 7所示.

如图 7(a)、(b)为不同速比时的水轮机推力的上、下包络线时历曲线对比图.图 7(c)、(d)为不同速比时的水轮机侧向力的上、下包络线时历曲线对比图.当水轮机在均匀流中不发生纵摇运动时，水轮机的推力和侧向力的包络线都为直线形式.从4幅图中可以看出水轮机发生纵摇运动时，推力和侧向力的包络曲线明显发生波动，而且波动形式相似.随着速比增加，上包络线的值随着增大，而下包络线正好相反，随着速比增大而减小.推力和侧向力在纵摇情况下的波动对水轮机的结构强度和疲劳效应都有不利影响.

如图 8为立轴水轮机在均匀流的推力、侧向力和波浪中纵摇运动的推力、侧向力的波动幅值在不同速比下的对比曲线.图 8中可以看出各个力随速比的变化规律相同，都是随着速比增大波动幅值也随之增大，但纵摇运动时两个力的波动幅值比均匀流的明显要大.例如在速比3.0时，纵摇运动下的推力波动幅值比均匀流的增大18%，侧向力波动幅值则增大了将近30%，这是因为水轮机在波浪中发生纵摇运动时，叶片附近的水流受纵摇影响变得复杂，此处水流的合速度增大，从而增大了叶轮受力，这样不利于电能的稳定输出及水轮机的结构安全.

3 载荷系数拟合分析 3.1 纵摇运动下推力系数及侧向力系数的拟合

纵摇运动下的水轮机推力和侧向力的时历曲线可以通过CFD方法得到，利用最小二乘法根据式(1)、(2)可以拟合得到推力和侧向力的阻尼项和附加质量项，图 9为推力和侧向力的时历曲线和拟合值的对比图，从图 9中可以看出拟合值和被拟合的曲线高度吻合，可以对其进行对比分析.

对不同速比下推力和侧向力的时历曲线进行最小二乘法拟合，归一化后得出C_FX和C_FY级数形式的各项系数，见表 2、3，表中给出了C_FX和C_FY展开各项系数中的常数项和一次项，对其进行对比分析.

表 2中C_xx⁰、C_xx¹分别为推力系数的均匀水动力展开项，都与水轮机的纵摇速度和加速度无关，而且随着速比增大而增大，但这两项相比数值上差别不大，这说明在非纵摇运动时推力系数的平均值和波动幅值相差不大，波动幅值稍微大于平均值.ψ_xx¹、ψ_ω′x¹、ψ_a′x¹分别为纵摇运动下水轮机各项和叶轮位置角的滞后和超前，ψ_a′x¹的绝对值随着速比增大而增大，其他两项随速比变化没有明显变化.阻尼项n_xx⁰、n_xx¹随着速比的增大而增大，而且和附加质量项m_xx⁰、m_xx¹相比，数值上大一个量级左右.因此纵摇运动时推力的波动主要是阻尼力造成的，附加质量力影响较小可以忽略不计，而且随着速比增大，阻尼力的平均值和波动幅值随之增大.

表 3中C_yy⁰、C_yy¹分别为均匀水动力的展开系数与水轮机的纵摇速度和加速度无关，C_yy⁰绝对值小于C_yy¹，说明非纵摇时侧向力系数平均值小，而波动幅值较大；随着速比增加，C_yy⁰、C_yy¹的绝对值随之增大，这说明非纵摇时水轮机载荷的平均值(绝对值)随着速比增大而增大，而且其波动值也随着速比增大而增大.侧向力系数展开项中的阻尼力常数项n_yy⁰比m_yy⁰大，说明阻尼力是水轮机运动中的主要成分，即在纵摇运动时侧向力曲线的波动主要是阻尼力造成的，附加质量力占较小的影响.

3.2 水轮机计算结果的应用

本文研究了速比从1.0~3.0的强迫纵摇运动下水轮机推力和侧向力，并通过最小二乘法拟合得出了附加质量和阻尼项随速比的变化情况.因此通过插值计算，在任一速比下的纵摇水轮机的推力和侧向力的附加质量和阻尼项都可以得到.然后将得到的附加质量和阻尼系数应用到浮式平台上来分析不同速比下纵摇运动对浮式潮流电站的运动响应.

图 10为浮式垂直轴潮流电站的运动响应研究流程图，本文主要研究的是第1步垂直轴潮流能水轮机的水动力系数分析.通过CFD方法得到不同速比强迫纵摇运动下的水动力系数时历曲线，然后最小二乘法拟合分析得到附加质量和阻尼系数.

4 结论

1) 强迫纵摇运动的立轴水轮机和均匀流中水轮机的效率曲线变化规律相似，都是在速比2.5左右达到峰值，而且纵摇情况下效率值略大于均匀流中，在速比2.5时尤其明显.

2) 随着速比增加，推力和侧向力上包络线的值随之增加，而下包络线正好相反，随着速比增加而减小.

3) 立轴水轮机发生纵摇运动时，推力和侧向力的波动幅值明显比均匀流中的要大，这对叶轮的结构强度和疲劳寿命都会产生不利影响.

4) 立轴叶轮的推力和侧向力系数通过最小二乘法分析得到纵摇运动下的阻尼项和附加质量项，其中阻尼项比附加质量项大，所以阻尼力是水轮机纵摇运动中的主要成分，附加质量力占的影响较小.