基于重叠网格方法的水轮机非定常水动力数值仿真

刘垚，蔡卫军，王明洲

西安精密机械研究所，陕西西安710077

0 引 言

随着我国海洋空间站、深海预置武器等各类海洋探测仪器设备的快速发展应用，其待机状态下长期的能源补给成为一个急需解决的问题。海流能作为海洋中一种分布广泛、有较强规律性和可预测性、能量大且稳定的可再生能源，可以为海洋设备提供充足的能源补给，因此海流发电技术对海洋工程的发展具有重要意义。我国大部分海域海流的流速普遍较低，海洋装备所处的海洋环境较为复杂，周围海流流速不断变化，对海流能水轮机在低速流下的各项水动力性能要求较高。因此，如何在低速、不稳定流条件下捕获海流能进行高效发电成为亟待解决的问题。

在研究海流能水轮机的水动力性能过程中，为了得到更精确、详细的水轮机性能参数和更详细的流场信息，在试验之前，利用计算流体动力学（Computational Fluid Dynamics，CFD）仿真方法进行数值求解是最有效、便捷的方法之一。目前，国内外对水轮机的仿真研究聚焦于利用滑移网格技术实现水轮机的定转速旋转过程，通过调节水轮机的额定转速改变叶片减速比，在假设来流流速不变的条件下，对水轮机进行离散化的水动力性能研究。Batten 等［1-2］和 Bahaj等［3］结合获能原理和试验分析，对水平轴海流能水轮机的水动力特性进行了系统研究；赵陈等［4］对海流能水轮机进行优化设计，在不同额定转速下对水轮机的水动力特性进行了离散化仿真，在最大输出功率条件下分析了水轮机的最佳转速；陈存福［5］运用多参考系模型（MRF模型）研究了来流流速、水轮机转速和叶片安装角对海流能水轮机功率的影响。

在对水轮机非定常水动力性能仿真研究中，常规的滑移网格方法需要对网格进行局部调整和再生成，这给研究带来了极大的不便。为此，本文拟研究利用重叠网格来实现水轮机在流场中的动态旋转过程，结合六自由度动态流体固态相互作用（DFBI）运动实现水轮机在来流作用下的被动旋转。

1 数学模型

1.1 控制方程

CFD仿真软件STAR-CCM+在数值仿真计算中的理论控制方程是质量守恒方程、动量守恒方程和能量守恒方程［6］，微分方程的守恒通式为

式中：ρ为流体密度；φ为通变量；w为待求解的速度；Γ为广义扩散系数；u为速度矢量；S为广义源项；v为运动粘性系数。

在海流能水轮机叶片旋转运动过程中，由于水为不可压缩流体，周围流场环境为常温常压状态，故可假设：1）流体为不可压缩流；2）粘性系数为常数；3）忽略体积力；4）不涉及求解能量守恒方程。流体的连续方程为

动量方程为

式中，p为压力项。

1.2 基于重叠网格技术的DFBI方法

目前，对水轮机水动力性能的CFD仿真研究主要是运用滑移网格方法，滑移网格是简化的计算模型，在一个独立区域内所有网格一起运动，插值计算只在滑移交界面上进行。与滑移网格相比，重叠网格是将复杂的流动区域分成多个几何边界比较简单的子区域，各子区域中的计算网格独立生成，彼此存在着重叠、嵌套或覆盖的关系，流场信息通过插值在重叠区域进行匹配和耦合［7-9］，图1所示为重叠网格示意图。对网格2固体边界进行挖洞处理，不参与计算，插值点是人为设置的边界，与网格1重叠部分的边界会相互传递计算流场信息，完成流场信息之间的链接。

在处理物体非定常运动数值仿真时，重叠网格可以解除物体与网格之间的拓扑结构约束，物体在计算域内自由运动，计算网格无需重新生成。六自由度运动模型DFBI可以模拟刚体在流体作用下的被动运动，刚体包括6个自由度，可以根据需要确定刚体的自由度。将重叠网格方法与DFBI结合，可以实现在流体作用下刚体运动特性的仿真模拟，对刚体运动过程的仿真更加接近真实运动情况。

在海流能水轮机的非定常水动力性能数值仿真中，使用基于重叠网格技术的DFBI方法不仅可以模拟水轮机叶片匀速旋转，还可以针对复杂海流环境条件下水轮机叶片被动旋转的启动过程和动态旋转进行仿真研究，精确监测水轮机因周围流场变化导致的加速和自适应转动等过程，实现对水轮机非定常水动力特性的研究。

1.3 计算模型验证

Batten 等［1-2］和 Bahaj等［3］研发了海流能水轮机，并进行了多次水槽试验。受实验条件的限制，在进行仿真计算前，为了验证本文选用仿真方法的准确性，利用本文的数值计算方法对Bahaj等的水轮机模型进行了仿真计算，本文设置的边界条件与Bahaj等的相似。在给定水轮机转速的条件下，得到了水轮机的获能系数CP与水轮机尖速比（TSR）之间的关系。将仿真结果与Bahaj等的水槽实验结果进行对比，结果如图2所示。图2中的黑色曲线为计算值的拟合曲线。由图2可以看出，水轮机获能系数的整体增长趋势基本相同，计算值与实验值的差值在15%以内，仿真结果与实验结果较为吻合，因此可以验证本文数值仿真方法的有效性。

2 计算域、网格及边界条件

2.1 海流能水轮机模型建立

本文选取NACA 63815翼型作为水轮机叶片翼型，水轮机平均半径为400 mm，用R表示。各点的水轮机半径为r，叶片数为3，轮毂直径为160 mm，叶片的安装角为25°，叶片弦长L和扭角β沿径向的分布情况如图3和图4所示，水轮机的三维模型如图5所示。叶片材料为铝质材料，查得水轮机质量为80 kg，绕x轴的旋转转动惯量为0.32 kg/m2。

2.2 计算区域及网格划分

在建立水轮机几何模型的基础上，对海流能水轮机计算区域进行划分。计算区域由旋转区域和流场区域组成，如图6所示。旋转区域为水轮机在流场内旋转的区域，将旋转区域设置为重叠网格区域以实现水轮机的真实旋转。流场区域为背景域，为了让来流充分发展，避免外流场域边界对水轮机水动力性能仿真计算产生影响，取直径为水轮机直径6倍的圆柱体作为外流场域边界。

利用CFD仿真软件的网格划分功能对水轮机计算域进行三维网格划分，网格划分截面如图7所示，模型的整体网格规模为1 220 000。旋转区域设置为重叠网格区域，该区域采用切割体网格和棱柱层网格，为了更加精确地监测叶片表面的各类参数，对叶片表面设置面控制进行加密处理，叶片表面棱柱层数为5层，棱柱层增长率为1.2。外流场域采用切割体网格，并在旋转区域周围设置网格加密过渡区域。为了保证加密区域和旋转区域之间重叠网格边界计算数据的准确传递，设加密区域与重叠网格接触面附近区域的网格大小一致。

2.3 边界条件

整个计算域的边界包含进口边界、出口边界、远场、叶片、旋转域与外流场交界面，这些边界的设定和处理方法如表1所示。

结合水轮机周围的海流环境，针对水轮机数值模拟文献中的模型［10-11］，通过理论分析和数值模拟的方法比较不同湍流模型，验证了标准SSTk-ω模型可以更高效地模拟海流能水轮机的水动力性能。因此本文选取标准SSTk-ω模型作为湍流模型。采用标准壁面函数法解决近壁面流体低雷诺数问题，使用SIMPLE算法对压力速度耦合方程进行求解，对流项采用二阶迎风格式，离散方程粘性项采用二阶中心差分格式。在来流作用下，水轮机随流被动旋转运动的数值模拟采用非稳态计算，在DFBI中设置六自由度体为水轮机，水轮机沿着x轴方向做被动旋转运动，设定水轮机的运动形式为单自由度旋转运动，设置为绕x轴旋转。

表1 边界条件设定Table 1 Boundary conditions

3 仿真计算结果与分析

3.1 水轮机特性参数定义

水轮机的获能功率P的表达式为

式中：J为水轮机绕轴旋转时的转动惯量；ω为水轮机转速。

水轮机在海流中获取的能量效率用获能系数Cp表示，Cp是评价水轮机水动力性能的重要参数，其表达式为

式中：V为来流流速；ρ为海水密度。

3.2 水轮机的启动过程分析

设计流速为1.2 m/s，对水轮机的启动过程进行非定常数值模拟，图8给出了水轮机在启动过程中角速度和转动力矩随时间的变化情况。由图可知，水轮机的启动过程可以分为3个阶段：起步、加速和稳定阶段。

在0～0.7 s时，在来流的冲击作用下，水轮机叶片开始旋转，启动力矩快速增长，叶片初始转速较慢；0.7～2 s为水轮机的加速阶段，水轮机力矩增加到峰值8.4 N·m后逐渐下降，转速则急剧增加，快速接近稳定值；2～7 s为水轮机稳定过程。由于本研究没有给水轮机附加阻尼力矩，最终水轮机受到的力矩波动逐渐变缓最终趋于0。

图9 为水轮机在0.5，1.5，4和7 s时的速度云图，图10为水轮机在0.5，1.5，4和7 s时叶片后端0.1 m处的压力云图。由图9和图10可以看出：在启动阶段，来流在水轮机的阻挡作用下会产生一个直径为1.6 m的低速尾流区域，结合叶片后端压力云图可知，在启动阶段该区域的压力梯度变化较大，流场紊乱，湍流强烈；在加速阶段，水轮机下游的尾流区域后移变长，下游流速影响范围扩大，在水轮机尾部产生了尾流涡，会影响水轮机的稳定性，叶片后端的压力梯度逐渐变缓，叶片后部流场开始稳定；在稳定阶段，水轮机尾流区域长度基本稳定，尾流区的湍流涡逐渐消失，整体流场趋于稳定。

3.3 水轮机在定常流速条件下的水动力性能分析

在来流流速0.5～2.0 m/s条件下，对水轮机被动旋转进行水动力性能仿真分析。图11给出了叶片在0.8，1.4和2.0 m/s流速条件下的涡量云图。在1.4 m/s之前，叶片表面产生的湍流涡较少，叶片周围流场较稳定；在2.0 m/s时叶片表面的流动工况明显恶化，边界层受到破坏，叶尖出现较大的涡流区，叶背尖端形成了大量湍流涡，导致叶片发生失速，叶片获能效率降低。

图12给出了水轮机的输出功率和获能系数随来流流速的变化情况。由图可知，在零阻尼的条件下，随着来流流速的增加，水轮机的输出功率一直呈增长趋势，在设计流速1.2 m/s之前水轮机的获能系数Cp稳定在0.23，在设计流速点附近效率最高达到0.24。随着流速继续增加，由于水轮机叶片发生失速，水轮机获能系数Cp开始下降，水轮机获能效率降低。

3.4 水轮机在非定常流速条件下的水动力性能分析

在真实的海洋环境中，海流流速是在一定范围内不断变换的。为了研究来流流速的不断变化对水轮机模型动态水动力性能的影响，在前文研究的基础上，设置计算域速度进口条件，利用周期T=20 s，流速为0.8～2 m/s范围的正弦函数模拟海流流速的变化规律，对水轮机水动力性能进行动态仿真分析研究。

为了研究随着来流流速的变化水轮机叶片表面压力和水轮机转速的变化情况，对水轮机在两个流速变化周期内水轮机转速和获能系数随时间的动态变化情况进行了监测，结果分别如图13和图14所示。图15为对水轮机周围不同时刻的速度流场进行监测的结果。

由图13和图14可见，水轮机转速和获能系数受流速变化范围影响较大，转速在140～430 r/min内周期性波动，获能系数Cp也在0.09～0.32间波动。

由图15可以看出，在T/4时，来流流速为0.8 m/s，来流流速较低，叶片周围流场速度梯度较小，叶尖涡基本不存在，在此流速附近水轮机的获能效率达到0.29；从T/4到3T/4，流速不断增加，水轮机周围流场速度梯度明显增加，叶尖涡范围不断扩大，水轮机在流速增加阶段获能系数快速下降，在较低值范围内波动；3T/4之后流速开始下降，叶尖涡逐渐减小，水轮机获能系数也开始缓慢上升。

在变流速条件下，计算水轮机的平均获能系数约为0.181，比设计流速1.2 m/s时的获能系数下降了约33%，由此可知水轮机获能效率受流速变化影响较大，有必要针对流速变化对水轮机的影响进行优化设计，以提高水轮机在动态海流环境中的获能效率。

4 结 论

本文利用基于重叠网格技术的DFBI方法对水轮机模型进行了非定常仿真研究，实现了对水轮机被动旋转启动过程和不同流速条件下水动力性能和流场特征的非定常仿真分析。研究发现：

1）应用基于重叠网格技术的DFBI方法，可以实现仿真模拟水轮机在来流作用下被动旋转的非定常过程，对实际海流情况下的水轮机被动旋转非定常水动力特性进行监测，这对实际工程设计具有较好的参考价值。

2）水轮机在来流作用下的启动过程会经历启动、加速和稳定阶段，水轮机受到的力矩峰值集中在启动阶段，此时水轮机周围流场变化紊乱，压力梯度较大；加速阶段水轮机转速迅速上升，并且在水轮机尾部产生尾流涡，影响水轮机的稳定性；在稳定阶段水轮机受到的力矩波动较小并逐渐趋于0。

3）水轮机在来流流速0.5～2.0 m/s作用下被动旋转，叶轮获能系数Cp在设计流速1.2 m/s附近达到峰值0.24，但在设计流速1.2 m/s后，随着流速的增加，叶尖流动分离和失速导致水轮机获能效率快速下降，因此水轮机叶片的设计流速应尽量符合水轮机工作环境的流速，避免因叶片失速而导致效率降低。

4）海流流速的不断变化对水轮机周围流场和水动力性能的影响较大。在周期性波动的变流速条件下，水轮机的平均获能系数约为0.181，比设计流速1.2 m/s时的获能系数下降了约33%。有必要针对流速变化对水轮机的影响进行优化设计，以提高水轮机在动态海流环境中的获能效率。

［1］BATTEN W M J，BAHAJ A S，MOLLAND A F，et al.Hydrodynamics of marine current turbines［J］.Renewable Energy，2006，31（2）：249-256.

［2］BATTEN W M J，BAHAJ A S，MOLLAND A F，et al.Experimentally validated numerical method for the hydrodynamic design of horizontal axis tidal turbines［J］.Ocean Engineering，2007，34（7）：1013-1020.

［3］BAHAJ A S，MOLLAND A F，CHAPLIN J R，et al.Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank［J］.Renewable Energy，2007，32（3）：407-426.

［4］赵陈，刘羽，陈正寿，等.潮流能水轮机叶片设计与计算［J］.浙江海洋学院学报（自然科学版），2015，34（1）：59-67.ZHAO C，LIU Y，CHEN Z S，et al.Design and computation about blades of tidal current turbine［J］.Journal of Zhejiang Ocean University（Natural Science），2015，34（1）：59-67（in Chinese）.

［5］陈存福.潮流能水平轴水轮机叶片优化及水动力性能研究［D］.青岛：中国海洋大学，2012.CHEN C F.Study on blade optimization and hydrodynamic performance of horizontal axis marine current turbine［D］.Qingdao：Ocean University of China，2012（in Chinese）.

［6］MENTER F R.Two-equation eddy-viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32（8）：1598-1605.

［7］CHAN W M，GOMEZ R J，ROGERS S E，et al.Best practices in overset grid generation［C］//Proceedings of the 32nd AIAA Fluid Dynamics Conference and Exhibition.Missouri：AIAA，2002.

［8］NOACK R W.A general capability for moving body overset gridassembly［C］//Proceedingsof the17th AIAA Computational Fluid Dynamics Conference.Toronto：AIAA，2005.

［9］CARRICA P M，WILSON R V，NOACK R W，et al.Ship motions using single-phase level set with dynamic overset grids［J］.Computers&Fluids，2007，36（9）：1415-1433.

［10］王福军.流体机械旋转湍流计算模型研究进展［J］.农业机械学报，2016，47（2）：1-14.WANG F J.Research progress of computational model for rotating turbulent flow in fluid machinery［J］.Transactions of the Chinese Society for Agricultural Machinery，2016，47（2）：1-14（in Chinese）.

［11］王树杰，盛传明，袁鹏，等.潮流能水平轴水轮机湍流模型研究初探［J］.中国海洋大学学报，2014，44（5）：95-100，113.WANG S J，SHENG C M，YUAN P，et al.A Study onturbulence models of horizontal axis tidal current turbines［J］.Periodical of Ocean University of China，2014，44（5）：95-100，113（in Chinese）.