漂浮式海上风机(floating offshore wind turbine, FOWT)在风、波浪和流等多重因素的共同作用下产生不规则的运动，其空气动力学特性及动态响应变得尤为复杂。FOWT所产生的不规则运动可视为不同幅值和频率的规则简谐运动的叠加。通过对简谐运动下FOWT气动特性的分析，可以建立其运动与气动载荷变化之间的对应关系，为研究复杂不规则运动下的动态响应奠定基础。
目前已有相关研究使用多种方法对FOWT的气动特性开展分析。Lijima等[1]将自行开发的水动力弹性时域分析程序SSODAC与美国国家可再生能源实验室(NREL)开发的耦合分析程序FAST结合，对风浪联合作用下的FOWT的耦合动力特性开展了分析。Manolas等[2]分别采用叶素动量理论(blade element momentum, BEM)方法和自由涡尾迹法对FOWT开展了气动-结构耦合动力分析，并指出了2种方法用于FOWT的气动分析存在的差异。Sebastian等[3]使用自行开发的基于自由涡尾迹方法的程序WInDS预测运动状态下的FOWT的气动载荷和尾流演变。Farrugia等[4-5]进一步开展了FOWT的气动试验，并应用WInDS程序对试验进行了数值仿真再现。Xu等[6]将三维旋转效应和动态失速效应补充到自由涡尾迹方法中，提升了非定常条件下气动性能的预测能力。Shen等[7]提出了一种基于自由涡尾迹方法的时间推进的升力面数值求解涡格法，利用该模型对纵荡下的FOWT的非定常气动性能进行了预测。对FOWT的非定常气动性能、释放涡稳定性和气动稳定性进行了详细的研究。Jeon等[8]采用涡格法分析纵摇运动下FOWT的非定常气动特性。Ding等[9]使用非定常叶素动量理论方法分析了不同运动幅值下的FOWT气动性能，通过计算得到六自由度平台推力、扭矩和功率的时域结果，总结六自由度平台运动下FOWT气动性能的变化规律。黄扬等[10]使用计算流体力学(computational fluid dynamics，CFD)方法开展了FOWT的气动-水动-气弹性耦合数值仿真。徐顺等[11]结合CFD方法和致动线模型，对复杂大气入流下海上风机力学特性开展了研究。李政[12]采用CFD方法对FOWT的气动特性进行了研究。
有学者对风机的不同气动分析方法进行了比较。万德成等[13]对目前主要的FOWT气动力计算方法和水动力计算方法进行了介绍。Sant等[14-15]应用数值仿真工具FAST和WInDS研究了张力腿式FOWT在不同偏航角和规则波浪条件下风机的推力和功率时变特征，且研究着重分析了在纵荡和纵摇运动下的结果。Thanh-Toan和Kim等[16-17]基于非定常叶素动量理论方法和CFD方法对纵摇和偏航运动下的非定常气动特性进行了研究。EIqatary等[18]使用CFD方法和BEM方法分别计算了变速和变桨工况内的风机功率，并进行了比较。结果表明变速工况内2种方法的计算结果一致，但变桨工况下的计算结果差异较大。Plaza等[19]分别使用CFD方法和BEM方法对不同风速下的MEXICO风机进行气动分析。结果表明低风速下BEM方法的准确性高于CFD方法，高风速下则相反。对比CFD方法，Tahani等[20]的研究表明BEM方法更适用于工程估算。
近年来，Wu等[21]对纵荡和垂荡运动下的双转子风机(dual-rotor wind turbines, DRWTs)的功率输出和尾流特性进行了研究。Yang等[22-23]将FOWT的推力波动等效为作用在转子上的力和力矩，基于推力波动表达式提出了一种频域解析模型计算FOWT的等效气动质量和阻尼，并进行了水池试验。Zhang等[24]结合CFD方法和有限元(finite element method, FEM)方法，对国际能源署(International Energy Agency, IEA)提出的15 MW的FOWT的叶片结构进行了单向流固耦合分析，研究中考虑了FOWT的结构运动。Zhang等[25]利用Reynolds平均Navier-Stokes(Reynoldsaveraged Navier-Stokes, RANS)湍流模型对DRWTs的功率和尾流特性进行了研究。Wang等[26]应用CFD方法研究了六自由度运动对FOWT的非定常气动特性和干扰效应的影响。此外，Wang等[27]和Luo等[28]通过数值模拟提出了新的FOWT气动相似准则用于风洞试验，为后续开展FOWT风洞试验提供了新的参考。Yang等[29]通过5 MW海上风机的水池试验对FOWT气动效应的影响进行了研究。研究发现风速增大导致FOWT运动诱导的气动阻尼减小，同时入流风的存在会降低塔筒的波频振动响应。
许多研究者采用BEM方法、自由涡尾迹法、涡格法或CFD方法对运动状态下的FOWT的气动特性进行了分析，并开展了FOWT的气动-结构耦合动力分析。但上述研究所使用的模型较为复杂，计算成本高，难以在实际工程中大规模应用。不同研究中使用的数据和实验条件存在差异，导致结果之间的对比和一致性存在一定问题。因此，本研究选取NREL 5 MW风机为研究对象，采用BEM方法对单自由度简谐运动下的FOWT气动特性进行研究。进一步在ExtPtfm模块中建立FOWT底部基础的超单元模型，并实现FOWT的单自由度简谐运动。最后，基于NREL开发的OpenFAST软件对六自由度运动下的FOWT开展全面的气动分析。
1 理论背景
1.1 气动载荷计算方法
本研究采用BEM方法分析简谐运动下风机所受到的气动载荷。在风机转子旋转平面内，以dr为宽度的圆环所受到的推力dT以及转矩dQ为[30]
其中：B为叶片数量，ρ为空气密度，Vtotal为叶素的总速度，Cl和Cd为叶片翼型的升力和阻力系数，φ为局部入流角度，c为弦长，r为局部半径。
结合动量理论和式(1)与(2)，每个转子圆环受相对速度的作用而产生的推力和扭矩为[30]
其中：U∞为平均风速，Ω为转子转速，a和at为轴向诱导系数和旋转诱导系数。
1.2 FOWT运动实现方法
ExtPtfm模块中所定义的平台质量需足够大, 以确保除自身重量外风机其他部分的重量对运动几乎不产生影响。该模块中所定义的平台的质量矩阵M的形式如下：
其中：矩阵中主对角线上的元素(统称为mDOF)表示单独每一运动自由度上平台的质量和惯性矩，独立于其他运动自由度；下标DOF表示平台的六个运动自由度。本研究不考虑各自由度运动之间的耦合作用，因此质量、刚度和阻尼矩阵的非对角线元素均为0[31]。
平台的刚度矩阵K和阻尼矩阵C的形式与质量矩阵一致。相应的kDOF和cDOF分别表示刚度矩阵和阻尼矩阵中主对角线上的元素，可由下式计算
考虑上述M、K和C，单自由度阻尼系统受简谐激励作用下的运动方程如下：
其中：x、$\dot{x}$ 和$\ddot{x}$ 分别为结构的位移、速度及加速度；FDOF为外部简谐激励的幅值，A为由外部激励产生的简谐运动的幅值，WDOF为简谐激励的圆频率，t为时间。
基于此，为使FOWT产生规定的简谐运动，在ExtPtfm中所需施加的外部时程激励可由下式计算：
其中：Fex为施加的简谐激励；xDOF为风机受到激励作用产生的位移。
1.3 风机一体化分析方法
本研究使用OpenFAST对进行简谐运动的FOWT开展耦合分析。OpenFAST基于多体动力学中的Kane方程对FOWT进行一体化分析[32]。一个具有n个自由度的简单完整约束系统的Kane运动方程如下：
其中：Fj为第j个广义主动力，Fj*为广义惯性力；该系统由w个刚体组成，每个刚体以参考系Ni表示，i为[1, w]。假设每个刚体Ni主动力作用在质心Xi位置，EvjXi和EωjXi分别为局部线速度和局部角速度。EajXi为质心Xi处的加速度。${ }^{\mathrm{E}} \dot{H}^{N_i}$ 为惯性参考系内每一刚体关于质心Xi的角动量对于时间的导数。上标E表示惯性参考系。
2 模型描述
2.1 NREL 5MW海上风机模型
2.2 浮式平台的超单元模型
为实现风机进行简谐运动，本研究在OpenFAST的ExtPtfm模块中将风机的浮式平台建立为超单元模型。在该模块中，对平台的质量、刚度及阻尼矩阵进行定义，具体定义方式见1.2节。为保证转子受到的风载荷不会影响风机的简谐运动，底部平台的质量被设定为非常大的值，并非参照某种实际平台的质量。此外，本研究的关注重点为简谐运动下FOWT的气动性能，因此这不会影响本研究的结果。本研究中浮式平台的M、C和K的具体取值如下：
2.3 FOWT一体化模型
本研究在OpenFAST中建立FOWT的耦合模型。浮式平台的超单元模型作为ExtPtfm模块的输入文件引入到OpenFAST中。平台在交界处的反作用力由ExtPtfm模块计算。超单元模型被整合到风机的耦合模型中，并对整体进行时域分析。基于输入文件，ExtPtfm模块接收交界位置的位移、速度和加速度。完整耦合模型如图 2所示。OpenFAST中的InflowWind模块为输入入流风载荷模块；AeroDyn模块为气动分析模块；ServoDyn模块为伺服动力学分析模块；ElastoDyn为结构动力学分析模块；ExtPtfm为风机底部基础结构动力学模块。
3 分析与讨论
3.1 工况设置
本研究考虑了浮式平台六自由度运动及其不同的幅值和频率。六自由度运动的幅值和频率设置如表 2所示。
表 2 单自由度简谐运动设置 |
| 运动 | 幅值 | 频率/Hz |
| 1 m | 0.100 0 | |
| 纵荡 | 1 m | 0.066 7 |
| 5 m | 0.100 0 | |
| 1 m | 0.100 0 | |
| 横荡 | 1 m | 0.066 7 |
| 5 m | 0.100 0 | |
| 垂荡 | 1 m | 0.100 0 |
| 5 m | 0.100 0 | |
| 1° | 0.100 0 | |
| 横摇 | 1° | 0.066 7 |
| 5° | 0.100 0 | |
| 1° | 0.100 0 | |
| 纵摇 | 1° | 0.066 7 |
| 5° | 0.100 0 | |
| 首摇 | 1° | 0.100 0 |
| 5° | 0.100 0 |
风机在0.1 Hz频率下的位移曲线如图 3所示。本研究中描述的平移运动包括纵荡、横荡和垂荡运动，旋转运动包括横摇、纵摇和首摇运动。图中的“A1”和“A5”分别对应幅值为1 m和5 m的简谐运动。由图可知，本研究中生成的风机的简谐运动良好。
3.2 单自由度简谐运动下FOWT的推力系数和转矩系数
本文中风机进行简谐运动的仿真总时长为3 600 s，时间步长设置为0.05 s。在仿真的初始阶段，数值结果表现出不稳定性。此外，风机还受到初始条件、自身固有频率及阻尼等因素影响，产生逐渐衰减的暂态响应。因此，本节分析不考虑仿真初始的600 s。
针对频率为0.1 Hz，幅值不同的六自由度运动下风机的推力系数CT和转矩系数CQ进行分析，统计结果如表 3—6所示。变异系数一定程度上反映了数据的离散程度。运动幅值为1 m的纵荡运动下，风机的CT和CQ的变异系数分别为0.023 71和0.035 28，分别是固定式风机的2倍和3倍(见表 3和表 5)；而在运动幅值为1°的纵摇运动下，则为固定式的3倍和6倍(见表 3和表 5)。其他自由度运动下变异系数与固定式基本相同。这说明在相同运动幅值和频率的情况下，六自由度运动中，进行纵荡和纵摇运动时风机的CT和CQ发生较大的波动。相较于纵荡和纵摇运动，其他四个自由度方向上CT和CQ的波动相对较小。运动幅值为5 m的纵荡运动下风机的CT和CQ增大至固定式风机的8倍和23倍(见表 4和表 6)；而在运动幅值为5°的纵摇运动下，则为20倍和44倍(见表 4和表 6)。这说明运动幅值增大导致风机的CT和CQ的波动更剧烈。
表 3 推力系数统计结果 |
| 运动模式 | CT平均值 | CT标准差 | 变异系数 |
| 固定式 | 0.713 20 | 0.007 81 | 0.010 95 |
| 纵荡 | 0.713 20 | 0.016 91 | 0.023 71 |
| 横荡 | 0.713 06 | 0.007 84 | 0.010 99 |
| 垂荡 | 0.713 05 | 0.007 99 | 0.011 21 |
| 横摇 | 0.712 90 | 0.008 48 | 0.011 90 |
| 纵摇 | 0.698 01 | 0.023 08 | 0.033 07 |
| 首摇 | 0.712 93 | 0.007 82 | 0.01097 |
注：幅值为1 m(或(°))、频率为0.1 Hz。 |
表 4 推力系数统计结果 |
| 运动模式 | CT平均值 | CT标准差 | 变异系数 |
| 固定式 | 0.713 20 | 0.007 81 | 0.010 95 |
| 纵荡 | 0.561 80 | 0.049 78 | 0.088 62 |
| 横荡 | 0.711 74 | 0.008 39 | 0.011 79 |
| 垂荡 | 0.712 33 | 0.011 68 | 0.016 39 |
| 横摇 | 0.683 24 | 0.043 21 | 0.063 25 |
| 纵摇 | 0.440 04 | 0.094 16 | 0.213 97 |
| 首摇 | 0.709 72 | 0.008 52 | 0.012 00 |
注：幅值为5 m(或(°))、频率为0.1 Hz。 |
表 5 转矩系数统计结果 |
| 运动模式 | CQ平均值 | CQ标准差 | 变异系数 |
| 固定式 | 0.065 12 | 0.000 76 | 0.011 61 |
| 纵荡 | 0.064 68 | 0.002 28 | 0.035 28 |
| 横荡 | 0.065 12 | 0.000 76 | 0.011 70 |
| 垂荡 | 0.065 12 | 0.000 78 | 0.011 96 |
| 横摇 | 0.065 11 | 0.000 88 | 0.013 48 |
| 纵摇 | 0.063 59 | 0.004 37 | 0.068 74 |
| 首摇 | 0.065 13 | 0.000 77 | 0.011 76 |
注：幅值为1 m(或(°))、频率为0.1 Hz。 |
表 6 转矩系数统计结果 |
| 运动模式 | CQ平均值 | CQ标准差 | 变异系数 |
| 固定式 | 0.065 12 | 0.000 76 | 0.011 61 |
| 纵荡 | 0.055 17 | 0.014 80 | 0.268 34 |
| 横荡 | 0.065 26 | 0.000 82 | 0.012 60 |
| 垂荡 | 0.065 22 | 0.001 27 | 0.019 43 |
| 横摇 | 0.065 26 | 0.000 82 | 0.012 60 |
| 纵摇 | 0.045 24 | 0.023 07 | 0.510 07 |
| 首摇 | 0.065 82 | 0.001 19 | 0.018 14 |
注：幅值为5 m(或(°))、频率为0.1 Hz。 |
3.3 简谐纵荡运动下FOWT的气动特性
3.3.1 纵荡运动下风机的轴向气动载荷
固定式风机和FOWT所受的轴向气动载荷时程曲线如图 5a所示，其中FOWT进行幅值为1 m、频率为0.1 Hz的纵荡简谐运动。在时域结果中，FOWT受到的轴向气动载荷的标准差为4.32×104 N，而固定式风机的为5.63×103 N。由此可见纵荡运动下风机受到的轴向气动载荷发生更为显著的波动。
进一步对叶片不同翼型截面位置的攻角进行分析，4处截面对应的翼型依次为DU35、DU25、DU21和NACA64。攻角变化曲线如图 7所示。由图可知，第3 017.5 s、3 019.5 s和3 021 s时刻，叶片3、叶片2和叶片1的攻角分别出现短暂的下降，这与气动载荷曲线出现短暂下降的位置一致。这说明叶片翼型的攻角的下降导致了风机气动载荷的减小。结合叶片方位角变化曲线，对应时刻FOWT的叶片2、叶片1和叶片3均位于180°方位，由此说明叶片经过塔筒是风机气动载荷出现短暂下降的原因。
3.3.2 纵荡运动下风机的气动转矩
风机的气动转矩的变化一定程度上反映了风机的发电功率的情况。由图 9可知，固定式风机的发电功率稳定，而进行纵荡下的风机的发电功率出现剧烈波动，且功率波动随运动幅值增加而愈加剧烈。此外，本文对固定式风机和进行不同幅值纵荡运动下的风机在所考虑的仿真时间内的发电量进行了计算。在时长为3 000 s的仿真内，固定式风机的发电量为3 823 kWh，纵荡幅值为1 m和5 m的风机发电量分别为3 834 kWh和3 835 kWh。相同时间段内进行纵荡风机的总发电量略有增加。
3.4 简谐纵摇运动下FOWT的轴向气动载荷和气动转矩
3.4.1 纵摇运动下风机的轴向气动载荷
3.4.2 纵摇运动下风机的气动转矩
不同幅值的纵摇运动下风机的发电功率变化曲线如图 12所示。与纵荡运动下的结果类似，纵摇运动同样造成风机发电功率出现明显波动，且运动幅值越大，发电功率的波动程度越剧烈。此外，纵摇幅值为1°和5°时，风机发电量分别为3 834 kWh和3 900 kWh。这说明与固定式风机相比，纵摇运动下风机的发电量略有增加。纵摇运动幅值的增加对风机总发电量的影响相对较小。
4 结论
本文对FOWT在六个自由度运动下的推力系数和转矩系数进行了统计分析，并重点研究了对气动特性影响显著的纵荡与纵摇运动。假设入流风向与纵荡、纵摇运动同向，得出以下结论：
1) FOWT的六个自由度运动中，纵荡和纵摇对推力系数与转矩系数影响最大。这2种运动与风向同向，叶片表面相对速度发生显著变化，导致轴向气动载荷和气动转矩随之波动。其他自由度运动对叶片表面相对速度的影响较小。
2) 纵荡和纵摇对风机气动特性的影响相似，均引入随运动周期性变化的载荷分量，该分量决定了气动载荷波动的幅值和频率。简谐运动引起的周期性相对速度变化是其根本原因。
3) 风机叶片经过塔筒时，轴向气动载荷与气动转矩出现短暂下降，归因于塔筒附近流速降低所致的叶片相对速度减小。
4) 纵荡与纵摇运动引起风机发电功率显著波动，但对总发电量影响较小。功率波动可能危害电网稳定性、加剧设备磨损并增加控制系统负担。因此，相较于发电量的小幅变化，更应关注功率波动的不利影响。
轴向气动载荷直接关系到叶片和塔筒的结构强度，运动引起的载荷波动可能增加叶片、塔筒及锚泊系统的疲劳破坏风险。未来研究应进一步关注纵荡和纵摇运动下FOWT的结构强度及疲劳破坏问题。
