随着化石能源的日益枯竭和全球经济的持续发展，寻找可持续新能源成为世界关注的热点，风能作为一种清洁、廉价而丰富的新能源受到广泛关注. 到2015年年底，全球风电装机容量较上一年增加17%，达到433 GW^[1]. 当前，国内外陆上风电的发展较成熟，但陆地空间有限，风力资源丰富的地区大多人口稀疏，远离耗电区域，陆上风电的开发潜力日渐降低. 海上风能资源丰富，世界上绝大部分海域风能密度大于200 W/m^2[2]，能够满足日益增长的能源需求，因此受到越来越多的关注. 考虑到经济性，在水深大于60 m的海域，固定式风电机组将不再适用,故深海海上风力发电多采用浮式平台基础^[3]. 目前，研究较多的风力机浮式平台基础包括驳船型、半潜型、TLP型和Spar型等.

各国学者利用各种数值模拟手段对Jonkman等^[4]提出的NREL 5 MW风机的各种浮式风力机概念进行广泛的研究和探讨. Jonkman等^[5]应用FAST，对ITI驳船型、TLP型、OC3-Hywinds和坐底式风力机进行模拟计算，在模型开发和载荷分析方面进行比较研究. Zhang等^[6]介绍海上风力机浮式基础设计的关键问题，如耦合动态分析模型和计算方法，建立有限元模型，利用SESAM软件计算不同工况下浮式风力机的总体运动响应和系泊力；结果表明，该半潜型浮式基础和系泊系统可以在有义波高小于4 m的海洋环境中安全运行. Cheng等^[7]利用SIMO-RIFLEX-DMS程序，对Spar型、半潜式和TLP型3种垂直轴浮式风力发电机进行全耦合非线性时域仿真分析比较；结果表明，与TLP型浮式基础相比，Spar和半潜型浮式基础由于采用悬链系泊系统，可以减轻对结构负载和系泊张力的2P气动载荷影响，代价是具有更大的平台运动. Tran等^[8]基于重叠网格，采用计算流体力学（CFD）方法对半潜式风力机进行耦合动态响应计算，采用VOF的耦合流体流动和多体动态分析方法，研究半潜式浮体的典型水动力响应. Zhao等^[9]基于势流理论和BEM理论，对5 MW多柱张力腿式风力机进行动态耦合分析，比较与陆地式风力机的载荷影响差异，分析机舱加速度、塔底弯矩和平台纵荡运动响应；结果表明，该TLP式风力机具有良好的工作和生存能力. Goupee等^[10]采用1/50的缩尺比，对TLP、半潜式和Spar型浮式风力机在不同的工况下进行实验比较；结果表明，不同的浮式基础各自具有优缺点，如TLP具有较小的机舱加速度和总体运动响应，但承受的系泊载荷远大于Spar和半潜式平台.

本文针对设计出的新型6 MW Spar型浮式风力机，考虑气动载荷和水动力载荷，利用气-液-固-弹性数值模拟软件FAST及Sesam软件，分析并验证该Spar型浮式基础在不同风浪组合下具有较好的总体运动性能，且首摇角在delta-line的控制下不会过大.

1 6 MW风力机设计参数

目前，国内外浮式风力机实型机发展仍未成熟，针对100 m水深条件设计出新型Spar浮式风力机.

1.1 风机和塔架

风机为传统的水平轴上风向三叶片中速永磁型风力发电机组，其中叶片采用NACA及DU系列改型的翼型. 塔架材料采用Q345钢，塔底位于静水面以上13 m处，塔身高度为83.806 m. 塔顶要承受机舱和风机的载荷，故塔顶厚度提高. 具体参数如表1所示，风机的弦长和扭角等外形数据如图1所示. 图中，L_c为弦长，T_d为扭角，R_t为相对厚度，T_pr为预弯，d_b为叶片展向距离.

在额定推力下，塔架底部弯曲应力为

$\begin{split} \sigma = &\frac{{Mh}}{W}= \displaystyle\frac{{891.7 \times 100}}{{\displaystyle\frac{\pi }{{32}} \times {8^3} \times \left[ {1 - {{\left( {\displaystyle\frac{{8 - 2 \times 0.04}}{8}} \right)}^4}} \right]}}=\\ &45.02\;{\rm MPa}. \end{split} $

该结果远远小于Q345钢的屈服强度，塔底外径与壁厚之比满足DNV规范^[11]关于屈曲的要求.

1.2 Spar型平台

平台作为浮式风力机的重要支撑结构，平台性能对浮式风力机的总体运动响应起到重要影响. Spar型平台具有小水线面和大排水量的特点，能够提供浮式风力机所需的浮力和静水回复力. 平台顶部外径为8 m，底部外径为14.6 m. 表2给出平台的几何参数. 图2、3给出该风力机整体示意图和浮式基础的结构尺寸外形二维图，坐标系和平台自由度如图4所示.

设计出的该浮体的静倾角为4.3°，符合设计初稳性纵倾角不超过5°的设计原则^[12]. 平台各自由度的固有周期都避开了波浪的能量集中范围. 通过进一步的计算表明，该平台的设计是可行的.

1.3 系泊系统

平台的垂荡、横摇和纵摇运动由静水回复力、力矩来控制，由于外力作用而导致的偏移在外力消失后能够回到原有的平衡位置；纵荡、横荡的偏移须由锚链提供回复力. 该系泊系统包括3根互相成120°的悬链线，如图5所示，锚链的直径和单位长度质量及破断力根据API^[13]的规范确定. 锚链设计的另一个重要方面在于控制平台的首摇运动，采用连接到导缆孔处的delta-line提供首摇回复力. 附加首摇刚度^①为

① NIELSEN F G. Tuning the geometric and mass properties of hywind. Statoil (Internal document), 2009.

${C_{66}} = \sum\limits_{i=1}^n {TR\cos \varphi \left( {1 + \frac{R}{{{L_{\rm H}}}}} \right)} .$

(1)

式中： ${C_{66}}$ 为附加首摇刚度， $T$ 为锚链顶端张力， $R$ 为有效锚链线连接点到平台中心的水平距离， $\varphi $ 为锚链顶部与水平面的夹角， ${L_{\rm H}}$ 为锚与有效锚链线连接点的水平距离. 表3给出系泊系统的具体参数.

2 数值模型 2.1 系统时域运动方程

耦合风机、塔架、平台和锚链运动的非线性时域运动方程为

$\begin{split}&\displaystyle\sum\limits_{j = 1}^6 {\left\{ {\left( {{M_{ij}} + {m_{ij}}} \right){{\ddot x}_j}\left( t \right) + \int\limits_{ - \infty }^t {{{\dot x}_j}\left( \tau \right){K_{ij}}\left( {t - \tau } \right){\rm d}\tau + } } \right.} \\&\left. {{b_{vij}}{{\dot x}_j}\left( t \right) + {C_{ij}}{x_j}\left( t \right)} \right\}={F_i}^{{\rm Aero}}\left( t \right) + {F_i}^{{\rm Hyd}}\left( t \right) + {F_i}^{{\rm Moor}}\left( t \right) .\end{split}$

(2)

式中： ${M_{ij}}$ 为耦合系统的质量矩阵， ${m_{ij}}$ 和 ${K_{ij}}$ 分别为时域内的附加质量和延迟函数， ${C_{ij}}$ 为静水回复刚度系数， ${F_i}^{{\rm Aero}}\left( t \right)$ 、 ${F_i}^{{\rm Hyd}}\left( t \right)$ 和 ${F_i}^{{\rm Moor}}\left( t \right)$ 分别为气动载荷、水动力载荷和锚链力在第i个运动自由度上的分量.

2.2 气动载荷

将基于叶素理论和动量理论的BEM理论应用于FAST程序中. 基本原理是将叶片沿展向离散成若干个独立的叶素，并在每个叶素中均定义一个叶切面；根据尾部气流状况，计算每个叶切面上的升阻力并投影到叶轮平面的法向和切向；然后将作用于所有叶切面上的力沿着叶片展向积分，可得作用于叶片的推力和扭矩，进而得到作用于整个风机的推力和输出功率. 作用于整个叶片的推力、扭矩和功率如下：

$T = \frac{1}{2}\left( {{C_{\rm L}}\cos \phi + {C_{\rm D}}\sin \phi } \right)\rho {v^2}{(1 - a)^2}cB\frac{{\Delta r}}{{{{\sin }^2}\phi }},$

(3)

$\begin{split}M =& \frac{1}{2}\left( {{C_{\rm L}}\sin \phi - {C_{\rm D}}\cos \phi } \right)\rho v\times\\ &(1 - a)\left(1 + a'\right) cB\frac{{\Delta r}}{{\sin \phi \cos \phi }}\varOmega r,\end{split}$

(4)

$P = M\varOmega r.$

(5)

式中： $a$ 和 $a'$ 分别为轴向和周向诱导因子， ${C_{\rm L}}$ 和 ${C_{\rm D}}$ 分别为升力、阻力系数， $\phi $ 为相对风向角， $v$ 为上风向风速， $\Delta r$ 为叶片截面的径向长度， $\varOmega $ 为角速度， $c$ 为翼型弦长， $B$ 为叶片数量. 从叶尖来看，叶素的横切面图面^[14]如图6所示.

对于浮式风力机，受到的风载荷叶片气动载荷包括塔架受到的风载荷. 塔架受到的风载荷可以根据中国船级社（CCS）^[15]规范计算得到：

$F = 0.613\sum\limits_{i = 1}^n {\left( {{C_{{\rm h}i}}{C_{{\rm S}i}}{A_i}{V_{10}}^2} \right)} .$

(6)

式中： ${A_i}$ 为沿塔架高度方向第 $i$ 个受风分段在受风方向上的投影面积， ${C_{{\rm h}i}}$ 为第 $i$ 个受风分段的高度系数， ${C_{{\rm S}i}}$ 为第 $i$ 个受风分段的形状系数， ${V_{10}}$ 为海平面以上10 m处的参考风速.

2.3 水动力载荷

与陆上风力机相比，海上浮式风力机还受到水动力载荷的作用. 频域内的波浪激励力计算公式为

${f_{{\rm w}j}} ={\rm i}\omega \rho \iint\limits_{{S_0}} {\left( {{\varphi _{\rm I}}\frac{{\partial {\varphi _j}}}{{\partial n}} - {\varphi _j}\frac{{\partial {\varphi _{\rm I}}}}{{\partial n}}} \right){\rm d}S}.$

式中： ${f_{{\rm w}j}}$ 为对应j运动模态的波浪激励力, $\omega $ 为波浪频率， $\rho $ 为海水密度， ${\varphi _{\rm I}}$ 为入射波速度势， ${\varphi _j}$ 为对应j运动模态的辐射速度势， ${S_0}$ 为湿表面面积.

入射势可以基于微幅波理论直接求解，辐射势通过求解具有自由面条件、物面条件、海底条件、无限远边界条件的拉普拉斯方程得到，进而求出附加质量、阻尼系数.

对于波浪力在时域内的计算，利用Cummins^[16]提出的时域中波浪力与频域中波浪力传递函数之间的卷积关系，通过对脉冲响应函数和波面升高的卷积积分形式来表示波浪力的时程：

${F_i}^{{\rm Hyd}}\left( t \right) = \int\limits_0^t {{H_j}\left( {t - \tau } \right)\eta \left( \tau \right){\rm d}\tau } .$

(7)

式中： ${H_j}\left( \tau \right)$ 为第j个运动自由度阶波浪力的脉冲响应函数，可以通过对频域波浪力传递函数进行Fourier变换得到； $\eta \left( t \right)$ 为平台处的波面瞬时高度.

对于本文的Spar浮式基础，通过势流理论计算得到的波浪力需要考虑由流体黏性引起的拖曳力进行修正. 根据Morison公式可知，由流体黏性在圆柱切片引起的水平拖曳力可以表示为

${\rm d}{F_i}^{{\rm viscous}}\left( {t,z} \right) = \frac{1}{2}{C_{\rm D}}{\rho _{\rm w}}\left( {D{\rm d}z} \right){u_i}\left( {t,z} \right)\left| {{u_i}\left( {t,z} \right)} \right|.$

(8)

式中： $D$ 为圆柱的直径； ${\rm d}z$ 为圆柱沿径向分成的切片长度； ${\rho _{\rm w}}$ 为海水密度； ${C_{\rm D}}$ 为黏性拖曳力系数，此处取0.7； ${u_i}$ 为流体水平速度.

2.4 锚链力

锚链力计算采用考虑惯性力和阻力的有限元（FEAM）模型，将锚链视为弹性杆模型^①. 建立流体中杆模型的控制方程如下：

① BAE Y H. Development of a dynamic mooring module FEAM for fast V8, 2014.

$\rho \mathop r\limits^{..} + {C_{\rm{a}}}{\rho _{\rm{w}}}\mathop {{r^n}}\limits^{..} + \left(EIr''\right)'' - \left(\overline \lambda r'\right)' =m + {\overline F ^{\rm{d}}},$

(9)

$\frac{1}{2}\left(r' \cdot r' - 1\right) = \frac{T}{{AE}} \approx \frac{\lambda}{{AE}}.\;$

(10)

式中： $q$ 为单位长度锚链线的力， $\rho $ 为单位长度质量， $F$ 和 $M$ 分别为沿着杆中心线的力和力矩， $T$ 为杆的张力， $m$ 为杆在水中的质量.

2.5 数值模拟工具和环境工况

采用2种数值模拟工具：挪威船级社（DNV）开发的Sesam软件和美国国家可再生能源实验室（NREL）开发的FAST软件.

通过Sesam软件中的Wadam模块进行频域水动力计算，得到关于Spar平台的水动力系数、静水回复力、单位波幅下的不同频率和方向的波浪激励力等一系列wamit文件. FAST基于BEM理论计算气动载荷，调用上述wamit文件计算水动力载荷，基于有限元法计算锚链力，从而实现针对该Spar浮式风力机的气-液-固-弹时域耦合数值模拟计算.

 

表4中，D_wind为风向，D_wave为浪向. 根据表4，定义4种类型工况（LC）：LC1为衰减实验工况；LC2只有波浪载荷，波浪周期为4~125.66 s；为了验证环境载荷方向改变时，设计出的浮式风力机能够满足使用有效性，取北海工作海况^[17]为LC3，其中包括只考虑风和风浪共存2种类型，考虑波浪作用时的波浪谱采用JONSWAP谱.

3 动态耦合响应数值模拟 3.1 固有周期

利用自由衰减法确定固有周期. 在测试中，风机处于停机状态且桨距角设置为90°，即风轮处于静止状态且风轮平面与总体坐标系的x轴平行. 表5给出该Spar浮式风力机各自由度的固有周期T_n.

3.2 不同浪向下平台的幅值响应算子（RAOs）

响应幅值算子表征浮式风力机系统的水动力性能. 环境条件LC2通过一系列规则波模拟，得到风力机总体运动响应的RAOs. 图7~10给出不同浪向下的各自由度RAOs仿真结果. 图中，T_m为周期，A_sur为纵荡幅值，A_hea为垂荡幅值，A_pit为纵摇幅值，A_Yaw为首摇幅值.

图7~10表明，在不同浪向情况下，同一自由度的RAOs具有相似的趋势，且各自由度RAOs的最大值均在该自由度对应的固有周期处. 对于纵荡、纵摇运动，各浪向均存在2个峰值，T_pit和T_sur分别为纵摇和纵荡的固有周期，表明纵荡和纵摇耦合作用较明显，且RAO随着浪向角的增大呈现下降趋势；对于垂荡运动，不同浪向对RAO的影响较小，在垂荡固有周期处存在明显的峰值，在纵荡固有周期处，运动幅值略微增大但不明显；对于首摇运动，当浪向为0°时，首摇运动几乎为0，其他各浪向均存在3个峰值，分别对应于首摇、纵摇和纵荡的固有周期，随着浪向角的增大，RAO在同一波浪周期处呈现上升趋势.

3.3 不同风浪向条件下平台运动响应研究

6 MW风机的额定工作风速为10.5 m/s，此时桨距角为0.7°；当风速增大时，桨距角增大. 在额定风速下，该风力机所受到的推力最大，不同风速下的推力如表6所示. 表中，v为风速，n为风轮转速，θ为桨距角，F_T为推力. 针对LC3工况中考虑波浪载荷时的计算结果表明，额定风速下平台的运动响应和锚链力最大，只给出10.5 m/s风速下的运动响应结果. 在工作海况下，风浪向的变化对平台运动响应产生较大的影响，共计算16个风浪共同作用的工作工况，如表7所示，并与4个只有风的工况进行比较.

3.3.1 统计结果分析

取不同风浪入射方向下Spar浮式风力机230~2 030 s时的运动结果进行统计分析，如图11~14所示. 如图11所示为x-y平面即水平面线运动响应统计值. 图中，A_h为水平运动幅值. 如图13所示为与z轴成一定夹角的角运动响应统计值. 图中，A_ang为角运动幅值.

结果表明，风浪同向时，Spar平台的水平位移幅值和角运动幅值最大，风浪夹角越大，水平位移、角运动越小；风向60°下，水平位移幅值较小，且各工况的水平运动平均值和最大值趋势相同；风浪向的变化对垂荡运动的最大值、均值及标准差影响都不大，可见风浪向的改变不会影响垂荡的运动响应；对于首摇运动，当风浪同方向时，平台首摇角最大值、均值以及标准差最小，当风浪夹角增大时，首摇角最大值、均值以及标准差都增大，这表明由于Spar平台的对称性，风浪同向不会产生较大的首摇响应，较大的风浪夹角会导致较大的首摇运动响应. 最大首摇角只有6.496°，这表明具有delta-line的6 MW Spar型风力机具有良好的首摇运动响应，避免锚链缠绕现象的发生. 对于各自由度，各风浪共同作用下的运动响应均值与只有稳定风作用时的运动幅值基本保持一致；可见，在额定工作条件下，平台运动响应主要取决于风载荷，波浪载荷会改变平台在平衡位置附近的振荡运动，即风引起定常效应，波浪会引起周期效应.

3.3.2 时域频域结果分析

纵摇和首摇运动响应的好坏对Spar浮式风力机的生存情况起关键作用. 取风向为0°，浪向为90°，Spar平台纵摇和首摇的时域结果如图15、16所示. 如图17、18所示为经过傅里叶变换后，风向为0°、浪向变化的纵摇和首摇功率谱S(ω). 图中，ω为角频率. 如图19、20所示为浪向为0°、风向变化的纵摇和首摇功率谱.

图17、18表明，对于纵摇，存在3个峰值点，分别对应风力机的纵荡、纵摇固有周期及波浪周期. 在纵荡固有周期处，峰较明显，说明纵摇和纵荡耦合作用较明显. 当风向不变，波浪入射角增大时，纵摇在波浪周期处的功率谱幅值明显下降，说明波浪方向对纵摇动态响应的影响较大；浪向不变，风浪同向时纵摇的功率谱幅值最大，但变化较小，可见风浪同向会产生最大的纵摇响应，但风向变化对纵摇的动态响应影响较小.

图19、20表明，对于首摇，仅存在一个峰值点，说明在工作海况下，首摇与其他自由度的耦合作用很小. 当风浪同向时，首摇在波浪周期处的功率谱幅值最低；当风浪夹角增大时，首摇在波浪周期处的功率谱幅值明显增大. 可见，当风浪夹角较大时，会激发大幅首摇动态响应，与时域计算的结果一致.

3.4 不同风浪向条件下系泊系统响应研究

仅讨论会出现最大张力的系缆3，图21给出系缆3导缆孔处的张力T₃在不同风浪方向条件下的统计结果，图22给出当风、浪入射角均为60°时的时历响应曲线，如图23、24所示为风向60°、浪向变化的张力频谱分析图和浪向60°、风向变化的张力频谱分析图.

图21表明，系缆3的张力在额定工况条件下的均值与只有稳定风作用时的幅值基本保持一致. 在额定工作条件下，系缆张力主要取决于风载荷，波浪载荷会改变系缆张力的动态响应.

图23、24表明，系缆3的张力存在2个峰值点，分别对应风力机的纵荡和纵摇固有周期. 在纵荡固有周期处，峰较明显，说明能量主要集中在低频运动. 当风浪向夹角均为30°时，张力功率谱的幅值变化不大；当风浪同向时，在波浪谱峰周期处具有最大的张力功率谱幅值；在不同的风浪入射角下，系缆的最大张力不超过1 800 kN，证明采用该系泊系统能够保证系统安全运行.

4 结　论

（1）Spar型浮式风力机纵荡和纵摇耦合作用较明显；当风浪同向时，会激发Spar型浮式风力机在该方向的大幅动态响应.

（2）风浪方向夹角越大，平台首摇运动响应越大；环境载荷方向对垂荡运动的影响不大.

（3）设计出的6 MW Spar浮式风力机在风浪的各种入射角下具有良好的运动响应，delta-line的存在能够避免首摇角过大而导致锚链互相缠绕的现象发生，系泊系统完全能够保证整个风力机的安全运行.