随着经济的发展, 人类对能源的需求量不断增长, 因此各国大力发展绿色可再生能源开发利用技术。潮流能^[1-2]作为海洋可再生能源, 因具有较好的可预测性、低流速发电利用率高、同时可提供海岸防护等特点日益受到人们的重视。潮流能水轮机^[3-4]是将潮流能动能转换为机械能的关键部件, 转子水动力载荷和轮机性能直接影响到潮流能发电站的可靠性和经济性。潮流水轮机水动力载荷研究是电站整体设计的基础, 与基于动量定理的流管法^[5-6]和基于旋涡理论的涡方法^[7-9]相比, 基于求解纳维-斯托克斯方程的计算流体力学方法^[10-13]因具有较为稳健的计算能力同时可以提供详细的流场信息日益得到广泛应用。

潮流能电站载体形式主要有坐底式和漂浮式, 漂浮式载体因安装便利、部署灵活、可回收性好等特点逐渐成为新增潮流能电站载体的主流。实际漂浮式潮流电站^[14]在工作环境中必然受到波浪的影响。王凯等^[15-16]。盛其虎等^[17]分别研究了水轮机艏摇频率对水动力系数的影响和纵荡运动对水载荷的影响。本文采用计算流体力学软件ANSYS-CFX对漂浮式竖轴水轮机进行横荡和艏摇耦合运动时的水动力载荷进行研究, 得到固定速比下同频和不同频横荡和艏摇耦合运动对水轮机水动力性能及载荷的影响规律。

1 计算流体力学数值仿真 1.1 模型建立

竖轴水轮机计算模型采用的叶片是直叶片, 由于叶片展弦比较大, 因此可采用二维模型模拟水轮机在均匀来流下耦合运动的水动力载荷。假设数值模拟中不考虑支撑结构对轮机的影响, 载体在波浪作用下艏摇和横荡运动均为简谐运动, 假设横荡运动Z_S和艏摇运动θ_Y分别为：

$ {Z_S} = {\xi _S}\sin \left( {{\omega _S}t} \right) $

(1)

$ {\theta _Y} = {A_Y}\sin \left( {{\omega _Y}t} \right) $

(2)

式中：ξ_S为横荡幅值；ω_S为横荡频率；A_Y为艏摇幅值；ω_Y为艏摇频率。计算模型参数为轮机直径D为0.8 m、叶片数N为2、叶片弦长C为0.12 m, 叶片翼型为NACA0018, 为了今后与模型实验数据对比, 计算来流速度为1 m/s。根据轮机耦合运动的组成方式与特点将尺度为28D×16D的计算域分为旋转域、艏摇域、横荡域、变形域和静止域, 速度入口距水轮机中心为8D, 各区域之间通过交界面相接, 计算域具体划分如图 1所示。在艏摇域中心建立大地坐标系(O-XY), 在旋转域中心建立叶片参考坐标系(o-xy), 如图 2所示, ξ_Y为艏摇角度, ω_T为叶轮旋转角速度。

1.2 控制方程和湍流模型

不可压缩流体雷诺时均方程(Reynolds-averaged navier-stokes equations, RANS)可写成笛卡尔哑指标求和形式：

$ \frac{{\partial {{\bar u}_i}}}{{\partial {x_i}}} = 0 $

(3)

$ \rho \frac{\partial \bar{u}_{i}}{\partial t}+\rho \bar{u}_{j} \frac{\partial \bar{u}_{i}}{\partial x_{j}}=\bar{f}_{i}-\frac{\partial \bar{p}_{i}}{\partial x_{i}}+\frac{\partial}{\partial x_{j}}\left(\mu \frac{\partial \bar{u}_{i}}{\partial x_{j}}-\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}\right) $

(4)

式中：ρ为流体密度；u_i、u_j为i、j方向的时间平均流动；p_i为压强；${ - \rho \overline {u_i^\prime u_j^\prime } }$为雷诺应力项。控制方程需要进一步补充湍流模型使之封闭, 由于SST k-ω模型可以较好模拟动态失速特性, 因此在计算中采用该湍流模型。

网格生成作为计算流体力学中的首要部分占据大量的人力处理时间。同时对数值模拟的精度和计算效率起着至关重要的作用。本文中整个计算域的网格划分全部采用结构化网格, 根据文献[18]中对网格尺寸的研究, 叶片表面网格第1层高度0.000 1 m, Y+值范围5.6~9.5, 网格总数接近10万, 图 3为旋转域和艏摇域的网格划分结果。

由于旋转域和艏摇域分别绕轴O旋转, 因此采用滑移网格对该域进行划分, 特点是不涉及网格的变形和重生；横荡域和变形域采用动网格划分, 用于模拟流场形状随时间改变, 只涉及Y向平动与拉伸。旋转、横荡和艏摇耦合运动的设置策略为：整体给定旋转域、横荡域及艏摇域横荡运动, 即整体沿竖水流方向来回平动；艏摇域和旋转域在发生横荡运动的同时, 又给定艏摇角速度；旋转域在发生横荡运动和艏摇运动的同时, 又整体给定绕自身主轴的旋转角速度；变形域与横荡域相邻, 由于横荡域的整体平动, 网格在Y方向上会产生拉伸或压缩, 但没有网格的扭曲；静止域网格恒定不变。

1.3 边界条件设置

设定大气压为参考气压；入口边界和上下边界设为速度入口, 给定入口来流速度U为1 m/s；出口边界采用压力出口, 设置出口相对压力为0；水轮机表面为无滑移壁面；计算域各分块之间通过交界面连接。模型中湍流模型采用SST k-ω模型, 求解器为隐式非定常, 计算时间步为水轮机旋转1°所需时间。

2 水轮机载荷分析

为了后续数据分析方便, 定义无量纲参数推力系数、侧向力系数、转矩系数及速比为：

$ {C_{Fx}} = \frac{{{F_X}}}{{0.5\rho {U^2}D}} $

(5)

$ C_{F y}=\frac{F_{Y}}{0.5 \rho U^{2} D} $

(6)

$ {C_{Mz}} = \frac{{{M_Z}}}{{0.5\rho {U^2}DR}} $

(7)

$ \lambda = \frac{{{\omega _T}R}}{U} $

(8)

式中：M_Z为水轮机绕旋转轴的力矩, Nm；ρ为流体密度, kg/m³；F_X为整体所受推力, N；F_Y为整体所受侧向力, N；ω_T为旋转角速度, rad/s；R为半径, m。

图 4~5为水轮机推力系数、侧向力系数和转矩系数在不同运行条件下的时历曲线。从图中可以看出, 在轮机仅以旋转频率ω_T=6.25 rad/s发生旋转运动时, 上述3个系数的时历曲线发生波动, 波动频率为12.56 rad/s, 约是旋转频率的2倍, 但波动幅值不发生变化, 其极值的上下包络线为直线。当水轮机发生旋转和其他运动(横荡、艏摇、横荡和艏摇)耦合时, 极值包络线也发生与运动频率相关的波动。因此, 水轮机的载荷波动由2部分组成：一部分是基于旋转频率的高频波动；另一部分是基于其他运动频率的低频运动。水轮机的载荷变化规律可通过极值的包络线来体现, 因此, 在后续的分析中只给出极值包络线的变化规律。

2.1 同频耦合运动

为了研究竖轴水轮机发生旋转及同频横荡和艏摇耦合运动时的水动力特性, 设置水轮机旋转在能量利用率最大的对应速比附近(λ=2.5), 艏摇幅值15°及横荡幅值0.2 m。选取3种不同频率进行计算, 并与水轮机只发生旋转和艏摇耦合运动及旋转和横荡耦合运动时的计算值进行比较。推力系数、侧向力系数及转矩系数计算值的对比曲线如图 6所示。其中的0.6、1.0及1.4表示水轮机发生旋转、艏摇及横荡耦合运动时的计算结果, S0.6、S1.0及S1.4表示水轮机只发生旋转和横荡耦合运动时的计算结果, Y0.6、Y1.0及Y1.4表示水轮机只发生旋转和艏摇耦合运动时的计算结果。

从图 6中可以看出, 当水轮机发生旋转、横荡和艏摇耦合运动时, 载荷系数包络线的变化规律与水轮机发生旋转和横荡耦合运动、旋转和艏摇耦合运动时的变化规律基本一致；随着耦合频率的增大载荷包络线(转矩系数下包络线除外)的波动幅值有了明显的增加。水轮机发生旋转和横荡耦合运动时, 载荷包络线(转矩系数下包络线除外)的波动幅值明显大于旋转和艏摇耦合运动时的波动幅值, 说明横荡运动比艏摇运动对水轮机载荷能产生更为激烈的影响, 同频耦合运动不利于水轮机电站的长期安全运行。

2.2 不同频耦合运动

取2种不同频率组合(艏摇频率0.6 rad/s、横荡频率1.4 rad/s；艏摇频率1.4 rad/s、横荡频率0.6 rad/s)进行计算, 并与水轮机只发生旋转和艏摇耦合运动及旋转和横荡耦合运动时的计算值进行比较。推力系数、侧向力系数及转矩系数计算值的对比曲线如图 7所示。

从图 7中可以看出, 在水轮机发生旋转、横荡和艏摇耦合运动过程中, 当横荡频率较大而艏摇频率较小时, 载荷包络线的幅值和变化规律基本上和水轮机发生旋转和与相同横荡运动频率耦合时的载荷包络线变化一致, 即横荡运动的影响大于艏摇运动的影响；当艏摇频率较大而横荡频率较小时, 载荷包络线的数值和变化规律明显与水轮机只发生旋转和相同频率的艏摇运动时的不同, 此时载荷包络线(转矩系数下包络线除外)的波动幅值明显增加。因此, 通过对比分析发现, 在水轮机发生旋转、横荡和艏摇耦合运动时, 横荡运动对水动力载荷的影响明显大于艏摇运动。所以, 针对横荡运动和艏摇运动而言, 因尽量减少浮式载体在波浪中的横荡运动的频率和幅值。

3 流场分析

为了进一步分析水轮机耦合运动下载荷波动的机理, 提取了第2节中同频运动和不同频运动下的尾流场的涡量分布, 选取时刻为流场稳定后一个周期, 用0、T/4、T/2、3/4T表示, 分别对应从平衡位置开始运动到Y向最大位置, 然后恢复到平衡位置最后到-Y向最大位置的过程, 如图 8为λ=2.5, ω_S=ω_Y=0 rad/s时仅旋转时涡量分布图, 图 9为ω_S=ω_Y=1.4 rad/s时同频耦合运动涡量分析, 图 10(a)为ω_S=1.4 rad/s, ω_Y=0.6 rad/s时不同频耦合运动涡量分析, 图 10(b)为ω_S=0.6 rad/s, ω_Y=1.4 rad/s时不同频耦合运动涡量分析。

从图 8可见, 水轮机仅旋转时尾流场涡量分布沿Y轴对称, 涡量大小随着旋涡扩散与耗散沿流向衰减。当水轮机旋转且与横荡或艏摇运动耦合时, 尾流场涡量分布呈现出非对称的特点, 水轮机上下叶片脱泄出的尾涡在耦合运动和涡系干扰下沿下游演化形成周期性交错非对称排列的卡门涡街, 这是载荷周期性波动的来源。综合比较图 7、10可得, 水轮机在同频横荡和艏摇时尾流场涡量分布非对称性最强, 体现在载荷包络线幅值最大。结合图 11可得, 横荡频率大的尾流场涡量分布非对称性比横荡频率小的要强, 对应上节对载荷的分析横荡的影响比艏摇大的结论。

4 结论

1) 水动力载荷系数随时间发生波动, 且波动由2部分组成：一部分是基于旋转频率的高频波动；另一部分是基于横荡和艏摇运动频率的低频运动。

2) 水轮机发生旋转、横荡和艏摇耦合运动过程中, 横荡运动对水动力载荷的影响大于艏摇运动。

3) 横荡和艏摇同频耦合运动时, 水动力载荷包络线波动频率与横荡频率一致, 波动幅值随横荡频率增加, 水轮机总体设计时应避免该工况发生。

4) 水轮机尾流涡量场分布在轮机只旋转时沿Y轴对称, 耦合运动导致尾流场涡周期性脱泄是载荷包络线波动来源, 其非对称性越强载荷幅值越大。