潮流能水平轴叶轮纵摇运动水动力分析

张亮，王树齐，马勇，张学伟，国威

（哈尔滨工程大学深海工程技术研究中心，黑龙江哈尔滨150001）

潮流能水平轴叶轮纵摇运动水动力分析

张亮，王树齐，马勇，张学伟，国威

（哈尔滨工程大学深海工程技术研究中心，黑龙江哈尔滨150001）

漂浮式潮流能装置载体的波浪运动响应使水平轴叶轮的水动力特性发生变化。采用滑移网格及动网格技术对无界均匀来流中旋转叶轮强迫纵摇时的三维水动力特性进行分析，研究不同纵摇频率、纵摇幅值、速比等参数对叶轮水动力的影响规律，以及根据叶轮纵摇运动时的水动力时历曲线，拟合得到叶轮的纵摇阻尼系数。研究表明：叶轮轴向载荷和能量利用率瞬时值产生波动，波动幅值随纵摇频率、纵摇幅值及速比的增大而增大；叶片表面压力分布时刻发生变化，纵摇频率、纵摇幅值及速比越大，变化越明显；纵摇阻尼系数与纵荡频率、振幅无关，而与叶轮旋转速度相关，旋转速度越大，阻尼系数越大。

潮流能；水平轴叶轮；强迫纵摇；水动力载荷；阻尼系数；计算流体力学

21世纪以来，面对常规能源日益枯竭、环境恶化的局面，越来越多的国家开始关注清洁的可再生能源。其中潮流能［1］没有明显的丰枯年以及季节的大幅度能量变化，具有可持续、能量密度大、可预测、同时不占有陆地面积等优势越来越受到人们的青睐。潮流能水轮机［2⁃4］是潮流能电站的关键部件，它的水动力性能直接关系到电站的效益。目前，在水平轴水轮机［5⁃11］水动力性能研究中，基本都是水轮机叶轮只绕固定轴做旋转运动。而在实际漂浮式［12］潮流电站中，漂浮式载体由于波浪的影响会产生运动，从而导致水轮机跟随漂浮式载体一起运动。例如，当漂浮式载体发生纵摇时，水平轴叶轮也会跟着载体发生纵摇运动，改变了叶轮旋转平面处的速度分布，进而影响叶轮的水动力载荷。因此研究叶轮在纵摇运动条件下的水动力载荷显得尤为重要。

本文基于CFX软件，研究纵摇频率、幅值、速比等参数对叶轮水动力载荷的影响规律；以及对叶轮纵摇运动时的水动力载荷时历曲线进行拟合来获得阻尼系数。为研究漂浮式载体的波浪运动响应、水轮机的结构设计及电能控制等提供相关数据。

1 数值模拟

1.1 理论基础

水平轴叶轮旋转中心处的纵摇速度为u，来流速度为U，水轮机直径为D，叶片数为N，旋转角速度为ωT，为了方便分析，定义无量纲参数如下：

式中：λ为叶尖速比，R为叶轮半径，u-为无因次运动速度，a-为无因次运动加速度，Cz为叶轮轴向载荷系数，Fz为叶轮轴向载荷，Cp为能量利用率，T为叶轮轴向转矩，ρ为介质密度。Cz对u-、u-2、a-和a-2取变分可得

式中：nzz为阻尼系数，mzz为附加值量系数。考虑到叶片数N的影响，当水平轴水轮机以频率ωT旋转时，水轮机整体水动力以N倍频变化，因此式（1）中的各项水动力系数都应是以为基频的周期函数。

假设水动力系数可写成级数形式：

由于动态失速以及尾流场的存在，水平轴水轮机的流场非常复杂，很难采用直接的方法来求解相关水动力系数。本文基于最小二乘法原理，对叶轮纵摇运动时的水动力载荷时历曲线进行拟合来获得各项纵摇系数。

1.2 CFX数值模拟

1.2.1 建立计算域

计算模型的叶片由哈尔滨工程大学海洋可再生能源研究所研发，除叶片根部外，展向各位置处的截面翼型为S809。叶轮直径0.7 m，叶片数为2，轮毂直径为0.1D，如图1所示。在CFX数值模拟时，需要同时模拟水轮机的旋转运动和纵摇运动（a＝Asin（ωt），A为纵摇幅值，ω为纵摇频率），因此将整个计算域分为2个区域：静止域和旋转域，旋转域为圆柱体，旋转域又分为纵摇旋转域和叶轮旋转域，叶轮旋转域在纵摇旋转域的内部，叶轮旋转轴距纵摇中心0.4 m，计算域如图1所示。

图1 模型及网格示意图Fig.1 Schematic of blade and mesh models

1.2.2 划分网格

网格化分是数值模拟技术的一个非常重要的组成部分，网格质量的好坏直接影响到数值计算的精度和效率。整个计算域的网格划分采用结构网格和非结构网格结合的方式，叶轮旋转域由于叶片外形复杂，采用非结构化网格形式，纵摇旋转域和静止域采用结构化网格形式。叶轮表面网格的第1层高度为0.000 5 m，网格总数为201×104个。如图1所示。

1.2.3 边界条件设置

设定大气压为参考气压；将进口边界设为速度进口，给定均匀来流速度U和湍流参数；出口边界采用压力出口边界条件；流体计算域的侧面设定为自由滑动墙面；叶轮表面为无滑移壁面；旋转域给定纵摇角速度（滑移网格形式），叶轮旋转域又以子域的形式给定旋转运动（动网格形式），即叶轮旋转域在旋转运动的同时又和纵摇旋转域一起产生纵摇运动，网格不发生变形；旋转域和静止域及叶轮旋转域和纵摇旋转域之间通过交界面连接。计算中湍流模型采用的是SST模型，求解器为瞬态求解器，时间步长为叶轮旋转3°所用的时间。

2 叶轮载荷分析

在同样的网格模型下，水平轴叶轮在无界均匀来流中仅做旋转运动时的能量利用率及轴向载荷系数曲线如图2所示。图中实验值是在哈尔滨工程大学船舶拖曳水池中进行测试的，由图2可知，计算值与实验值吻合很好，验证了此网格模型模拟水平轴叶轮的准确性，下面对纵摇进行分析。

图2 计算值与实验值对比曲线Fig.2 Comparison between calculated and experimen⁃tal values

2.1 纵摇频率对叶轮载荷的影响

选取纵摇频率在常见波浪频率范围内，即取纵摇频率为0.4、0.8、1.2、1.6、2.0 rad／s。固定叶轮转速230 r／min，纵摇幅值取15°，入口速度1.5 m／s。图3给出了不同纵摇频率下叶轮轴向载荷系数及能量利用率随时间的变化曲线，ω＝0表示无纵摇的情况。由图3可知，当不发生纵摇时，叶轮的轴向载荷系数及能量利用率不随时间推移而改变。而当叶轮发生纵摇运动时，叶轮轴向载荷系数及能量利用率随着叶轮的纵摇而产生明显波动，纵摇频率越大，波动幅值越大。当纵摇频率为2.0 rad／s时，叶轮轴向载荷系数和能量利用率波动幅度分别达到叶轮不发生纵摇时的45.7%和103.2%左右。原因是纵摇改变了叶轮旋转平面处的速度分布，相当于来流速度的大小及方向时刻发生变化，从而使叶轮载荷产生波动。从图中还可以看出，在叶轮旋转过程中，叶轮的轴向载荷系数及能量利用率会发生微幅波动，这是由于叶轮旋转平面不是正对来流，最大存在15°的夹角。因此，在实际海洋环境中，当遭遇高频波浪时，叶轮载荷将发生高频率、高幅值波动，这对水轮机的电能输出的稳定性及结构安全可靠性来说是非常不利的。

图3 不同纵摇频率瞬时值对比Fig.3 Instantaneous value of different oscillation fre⁃quency

2.2 纵摇幅值对叶轮载荷的影响

取水平轴叶轮转速230 r／min，水轮机纵摇频率为1.2 rad／s，取3种不同纵摇幅值A（10°、15°和20°）进行计算，轴向载荷系数及能量利用率随时间的变化曲线如图4所示，从图4可以明显看出，纵摇幅值越大，轴向载荷系数及能量利用率瞬时值的波动幅值越大，这是由相对来流方向的改变引起的。在实际潮流电站中，水轮机运动与浮体平台的运动是相互耦合的，因此在潮流电站的设计中应尽可能减小浮体平台的运动幅度，以提高水轮机的结构安全性。

图4 不同纵摇振幅瞬时值对比Fig.4 Instantaneous values of different oscillation am⁃plitude

2.3 速比对叶轮载荷的影响

考虑叶轮在纵摇条件下不同速比时叶轮载荷的变化规律。取纵摇频率为1.2 rad／s，纵摇幅值为15°，来流速度1.5 m／s，取7种不同转速进行计算，轴向载荷系数及能量利用率瞬时值随时间的变化曲线如图5所示，从图5可以看出，在同一纵摇频率及幅值条件下，随着叶轮转速的增加，叶轮轴向载荷系数及能量利用率波动频率保持不变，但波动幅值增大，这样不利于电能的稳定输出及水轮机的结构安全。

图5 不同速比瞬时值对比Fig.5 Instantaneous values of different tip speed ratios

2.4 翼型表面压力分析

为了研究纵摇运动时载荷发生波动的机理，对叶片表面压力分布进行分析，对不同纵荡幅值，分别提取2个纵摇位置（1＃位置是纵摇平衡位置，纵摇速度与来流速度相反，叶轮处的实际流速是来流速度与纵摇速度之和；2＃位置是纵摇平衡位置，纵摇速度与来流速度相同，叶轮处的实际流速是来流速度与纵摇速度之差）时的叶片表面压力，图6给出了叶片0.4R处时的翼型表面压力分布。图中l是截面翼型弦长坐标归一化，Cpre表示表面相对压力系数，表达式为

式中：p为绝对压力，p0为参考大气压。

图6中，1＃、2＃分别表示1＃位置、2＃位置，A表示纵摇振幅，A＝0表示纵摇振幅为0，即无纵摇。从图6中可以看出，当叶轮在1＃位置时，随着纵荡幅值的增加，叶片压力面的压力逐渐增大，吸力面的压力逐渐减小，从而使叶片表面压力差随着幅值的增大而增大，而当叶轮在2＃位置时，随着纵荡幅值的增加，叶片压力面的压力逐渐减小，吸力面的压力逐渐增大，从而使叶片表面压力差随着幅值的增大而减小。同理，当叶轮在1＃位置时，叶片表面压力差随着纵摇频率的增大而增大，而当叶轮在2＃位置时，叶片表面压力差随着频率的增大而减小。因此当叶轮发生纵摇运动时，叶轮水动力载荷会发生波动，且波动幅值的大小随频率及幅值的增加而增加。

图6 表面压力系数分布Fig.6 Distribution of the airfoil surface pressure

3 纵摇阻尼系数拟合

3.1 纵摇频率对阻尼系数的影响

对不同纵摇频率下轴向载荷的时历曲线进行最小二乘法拟合，得出Cz级数形式的各项系数，如表1，表中只给出了Cz展开中前3项的常数项和一次项。

表1 Cz展开系数Table1 Czexpansion coefficients

从表中可以看出，Cz展开的各项系数中，前2项的常数项远大于一次项系数，且随着纵荡频率的增加，常数项基本不发生变化，因此，纵摇频率对阻尼系数没有影响。

3.2 纵摇幅值对阻尼系数的影响

对不同纵摇振幅下轴向载荷的时历曲线进行最小二乘法拟合，得出Cz级数形式的各项系数，如表2，表中只给出了Cz展开中前3项的常数项和一次项。

表2 Cz展开系数Table2 Czexpansion coefficients

从表中可以得到与3.1节相同的结论，即Cz展开的各项系数中，前2项的常数项远大于一次项系数，随着纵摇振幅的增加，常数项基本不发生变化，因此，纵摇振幅对阻尼系数没有影响。

3.3 速比对阻尼系数的影响

对不同转速（速比）下轴向载荷的时历曲线进行最小二乘法拟合，得出Cz级数形式的各项系数，如表3，表中只给出了Cz展开中前3项的常数项和一次项。

表3 Cz展开系数Table3 Czexpansion coefficients

从表中可以看出，Cz展开的各项系数中，前两项的常数项远大于一次项系数，随着转速（速比）的增加，阻尼系数常数项（绝对值）逐渐增加，因此，转速对阻尼系数影响明显，即转速越大，阻尼系数越大。

4 结论

采用CFX对三维水平轴叶轮在均匀来流中强迫纵摇时的水动力特性进行了计算分析，研究结果表明：

1）叶轮纵摇时，叶轮轴向载荷及能量利用率产生明显波动，波动幅值随着纵摇频率、纵摇幅值及速比的增加而增加，这对叶轮的结构安全、疲劳寿命及潮流电站电能的稳定输出产生影响，在结构和电控系统设计中需要特别考虑。

2）叶轮纵摇时，叶轮旋转平面处的速度大小及方向时刻发生变化，从而使叶片表面压力分布发生明显改变，这是导致叶轮载荷及能量利用率产生波动的原因所在。

3）根据最小二乘法原理，对叶轮纵摇时的轴向载荷进行数值拟合，得出纵摇时的阻尼系数，阻尼系数的大小与纵摇频率、幅值无关，而与叶轮旋转速度相关，旋转速度越大，阻尼系数越大。这为研究漂浮式载体的运动响应提供了重要的数据。

［1］吕忻，郭佩芳.我国潮流能资源开发评述［J］.海洋湖沼通报，2011（1）：26⁃30.LYU Xin，GUO Peifang.Review of China tidal current ener⁃gy resource development［J］.Transactions of Oceanology and Limnology，2011（1）：26⁃30.

［2］戴军，单忠德，王西峰，等.潮流水轮机的研究进展［J］.可再生能源，2010，28（4）：130⁃133.DAI Jun，SHAN Zhongde，WANG Xifeng，et al.Current re⁃search progress of water turbine［J］.Renewable Energy.2010，28（4）：130⁃133.

［3］陈展，马勇，张亮，等.矩形潮流能水轮机性能研究［J］.华中科技大学学报：自然科学版，2013，6（41）：128⁃132.CHEN Zhan，MA Yong，ZHANG Liang，et al.Study on properties of rectangular tidal turbine［J］.Journal of Hua⁃zhong University of Science and Technology：Natural Science Edition，2013，6（41）：128⁃132.

［4］LEE J H，PARK S，KIM D H，et al.Computational meth⁃ods for performance analysis of horizontal axis tidal stream turbines［J］.Applied Energy，2012，98：512⁃523.

［5］BATTEN W，BAHAJ A S，MOLLAND A F，et al.The pre⁃diction of the hydrodynamic performance of marine current turbines［J］.Renewable Energy，2008，33（5）：1085⁃1096.［6］LI Y，CAHSAL S M.A discrete vortex method for simulating a stand⁃alone tidal⁃current turbine：modeling and validation［J］.Journal of Offshore Mechanics and Arctic Engineering，2010，132（3）：1⁃9.

［7］CHOI H J，ZULLAH M A，ROH H W，et al.CFD validation of performance improvement of a 500 kW Francis turbine［J］.Renewable Energy，2013，54（SI）：111⁃123.

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

［9］MYERS L E，BAHAJ A S.Experimental analysis of the flow field around horizontal axis tidal turbines by use of scale mesh disk rotor simulators［J］.Ocean Engineering，2010，37：218⁃227.

［10］沈云，李龙，朱多彪.水平轴潮流能水轮机转轮尾流特性数值分析［J］.水电能源科学，2013，31（10）：149⁃151.SHEN Yun，LI Long，ZHU Duobiao.Numerical analysis of wake performance of horizontal⁃axis tidal current turbine［J］.Water Resources and Power，2013，31（10）：149⁃151.

［11］王晓航.潮流能水平轴叶轮设计与性能预测方法研究［D］.哈尔滨：哈尔滨工程大学，2011：3⁃8.WANG Xiaohang.Study on the design and performance pre⁃diction of tidal current horizontal axis turbine［D］.Harbin：Harbin Engineering University，2011：3⁃8.

［12］马勇，张亮，盛其虎，等.漂浮式潮流能发电装置水动力特性试验研究［J］.华中科技大学学报：自然科学版，2012，10（40）：123⁃127.MA Yong，ZHANG Liang，SHENG Qihu，et al.The test study on hydrodynamic characteristics of floating tidal power generation devices［J］.Journal of Huazhong University of Science and Technology：Natural Science Edition，2012，10（40）：123⁃127.

The pitch hydrodynamic analysis of tidal current energy horizontal axis impeller

ZHANG Liang，WANG Shuqi，MA Yong，ZHANG Xuewei，GUO Wei
（Deepwater Engineering and Technology Research Center，Harbin Engineering University，Harbin 150001，China）

Under the condition of actual sea state，hydrodynamic characteristic of floating horizontal axis impeller is related to wave characteristics and floating carrier wave motion response.This paper uses slipping mesh and dynam⁃ic mesh to analyze the hydrodynamic characteristic in uniform unbounded stream when the rotating impeller is forced to pitch，and studies the influence of different pitch frequencies，pitch amplitudes and tip speed ratios on turbine hydrodynamics.Pitch damping coefficient can be got by hydrodynamic time⁃history curve of pitching impeller.The results showed that the impeller axial load and energy utilization ratio will make a fluctuation in pitching and the wave amplitude of momentary value increases with the increasing of pitch frequency，pitch amplitude，and tip speed ratio.The blade surface pressure distribution changes when the impeller is pitching，and the larger the pitch fre⁃quency，pitch amplitude and tip speed ratio，the more obvious the changes.The frequency and amplitude of the pitch have little impact on damping coefficient，but this coefficient is positively proportional to the rotational speed of the impeller.

tidal current energy；horizontal axis impeller；forced pitching；hydrodynamic load；damping coeffi⁃cient；computational fluid dynamics（CFD）

10.3969／j.issn.1006⁃7043.201403015

http：／／www.cnki.net／kcms／detail／23.1390.U.20150109.1505.006.html

TK730；O352

A

1006⁃7043（2015）03⁃0307⁃05

2014⁃03⁃06.网络出版时间：2015⁃01⁃09.

国家自然科学基金资助项目（51309069，51309068，51106034）；海洋能专项基金资助项目（GHME2010GC02）；哈尔滨市科技创新人才研究专项资金资助项目（RC2014QN001008，2012RFQXG086）；中国博士后科学基金面上资助项目（2014M561334）；黑龙江省博士后基金资助项目（LBH⁃Z14060）.

张亮（1959⁃），男，教授，博士生导师；马勇（1980⁃），男，讲师，博士.

马勇，E⁃mail：mayong02＠hrbeu.edu.cn.