海上风力发电机组环境载荷比陆上机组更加复杂，设计前期必须详细评估海洋环境参数。海冰会对风力机性能产生重大的影响，因此，包含冰荷载的计算工况成为风机设计的重要控制工况，而冰对塔筒负载将是一个重要的考虑因素。

国内外学者对风力机动力学特性进行了广泛研究。Yao等^[1]根据疲劳载荷进行了塔筒疲劳分析，采用定量可靠性理论计算塔筒平均故障时间；Mitsumasa等^[2]对比了不同控制策略对海上风力机的影响；金鑫等^[3-4]基于动力学模型研究了独立变桨控制策略对风电机组的影响; Bazeos等^[5]把风力机塔筒看作悬臂梁并进行了地震波作用下的动力学分析及稳定性研究；Murtagh等^[6]研究了风力机基础在风载荷作用下的动力反应以及叶片刚化效应的影响；Arka-diusa等^[7]研究了风力机在冰载作用下的整体动力响应；何玉林^[8]证明了风机叶片覆冰时，叶片气动性能将会受到影响；刘强等^[9]分析了风浪作用对机组的运动特性和发电性能的影响；岳前进等^[10]对锥体结构上的冰力作用形式进行观测分析，提出了随机冰力时程的概念。

笔者运用多体动力学方法建立了详细的风力发电机动力学模型，并借鉴随机冰力函数模型，将冰载荷添加到风力发电机模型中，建立了海冰作用下的风力机动力学模型并进行仿真，仿真结论对后期风电机组塔基载荷优化控制研究具有重要意义。

1 海上风力机塔筒受载描述

海上风电机组的外界载荷比较复杂，与陆上风机最大的区别是要考虑海上载荷条件。海上风力机塔筒可能受到的载荷有：风力发电机组的重力、气动载荷、波浪载荷、冰载荷等。

2 海上风力机多柔体动力学模型

在柔性多体系统动力学建模中，柔性体的位移通过使用动参考系和弹性变形两组广义坐标来描述。因此，基于拉格朗日方程的柔性多体系统控制方程为

$ \left\{ \begin{array}{l} \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial L}}{{\partial \xi }}} \right)-\frac{{\partial L}}{{\partial \xi }} + \frac{{\partial F}}{{\partial \xi }} + {\left( {\frac{{\partial C}}{{\partial \xi }}} \right)^{\rm{T}}}\lambda = Q, \\ C\left( \xi \right) = 0; \end{array} \right. $

(1)

式中：L为拉格朗日量；F为能量耗散方程；λ为拉格朗日算子；Q为广义作用力；C(ξ, t)为约束方程。

柔性体叶片采用铁木辛柯梁建模并综合考虑了叶片的刚度硬汉效应。根据叶片长度将整个叶片分成若干段，每段采用柔性梁建模，通过每段单元的预弯、面内外刚度、扭角等参数生成叶片模型。

轮毂、主轴、主机架等复杂柔性体零部件的建模思路为：根据部件三维模型划分网格，采用模态综合法进行子结构的缩减，完成模态计算，生成符合要求的输入文件，输入文件包含了有限元模型节点单元材料等信息和子结构分析后生成的质量矩阵、刚度矩阵等信息，最后将输入文件导入多体动力学软件中生成柔性体文件，完成柔性体的建模。

齿轮建模通过参数化方式进行，齿轮轴考虑为柔性体轴建模。齿轮啮合力元根据DIN3990^[11]标准完成模拟，该标准可以计算齿轮啮合的切向力、径向力和轴向力等参数。发电机建模类型为刚体，发电机箱体通过力元模拟4点弹性支撑。塔筒建模过程直接视为梁单元处理，通过建立一系列标记点，添加各截面的属性和参数，自动生成柔性体文件即可。

依据风力发电机认证规范要求，表 1所示为风力机系统各主要零部件刚柔属性要求及建模自由度。

图 1所示为风电机组整机多柔体动力学系统拓扑图，包含了风力发电机组的主要零部件，零部件之间载荷传递路线以及载荷传递实现方法等。风力机通过气动力元添加风载荷驱动，发电机定子和转子之间的反馈力矩由力元实现。控制系统主要由变桨系统、偏航系统和发电机系统组成。

风电机组由于自身机械结构强度的限制，当风速超过额定风速后，需要进行变桨距控制来降低风能利用率，使机组稳定在正常工作范围之内，不至于使风力发电机过载运行，以此保护发电机组。变桨距控制策略是通过将发电机转速反馈信号作为控制系统的输入变量来计算桨距角实现变桨控制。图 2所示为在Simulink中搭建的控制系统模型并与多体动力学模型模块实现了接口连接。

3 海洋环境载荷模型 3.1 波浪载荷模型

利用线性波浪理论求解水质点速度和加速度时，水平方向水质点的速度v_m为

$ {v_{\rm{m}}} = {A_{\rm{m}}}\omega \frac{{\cosh \left[{k\left( {z + h} \right)} \right]}}{{\sin \left( {kh} \right)}}\sin \left( {\omega t} \right), $

(2)

式中：A_m为波浪幅值；ω为角速度；h为海水深度；k为波浪数目；z为水面高度；t为波浪周期。

波浪大小通常由波浪力水平分量决定，对于小直径圆柱塔筒，柱体对波浪运动无显著影响，常采用莫里森(Morison)公式^[12]计算。Morison方程中，作用在柱体上的水平力是水平速度和加速度的函数，作用在风力机塔筒的水平波浪力包括惯性力和阻力两部分，则塔筒任意高度处的水平波浪力为

$ {\rm{d}}F = {C_{\rm{m}}}{\rho _{\rm{w}}}\frac{{{\rm{\pi }}{D^{2 \cdot }}}}{{\rm{4}}}v{\rm{d}}z + {C_{\rm{d}}}{\rho _{\rm{w}}}\frac{D}{2}\left| {{v_{\rm{m}}}} \right|{v_{\rm{m}}}{\rm{d}}z, $

(3)

式中：C_m为惯性力系数；C_d为阻力系数；D为波浪作用处直径；ρ_w为海水密度。

3.2 冰载荷模型

在结冰环境中，作用于海上结构物的冰载荷远大于其他载荷，因此，可将冰载荷当作结构物的控制荷载。根据海冰作用于结构的形式不同，海冰可能发生挤压、弯曲和屈曲形式破坏。当冰排作用于垂直结构时，冰以挤压的方式破碎；当冰排与锥体结构作用时，冰排主要发生弯曲破碎，使得冰排对结构的作用比挤压冰破坏作用要小得多。

海冰以弯曲这种形式的破坏，研究者对原型冰锥进行了大量测量，并分析锥体与冰作用过程中冰力大小变化及冰排破坏形式后得出结论：冰与锥体作用的冰力时间历程实际上是若干个脉冲函数，进一步对冰力时程做简化得到了确定性冰力函数。后来发现，冰与锥体结构实际破坏过程是一个随机过程，振动也是随机振动。Qu Yan等^[13]对实测的多条冰力时程分析中引入随机冰力函数的概念，并统计发现幅值和周期服从正态分布，将随机过程分为若干周期的组合，并给出了任意周期i内的冰力函数表达式：

$ {f_i}\left( t \right) = \left\{ \begin{array}{l} \frac{{6F_i^m}}{{{T_i}}}, 0 < t < \frac{{{T_i}}}{6}, \\ 2F_i^m -\frac{{6F_i^m}}{{{T_i}}}t, \frac{{{T_i}}}{6} < t < \frac{{{T_i}}}{3}, \\ 0, \frac{{{T_i}}}{3} < t < {T_i}; \end{array} \right. $

(4)

式中：F_i^m为第i个周期内的冰力幅值；T_i第i个冰力周期。

一般认为作用在锥体上的最大静冰力即为动冰力幅值的大小，文献[14]对锥体的冰力模型进行了总结：锥体水平静冰力主要由冰弯曲强度和推动碎冰上爬力两部分决定，对于窄锥体，一般不考虑上爬力，取Hirayama-Obara静冰力计算公式^[15]，冰力幅值为

$ {F^m} = K{\sigma _{\rm{f}}}{h^2}{\left( {\frac{d}{{{L_{\rm{b}}}}}} \right)^{0.34}}, $

(5)

式中：K为待定系数；σ_f为冰的弯曲强度；h为冰厚；d为冰作用处锥体直径；L_b为冰的破碎长度。

冰排弯曲破坏时，冰力周期由破碎长度和冰速决定，此时冰力周期表达式为

$ T = \frac{{{L_{\rm{b}}}}}{{{v_{{\rm{ice}}}}}}, $

(6)

式中：v_ice为冰速，短时间内为恒定值。因此，冰力周期主要由破碎长度决定，破碎长度是锥体几何形状、冰厚等参数的函数，当结构尺寸一定时，破碎长度主要由冰厚决定。根据冰力时程曲线统计结果，冰与锥体结构作用的时间大约占整个周期的1/3。

4 仿真实例分析

采用多体动力学软件进行风力发电机系统建模。使用NREL提供的5 MW近海水平轴单桩固定式风力机为研究对象进行仿真分析，塔架结构为圆筒型。

根据风力机认证规范IEC61400-3设计典型载荷工况，如表 2所示。波浪载荷通过水动力力元添加，完成多体动力学与水动力学耦合；海冰载荷通过力元实现，将载荷时间历程输入到力元中。把波浪载荷和海冰载荷作为模型输入控制量，图 3为风力发电机仿真过程图。冰厚为正常冰厚，取8 cm；参考文献[16]，冰速取值为0.3 m/s，图 4为在Matlab中生成的冰载荷时间历程曲线。

4.1 风力发电机固有频率验证

整机模态分析过程按质量分布和刚度参数进行，由于风力机主要工况处于低频率状态，对风力发电机影响较大的是前几阶模态，表 3分别是整机模型在SIMPACK和Bladed中计算的前3阶固有频率，由表可以看出两者计算结果很接近，相对误差最大为2.21%，小于5%。根据GL2010规范要求，证明该模型合理。

4.2 塔筒顶部动力特性响应

图 5分别为风载荷、波浪载荷和海冰载荷单独作用下塔顶位移响应曲线。由图 5可知，风荷载作用下塔顶位移均值为0.21 m，开始阶段波动幅度比较大，运行稳定后，波动幅值逐渐变小，这与实际情况吻合。波浪载荷和海冰载荷作用下塔顶的位移幅值分别为0.047 0 m和0.003 9 m，两者差别较大，但均远小于风荷载作用对顶部位移的影响。因此, 风波、风冰联合作用时，由于风力作用对顶部位移影响更大，因此，计算塔顶位移主要由风载荷决定, 波浪和海冰作用可忽略不计。

4.3 塔基载荷计算比较

将无海冰情况风力发电机的计算结果纳入比较，分析塔基载荷变化情况。图 6(a)-(f)分别为塔基受力及力矩的时域响应。在海冰作用下，塔基各个方向的振动特性均有明显的提高，其中F_X的幅值增大最多，M_Z的载荷波动变化最为显著，分析可知载荷波动突变的原因是由于海冰对塔基的正面持续撞击作用所导致，撞击作用会使塔基在机舱前后、垂直方向的力和机舱左右方向的力矩普遍大于无海冰作用时的塔筒受载。

塔基等效疲劳载荷结果如图 7所示。由图可知，联合作用下的疲劳载荷介于各自单独作用下的疲劳载荷之间，因此，联合作用在一定程度上可降低疲劳载荷。分析图 7可知，风冰联合作用下的疲劳载荷大于风波联合作用下的疲劳载荷，说明海冰对风力机塔筒疲劳载荷的影响大于波浪的影响。单独冰载荷作用对塔基的疲劳载荷最大，海冰联合作用下载荷也比较大，可见海冰对塔筒的疲劳影响不容忽视。

5 结语

将波浪载荷和冰载荷施加到风力机多体动力学模型中，建立海上风力机多体动力学仿真模型，并对部分典型工况进行仿真对比，结果表明：波浪和海冰对风力机塔顶的位移波动相对于风的影响比较小，一般忽略不计，塔顶位移主要受风载作用；但海冰对风电机组塔基载荷性能的影响显著，海冰的撞击作用会使塔基在机舱前后、垂直方向的力和机舱左右方向的力矩普遍大于无海冰作用时的塔筒受载。本文结论将为风力机塔架结构的改进和载荷优化提供技术参考，尤其是在破冰期，海冰对塔筒基础的破坏作用更为显著。