一种半潜式商用6 MW级风力机载荷数值模拟研究

方龙, 翟恩地, 李荣富, 宁巧珍, 赵斌, 李晔, 周雅容, 章丽骏

(1.天津大学 建筑工程学院, 天津 300350; 2.北京金风科创风电设备有限公司, 北京 100176; 3.南方海洋科学与工程广东省实验室(湛江), 广东 湛江 524002; 4.上海交通大学 船舶海洋与建筑工程学院, 上海 200240)

海上风电具有风资源储量高、远离居民区、选址方案限制低等诸多优势,近年来发展迅猛。我国近海海域的海上可开发风资源丰富,海上年均风速高,尤其是南海的部分海域,其年均风速值在5.5～10.7 m/s区间[1],适合安装大型海上风力机装置。根据支撑基础的形式,海上风力机可分为固定式和漂浮式。随着海上风电产业逐渐走向深远海,固定式风力机在其结构可靠性及经济性方面均难以满足需求。漂浮式风力机由于其在深远海海域的极强的适用性,已经成为全球海上风电领域新的研究热点[2]。我国也在浮式风电领域积极探索,国内首台漂浮式风电机组“三峡引领号”[3]于2021年12月在广东阳江海上风电场成功并网发电。2022年5月,国内首台深远海浮式风力机“扶摇号”[4]在广东湛江罗斗沙海域进行示范应用。这标志着我国在全球率先具备大容量抗台风型深远海漂浮式海上风电机组自主研发、制造、安装及运营能力。

浮式风力机工作在复杂的风浪流条件下,风力机、浮式平台以及系泊系统之间的作用相互耦合[5]。作用在叶轮上的力对系统的载荷有贡献,浮式平台的六自由度运动又改变了叶轮的气动性能,系泊系统进一步复杂化了整个漂浮式风力机系统的载荷与运动响应,漂浮式风力机的载荷分析是典型的多场-多体耦合问题[6]。同时,深远海环境的极端天气给漂浮式风力机的安全运行带来了巨大的挑战,因此,系统地评估漂浮式风力机在各种风浪流条件下的运动响应和载荷变化至关重要。由于漂浮式风力机成本昂贵,目前针对漂浮式风力机的研究以模型实验和数值分析居多。

在实验研究方面,国内外一些学者对不同形式的漂浮式风力机在不同环境下的动力响应进行了大量的实验研究。Koo等[7-8]对比了Spar和半潜式等不同形式的漂浮式风力机在不同波浪条件下的水动力响应,研究了不同浮式基础对塔筒固有频率的影响;Shin[9]研究了一种新型张力腿浮式平台在不同波浪条件下的水动力响应;Duan[10]评估了风载对浮式风力机纵荡和纵摇的影响,以及风载对系泊系统的动态响应的影响,揭示耦合响应的本质;Chen[11]针对模型浮式风力机的推力远小于实验弗劳德缩放值的问题,发展了一种叶片的高阶再设计方法,可以显著提高转子推力,提高气动力相似水平。

在数值分析方面,美国国家可再生能源实验室还与DNV GL、丹麦可持续能源国家实验室和挪威理工大学合作开发了数个数值模拟模块,如FAST[12]、HAOC2[13]等,并进行了代码比较协作[14-15]。张立等[16-17]结合叶素-动量理论、辐射/绕射理论及有限元方法,研究了不同风浪条件下对各种浮体形式的漂浮式风力机平台运动响应的影响;丁勤卫[18]基于辐射/绕射理论对比研究了垂荡板及其安装位置对漂浮式风力机Spar平台动态响应特性的影响。这些方法通常使用莫里森方程结合势流理论方法来计算水动力载荷,该方法一般忽略粘性项,引入阻力系数弥补粘性损失。此外,势流理论中的阻尼模型不能考虑与涡旋脱落有关的横向力或升力,这将影响了浮体在横向方向上的预测运动响应的精度[19]。蔡恒等[20]结合势流理论和叶素动量理论对比研究了单纯波浪以及风浪联合作用对浮式风力机系泊力的影响。姜劲等[21]对一种5 MW半潜式垂直轴风力机进行了气动、水动及系泊载荷的全耦合数值模拟,探究了速比对平台纵荡运动响应以及叶轮气动载荷的影响规律。

本文以一种商用6 MW级漂浮式风力机为研究对象,在荷兰MARIN实验室对其缩尺模型进行了风浪耦合试验研究。同时,基于Spalart-Allmaras湍流模型、流体-刚体动态相互作用模型和重叠网格技术,对该模型进行了气动力-水动力-系泊力的全耦合数值模拟。研究了波浪状况对浮式风力机平台运动状态以及叶轮气动性能的影响,并与实验值进行了对比,分析了引起浮式风力机气动载荷波动的原因。

1 研究对象

本文所述仿真模型为一种商用6 MW级浮式风力机的水池试验缩比模型,该机组根据我国南海的实际海况设计,缩比后的模型如图1所示。在该漂浮式风力机的水池实验[22]中,试验模型整体满足弗劳德数相似准则,该相似准则可以尽可能地满足模型机组与原机组之间的几何相似、动力相似和刚度相似。在满足弗劳德数相似的条件下,叶轮空气动力的雷诺数相似难以实现,但可以通过调整叶片的外形参数,达到模型机组与原机组在推力系数随叶尖速比变化趋势上的一致。限于实验条件,确定缩尺比例为1∶55,缩尺后的模型机组的基本尺寸如表1所示。

表1 实验模型浮式风力机基本物理量参数Table 1 Basic physical parameters of the scaled test model

图1 水池实验缩比模型Fig.1 Scaled model of the test in offshore basin

2 数值模型

本文基于非稳态的Spalart-Allmaras 雷诺平均湍流模型,利用重叠网格技术和滑移网格技术,以及DFBI六自由度运动模块,借助Star-CCM+软件模拟漂浮式风力机在风浪作用下的运动响应与载荷变化。SA(Spalart-Allmaras)湍流模型是工程上广泛应用的湍流模型之一,具有计算收敛性好和可靠性高等优点。对于漂浮式风力机,其运动响应受到水动力载荷和空气动力载荷的影响,同时,其运动状态也同样会影响作用在机组上的载荷。因此,漂浮式风力机载荷和运动响应是一个动态的相互作用过程。流体-刚体动态相互作用模型(dynamic fluid body interaction, DFBI)可以通过模拟刚体在流体域中的压力和剪切力作用下的运动,以及系泊线的恢复力,该模型可以很好地模拟浮式风力机与流体的相互作用。

2.1 计算域设置

仿真计算域的设置参考实验条件,图2给出了计算域示意图,入口边界和出口边界到浮式风力机的距离分别为5.2 m和12 m,两侧边界到浮式风力机的距离均为5 m,水深为2.46 m,顶部边界到水面的距离为7 m。为了尽量减少网格对于风速及波浪的耗散,在入口边界(图2中左侧边界)、顶部边界和两侧边界均设为速度入口。出口边界(图2中右侧边界)设为压力出口,以保证计算域内外压力平衡。在本文中,不考虑海床的粗糙度,将海底视作一个光滑的壁面,因此计算域的底部边界设为壁面。

图2 计算域尺寸定义及边界设置Fig.2 Size definition and boundary setting of the fluid domain

2.2 网格划分

本文选用切割体网格对计算域进行离散,如图3,利用重叠网格和滑移网格技术,将计算域划分为背景网格域、重叠网格域以及滑移网格域。重叠网格域将浮式风力机包络在其中,用来模拟机组的六自由度运动,滑移网格域为叶轮所在的区域,用该域网格的旋转来模拟叶轮的旋转。为了提高模拟精度,对叶轮的尾流区和水面线附近区域进行加密。叶片近壁面第1层网格高度为5×10-5m,经过验证,叶片大部分位置的y+值接近1,满足SA湍流模型的近壁面网格计算要求。

图3 网格划分示意Fig.3 Diagram of the mesh distribution

3 数值验证

仿真采用隐式非稳态计算,需要设置时间步长来实现时间上的离散。一般情况下,时间步长越小,模拟精度越高,计算量越大。反之,时间步长越大,计算量越低,模拟精度也越差。对于漂浮式风力机,空气动力和水动力的模拟对时间步长的要求不同,为了在控制计算成本时不影响模拟精度,本节分别从叶轮推力和浮体自由衰减2个方面来讨论时间步长对仿真结果的影响。

3.1 时间步长无关性验证

为了验证空气动力模拟所需的时间步长要求,本文在保持浮体固定的情况下模拟了叶轮的气动载荷,分别采用了0.05、0.03、0.02、0.01、0.005和0.003 s 6个时间步长。图4为叶轮推力系数随时间步长变化的结果,可以看出,时间步长越小叶轮的推力越大,并逐渐趋于平稳,时间步长为0.005 s和0.003 s的推力系数值的相对误差为0.17%,说明当时间步长满足Δt≤0.005 s时,可以达到空气动力的模拟精度要求。图5给出了4个时间步长下浮式平台在垂荡方向上的自由衰减模拟值,时间步长越大,垂荡自由衰减得越快,随着时间步长的减小,平台垂荡的自由衰减值趋向于收敛,时间步长为0.01 s和0.005 s的衰减周期的相对误差为2.05%,即当时间步长满足Δt≤0.01 s时,可以达到水动力计算的精度要求。

图4 叶轮推力系数时间步长验证Fig.4 Time step validation of the thrust coefficient

图5 自由衰减时间步长验证Fig.5 Time step validation for free decay

综合考虑水动力和空气动力的模拟精度要求,即当时间步长满足Δt≤0.005 s时,可满足模拟漂浮式风力机的精度要求,故后文模拟均采用0.005 s的时间步长。

3.2 土地利用动态度。在人和自然的共同作用下，区域内各种土地利用类型的数量在不同时间段变化的速度是不同的。土地利用类型的变化率与土地利用动态度可定量描述区域土地利用变化的速度。单一土地利用类型动态度可运用土地利用动态度模型，分析土地利用类型的动态变化，进而真实反映区域土地利用类型的变化速度和剧烈程度。公式为[13]：

3.2 浮体自由衰减对比

自由衰减通常用于确定浮体的自然周期,设置波浪条件为静水,给定浮体一个初始位移,然后释放,使其从初始位置自由移动,为了验证浮式风力机水动力模拟的准确性,对浮体在纵摇、横揺、纵荡和横荡上的自由衰减进行了模拟,并与实验测量值进行对比。图6和图7分别给出了浮式平台在纵摇和横揺上的自由衰减模拟值与实验值的对比。结合表2可以看出,模拟值在纵摇、横揺和垂荡上与实验值吻合程度较高,纵荡和横荡的模拟值在衰减周期上大于实验值,艏摇衰减周期的模拟值低于实验值,其原因为数值模拟中采用悬链线模型对系泊系统进行了简化,该模型计算得出的系泊系统作用在浮体上的载荷与实际值存在一定误差。因此,全耦合数值模拟的方法对浮式风力机水动力的模拟较为准确。

表2 浮式风力机自然周期Table 2 Natural period of floating wind turbine s

图6 纵摇自由衰减对比Fig.6 Comparison of free decay results of pitching motion

图7 横揺自由衰减对比Fig.7 Comparison of free decay results of rolling motion

4 浮式风力机的运动响应与载荷分析

本文分别对静水、规则波和非规则波作用下的漂浮式风力机进行了模拟,研究漂浮式风力机在不同外部条件作用下的载荷与运动响应,并在时域和频域上分析了漂浮式风力机的载荷波动与运动响应的关系。

4.1 静水工况

本文模拟了静水中的浮式风力机在来流风作用下的载荷与运动响应,来流风速设为U0=1.53 m/s,叶轮转速设为ω=9.32 rad/s,对应风力机叶尖速比为RTS=9.47。模拟过程中监测了浮式风力机的六自由度运动,通过转换可得到浮式风力机机舱的位移、速度和加速度。本文仅对来流风作用下机舱沿X轴方向(假设X轴正方向与风的来流方向相同)的运动进行分析。风力机机舱的沿X轴方向的位移xhub、速度vhub和加速度ahub的表达式为:

(1)

(2)

(3)

式中:xsurge代表漂浮式风力机平台沿X方向的线位移;h为机舱中心距离风力机重心的高度;θpitch代表漂浮式风力机绕Y轴的角位移。机舱相对速度指机舱运动速度与来流风速的比值,其表达式为:

(4)

图8给出了浮式风力机叶轮推力系数与机舱加速度波动的时历曲线,可以看出机舱加速度和叶轮推力呈正相关,当叶轮推力增大时,浮式风力机的倾覆力矩越大,机舱沿来流方向的加速度也越大;反之,当叶轮推力减小时,浮体回复力矩和系泊力的作用大于叶轮的推力作用,机舱加速度减小为负值。

图8 定常风静水工况下机舱加速度与叶轮推力系数Fig.8 Nacelle acceleration and turbine thrust coefficient curves under the condition of steady wind and still water

为了分析浮式风力机叶轮载荷波动的原因,图9给出了叶轮推力系数的频域分布,可以看出,推力系数的峰值①的频率为4.250 Hz,接近浮式风力机叶片的通过频率4.449 Hz,说明在纯风作用下,浮式风力机气动载荷的波动主要受叶轮转动的影响,一般情况下,叶片通过频率附近气动载荷的影响由风剪切和塔影效应造成,由于模拟中的采用的是均匀来流而非剪切来流,因此排除风剪切的影响,说明叶轮推力系数的波动是由塔影效应造成的。

图9 定常风静水工况下叶轮推力系数频域Fig.9 Frequency domain analysis of rotor thrust coefficient under the condition of steady wind and still water

4.2 规则波工况

对规则波作用下漂浮式风力机进行模拟,模拟与实验的来流、波浪以及风力机的运行参数一致,具体参数如下:来流风速为1.53 m/s,波高0.091 m,波周期1.079 s,浪向和来流风向一致,叶轮转速为9.32 rad/s,对应风力机叶尖速比为9.47。

由于纵摇和纵荡直接影响风力机的实际来流速度,推测这2种运动对风力机气动载荷的影响最大,通过合成纵摇和纵荡在机舱高度处沿来流方向的运动研究机舱运动与叶轮推力之间的关系。图10给出了机舱沿X轴方向相对速度与推力系数波动的时历曲线,模拟推力系数的平均值为0.703,实验测量的平均推力系数为0.671,相对误差为4.77%,说明全耦合数值模拟的方法较好地模拟了漂浮式风力机在风浪作用下的气动载荷。规则波作用下机舱运动的波动十分规律,呈现出正弦曲线的变化趋势。漂浮式风力机缩尺模型的纵摇运动时历曲线的数值模拟和水池模型试验结果如图11所示,可见两者在波动频率上几乎一致,而水池模型试验的纵摇幅值略大于模型试验,这主要是由于模型实验中的物面粗糙度和系泊结构的简化处理以及数值模拟时没有环境噪声的干扰。因此可以认为,在风浪耦合工况下,本文所述研究方法对于半潜式风力机的数值模拟结果可靠。此外,推力系数和机舱相对速度在1.078 s的周期上均有较大的波动,且波动峰值出现的时间节点接近,当机舱位移速度最小时,速度值为负,即机舱向来流的反方向运动,叶轮与空气的相对速度最大,此时叶轮的推力系数也达到了最大。相反,当机舱位移速度最大时,速度值为正,即机舱和来流同方向运动,叶轮与空气的相对速度最小,此时的推力系数也达到最小值,即叶轮推力系数变化和机舱沿来流方向的速度呈负相关。说明浮式风力机纵摇和纵荡运动对风力机的气动载荷影响较大,其原因为纵荡和纵摇运动使得机舱在来流方向上前后运动,进而影响了叶轮平面的实际来流速度。

图10 规则波下浮式机组机舱速度与叶轮推力系数Fig.10 Nacelle velocity and rotor thrust coefficient under regular wave

图11 漂浮式风力机纵摇运动曲线Fig.11 Pitch motion curves of the floating wind turbine

图12给出了浮式风力机推力系数的频域分布,可以看出,推力系数峰值①的频率为0.924 Hz,该频率接近波浪频率0.927 Hz,说明规则波工况下,浮式风力机气动载荷的波动主要受波浪的影响。进一步证明在波浪载荷的影响下,浮式风力机发生纵摇和纵荡引起叶轮实际来流速度的波动,从而导致了叶轮气动载荷的波动。推力系数的峰值①频率为4.276 Hz,接近叶片的通过频率4.449 Hz。根据4.1节的分析,该频率气动载荷的波动是由塔影效应造成的。图13和图14分别给出了规则波下浮式风力机纵荡和纵摇的频域分布,可以看出,浮式风力机纵荡和纵摇的峰值频率均出现在波浪频率0.927 Hz附近,说明规则波工况下,浮式风力机的运动的主要受波浪载荷的影响,且纵荡和纵摇运动的频率与波浪频率一致。图13与图14中在②号峰值频率左侧分别出现了一个低频频峰,该低频运动为平台的固有运动频率导致的。在图13中,②号峰值频率右侧分别出现了多个高频频峰,其分别为②号峰值频率的2、3和4倍频。由于叶轮的气动载荷相较于浮体的水动力载荷小得多,浮式风力机的纵荡和纵摇运动在叶片通过频率处的峰值不明显,说明波浪载荷对浮式风力机运动的影响比气动载荷大得多。

图12 规则波下叶轮推力系数频域图Fig.12 Frequency domain of rotor thrust coefficient under the regular wave

图13 规则波下浮式风力机纵荡运动频域图Fig.13 Frequency domain of the floating wind turbine′s surging motion under the regular wave

图14 规则波下浮式风力机纵摇运动频域图Fig.14 Frequency domain of the floating wind turbine′s pitching motion under the regular wave

4.3 非规则波工况

对非规则波作用下的漂浮式风力机进行模拟,来流风速为3.371 m/s,有义波高0.12 m,谱峰周期1.48 s,波谱采用JONSWAP谱,浪向与来流同向,叶轮转速9.32 rad/s,对应风力机的叶尖速比为4.29。

图15给出了非规则波下机舱X轴方向相对速度与推力系数波动的时程曲线,相较于规则波工况,非规则波工况下浮式风力机机舱的运动速度以及推力系数的波动更加无序。风力机推力系数和机舱运动速度直接相关,当机舱速度为正,即机舱向来流同方向运动时,叶轮与空气的相对来流速度减小,推力系数整体上低于平均值;反之,当机舱速度为负,即机舱向来流的反方向运动时,叶轮与空气的相对来流速度增加,推力系数整体上高于平均值,即叶轮推力系数变化和机舱沿来流方向的速度呈负相关,浮式风力机纵摇和纵荡运动对风力机的气动载荷影响较大。

图15 非规则波下浮式机组机舱速度与推力系数Fig.15 Nacelle velocity and thrust coefficient under the irregular wave

对非规则波工况下浮式风力机的推力系数进行频域的分析,如图16,推力系数频域峰值①的频率为0.623 Hz,接近非规则波谱峰周期对应的频率0.674 Hz,说明非规则波工况下,浮式风力机气动载荷的波动主要受波浪载荷的影响。推力系数频域峰值②的频率为4.468 Hz,接近风力机的叶片通过频率4.449 Hz,同理,该频率下风力机气动载荷的波动是由塔影效应造成的。图17和图18给出了浮式风力机在非规则波下纵荡和纵摇的频域分布,纵荡和纵摇峰值的频率均集中在0.674 Hz附近,而在叶片通过频率和叶轮转频上未出现峰值,同样说明在非规则波工况下,浮式风力机纵摇和纵荡主要受波浪的影响,叶轮的气动推力的波动对其运动的影响可以忽略不计。

图16 非规则波下叶轮推力系数频域图Fig.16 Frequency domain of rotor thrust coefficient under the irregular wave

图17 非规则波下浮式风机纵荡运动频域图Fig.17 Frequency domain of the floating wind turbine′s surging motion under the irregular wave

图18 非规则波下浮式风机纵摇运动频域图Fig.18 Frequency domain of the floating wind turbine′s pitching motion under the irregular wave

5 结论

1)模拟值在纵摇、横揺和垂荡上与实验值吻合程度较高,纵荡和横荡的模拟值在衰减周期上略大于实验值,艏摇衰减周期的模拟值低于实验值,其原因为本文所述研究方法中采用悬链线模型对系泊系统进行了简化,难以完全匹配浮式平台复杂的六自由度运动状态,计算得出的载荷与试验测得的载荷数值存在一定误差。

2)在静水工况下,浮式风力机仅受定常风作用,叶轮气动载荷的波动主要由塔影效应产生,机舱加速度的波动与叶轮推力系数变化呈正相关。

3)在规则波和非规则波工况下,叶轮推力的波动受浮式基础纵荡和纵摇的影响较大,其主要原因为这2种运动使得机舱沿来流方向前后运动,进而导致了叶轮的实际来流的速度不断发生变化,叶轮推力系数变化和机舱速度呈负相关。另外气动载荷的波动还与叶片的转动相关,这主要是塔影效应造成的。

4)在规则波和非规则波工况下,浮式风力机的运动主要受波浪载荷的影响,气动载荷对其影响相对较小。