水轮机转轮叶片是水轮机运行时的关键部件,叶片长时间受水流冲击,将引发结构的水弹性振动,可能会出现不同程度的裂纹破坏[1]。大量研究表明,激振力引发的结构共振所产生的动载荷作用是导致叶片裂纹的主要原因之一[2,3]。到目前为止,理论计算对于复杂水轮机周围的流场流动与水轮机自身的结构研究方面还不完善,只有实验研究与数值计算两种研究方式比较理想。数值计算利用CFD和CSD相结合的方式对水轮机进行流固耦合研究,是分析水轮机压力脉动的有效方法。潮流能水轮机的流固耦合问题主要指水轮机在水下运转时,其整体结构(特别是叶轮部分)会持续受到水流冲击作用,激振力导致水轮机结构产生变形,变形严重时会产生结构破坏,同时旋转的叶轮带动其周围流体运动,导致周围流体质量分布的改变,也会影响水轮机的正常工作。肖若富等[4]对水轮机转轮研究后发现,基于单、双向流固耦合方法分析下,两者最大应力点的位置与大小基本相同,应力与位移分布趋势大致相同。张立翔等[5]利用双向流固耦合方法在非定常流中对混流式水轮机叶片进行仿真计算,发现叶片的流激振动与实验结果吻合。金连根等[6]提出流体域、固体域在相同节点处双向流固耦合都大于单向流固耦合计算值。

本文以双转子悬浮式水平轴水轮机的叶轮为研究对象,基于单、双向流固耦合方法,对水平轴潮流能水轮机的叶轮实施三维数值模拟,以研究水轮机周围流场的非定常流动,得出在三种不同工况下水轮机叶轮的应变特征、各监测点的水压力脉动特性和频域振动特性。悬浮式潮流能水轮机运行时其叶轮上的压力脉动较为复杂,叶片主要受水流的激振力作用,非稳态的水流冲击力是激振力的主要来源,叶片在其作用下会发生振动,当振动的频率与叶片的固有频率相近时,会产生结构共振,诱发叶片裂纹破坏,影响水轮机的正常工作。

1、水轮机模型的建立

本文相关的数值仿真工作是基于本文作者提出的一种悬浮式水平轴潮流能水轮机,其结构分为三个部分:水轮机、锚泊系统与浮体[7,8]。如图1所示,该水轮机是双转子水平轴水轮机结构,两转子旋转方向不同,叶片数为5,采用NACA63215翼型,叶片采用Q235B结构钢[9,10],同时此材料与初期实验模型所选材料相同;上、下有一大一小两平衡翼结构,与导流罩一样也采用中空薄壁结构,这样不仅可以增强水轮机的抗倾覆能力,还可与导流罩一起为水轮机提供浮力,两结构提供的浮力可以使水轮机悬浮于水中。经模型试验与数值仿真发现,该水轮机设计方案在静水到设计最大流速,以至更高流速范围内均可以工作姿态悬浮于水中,由于其特殊的空心导流罩和上下水翼结构,可保证其在各种流速下均不发生倾覆,并可随海流方向自适应调整姿态;4根钢质锚系缆绳分别系于上、下水翼支架的前缘,另一端系于固定在海底的万向接头;Doherty等[11]对水轮机发电机组纵、横向布设优化方案的研究发现,水平轴水轮机横向间距约为2D时比间距更小的机组更具有优势。经本文作者前期的数值仿真研究发现,由于每个叶轮采用独立的导流罩设计方案,相邻叶轮之间的尾流干扰会更小。实际模型加工时为了保证足够的结构强度安全余量,两导流罩中心轴的间距采用1.8D,叶轮直径D=0.6m。

图1悬浮式水平轴潮流能水轮机结构示意图

本文以舟山海域实际流速作为参考,流速取1.0~3.0m/s进行研究[12],并通过定义水轮机的尖速比(λ=Dω/V)来确定流速与角速度的关系,水轮机的整体参数如表1所示。

表1水轮机整体参数

2、计算模型建立

2.1强度计算方法

叶轮结构静力学有限元方程为

公式1

式中:[K]为整体刚度矩阵;{u}为位移矢量;{F}为静力载荷。

叶轮结构动力学有限元方程为

公式2

式中:[Ms]、[Cs]、[Ks]分别为结构的质量矩阵、阻尼矩阵、刚度矩阵;{u}为节点的位移矢量;{F(t)}为激振力[13]。

该水轮机材料为Q235B结构钢,弹性模量为200GP,泊松比为0.3,Workbench中对于实体的默认单元类型是10节点四面体结构单元SOLID187。

2.2网格建立与边界条件的设定

选取悬浮式水平轴潮流能水轮机的单叶轮作为研究对象,流体为黏性不可压缩流体,流体域被划分为两部分,外部流体计算域的长、宽、高分别为10D(D为叶轮直径)、7D、5D,内部旋转流体域是直径为1.03D,长为0.33D的圆柱体区域。

如图2所示,网格划分时,外部流体域划分为结构网格,内部旋转流体域划分为非结构网格,并对叶轮近壁面区域进行网格加密,设置边界层,采用分层拉伸处理方法,对远离叶轮的流体区域的尺寸适当放宽,以减少网格数量。网格数的增加会使求解的误差逐渐缩小,直到误差消失[14]。y+值对水轮机计算结果的精度有很大的影响[15]。叶轮近壁面第一层网格高度取4×10-1mm,网格边界层数为6,y+此处取值范围为15.8~39.5,与文献[16]推荐值吻合,整个流体域的网格单元总数为296万个。

图2悬浮式水平轴潮流能水轮机叶轮的流体区域计算网格

2.3计算求解

本文使用ANSYS-CFX软件完成相关的数值仿真工作,采用基于压力求解的耦合算法,其优点是对网格质量依赖较低,具有良好的收敛性,而且在流场计算中,N-S方程、连续性方程两者可以一起耦合,实现动量守恒方程与质量守恒方程联立后的同时求解。瞬态模拟采用隐式时间积分方式,非定常迭代计算,采用二阶迎风差分格式,选用SST湍流模型进行方程封闭,SST模型能有效适应流体压力梯度的变化,而且可以利用壁函数精确模拟边界层现象,同时可以通过限制涡黏性系数来考虑湍流剪应力的传播,在近壁自由流中有更好的精度[17,18,19]。由尖速比确定叶轮的转速,采用MRF(多重参考坐标系法)定义叶轮模型的旋转速度和旋转域,MRF法可将流体的旋转域单独进行划分,并用与旋转域相同转速的相对坐标对旋转域进行定义[20]。外部流体域边界设置为无摩擦滑移壁面;叶轮表面设置为无滑移壁面[21,22],近壁面区域采用标准壁面函数来求解此区域内的低雷诺数流动,同时根据实际情况设定壁面的粗糙度。

流固耦合计算以稳态时的计算结果作为非定常流动计算的初始条件,计算精度精确到10-5。为能更加精确的了解水轮机叶轮附近的流场流动非定常信息,时间步长设置为5.83×10-4s(选取工况3叶轮转过1°时所用时间),计算总时长为16s。

水轮机叶轮的压力脉动频域分析采用快速傅里叶变换方法,选用矩形窗函数来截取分析数据。矩形窗函数在捕捉主频脉动频率,保持频率幅值不失真方面具有一定的优势[23]。

3、计算结果分析

基于单、双向流固耦合方法对悬浮式水平轴潮流能水轮机的叶轮进行迭代计算,耦合场收敛后得到流体域与固体域的计算结果。在三种工况下,水轮机的流速与转速各不相同,得到的计算结果也有差异。根据本项目之前的数值模拟及实验数据对比,可知此次仿真模拟的准确性。

3.1应变结果分析

由仿真计算结果可知三种工况下,叶片最大应变处都在叶尖部位,这是因为叶片尖端的弦长和厚度都为整个叶片的最小值,其在激振力作用下极易产生变形。以工况3为例,如图3所示,采用单、双向流固耦合计算方法,叶片最大变形量分别为5.389×10-1mm,6.949×10-1mm;叶片总长为202mm,叶片长度的5%为10mm,即最大允许变形量为10mm,对比后发现水轮机叶片满足叶片设计的变形要求。

图3工况3时叶轮的变形图

3.2应力结果分析

基于不同流固耦合方法计算时叶片叶根处应变分布图略有差异,但在三种工况下,叶片最大应力都位于叶片中部靠近叶根处。以工况3为例,如图4所示,采用单、双向流固耦合计算方法,叶片最大应力值分别为28.66MPa、29.09MPa,远低于Q235B钢的屈服强度σs=235MPa,属塑性材料,取安全系数为2.0,叶片的最大应力值也远小于叶片材料的许用应力127.5MPa,故叶片结构强度满足设计要求。

图4工况3时叶轮的应力云图

图5为根据工况3下叶轮最大应力处的动应力变化情况绘制的动应力脉动时域图,取样时间为11~16s,旋转周期为209.8ms,此时间段的开始周期为动应力波动逐渐平稳后的第52个周期,一周期内出现5次较为明显的波动,因所选监测点只在其中单一叶片,所以每次波动的峰值不尽相同,脉动主频与叶频相同,均为23.83Hz[24]。

图5叶轮最大应力处动应力脉动时域变化曲线(工况3)

表2、表3分别为基于单、双向流固耦合计算的叶轮应变量与应力大小的对比表。Do、Dt分别为单、双向流固耦合条件下的叶轮应变量,σo、σt分别为单、双向流固耦合条件下的叶轮应力值。水轮机叶轮在结构强度要求范围内,随着水流速度的不断增加,叶轮叶片的应变量不断加大,但增幅开始降低;与此同时,流速增加,作用在叶轮上的有效应力不断增大,但其增幅与应变量一样也开始降低。

表2基于单、双向流固耦合计算的叶轮应变量比较

表3基于单、双向流固耦合的叶轮应力大小比较

通过对比基于单、双向流固耦合方法的结构强度计算结果可以发现,相同工况,相同监测点下,与单向流固耦合的计算结果相比,双固耦合计算得到的叶片的变形量与应力值均略大,随着流速增加,两者的差值逐渐变小。考虑到双向流固耦合计算时固体结构与流体之间的相互作用及基于双向流固耦合方法下叶片的变形量与应力值比单向流固耦合计算时都略大,故从结构强度校核角度考虑,工程应用中使用双向流固耦合的计算结果更保守、更安全一些。金连根等也得出流体域、固体域在相同节点处双向流固耦合都大于单向流固耦合计算值的相似结论。

3.3压力脉动分析

为了度量水轮机的压力脉动,特引入压力脉动参数[25]

公式3

式中:Δp为压力与其平均值之差;ρ为流体密度;u为入流速度。

表4为不同工况下水轮机的叶频,叶片在振动过程中其表面不同位置的压力不断变化,为监测这种变化规律,先后在叶片的正面与背面分别设置了4个监测点,记录压力脉动情况(见图6)。bp-1、bs-1如图所示分别在叶片的尖端中点,bp-2、bs-2在整个叶片的1/3处,bp-3、bs-3在1/2处,bp-4、bs-4分别为正、背面应力最大处。

表4不同工况下水轮机的叶频

图6叶片正、背面压力监测点分布

3.3.1叶片正面的压力脉动分析

采用上述流场计算方法,在三种工况下分别进行了16s的(工况1、2、3分别完成了25、50、76个周期计算)仿真计算,并提取计算收敛后的4个稳定周期进行分析。

由图7~图12可知压力脉动系数随入流速度的增大而减小。由压力脉动系数的定义可知,入流速度对压力脉动系数影响的权重较大,水流激振力增大的压力值不足以抵消速度对压力脉动系数的影响,所以随着入流速度的增大会出现压力脉动系数减小的情况。在相同工况下在bp-1点处的压力脉动系数最大值相对于其他三点偏小,其原因主要有以下三点:(1)bp-1点选在叶片的边缘,此处提取压力数值时可能会产生误差;(2)叶片边缘的变形较大,使bp-1点周围的流场变得十分复杂,湍流、脱流、回流、二次流等复杂流动都可能出现,这些现象与叶轮之间的相互作用都会影响叶片尖端的压力脉动;(3)叶轮在转动过程中由于底部叶片对水流的阻挡,会使作用在叶片上的部分水流流向中的垂直分量转化为切向分量,从而影响位于叶片尖端bp-1点的压力脉动。

图7工况1下bp-1~bp-4点的压力脉动时域

图8工况1下bp-1~bp-4点的压力脉动频域

图9工况2下的bp-1~bp-4点的压力脉动时域

图10工况2下bp-1~bp-4点的压力脉动频域

3.3.2叶片背面的压力脉动分析

由图13~图18知,与叶片正面的压力脉动系数的变化趋势相同,在相同工况下在bs-1点处的压力脉动系数相对于其他三点偏小,其原因主要有以下三点:(1)bs-1点选在叶片的边缘,此处提取压力数值时可能会产生误差;(2)叶片边缘的变形较大,使bs-1点周围的流场变得十分复杂,湍流、脱流、回流、二次流等复杂流动都可能出现,这些现象与叶轮之间的相互作用都会影响叶片尖端的压力脉动;(3)叶轮尖端变形较大,改变叶片背面水流的流向,从而影响bs-1点的压力脉动。

图11工况3下bp-1~bp-4点的压力脉动时域

图12工况3下bp-1~bp-4点的压力脉动频域

图13工况1下bs-1~bs-4点的压力脉动时域

由表5可知,在三个工况下,监测点bp-2、bp-3、bp-4的压力脉动系数,其脉动系数幅值依次减小,以工况1为例,bp-2点压力脉动系数最大幅值是bp-3的1.24倍,bp-3是bp-4的1.20倍。监测点bs-2、bs-3、bs-4的压力脉动系数,其脉动系数幅值依次减小,以工况1为例,bs-2点压力脉动系数最大幅值是bs-3的1.25倍,bs-3是bs-4的1.16倍。

图14工况1下bs-1~bs-4点的压力脉动频域

图15工况2下bs-1~bs-4点的压力脉动时域

图16工况2下bs-1~bs-4点的压力脉动频域

图17工况3下bs-1~bs-4点的压力脉动频域

图18工况3下bs-1~bs-4点的压力脉动频域

正如3.3.1所述,由于压力脉动系数随流速增加而降低,所以工况1、工况2时压力脉动变化明显。由图19知,工况1、工况2中,从监测点1~4中能明显看出水轮机叶片背面的压力脉动比叶片正面变化幅度大,工况3下压力脉动系数相对较小,变化幅度没有工况1、工况2明显。对于工况1、工况2下水轮机叶片背面的压力脉动与正面相比变化幅度要大的原因主要是因为水轮机背面更容易产生涡流,涡流会加剧其波动,这一点可以从图20中的流线变化情况中可明显看出。

表5不同工况下各监测点的压力脉动系数的最大值

由表6可知,虽然在三种工况下水轮机的叶频随转频的变化而变化,但水轮机叶片正、背面的压力脉动的主频率集中在相应工况的叶频附近。

图19各监测点不同工况压力脉动系数最大值的比较

图20水轮机叶片周围的流线分布

表6不同工况下各监测点压力脉主频与相应幅值对比表

4、结论

本文通过分析潮流能水轮机周围流场的非定常流动,得到了在三种不同工况下叶轮的应变特征、各监测点的水压力脉动和频域振动特性,得出如下结论:

(1)与单向流固耦合结果相比,双向流固耦合计算得到的叶片的最大变形量与最大应力值均略大,随着水流速度的增大,叶轮的叶片的应变量与有效应力增幅逐渐降低。

(2)水轮机的压力脉动系数的幅值由叶尖到叶根基本呈递减趋势。

(3)水轮机的压力脉动的主频率集中在相应工况的叶频附近,且水轮机叶片背面的压力脉动比正面的波动要剧烈。

参考文献：

[1]王向红,朱昌明,毛汉领,等.基于核主成分分析及支持向量机的水轮机叶片裂纹源定位[J].振动与冲击,2010,29(11):226-229.

[2]丛国辉,王福军.双吸离心泵隔舌区压力脉动特性分析[J].农业机械学报,2008,39(6):60-63.

[3]朱文龙,周建中,夏鑫,等.基于水电机组运行工况的水轮机压力脉动诊断策略[J].振动与冲击,2015,34(8):26-30.

[4]肖若富,朱文若,杨魏,等.基于双向流固耦合水轮机转轮应力特性分析[J].排灌机械工程学报,2013,31(10):862-866.

[5]张立翔,王文全,姚激.混流式水轮机转轮叶片流激振动分析[J].工程力学,2007,24(8):143-150.

[6]金连根,毛建生,方兵.混流式水轮机转轮流场单向､双向流固耦合数值的分析比较研究[J].机电工程,2014,31(12):1564-1568.

[7]陈正寿,刘羽,赵陈,等.水平轴潮流能水轮机尾流场数值模拟[J].水力发电学报,2015,34(10):130-137.

[8]陈正寿,张国辉,刘羽,等.轴流式水轮机导流罩与叶轮尾流流场仿真研究[J].水动力学研究与进展,2016(1):90-98.

[10]张亮,何环宇,张学伟,等.垂直轴水轮机单向流固耦合数值研究[J].华中科技大学学报(自然科学版),2014(5):80-84.

[12]侯放,于华明,鲍献文.舟山群岛海域潮流能数值估计与分析[J].太阳能学报,2014,35(1):125-133.

[13]叶正寅,张伟伟,史爱明,等.流固耦合力学基础及其应用[M].哈尔滨:哈尔滨工业大学出版社,2010:1-26.

[16]纪兵兵.ANSYSICEMCFD网格划分技术实例详解[M].北京:中国水利水电出版社,2012.

[17]吴军,谷正气,钟志华.SST湍流模型在汽车绕流仿真中的应用[J].汽车工程,2003,25(4):326-329.

[19]荆丰梅,肖钢,熊志民.潮流能水轮机单向流固耦合计算方法[J].振动与冲击,2013,32(8):91-95.

[20]王企鲲.导流型垂直轴风力机内部流场数值模拟中若干问题的探讨[J].机械工程学报,2011,47(4):147-154.

[21]王松林,谭磊,王玉川.离心泵瞬态空化流动及压力脉动特性[J].振动与冲击,2013,32(22):168-172.

[23]姚志峰,王福军,肖若富,等.离心泵压力脉动测试关键问题分析[J].排灌机械工程学报,2010,28(3):219-223.

[24]谭磊,王玉川,曹树良,等.离心泵蜗舌区非定常流动特性[J].北京理工大学学报,2014,34(7):670-675.

[25]王洋,代翠.离心泵内部不稳定流场压力脉动特性分析[J].农业机械学报,2010,41(3):91-95.

陈正寿,张国辉,赵宗文,程枳宁,郑武.潮流能水轮机叶轮压力脉动特性分析[J].振动与冲击,2017,36(19):98-105+138.

基金:国家自然科学基金(41476078);浙江省公益技术应用研究计划项目(2015C34013);大连理工大学海岸和近海工程国家重点实验室开放课题基金资助项目(LP1703)