海上浮式风机气动性能数值模拟

吴 俊，丁金鸿，何炎平，赵永生，刘亚东

(上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室，上海 200240)

海上浮式风机气动性能数值模拟

吴 俊，丁金鸿，何炎平，赵永生，刘亚东

(上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室，上海 200240)

采用计算流体动力学(CFD)方法，基于RANS方程和SST k-ω湍流模型，对OC3-Hywind Spar浮式平台支撑的NREL 5 MW风机进行气动性能模拟。对固定式风机的数值模型做网格无关性验证，同时考虑垂直风切变的影响，并将数值结果与NREL设计数据进行对比以验证模型的有效性。在FLUENT软件中,设定嵌入式滑移网格和用户定义程序(UDF)来模拟风机叶轮随平台的周期运动，分别研究浮式平台的纵荡、纵摇和首摇运动对风机气动荷载的影响。数值结果表明平台的纵荡和纵摇运动对输出功率影响较大，且平台运动幅值越大周期越低，其气动荷载变化越剧烈。合理控制平台的运动幅值对提高浮式风机的发电性能和疲劳强度有很大作用。

海上浮式风机；计算流体动力学；气动性能；六自由度运动

Abstract:Computational fluid dynamic (CFD) method is used to simulate the aerodynamic performance of floating offshore wind turbine (FOWT) experiencing platform motions.Based on the RANS equations and the k-ω shear-stress transport (SST) turbulence model,aerodynamic simulations for NREL 5MW turbine combined with OC3-Hywind spar buoy are studied.Grid independence studies are performed in fixed type offshore wind turbine,considering the effect of vertical wind shear,and the simulation results are compared with NREL published data to verify the numerical model.Furthermore,embedded moving mesh technique and user-defined functions (UDF) are adopted in FLUENT software to model the periodical motions of the rotor,and the effects of typical surge,pitch and yaw motions of the platform on turbine aerodynamic loads are investigated.The simulation results reveal that the surge and pitch motions have a greater impact on the fluctuation of the output power,and more significant variation of the aerodynamic loads would be introduced by larger amplitude and shorter period of the motion.The aerodynamic performance of the turbine and the fatigue strength of the structures would be significantly improved by reasonable controlling of the platform motions.

Keywords:floating offshore wind turbine (FOWT); computational fluid dynamics; aerodynamic performance; six-degrees-of-freedom motions

海上浮式风机在风浪流的联合作用下产生六自由度运动，使风机周围的流场发生剧烈变化，其复杂程度远远超出了陆上和离岸固定式风机的流场[1]。其中平台的纵摇和首摇运动将导致叶轮附近出现非轴向或扭曲的流场，且风机叶片也会与尾流相互作用，使叶片的气动荷载计算变得更加复杂。同时平台的纵摇运动使叶轮产生有效风剪切，故叶轮表面均匀风速的假设不再适用[2]。纵荡和纵摇运动引起的附加风速与入流风速叠加将改变叶片的相对入流速度以及叶轮-尾流之间的耦合作用。

目前大型风机气动性能计算方法主要有：叶素动量理论；涡格法(VLM)和CFD方法。叶素动量理论将叶素切片理论与简单动量方法相结合，并引入一些修正因子，包括动态失速效应，Glauert修正系数以及Prandtl叶尖损失因子等，该方法广泛地用到叶片设计和气动性能计算中。但对于浮式风机，由于高耸于平台上部的风机随平台产生剧烈运动，使风机处于风车和螺旋桨的瞬态变化状态中[3]，这种状态无法用传统的BEM理论来模拟。特别当风机往下风向纵摇时，叶轮与尾部湍流相互作用将产生垂直于叶轮的螺旋形回流，Sebastian[3]将这种瞬态气动效应定义为涡环状态 (VRS)。当风机处于涡环状态时，BEM理论中的动量方程将不再守恒[2]。还有一种方法称为涡格法，其是基于势流理论提出的，并在直升机旋翼和螺旋桨分析中得到广泛应用，Jeon等[4]曾使用该方法研究浮式风机在平台纵摇运动下的气动性能响应，但目前该方法在风机气动性能模拟方面应用较少。而CFD作为风机气动性能分析的新一代工具，能精确模拟浮式风机由于平台运动引起的复杂流场。目前，如Chow[5],Bazilevs[6]和刘强等[7]成功运用CFD方法研究了NREL 5 MW大型风机的气动性能和尾流特性。Tran等[8]用UBEM、CFD和FAST三种方法模拟了浮式平台的纵摇运动对风机气动荷载的影响并对三种方法进行对比。同时，斯图加特大学的研究团队利用FLOWer软件和重叠网格技术研究NREL 5 MW浮式风机在指定的纵摇运动下的非稳态气动响应[9]。任年鑫等[10]基于Fluent软件中的滑移网格技术，实现了浮式风机叶片在平台纵荡下与周围流场的耦合模拟；同时利用VOF方法和滑移网格，Ren等[11]成功建立了浮式风机与风-浪耦合流场之间的流固耦合模型。

这里利用嵌入式滑移网格和UDF模拟NREL 5 MW风机的旋转及随浮式平台的典型周期性运动。首先对模型做网格无关性验证，考虑垂直风切变的影响，并将无平台运动的数值结果与NREL数据对比以验证数值模型。分别模拟浮式风机在平台纵荡、纵摇和首摇运动下的气动荷载响应，并探讨不同运动幅值和周期对风机发电量和推力的影响，得出了平台运动对浮式风机气动性能的影响规律。

1 数值计算方法

1.1控制方程

基于雷诺时均的三维黏性不可压Navier-Stokes方程的连续性方程和动量方程分别为：

1.2湍流模型

采用SST k-ω湍流模型，该模型结合了k-ω模型在近壁区域计算的鲁棒性和精确性等优点以及k-ε模型在远场计算的优点，其输送方程可写为：

2 计算模型与网格

2.1浮式风机模型

研究对象选用美国可再生能源实验室(NREL)所开发的5 MW离岸风机[12]，具体几何参数如表1所示。沿叶片展向，截面的翼型不同叶厚的DU翼型以及NACA-64618翼型。下部支撑平台为Offshore Code Comparison Collaboration (OC3)研究计划所选用的OC3-Hywind Spar式浮式平台。

表1 NREL 5 MW海上风机总体参数Tab.1 The overall parameters of the NREL 5 MW offshore wind turbine

采用实尺度模型计算风机的气动性能，为提高计算效率，模型不包括塔架和浮式平台。整个计算域划分为外流场域、球形旋转域和叶片旋转域，如图1所示。其中球形旋转域用于模拟在平台转动下整个叶片旋转域的转动，而叶片旋转域用来定义叶轮相对球形旋转域的自转。每一时间步迭代前，都需用UDF定义旋转域的运动速度，由于叶片旋转域存在嵌套转动，故需要使用嵌入式滑移网格以定义其运动。两个区域彼此相互接触的面都设为交界面，交界面之间存在相对滑移并传递流场信息。叶片表面设为无滑移壁面，外部流场分别设置为速度入口和压力出口，其他四周均为对称壁面。

图1 计算模型示意Fig.1 Sketch of the whole computational domain

2.2网格划分

采用结构化网格，加密旋转域和叶片附近的网格，且对叶片赋予一定厚度的边界层。图2所示为叶片表面与周围流场网格分布以及叶片展向2/3R处的网格截面图。

图2 网格划分示意Fig.2 Illustration of the mesh

3 模拟结果与分析

采用FLUENT商业软件，选取SST k-ω湍流模型，采用压力基求解器，压力速度的耦合选取PISO算法，压力项采用PRESTO！格式进行离散。扩散项为中心差分格式，对流项、湍流动能方程以及湍流耗散率均采用二阶迎风格式。

3.1网格无关性验证

为确保数值模拟的精确性，需进行网格无关性验证。其中计算域宽度与高度均取10倍叶轮半径，上风向长度为6倍半径，为满足压力出口条件，下风向取20倍半径。对三种不同网格方案模型，在额定风速下进行气动荷载模拟，并将结果与设计值对比，如表2所示。由表可知方案1的网格数量最少，但模拟结果偏差较大；而方案2的网格数量较适中，且功率和推力的偏差与方案3较接近，为了提高计算效率，取方案2，最小单元尺寸为2 cm，网格总量为832万。

表2 不同网格精度的对比Tab.2 The comparison of different grid sizes

3.2考虑垂直风切变的影响

5 MW风机尺度巨大，轮毂高度距海平面90 m，沿高度方向的风速存在较大梯度，故需考虑垂直风切变的影响。根据API-RP-2A-WSD规范，采用对数风剖面模拟海上风速沿高度的变化：

其中，U0为海平面10 m高度处的1小时平均风速(m/s)。

将风速梯度编制成UDF，初始化流场速度入口，以模拟在额定风速下风机气动荷载的变化。由于越往高处，风速越大，转至高处的叶片受到的推力与扭矩比处于低处的叶片要大，而我们所关心的风机总体气动荷载是对三枚叶片求和得到的。如图3所示，虽然叶轮输出功率和推力存在一定波动，以叶轮转1/3圈为周期类似简谐曲线变化，但总体波动很小，在0.2%～0.4%之间，其中平均功率相对均匀风的工况降低了1.2%。总体来说，垂直风切变对风机总体气动性能影响较小，为了简化，以下工况均取速度入口为均匀风。

图3 垂直风切变对风机总体气动性能的影响Fig.3 The effect of vertical wind shear on overall aerodynamic performance of wind turbine

图4 风机后方不同位置截面的流场云图Fig.4 Contour of velocity from different sections behind the rotor

3.3数值模型验证

根据NREL 5 MW风机转速随来流风速的变化情况，分别模拟三个低风速工况、一个额定风速以及两个高风速工况。首先采用旋转参考坐标系法(MRF)进行稳态模拟，当收敛精度达到10-4的且气动荷载变化较小后，即认为流场稳定。而瞬态计算是在稳态计算的基础上采用滑移网格法进行模拟，当叶轮旋转4圈后，推力和扭矩基本稳定。图4显示了风机后方稳定流场截面的速度云图，当来流刚经过风机后，存在明显的涡流，随着流动距离不断加大，漩涡逐渐消失。将以上的数值结果与NREL的设计数据进行对比，由图5可知，CFD的模拟结果略小于NREL设计数据。当风速较小时，CFD结果与设计值相差很小，当风速高于额定风速时，相差越大，其中当风速达到15 m/s时，水平推力相差9.4%，而输出功率相差9.1%。这主要是由于当风速较大时，叶片中部翼型截面处的流场逐渐进入动态失速状态，故在有限的网格数量和所选湍流模型下很难精确模拟出该复杂流场[10]。鉴于模拟结果的误差在可接受范围内，从一定程度上也验证了该数值模型模拟海上风机气动荷载的有效性。

图5 风机推力和输出功率的对比Fig.5 Comparisons of fixed wind turbine thrust and rotor power

3.4平台纵荡的影响

当浮式平台发生纵荡时，风机与平台以相同的速度运动。为了简单起见，这里将纵荡运动简化为一定周期和幅值的简谐运动，如式(7)所示。对上式求导即为纵荡速度，当运动开始时，该速度从零开始增加，故流场不会因为平台的运动而发生突变。在开始加载纵荡运动前，先计算600步左右的稳态模拟和4个叶片旋转周期的瞬态模拟，以确保流场已充分发展并消除瞬态波动。模拟工况选取额定风速为11.4 m/s，叶片旋转速度为12.0 rpm，浮式平台运动周期分别为12 s和20 s，运动幅值为2 m和4 m，叶轮推力及输出功率的变化如图6和图7所示。

图6 不同运动幅值下浮式风机气动性能比较Fig.6 Comparison of blade aerodynamic performance under different surge amplitudes

图7 不同运动周期下浮式风机气动性能比较Fig.7 Comparison of blade aerodynamic performance under different surge periods

由图6和图7可知，风机气动荷载的响应曲线也近似为简谐曲线，与平台的运动规律基本一致。由于平台的纵荡运动产生的附加速度，使风机与入流风之间的相对速度改变，从而导致叶片各截面翼型的攻角、升力系数和阻力系数发生周期性变化，最终反映到风机推力和功率的周期性变化。在纵荡运动下，推力和输出功率的平均值与固定式风机相比降低1%左右，主要原因是气动荷载的简谐曲线关于中心轴不是对称的，即风机往下风向运动所减小的气动荷载值要略高于往上风向运动所增加的值。对于周期相同而幅值为2 m和4 m的两种纵荡工况，风机推力变化幅度分别为10%和18%，而功率变化为22%和68%。该结果表明在纵荡运动下，平台运动幅值越大，周期越小，气动荷载响应的变化越剧烈，且输出功率变化较为敏感，故合理控制平台运动幅值对保证风力发电的稳定以及风机结构包括塔架、叶片和减速齿轮的疲劳寿命都有重要意义。

下一步分析纵荡运动对每个叶片的影响。考虑到叶轮的对称性，只取一个叶片作为研究对象，平台运动选取周期为12 s而幅值为4 m的工况。图8给出了从纵荡运动开始到气动荷载达到最大值期间，推力和转矩沿叶片展向的分布情况。当t=20 s时，此时纵荡速度为零，故气动性能的分布曲线处于中间位置；而当t=23 s时，即风机的下风向速度达到最大，此时转矩和推力的分布都处于最低位置；当t=26 s时，风机的纵荡速度又回到零，故与t=20 s的曲线几乎重合；而当t=29 s时，风机的上风向速度达到最大，故叶片气动性能分布曲线也处于最高位置。由图8可知，叶尖和靠近轮毂部分的气动荷载变化较小，而变化最剧烈的范围集中在0.5R～0.9R之间，故此处需着重考虑叶片的结构和疲劳强度。

图8 在纵荡运动下气动荷载沿叶片展向的分布Fig.8 The spanwise distributions of aerodynamic loads along the blade experiencing surge motion

3.5平台纵摇的影响

将平台的纵摇角位移简化为简谐运动，根据OC3-Hywind spar结构，取轮毂以下90 m作为纵摇中心。对于锚泊式浮式风机，其典型纵摇运动的频率接近入射波浪的频率，故这里选取周期为10 s且角位移幅值分别为2°和4°的两种纵摇工况进行分析。图9显示了浮式风机纵摇过程的尾流变化情况，当风机纵摇处于平衡位置时，叶片与后方流场相互干扰较小，故尾流相对稳定。而当风机纵摇至上风向和下风向时，由于叶片的一部分受到尾流的影响，同时由于叶片的存在对流场产生一定扰动，故叶轮与尾流之间的耦合作用使尾流变化更剧烈。图10给出了在两种工况下风机气动荷载的变化。由图10可知，纵摇角度越大所引起气动荷载的波动也越剧烈。与无平台运动的情况相比，在幅值为2°和4°的纵摇运动下推力的平均值分别降低2.7%和5.2%，而两者的功率平均值基本一致，均降低4.2%。其原因可能是纵摇运动使风机在入流方向的投影面减小，纵摇幅值越大投影面越小，从而受到的水平推力也就越低。对于纵摇角度为2°和4°的运动，水平推力的变化幅度分别为17.5%和35.3%，而输出功率变化幅度分别为23.8%和65.9%。与纵荡运动类似，纵摇运动产生的气动荷载的剧烈波动也会对浮式风机的机舱控制系统和结构疲劳造成一定影响。

当风机发生纵摇时，风轮盘面存在的附加速度差会引起有效风剪切[3]，如图11所示。为了研究这种效应对风机气动荷载的影响，分别计算每个叶片上的荷载响应。图12给出了三个叶片在纵摇运动下的水平推力的响应曲线，其中叶片#1的初始位置在Z轴正方向，叶片#2和#3与Z轴分别成±120°。由于叶轮的旋转作用，三个叶片的气动荷载曲线存在相位差，其中叶片#1比#3超前2/3π相位，比#2落后2/3π相位。由于叶片#1的截面到纵摇中心的垂直距离最长，相应其附加速度的变化也是最大的，故如图12所示，作用在叶片#1上的推力的变化幅值比其他叶片要大。然而这种差别低于5%，故由于纵摇引起的风剪切效应对风机叶片气动荷载的影响并不大。

图9 浮式风机纵摇至不同位置的尾流速度图Fig.9 Velocity contours of wake flow during different pitch positions of floating wind turbine

图10 不同纵摇幅值下风机的气动性能对比Fig.10 Comparison of blade aerodynamic performance under different pitch amplitudes

图11 纵摇运动所引起的风轮盘面的有效风剪切Fig.11 Effective wind shear across the rotor due to pitching motion

图12 作用在每个叶片上的水平推力Fig.12 The aerodynamic thrust force for each blade

3.6平台首摇的影响

与纵摇类似，平台首摇运动也会在叶轮盘面处产生风剪切，并改变相对入流速度，同时也造成叶片桨距角的改变，即叶片不同截面翼型处的攻角也相应发生变化。在附加速度的影响以及叶片桨距角的改变的共同作用，分析风机在首摇运动下的气动荷载是非常复杂的。文中研究了风机在周期为10 s、角幅值为2°的首摇运动下的气动性能，图13给出了风机在该运动下的总体气动荷载变化以及叶片#1的推力和扭矩的变化。由图可知，除了初始阶段的数值波动外，风机的总体性能相对较稳定，而单独每个叶片的气动荷载存在较大的波动。为了研究叶片翼型截面在首摇作用下的气动性能，借助BEM理论进行定性分析，如图14所示。以叶片#1为例，当首摇运动开始时(t=10 s)，叶片绕Z轴的旋转力臂为零，即桨距角和附加速度都从零开始增加，故此时叶片的每个截面的攻角都没有改变。显然，附加速度与首摇角速度和转动力臂成正比，附加速度增加会导致翼型攻角和气动荷载的增大，而增大逆时针方向的桨距角会降低翼型攻角和气动荷载。当t处于10～11 s之间时，虽然由于桨距角的增大会减小攻角，但附加速度的增加占主导地位，故此阶段推力和转矩都增大。过了11 s后，由于首摇角速度减小使附加速度的增幅放慢，而此时桨距角的增大占主导地位，故叶片气动荷载下降。到了1/4T(t=12.5 s)时，首摇角速度减为零，即附加速度为零，同时逆时针方向的桨距角也达到了最大值，故此时叶片扭矩降到最低值，而由于升力在推力方向上的分量增加所带来的补偿，故水平推力并没有降到最低值而是略低于正常工况。

图13 在首摇运动下风机的气动性能Fig.13 The aerodynamic performance of the rotor experiencing yaw motion

图14 利用BEM理论分析平台首摇的影响Fig.14 Analysis of the effect of yaw motion using BEM theory

4 结 语

1)在模型的网格无关性验证基础上，模拟垂直风切变的工况，结果表明风机总体气动性能变化很小可以忽略。同时不同风速工况的模拟结果表明CFD的模拟结果略小于NREL设计数据。在低风速下，两者相差很小，而在高风速下，相差略大，然而误差在可接受范围内，一定程度上验证了该数值模型的有效性。

2)在平台纵荡运动下，浮式风机的推力和输出功率的平均值相比固定式减少了1%。平台运动幅值越大，周期越小，气动荷载变化越剧烈，同时功率的变化较推力值更为敏感。当平台纵荡时，叶片气动荷载变化最剧烈的部位集中在0.5R～0.9R之间。

3)平台纵摇角度越大，其气动荷载波动也越剧烈，且输出功率的变化比推力更敏感，同时由于纵摇引起的风剪切效应可忽略不计。在纵摇运动下，浮式风机的推力和输出功率的均值较固定式风机要低，且运动幅值越大推力的均值越低。

4)在平台首摇运动下，风机的总体性能相对稳定，而每个叶片的气动荷载波动较大。由于首摇角速度产生的附加速度以及叶片桨距角的改变，两者的共同作用使风机的气动性能预报变得更加复杂。

[1] SEBASTIAN T,LACKNER M.Analysis of the induction and wake evolution of an offshore floating wind turbine[J].Energies,2012,5(12):968-1000.

[2] SEBASTIAN T,LACKNER M A.Characterization of the unsteady aerodynamics of offshore floating wind turbines[J].Wind Energy,2013,16(3):339-352.

[3] SEBASTIAN T,LACKNER M.A comparison of First-order aerodynamic analysis methods for floating wind turbines[C]//Proceedings of the 48th AIAA Aerospace Sciences Meeting.2010.

[4] JEON M,LEE S,LEE S.Unsteady aerodynamics of offshore floating wind turbines in platform pitching motion using vortex lattice method[J].Renewable Energy,2014,65:207-212.

[5] CHOW R,VAN DAM C P.Verification of computational simulations of the NREL 5 MW rotor with a focus on inboard flow separation[J].Wind Energy,2012,15(8):967-981.

[6] BAZILEVS Y,HSU M C,AKKERMAN I,et al.3D simulation of wind turbine rotors at full scale.Part I:Geometry modeling and aerodynamics[J].International Journal for Numerical Methods in Fluids,2011,65(1-3):207-235.

[7] 刘强,杨科,黄宸武,等.5 MW 大型风力机气动特性计算及分析[J].工程热物理学报,2012,33(7):1155-1159.(LIU Qiang,YANG Ke,HUANG Chengwu,et al.Simulation and analysis of the aerodynamic characteristics of a 5 MW wind turbine[J].Journal of Engineering Thermophysics,2012,33(7):1155-1159.(in Chinese) )

[8] TRAN T,KIM D,SONG J.Computational fluid dynamic analysis of a floating offshore wind turbine experiencing platform pitching motion[J].Energies,2014,7(8):5011-5026.

[9] MATHA D,SCHLIPF M,CORDLE A,et al.Challenges in simulation of aerodynamics,hydrodynamics,and mooring-line dynamics of floating offshore wind turbines[R].National Renewable Energy Laboratory,US Department of Energy,Office of Energy Efficiency and Renewable Energy,2011.

[10] 任年鑫,李玉刚,欧进萍.浮式海上风力机叶片气动性能的流固耦合分析[J].计算力学学报,2014,31(1):91-95.(REN Nianxin,LI Yugang,OU Jinping.The fluid-structure interaction analysis of aerodynamic performance of floating offshore wind turbine blade[J].Chinese Journal of Computational Mechanics,2014,31(1):91-95.(in Chinese) )

[11] REN N X,LI Y G,OU J P.Coupled wind-wave time domain analysis of floating offshore wind turbine based on computational fluid dynamics method[J].Journal of Renewable and Sustainable Energy,2014,6(2):23106.

[12] JONKMAN J,BUTTERFIELD S,MUSIAL W,et al.Definition of a 5 MW reference wind turbine for offshore system development[R].National Renewable Energy Laboratory Golden,CO,2009.

Numerical analysis of aerodynamic performance of floating offshore wind turbine

WU Jun,DING Jinhong,HE Yanping,ZHAO Yongsheng,LIU Yadong

(State Key Laboratory of Ocean Engineering,School of Naval Architecture,Ocean and Civil Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

P753; TM614

A

10.16483/j.issn.1005-9865.2016.03.005

1005-9865(2016)03-0038-09

2015-04-30

国家重点基础研究发展计划(973 计划)项目(2014CB046200)；高等学校博士学科点专项科研基金新教师类资助课题(20120073120014)

吴 俊(1990-)，男，浙江湖州人，硕士，主要从事海上风机的气动荷载模拟。E-mail:wujun2009110@sjtu.edu.cn

丁金鸿。E-mail：ahaha@sjtu.edu.cn