潮流能水轮机尾流场及涡特性DDES模拟

何 聪,史宏达,张 嶔,赵昌宇,翟仁彬

(中国海洋大学 工程学院,山东 青岛266100)

潮流能因可预测性高、储量丰富而被广泛研究。水平轴潮流能水轮机能很好地转换利用潮流能,但高昂的发电成本影响了潮流能的商业发展。采用多机组潮流能水轮机阵列能有效节省发电成本、增加供电量,通过对尾流结构研究可以确定相邻水轮机之间的最佳距离,从而获得更高的转换效率[1]。同时,水轮机尾流对海床冲刷、积淤等也有重要影响[2]。因此研究潮流能水轮机的尾流结构有重要意义。

在潮流能尾流场研究中,Gaurier等[3]利用实验对比了不同湍流下的水轮机性能,结果表明低流速下水轮机荷载会受内部湍流结构的影响。张玉全等[4]讨论了水轮机安装高度对尾流扩散的影响,提出增加水轮机与池底的距离可以加快尾流恢复速度,减小尾流长度。但尾流复杂的物理规律没有被完全揭示,随着科学技术进步,数值模拟因成本低、效率高逐渐应用于潮流能水轮机尾流的空间演变过程研究。Leroux等[5]使用RANS-SST模型对潮流能水轮机进行了数值模拟研究,结果表明稳态和瞬态的方法在潮流能水轮机性能模拟方面精确度相似,都能较好预测推力系数和功率系数,但是在对尾流的数值模拟中,瞬态模拟的结果更加符合实验结果。Afgan等[6]对比了基于RANS和LES模型对潮流能水轮机尾流场的数值模拟差异,结果表明LES模型可以捕捉到叶片尖端旋涡及解释支撑塔架对尾流的相互影响,而RANS未捕捉到这2种旋涡。相比RANS模型,LES模型对流场的描述更加准确详细。但需要昂贵的计算成本,且受限于计算设备。为了减少求解边界层湍流结构的计算量,DES模型被提出并后续优化为DDES模型。

本文基于DDES方法对潮流能水轮机尾流场进行数值模拟,通过与实验结果对比,验证数值模拟的计算精度和可靠性,叙述数值模拟的网格类型、模型参数。重点分析水轮机尾流场的典型物理量的分布特征,以及研究水轮机的涡结构发展过程。

1 几何模型与计算区域

潮流能水轮机几何模型采用Payne等[7]设计的三叶片水平轴水轮机模型(图1)。转盘半径(R)为0.6 m,叶片翼型为NACA 63-8XX。水轮机尾舱由两部分组成,直径(D)由0.12 m增加至0.16 m。数值模拟的计算域大小与实验水槽尺寸相符,长11D,宽3.3D,高1.7D(图2)。计算域原点位于水轮机转子顶部中心。

图1 水轮机模型Fig.1 Model of the water turbine

图2 计算域尺寸Fig.2 Di mension of the computation domain

对入口的边界条件设定为速度入口,(x,y,z)向流速分别为(0.81,0,0)m·s-1,湍流强度为12%。顶部、两侧和底部的边界类型为对称边界条件“sy mmetry”,出口处设置为压力出口,水轮机的叶片和尾舱部分设置为无滑移壁面边界类型。

2 数值模型建立

2.1 湍流模型

DDES模型较DES模型[8]引入一个延迟函数,重新构造了DDES的长度尺度,同时考虑了网格尺度和涡黏场,避免了从RANS到LES的切换太过靠近壁面,也防止模型预测过早分离。本文基于SST模型的DDES-SST方法[9],其控制方程为

式中:ρ为流体密度;U为速度;k为瞬时湍流动能;lDDES为DDES长度尺度;ω为比耗散率;σk=0.85、σω=0.5、σω2=0.856、a1=0.31、α=1、β=0.09,均为模型系数;S为应变力张量;Pk、Pw为结果项;P为压力;μt为湍流涡黏度;F1、F2为SST的混合函数。式(1)中DDES的长度尺度lDDES为

式中:ILES、IRANS分别为LES、RANS的长度尺度。

DDES的长度尺度l 和涡黏性场、时间有关,其作用能使RANS湍流模型具有自我延续功能。如果数学函数fd表明某点位于边界层内,则使用长度尺度lDDES的模型可以判断转变为LES模式的时机。当出现大范围的分离流动时,fd值就会从0转变为到LES模型。流体发生分离后从RANS到LES的转换比DES方法更迅速。DES结合了RANS和LES的优势同时相比LES计算量减少,而DDES改善了DES存在的灰色区域和模化应力不足的问题。

2.2 网格划分

结构化网格能够对尾流场信息捕捉效果更好,网格较为规整。非结构化网格对复杂区域处理更为高效,同时也能保证数值模拟的精确性[10]。因此本文利用Point wise18网格划分软件,采用混合网格类型对潮流能水轮机进行网格生成。即水轮机叶片面网格、叶片周围的体网格以及水轮机尾流区域采用结构化网格,转子顶部、尾舱表面网格使用非结构化三角形网格(图3)。借助软件中的“Voexl”算法完成计算域中体网格的填充,该算法以四面体、金字塔、棱柱等类型填充过渡区域,在尾流区域填充六面体结构化网格(图4)。网格拓扑结构分为动域与静域两部分,动域包含水轮机叶片与转子,静域为流场其余部分(图5)。此拓扑结构有利于水轮机进行定常的多重参考系、非定常的滑移网格方法计算。

图3 水轮机面网格划分Fig.3 Surface grid divisions of the water turbine

图4 Voexl算法对叶片附近填充的体网格Fig.4 The volu me grids filled near the blades based on the Voexl algorith m

图5 计算域网格Fig.5 Grids of the co mputation do main

2.3 网格无关性分析

以粗网格、中网格、细网格类型,设计了3种密度的水轮机网格数量,详见表1。在加密动域的同时,静域网格数也在增加,使得动域、静域交界面网格尺寸大小相近,以保证网格连续性和良好的过渡性,提高数值模拟计算的精确性。

表1 5种网格的静域、动域网格数Table 1 The grid nu mber of the static and the dyna mic do mains of 5 kinds of grids

采用湍流“Spalart-All maras”模型进行稳态计算,停止标准为迭代5 000次,在迭代3 000次时水轮机推力值趋于稳定。水轮机网格在计算得出壁面的最大Y+值在5以内,出现在叶尖处,最小首层网格厚度为1×10-7m。随着水轮机网格的不断加密,推力系数也在增加。当网格加密至4 182 W时,推力系数监测值变化幅度减小,而粗网格和细网格的推力差距较大。综合计算精度、计算资源成本,本文最终选择细网格作为数值模拟的最终网格。

3 计算结果与讨论

3.1 水轮机性能结果验证

水轮机性能与流场参数定义如下:

式中:RTS为转速比;CT,CP分别为水轮机的推力系数和功率系数;P为功率;ETK为湍流动能(m·s-1)2;ω为水轮机旋转速度(rad·s-1);R为水轮机半径(m);U0为入口流速(m·s-1);T为转子推力(N);ρ为水的密度(kg·m-3);A为转盘面积(m2);Q为水轮机转矩(N·s);Ux、Uy、Uz分别为x、y、z流向的瞬时速度与平均速度差值(m·s-1)。

以细网格进行定常计算,设置水轮机RTS分别为5.37,5.52,5.75,5.98,6.52,6.97和7.51。通过功率系数、推力系数式(6)～式(7)计算得出水轮机在不同转速下的推力系数与功率系数。由图6可见,数值模拟结果低于实验结果的7%,但是仍处于实验采集数据的相对误差10%范围内。由于本文主要对尾流场发展规律进行研究,未考虑实验中的支撑结构,该部分误差可能来源于支撑结构的影响[11]。

图6 实验与数值模拟的功率系数和推力系数对比Fig.6 Co mparison of power coefficient and thrust coefficient bet ween the experi ments and the nu merical si mulations

3.2 水轮机尾流结果验证

在入口x向流速为0.81 m·s-1、转速为9.01 rad·s-1工况时,非定常模型计算的时间步长为9.69×10-4s,每步最大迭代次数8次。达到停止时间32 s后,水轮机在旋转10圈后流场达到稳定状态。取数值模拟结果中水轮机旋转最后5圈数据进行相位平均处理。由此得出水轮机尾流场中X=0.1D,0.5D和1.0D截面上监测点的平均速度,并与实验监测数据[7]进行对比。在Y=0处,因为尾舱的存在,X=0.1D,0.5D截面下,尾舱直径的部分实验与数值数据均无法采集。由图7可见,数值模拟的尾流速度与实验尾流的截面数据吻合良好,下游数值结果中出现的差异主要与数值模拟中使用的简化有关,其中没有对支撑结构进行建模,而尾舱被短圆柱代替。上述结果说明数值计算模型的计算精度较高,对后续尾流场及涡特性分析具备一定的可信度。

图7 3个截面下数值模拟与实验尾流监测点速度对比Fig.7 Comparison of the wake speeds at the monitoring points of three sections bet ween the experi ments and the numerical si mulations

3.3 水轮机尾流特性分析

本文以不同转速、不同流速设计了4种工况,以此进行尾流场结果的对比分析,并研究尾流场的发展规律。同时借助高性能计算集群,以192核的计算资源对潮流能水轮机单个工况并行计算,各算例详细信息及计算所耗时间见表2。

表2 数值模拟中各工况详细信息及计算耗时Table 2 The detailed infor mation and co mputation ti me of each working conditions in the nu merical si mulations

湍流是非常复杂的随机过程,对湍流动能的分析可以了解流场湍流随时间和空间变化的情况。为贴近真实的海洋流场情况,入口边界至X=3D处的湍流动能值较高。对湍流动能数值高于10-4的区域进行遮挡,以便对尾流场湍流动能观察更为清晰。由图8可见,4种工况中流场湍流动能均沿Z=0对称分布,低流速流场在X=5D处湍流动能开始减弱。此后,流场的发展以自由剪切流与尾流混合发展为主。湍流动能在X=8.5D处向中心聚集。高流速流场中,尾舱后的湍流动能数值为低流速流场3倍,高转速下的流场,湍流动能在发展至X=6D时,湍流动能有增加的趋势,在X=9.5D处湍流动能明显增加。这说明在高转速情况下,流场与水轮机转速后的尾流与自由流存在混合现象,增强了尾流的湍流效应,从而导致湍流动能增加。低转速下,湍流动能从X=5.5D开始减弱,但流场发展仍以水轮机旋转产生的尾流为主。

图8 4种工况下Y=0截面的湍流动能云图Fig.8 Cloud map of t he tur bulent kinetic energy at Y=0 section under the four wor king conditions

图9 中给出了4种工况下的涡量云图。涡量是描述流体的旋转运动特性的物理量,能表征涡特性及涡分离、脱落演化过程,由对速度的旋度计算得来。x,y,z方向涡量的计算公式为

式中:u、v、w分别为x、y、z方向的速度。涡量(单位为S-1)值大小的计算方法为‖ω‖=

由图9可见,从Y=0截面观察到水轮机叶尖后产生2个主涡带,转轴后为2个次涡带,轮毂涡在经过尾舱后发展成为一个涡带。由于叶片的高速旋转,而叶尖处形成叶尖涡,并逐步脱落形成叶尖脱落涡。4种工况下的叶尖脱落涡以规则圆形向后发展。并且沿流场方向直径增加,涡结构失稳直至消散。在工况1与工况2中能清晰观察到绕尾舱旋转的轮毂涡。由水轮机叶片根部旋转与尾舱的相互作用产生,并从尾舱末端脱落。高流速流场叶尖脱落涡在X=7D处开始减弱、消散,而低流速流场叶尖脱落涡在X=5D处开始消散。

图9 4种工况下的涡量云图Fig.9 Cloud map of the vorticity under the f our working conditions

图10 为4种工况下的尾流场涡量等值面图,通过对水轮机的三维涡量等值面图分析,可以观察到水轮机中涡结构的详细发展过程。为了更清晰地观察流场涡结构的变化,以不同涡量值形成等值面,设置工况1～工况4的显示涡量值分别为3.0、3.8、3.5和5.0。图中展现的叶尖脱落涡最终都是以椭圆形状脱落,但是脱落过程不同,低流速流场较高流速流场的脱落过程往往更加迅速,叶尖脱落涡形成的涡环连续性更强,涡环之间的距离减少0.1D。高流速流场中,在叶片旋转与流场作用下,转盘中间形成的涡数量更多。轮毂涡在低流速流场下,绕尾舱旋转发展的迹象更突出,能观察到轮毂涡旋转时形成的涡环。轮毂涡在经过直径变大的尾舱时,部分轮毂涡发生破碎。流场速度增加时,叶片与尾舱、流场之间的相互作用增强,轮毂涡的发展情况更为复杂,涡破碎现象更为明显。当流速一定时,转速增加,叶尖涡脱落形成的涡环直径增加。具体变化为工况1涡环直径由0.18D增加至与工况3的0.2D,工况2涡环直径由0.19D增加至与工况4的0.22D。

4 结 论

本文采用DDES方法对水轮机尾流场进行数值模拟研究,通过对比不同入流与水轮机转速尾流场的数值模拟结果,得出以下结论:

①4种工况下湍流动能均沿Z=0轴对称分布,转速一定时,流速越大湍流动能越高。流速一定时,转速越大流场的自由流与尾流混合现象越明显。

②水轮机尾流场涡结构主要由叶尖涡、叶尖脱落涡和轮毂涡组成。流场的入流速度与转速会影响尾流场中涡结构发展距离与涡量形成大小。转速一定时,流场流速越大,叶尖脱落涡、轮毂涡的发展距离越远。流速一定时,转速越大,涡的发展距离越短。

本文对比了不同工况潮流能水轮机尾流涡结构发展变化规律,研究结论可为潮流能水轮机的阵列布局提供技术支持,对未来环境保护、海床冲刷研究奠定基础。同时,研究结论可为潮流能水轮机的阵列布局提供技术支持,对未来环境保护、海床冲刷研究奠定基础。