随着陆上风电场可开发资源的减少，海上风电开发已成必然。海上风力机与陆上风力机工作环境迥然不同，其主要区别在于，海上风力机除了气动载荷之外，还需要进一步考虑波浪等海洋环境的影响^[1-3]。漂浮平台可分为单柱平台、张力腿平台、驳船平台以及三者的组合^[4-5]。本文选取张力腿平台进行研究。

气动计算方法包括叶素动量理论、涡尾迹方法和计算流体力学方法(Computational fluid dynamics，CFD)。叶素动量理论简单快捷，工程应用广泛，但准确度并不高。CFD方法^[6]虽然能够模拟叶片上的复杂流动，但是计算效率低，消耗资源巨大，且受如网格划分、湍流模拟和转捩模拟等因素影响，限制了其在风力机计算中的应用。涡尾迹方法本质上具有旋涡特性，计算成本与CFD方法相比更小，比叶素动量理论更精确，且物理上比叶素动量理论更准确地考虑涡尾流的诱导作用，是模拟风力机较为灵活的数值工具。预定涡尾迹^[7-8]是根据大量的尾流场实验数据，建立尾迹形状的经验描述函数。与其相比，自由涡尾迹^[9-12]不需要涡元位置的实验数据，尾迹涡元允许在当地速度场的影响下自由变形，能够计算尾迹的畸变和卷起。对于气动、结构和水动的耦合计算，以及风力机的复杂运动和动态响应，自由涡尾迹有着更大的优势。

波浪模拟方法^[13]主要包括线性波理论、Stokes波理论和流函数方法。线性波理论假定波浪振幅与波长相比很小，水质点运动速度缓慢，因此可对势函数、边界条件以及连续条件作线性化处理，即对波浪运动一阶近似。Stokes波理论则保留波浪运动的三阶或者五阶项。流函数理论采用流函数描述波浪运动。通常采用频谱函数的方法，通过不同相位角的线性波叠加描述真实海况。得到波浪运动特性之后，采用Morison方程^[14]对波浪载荷进行计算。

本文构建张力腿漂浮平台动力学方程，以气动载荷和波浪载荷为激励，计算其纵荡、纵摇、横荡和横摇4个自由度响应，并反馈至气动计算，从而完成气动、结构和水动的耦合计算。最后，以大型风力机NH1500和张力腿漂浮平台为例，对漂浮式风力机进行仿真模拟，对比分析来流风速、浪高和波浪入流角对漂浮平台运动和风力机气动性能的影响，采用PM频谱模拟真实海况并计算漂浮式风力机动态响应，为漂浮式风力机设计提供了依据。

1 气动载荷 1.1 自由涡尾迹方法 1.1.1 叶片升力面

在Weissinger-L升力面模型中，叶片的升力特性用一条附着涡线模拟，附着涡线位于1/4弦线。叶片沿展向分成若干段叶素，每段叶素对应一个控制点，位于叶素中间剖面的3/4弦长处。叶片离散示意图如图 1所示。

1.1.2 Biot-Savart定律

涡元求解分析示意图如图 2所示。由Biot-Sarvart定律确定的涡元对点P的诱导速度为

$dV = {\Gamma \over {4\pi }}{{dl \times r} \over {{r^3}}}$

(1)

式中:V为诱导速度；Γ为环量；dl为涡元长度；r为点到涡元距离。积分可得

$V = {\Gamma \over {4\pi h}}\left( {\cos {\theta _A} - \cos {\theta _B}} \right)$

(2)

1.1.3 尾迹控制方程

涡线示意图如图 3所示。涡线的控制方程为

${{dr\left( {\psi ,\zeta } \right)} \over {dt}} = {V_{loc}}\left( {r\left( {\psi ,\zeta } \right),t} \right)$

(3)

式中:ψ为叶片方位角；ζ为尾迹寿命角；V_loc为当地速度，包括自由来流速度V₀和诱导速度V_ind，即

${V_{loc}}\left( {r\left( {\psi ,\zeta } \right),t} \right) = {V_0} + {V_{ind}}\left( {r\left( {\psi ,\zeta } \right),t} \right)$

(4)

本文采用时间步进法对涡线形状进行迭代求解，并使用预估-校正方法^[11]提高准确性。

1.2 塔架载荷

对于上风向风力机，塔架所受气动载荷影响较小，因此仅用阻力计算公式对其做简单计算

${F_d} = {1 \over 2}{C_d}\rho {v^2}S$

(5)

式中:C_d为阻力系数，本文取0.6；ρ为空气密度；v为来流速度；S为塔架横截面积。

不同高度处的风速采用风剪切模型计算

$v\left( h \right) = v\left( {{h_0}} \right){\left( {{h \over {{h_0}}}} \right)^\alpha }$

(6)

式中:v(h₀)为参考高度h₀处的风速；α为风剪切系数，与地形地貌相关，本文取0.12。

2 波浪载荷 2.1 波浪模拟

本文选取PM能量频谱进行波浪模拟，PM谱的峰频为

${\omega _p} = 0.877g/U$

(7)

$S\left( \omega \right) = \alpha {{{g^2}} \over {{\omega ^5}}}\exp \left\{ { - 1.25{{\left( {{{{\omega _p}} \over \omega }} \right)}^4}} \right\}$

(8)

式中:α=8.1×10^－3；g为重力加速度，取9.8 m/s²；U为海面以上19.5 m高处的风速，轮毂风速为10 m/s时，通过上述风剪切模型计算得到U=8.6 m/s。

对波浪频谱PM谱进行频率分割离散求出单个波的波幅为

${A_n} = \matrix{ {\sqrt {2S({\omega _n})\Delta \omega } } & {{\omega _n} = n\Delta \omega } \cr } $

(9)

浪高函数为

$\eta \left( {x,y,t} \right) = \sum\limits_{n = 1}^m {{A_n}} cos{\varphi _n}$

(10)

${\varphi _n} = - {\omega _n}t + {\theta _n}$

(11)

式中：θ_n为随机相位，在[0,2π]内均匀分布，以确保η(x,y,t)的随机性。确定叠加谐波阶数m的标准是随机波的谱密度函数S_ηη(ω)的积分面积至少要达到整个波普面积的95%；谐波数 k_n可以由频散方程式求得

$\left\{ {\matrix{ {\omega _n^2 = g{k_n}\tan h{k_n}d} & {浅水区\left( {d < {\lambda \over 2}} \right)} \cr {\omega _n^2 = g{k_n}} & {深水区\left( {d > {\lambda \over 2}} \right)} \cr } } \right.$

(12)

PM频谱与浪高如图 4所示。

得到浪高函数后，可参考文献^[14]得到水质点运动速度和加速度。

2.2 波浪载荷

波浪力的求解应用Morison方程。在深度z，水平波浪力由两部分构成：阻力和惯性力。波浪入流角定义为波浪方向与来流风向的夹角。作用在结构上的波浪载荷表达式为

$\matrix{ {{F_{wx}} = {C_m}\rho \pi {D^2}4\left( {{a_{wx}} - {a_{sx}}} \right) + \rho A{a_{wx}} + } \hfill \cr {{C_d}\rho D4\left( {{u_{wx}} - {u_{sx}}} \right)\left| {{u_{wx}} - {u_{sx}}} \right|} \hfill \cr } $

(13)

$\matrix{ {{F_{wy}} = {C_m}\rho \pi {D^2}4({a_{wy}} - {a_{sy}}) + \rho A{a_{wy}} + } \hfill \cr {{C_d}\rho D4\left( {{u_{wy}} - {u_{sy}}} \right)\left| {{u_{wy}} - {u_{sy}}} \right|} \hfill \cr } $

(14)

式中:ρ为海水密度；D为垂直圆筒塔架的直径；a_wx,a_sx分别为水平方向上水质点和塔架的加速度；A为垂直圆筒塔架的横截面积；u_wx,u_sx分别为水平方向上水质点和塔架的速度;C_m为质量系数(或惯性系数)，与被绕流物体的形状、绕流的流态以及物体剖面大小与波长的比值等因素有关。对于圆柱，当L$ \gg $D时，C_m可取值1.0。关于阻力系数C_d,除在个别情况下有理论解答外，大多数情况下要依据试验确定。对于圆柱，当L$ \gg $D时，如果考虑疲劳载荷，C_d=1.25。

3 张力腿平台动力学方程

张力腿平台分为纵荡、纵摇、横荡、横摇、垂荡和艏摇6个自由度，如图 5所示。由于垂荡和艏摇方向载荷和响应较小，因此只考虑前4个自由度。动力学方程为

$M\ddot x + C\dot x + Kx = F$

(15)

式中:M为质量矩阵;C为阻尼矩阵;K为刚度矩阵，具体形式可参考文献^[13];F为外载荷，包括气动载荷与波浪载荷。首先通过气动与波浪载荷计算得到外部激励，然后通过Newmark方法进行时间步长推进，得到张力腿平台4个自由度响应，并将其反馈至气动和波浪载荷计算，从而完成气动、结构和水动的耦合计算。

4 方法验证与算例计算分析 4.1 计算方法验证

对大型风力机NH1500的功率系数C_p进行计算，并与实验数据对比，如图 6所示。从图 6可以看出，自由涡尾迹方法计算所得功率系数与实验值相比偏大，是因为自由涡尾迹方法只在涡核处考虑了粘性，而没有考虑流场中的粘性。总体来说，自由涡尾迹方法计算结果与实验值比较吻合，验证了该方法的准确性与可靠性。NH1500叶片数目为3，叶片半径为40.5 m，额定转速为17.2 r/min，额定功率为1 500 kW，额定风速为10.4 m/s，功率控制为变桨变速，塔架高度为63.0 m。

4.2 算例结果与分析

张力腿平台参数如表 1所示。首先研究了风速大小对漂浮式风力机的影响，选取单一线性波模拟海况，波浪入流角为0°，浪高2 m，波浪周期10 s，来流风速取0,6,8,10和12 m/s。风速为0 m/s时，即考察波浪载荷对漂浮平台运动的影响。漂浮平台纵荡和纵摇响应，以及风速10 m/s时的波浪载荷与风轮推力如图 7所示。

通过对比不同风速时漂浮平台纵荡运动，可以看出，风速大小对平台纵荡运动的影响并不显著，波浪载荷是平台往复运动的主要因素。对于纵摇运动，风速大小的影响则比较明显，这是由于风轮推力的力臂较长，产生的力矩较大，从而放大了风轮推力的影响。通过图 7(c)可以看出，风轮推力相对于波浪载荷较小，这与上述结果和分析是一致的。

然后计算不同浪高对平台运动和风力机气动性能的影响。来流风速为10 m/s，选取单一线性波模拟海况，波浪入流角为0°，波浪周期10 s，浪高为0,1,2,3 m。浪高为0 m时，即考察气动载荷对平台运动的影响。计算结果如图 8所示。

通过图 8可以看出，浪高越大，水质点动能越大，波浪载荷也就越大，因此平台运动幅度也就越大，对功率系数造成的波动也就越大。但无论是平台运动，还是风力机功率系数响应，其平衡位置都是与无波浪入流时的响应基本一致。并且没有波浪入流时，只有气动载荷作用下的平台运动幅值较小，功率也比较稳定，因此波浪载荷是造成响应波动的主要原因。

再然后计算波浪入流方向对平台运动和风力机气动性能的影响。来流风速为10 m/s，选取单一线性波模拟海况，浪高2 m，波浪周期10 s，波浪入流角分别为0°,30°,60°和90°。计算结果如图 9所示。

随着波浪入流角从0°增大到90°，波浪载荷在纵向的分量逐渐减小，横向分量逐渐增大。因此漂浮平台的纵荡、纵摇逐渐较小，而横荡、横摇则逐渐增强。平台纵摇的平衡位置约为0.008°，而横摇的平衡位置为0°左右，是由于风轮推力方向与平台纵向运动一致导致的，由于波浪载荷较大，风轮推力在纵荡横荡中作用不明显。风力机功率系数平衡位置基本不变，波动幅值逐渐减小，说明平台纵向运动的减小，导致功率系数波动的减小。虽然入流角90°时的横向运动和0°时的纵向运动幅值基本相等，但是横向运动造成的功率系数波动明显不如纵向运动影响大，这同样是由于平台纵向与来流风向一致，其波动对风力机气动性能的影响更为直接显著。

最后应用PM频谱模拟波浪运动并进行计算，来流风速10 m/s，波浪入流角45°，波浪特性由频谱函数模拟得到。风力机功率系数、平台运动和载荷如图 10所示。

由于波浪模拟采取随机相位角，波浪载荷呈不规则波动。PM频谱选取的参考风速较小(8.6 m/s)，因此由图 4可以看出，各频率对应的浪高较小，导致与风轮推力相比，波浪载荷的幅值较小，所以此时平台运动的主要影响因素为风轮推力，但是波浪载荷的波动也对平台运动产生一定的影响，进一步导致风力机气动性能有一定波动。平台纵摇和横摇的幅值较小，纵荡和横荡的运动比较明显。由于风轮推力方向和平台纵向一致，导致纵荡方向平衡位置明显高于横荡方向，但是二者的波动幅值接近，说明波浪载荷是平台往复运动的主要原因，气动载荷则影响平台运动的平衡位置。同样，纵摇角度平衡位置明显大于横摇角度，而二者波动幅值较为接近。纵荡和横荡运动幅值均在0.5 m以内，纵摇和横摇运动幅值均在0.1°以内，表明了本文张力腿平台具有良好的稳定性。尽管平台运动幅度较小，但仍对风力机功率系数产生较大影响，因此对于漂浮式风力机，不仅需要严格控制平台运动，更需要针对其建立灵活的控制策略。

5 结论

本文采用自由涡尾迹方法对风力机气动性能进行了计算，并通过算例验证了该方法的准确性和可靠性。对漂浮式风力机气动性能和平台响应进行了多种工况计算分析，得到了以下结论：

(1) 平台运动受波浪载荷和气动载荷的综合作用，深海浪高超过2 m时，波浪载荷起主要作用。风轮推力决定平台运动的平衡位置，波浪载荷决定平台运动的幅值。

(2) 浪高越大，造成的平台运动和功率波动幅值越大，但平衡位置基本不变。

(3) 风力机的气动性能受平台纵向运动的影响较大，横向运动产生的影响较小。

(4) 波浪载荷对风力机气动性能影响显著，因此不仅需要在结构上减小漂浮平台运动，更要建立灵活的控制策略及时调整风力机运行状态。

本文方法可以较好地模拟漂浮式风力机气动性能和平台动态响应，为漂浮式风力机设计提供依据。