漂浮式风电机组平台运动对叶片应力特性的影响

2021-06-04陈子文王晓东

华北电力大学学报(自然科学版) 2021年3期

陈子文，王晓东，丁坤，康顺

(1.华北电力大学电站能量传递转化与系统教育部重点实验室，北京 102206；2.华能新能源股份有限公司，北京 100036)

0 引言

现代社会发展对于能源需求逐渐增加，可再生能源尤其是风能得到迅速发展，海上优质的风资源受到关注。近年来海上风力发电发展迅速，风电机组叶片尺寸不断增大，同时风电机组在海上风况作用下也承受着较大载荷。与其他类型的旋转机械相比，海上风电机组风轮以及传动系统会受到10倍以上的疲劳负载[1]。风电机组在运行过程中，风载荷的施加会导致叶片发生弯掠、扭曲变形甚至断裂，对其结构造成严重影响。

国内外诸多学者对风电机组叶片的结构特性已经进行了广泛研究。在模型建立方面，吕品等[2]基于几何精确梁理论建立了叶片的非线性结构模型，并通过与线性模型计算结果对比验证了结构模型的正确性。喻涛涛等[3]基于流固耦合的方法对两种不同材料的风电机组叶片进行了数值模拟，分析了不同材料对叶片产生形变的影响，得到了实木材料的叶片随风速的增大变形更为明显的结论。张立茹等[4]基于风力机偏航的Glauert动量定理与流固耦合方法，将理论推导计算与轴流过程相结合，对风力机在偏航工况下的输出功率进行了分析，揭示了耦合与未耦合数值方法的差异性，研究表明耦合方法对于轴向诱导因子的计算值大于未耦合，从而导致输出功率存在差异。流固耦合不仅应用于对来流风况的研究，在叶片安全性研究方面也有所体现。周勃等[5]针对风电机组叶片裂纹扩展的复杂性进行了仿真分析，在不同风速条件下计算得到裂纹应力强度因子，结果表明裂纹的初始尺寸与风力载荷的大小均对叶片的裂纹扩展速率与方向产生影响,且在高风速下裂纹发生失稳扩展的几率增大。白叶飞等[6，7]基于应力无线遥测技术对风力机叶片应力应变进行了实验研究，得到风力机叶片翼型截面弦向最大主应力值大于剪切应力值，且随叶尖速比的升高而增大的结论，为叶片的动态应变安全检测提供了实验数据参考。在漂浮式风电机组叶片气动载荷方面，也有学者进行了一定研究，Wang等[8，9]提出多轴角运动模型，以大型海上漂浮式风电机组为研究对象，对风电机组平台运动下叶片的整体和展向气动载荷进行了研究。研究表明，平台运动会引起叶片沿展向的载荷波动，对漂浮式风电机组的运维安全提供气动性能参考。

对于漂浮式风电机组的研究大多围绕气动特性开展[10，11]，而平台多自由度运动下风电机组载荷及结构响应的研究开展较少。海上风、浪、流载荷的共同作用会使风电机组运动更为复杂，对叶片结构也会造成显著影响，且考虑叶片尺寸较大，所以对漂浮式风电机组叶片形变与应力分析进行研究。本文以NREL 5 MW风电机组为研究对象，分别对流体域和固体域建立计算模型，首先进行总体性能及叶片截面气动载荷分析，并沿叶片展向进行载荷性能研究。再基于漂浮式风电机组气动性能的模拟计算结果，并对所建立的叶片固体结构进行模态校核，并将气动结果加载至固体模型，得到沿叶片展向的气动载荷变化以及叶片挥舞形变；最后分析等效应力，讨论叶片在不同平台运动诱导下的应力分布。

1 数值方法

1.1 计算模型

NREL 5 MW风电机组为三叶片、上风向，采用变速变桨运行控制方式，适用于海上IEC风区1A级的风力气候条件。风轮直径为126 m，轮毂高度90 m。叶片由DU系列翼型与NACA64系列翼型构成。切入风速为3 m/s，切出风速为25 m/s，额定风速11.4 m/s，额定转速12.1rpm，转速变化范围为6.9 r/min～12.1 r/min，叶片展向不同位置处的翼型参数以及刚度、质量分布等参数参考文献[12]。本文所研究流体和固体模型忽略了机舱以及塔架，固体模型在轮毂与叶片交界面处施加约束，如图1和图2所示。

图1 叶片流体计算模型Fig. 1 Computational model of fluid domain

图2 固体叶片结构模型Fig. 2 Solid structure model of blade

其中，固体模型部分将叶片划分为多个叶元段，通过叶片各截面的刚度、质量分布对每个叶元段计算出面积惯性矩和杨氏模量等参数，进行等效结构建模，搭建了叶片壳模型。

1.2 流体计算域

如图3所示，计算域由四个区域构成：远场域、平台平动域、平台转动域和旋转域。外部长方体区域为一个长75D、宽高均为50D的远场区，其中D为旋转直径。考虑到海平面较为广阔，所以边界尺寸设置的足够大，以减小边界对计算精度的影响。区域入口风速设置为11.4 m/s的均匀风速，湍流强度为5%，远场域的上下平面分别设置为自由滑移壁面和无滑移壁面，叶片表面为无滑移壁面。出口边界为大气压力。平台平动域采用动态边界层法的动网格进行控制，平台转动域及风轮旋转域界面匹配连接。

图3 流体计算域Fig. 3 Computational domains

1.3 网格划分

流体域采用ICEM和AutoGrid5分别划分外场和风轮旋转域流场三维结构化网格，网格总数约为1 380万，如图4所示，控制平台运动区域的网格总数约为730万；风轮域网格数为396万，采用Autogrid软件绘制，并对叶片壁面附近进行了加密；叶片表面进行了细化处理，壁面y+<5，满足T-SST湍流模型计算要求。固体叶片壳模型采用Meshing划分，数量约为4.9万，沿叶展方向呈均匀分布，如图5所示。

图4 流体网格Fig. 4 Fluid mesh

图5 固体网格Fig. 5 Solid mesh

表1展示了三种不同网格数对应的风轮扭矩，沿展向截面翼型周向进行了不同程度的加密处理，分别为193万、396万和613万。结果表明，中等数量的网格可以较好的满足研究需要，所以本研究采用396万网格，相对误差最小，计算耗时适中。

表1 不同风轮网格数量对应的扭矩计算Tab.1 Torque of different gird sizes

1.4 计算方法

流体域采用Fluent软件求解非定常URANS方程，采用SIMPLE算法，湍流模型为带转捩的SST模型。以达到稳态的定常结果作为初始流场，进行非定常计算。非定常计算采用双时间步法，单位步长对应的相位角为2.5度。控制平台纵荡和纵摇运动的运动规律分别为

Vs=2πfAscos(2πft)

(1)

ωp=2πfApcos(2πft)

(2)

式中：A表示纵荡或纵摇运动的振幅；Vs为平台纵荡运动速度变化；ωp为平台纵摇运动角速度变化；f为其运动频率。分别对振幅As=2 m、Ap=4°对应下的纵荡运动和纵摇运动，以及平台纵荡-纵摇耦合运动进行研究分析，频率为f=0.1 Hz。其中平台纵荡运动为速度控制方程，平台纵摇运动为转动控制方程。当尾迹分布达到稳定状态后，取最后两个平台运动周期结果做研究分析。流体域计算收敛准则为：各区域残差下降6个数量级以上，总体性能参数随时间呈周期性变化，达到收敛准则。

流体模块与固体模块通过Ansys workbench进行连接，将流体模块计算数据传递至固体模块，叶片表面为流固耦合面，在叶片与轮毂界面处添加固定约束，在流固耦合面处进行压力加载，对固体模块进行模态校核及静力学分析。结构力学分析中，节点位移由应力通过刚度矩阵进行计算求解，如式(3)所示，

{F}={K}{d}

(3)

式中:{F}为节点应力载荷向量，{K}为风电机组刚度矩阵，{d}为结点位移向量。流固耦合计算遵循守恒原则，在流固耦合交界面处，满足应力相等(τ)、位移(d)、热流量(q)和温度(T)守恒:

(4)

其中，等式左侧为流体域变量，右侧为固体域变量，满足守恒条件进行载荷和节点位移的数据传递。

2 结果与分析

2.1 总体性能

在研究平台运动对叶片展向截面气动载荷分布之前，首先依据风轮在平台运动诱导下所处的不同位置分为四个不同时刻进行讨论，如图6。位置1表示平台向后运动的起始时刻，对应t=0；位置2表示平台运动至向后极限位置处，对应t=1/4T，T为平台运动周期；位置3表示平台回到最初风轮旋转平面并向前运动时，对应t=1/2T；位置4表示平台运动至前向极限位置处，对应3/4T。选取位置1和位置3，分别对应风轮所受平台诱导速度最小和最大处，进行叶片展向截面气动载荷分析。

图6 平台运动下风轮位置随时间变化Fig. 6 Variation of platform position

由于缺乏可用的实验数据，选择文献[13]与文献[14]的计算数据与计算结果进行比较，以验证平台纵荡计算数值方法的准确性。此外选取文献[15-16]与文献[17]的数据验证平台纵摇运动计算结果，流体计算验证及分析参考文献[9]。

首先对平台运动下风电机组的总体气动性能进行分析。图7和图8分别为三种平台运动下总体性能随相位角的变化。功率输出和推力展现出波动性，且最大功率大于额定功率。平台纵摇运动所引起的波动幅度大于纵荡运动，但二者均小于耦合运动，原因在于两种平台运动的正相关叠加使得风轮平面的风载增大，输出功率以及推力载荷展现了动态效应。总体较为剧烈的波动也反映出叶片气动载荷与结构载荷受到影响，后文将分别对其进行分析。

图7 功率随相位角变化Fig. 7 Power variation

图8 推力随相位角变化Fig. 8 Thrust variation

2.2 截面受力

为进一步研究风电机组叶片在展向的载荷分布，如图9，将扭矩以及轴向力分布沿叶片展向进行积分，做出在平台静止下风电机组叶片展向截面的扭矩输出以及轴向载荷分布。可以看出，对风电机组只进行非定常模拟时，三个叶片表现出相同的气动载荷分布，且扭矩输出和轴向载荷都集中于80%至90%叶展处。

图9 平台静止工况下截面气动载荷分布Fig. 9 Aerodynamic load distribution with platform stationary

图10为平台前向运动时，对应图6位置3处，叶片截面的气动载荷分布。分别对三个叶片以及风轮整体做研究分析，可以看到，与图7相一致，扭矩输出最大值为叶片外叶展80%至叶尖90%叶展处，表明该叶元段承担叶片整体最大做功量。所不同的是，平台纵摇与耦合运动下，三个叶片表现出不同的载荷特性，且均大于平台静止时叶片展向载荷。原因是在平台前向运动的时刻，叶片1处在风轮上半部分接近于竖直方向，所以具有较大的平台诱导速度，使其扭矩输出于轴向载荷也较高于叶片2和叶片3。而平台纵荡运动下，三个叶片的载荷特性保持一致，因其诱导速度与竖直平面时刻保持垂直，平台平动下三个叶片呈现出一致的扭矩分布与轴向载荷。平台纵荡、纵摇与耦合运动下，风轮的整体载荷依次增大。

图10 平台前向运动截面气动载荷分布Fig. 10 Aerodynamic load distribution with platform forward

图11为平台后向运动时，即图6所示位置1所对应时刻，风轮以及叶片截面的载荷特性。可以看到，因诱导速度与来流方向相反，平台纵荡运动诱导的沿叶片展向扭矩输出和轴向力分布均大于平台纵摇、耦合运动，但均小于平台静止时叶片的展向载荷。平台纵摇运动诱导下，叶片1处于反向诱导速度最大的位置，所以扭矩输出处于较低水平，且靠近叶尖的展向位置处，输出扭矩基本为零。参照轴向载荷沿叶展的分布，该振幅与频率所诱导的平台纵摇运动对于风电机组的影响为：降低扭矩输出，加剧叶片载荷波动。如图14(a)和(b)的耦合运动所示，纵荡运动的参与使得风电机组在平台诱导下展现出较为强烈的载荷波动，平台多自由度运动使机组平均输出功率降低，沿叶片展向载荷分布幅值增大。

图11 平台后向运动截面气动载荷分布Fig. 11 Aerodynamic load distribution with platform backward

2.3 叶片形变量

漂浮式风电机组运行过程中，叶片在风载与波浪的共同作用下会发生形变，其中以挥舞和摆振效应最为明显。风电机组叶片采用铺层设计，内部结构复杂，本文依据NREL 5 MW风电机组叶片各截面的刚度、质量分布以及面积惯性矩等参数进行等效结构建模，将叶片划分为多个截面，建立了有限元模型，在图2中进行了展示。等效建模过程中选取挥舞和摆振方向的弹性模量为基础，进行模态对比分析。结果如表2所示，图12为振型示意图。其中Mechanical为本文所采用计算方法，与FAST和Adams[18]所计算得到的结果进行对比。其中，一阶挥舞频率计算结果吻合较好，一阶摆振频率计算结果Mechanical与FAST有着接近8%的差别，原因在于本文以挥舞方向的弹性模量为基准进行建模，细长梁模型的摆振效应有所减弱[20]，侧重考虑风电机组叶片的挥舞变形。但总体误差较小，验证了所建固体模型的准确性。

表2 叶片固有振动频率Tab.2 Natural frequency of blade

图12 模态振型示意图Fig. 12 Modal diagram

图13所示为平台静止和平台运动下，叶片挥舞形变量(沿z轴方向)随相位角变化，其中0°对应叶片前向运动时刻(位置3)。可以看出，叶片沿展向形变量呈非线性变化，叶尖处具有最大位移量，叶根变形最小，图中形变增速区对应叶片变形梯度最大位置。平台静止时，叶尖位移量为2.18 m，沿周向保持不变。平台纵摇运动诱导下，叶尖在z方向的变形量范围为0.388～4.85 m，平台纵荡和耦合运动诱导下叶尖的变形变量在1.625 1～2.99 m、0.349 2～7.57 m范围内变化。与图10中平台耦合运动诱导出较大的轴向力载荷相对应，垂直于风轮平面诱导速度较大，加载至固体模型后，在挥舞方向引起最大的叶尖位移量。

图13 平台运动诱导下叶片挥舞形变随相位角变化Fig. 13 Variation of flapwise deflection respecting the azimuth angle under platform motion

平台向后运动时，对应图6中180°～360°相位角，可以看到，平台纵荡运动诱导下，叶片仍具有较为明显的形变量，但平台纵摇与平台耦合运动的诱导却使叶片形变量较小。这是由于较大的反向速度叠加使得耦合运动下叶片在挥舞方向产生的形变量比较小。通过两个运动过程中叶片挥舞形变的对比，可知平台运动向前以及向后运动过程中，叶片形变量随相位角也发生显著变化，这会增加叶片疲劳负载，表明叶片需要一定的强度来抵抗复杂平台运动所带来的应力形变。

2.4 等效应力

载荷施加会伴随应力集中，图14展示了平台静止和平台运动诱导下叶片等效应力分布随相位角的变化。平台静止时，叶片最大与最小等效应力分别为7 MPa和0.94 kPa；而平台耦合运动下，叶片等效应力最值分别为31.5 MPa、0.22 kPa。总体上看，叶片中部出现应力集中，与图13中的叶片挥舞形变量增速区相对应。在0～180°相位角范围内，平台运动诱导下，叶片所承载应力与静止状态相比较大，因为风轮前向运动过程中叶片载荷与来流风载相叠加；而在相位角180°～360°变化过程中，平台静止状态下风电机组叶片应力集中程度更为严重。平台耦合运动诱导下，在靠近叶片前缘处应力集中呈多核状分布。平台前向运动时，对应图6中的位置3处，压力面应力集中的现象较为严重。后向运动时，即图6中位置1对应的时刻，应力集中程度较弱，梯度较小，与叶片的形变量呈正相关。表明漂浮式风电机组在平台运动较为复杂的情况下，中叶段40%～70%内应力集中程度较为严重且随相位角变化较为急剧，所以需要对叶片中部应力集中部分进行材料或者结构强化。