基于整机刚柔耦合的风力机叶片气弹稳定性分析

郭小锋，黄鑫祥，车江宁，张 忍

（中原工学院 机电学院，河南 郑州450007）

0 引言

为了提高风轮对风能的吸纳能力，降低发电成本，现代风力机叶片的长度、重量和极限载荷较以往有明显的增加。早期的2MW Tjæreborg风力机叶片长度仅为29.1 m，重量为6.8 t[1]。现代的2 MW风力机叶片长度一般在50～60m，重量为10～15 t[2]。风轮叶片重量及长度的增加使叶片的气弹扭转变形分析及气弹失稳问题成为风力机设计中的关键。

对于风力机的气弹动力学问题，传统的分析方法包含模态分析法、有限元法（FEM）和多体动力学方法[3]～[5]。模态分析的基本原理是将塔筒和风轮叶片的变形描述为若干形函数的线性组合，然后采用虚功法建立整个风力机的动力学方程。有限元方法一般采用Timoshenko梁或非线性梁理论将叶片离散成多个单元，根据结构动力学理论和哈密顿原理建立叶片的动力学方程。多体动力学方法是将风力机看作由若干刚性或柔性部件的组合，这些部件通过力元素（如弹簧）或运动学约束（如关节）相互连接，然后采用拉格朗日方程或牛顿-欧拉方程构建整机的动力学模型。

风力机的来流风况具有复杂多变性，风轮叶片在气动力、弹性力及变速变桨控制的耦合作用下呈现出复杂的动态响应。现有的分析方法和工具主要侧重于叶片变形及载荷的研究，针对当前风力机叶片大型化、柔性化的发展趋势，研究风力机叶片在复杂工况下气弹稳定性分析的新方法，具有重要的意义。

1 风力机多体动力学模型

在工作状态下，风力机叶片和塔筒存在明显的弹性变形，所以叶片和塔筒需作为柔性体处理，而机舱、主轴及轮毂等部件在工作时的弹性变形很小，可以看作为刚体。因此，风力机整机模型是典型的刚柔耦合多体系统[6]。

1.1 柔性体扭转变形的表达方法

基于多体系统（MBS）理论，柔性体的叶片被看作是若干个截面的组合，在叶片的某截面A处建立当地局部参考坐标系Oc-XcYcZc，当叶片发生弹性变形时，截面A相对于其初始状态会发生三维空间上的弹性移动和扭转变形，在某一时刻其位置和姿态如截面A＇所示，截面A上的参考坐标系Oc-XcYcZc转变为Oc＇-Xc＇Yc＇Zc＇（图1）。

图1 柔性叶片气弹变形示意图Fig.1 Generalmodel of a flexible body

截面A＇相对于A的弹性扭转变形可以用三维向量θc表示，并采用模态叠加法求得。

在叶片的气弹变形分析中，截面的气弹扭转变形是通过欧拉变换进行表达的。考虑截面A相对于叶根坐标系存在气动扭角β的变换，因此，将O-XYZ坐标系平移至Oc点，绕Z轴旋转β，即得到Oc-XcYcZc坐标系。同样，将Oc-XcYcZc坐标系平移至Oc＇点，先绕Zc轴旋转ψ，再绕新的Yc轴旋转θ，最后绕新的Zc轴旋转φ，得到Oc＇-Xc＇Yc＇Zc＇。

综上所示，ψ+β为截面A相对于叶根坐标系的角扭转角位移，θ为截面A在挥舞方向的扭转角，而φ为截面A在摆振方向的扭转角。

截面A＇相对于A的弹性位移变形可以用三维向量uc表示，同样采用模态叠加法求得[7]。

因此，截面A在某一时刻下的位置和形态为

1.2 整机动力学模型

本文将风力机模型进行合理简化，只包括叶片、轮毂、塔筒、机舱等主要部件，而不考虑齿轮箱、发电机等部件。偏航法兰与塔筒顶端固接，用于连接机舱，虽然机舱中不创建主机架、齿轮箱、发电机等模型，但是，需要把忽略部件的质量和转动惯量算入机舱中。低速轴与机舱为转动连接，轮毂与低速轴进行固接，叶片与轮毂由变桨法兰连接。所构建的风力机整机模型如图2所示。

图2 风力机整机动力学模型Fig.2 Dynamicmodel ofwind turbine

图2中：塔筒底部定义整机基础坐标系eb；机舱质心处定义en坐标系；主轴质心和轮毂中心分别定义了es和eh坐标系；在叶片根部定义了ebi坐标系。

在完成整机坐标系的构建后，风力机多体动力学系统的自由度参数可表述为

建立好风力机整机的刚柔耦合多体动力学模型后，采用凯恩方法对风力机动力学方程进行求解。

2 风力机多体动力学仿真方法

本文的研究对象为某内陆商用2MW风力机，叶片长度为54m，塔筒高度为85m，额定风速为9.5 m/s。基于多体动力学软件SIMPACK，本文所构建的风力机气-弹-控联合仿真模型及求解流程如图3所示。

图3 风力机气弹动力学模型与分析流程图Fig.3 Aerodynamicmodel and analysis flow chartof wind turbine

2.1 叶片参数化表达

首先采用EXCEL表格对叶片的外形、铺层结构进行参数化设计；其次编写MATLAB程序读取处理EXCEL表格中的叶片参数数据，调用PreComp软件计算得到SIMPACK叶片建模中所需叶片各截面的属性参数（如截面质心、弹性中心、剪切中心、扭转刚度、线密度等）；最后对上一步输出的数据进行整理，生成Rotorblade Generation模块的输入文件，完成风轮叶片参数化模型的构建。

2.2 塔筒部件

本文采用SIMPACK中的Flexible（model）模块进行柔性塔筒建模。在已知塔筒内外径、材料等参数的情况下，借助ANSYS自身的参数化语言APDL实现塔筒有限元模型的构建。由于塔筒有限元模型的网格规模庞大、自由度多，为了节省SIMPACK多柔体数值计算的仿真时间，须要对有限元塔筒模型进行模型缩减，即对模型进行子结构分析。最后由ANSYS计算输出包含有限元塔筒模型节点、单元、材料等信息的.cdb模型文件以及包含塔筒缩减模型质量矩阵和刚度矩阵等信息的.sub文件。Flexible模块调用这些文件，最后生成塔筒模型。

2.3 空气动力学模块

为了实现对风力机在复杂工况下的气弹动力学分析，本文采用SIMPACK与AeroDyn联合仿真的方法。由ICwind或Turbsim[8]软件生成相应的风模型文件，SIMPACK多体动力学求解模块通过241号力元接口调用AeroDyn计算当前时刻的气动数据，根据气动数据对风力机进行动力学方程求解，计算得到各部件的变形量和速度变化。

2.4 控制模块

风电机组采用变速变桨控制策略，即低于额定风速时采用转矩控制方式，以发电机转矩需求为输入，以发电机测量转速为输出建立风力发电机动态数学模型。高于额定风速时采用桨距控制方式，以叶片桨距需求为输入，以发电机测量转速为输出建立风力发电机动态数学模型。本文采用Fortran基于PI控制原理编写风力机转速及桨距角控制的动态链接库文件[9]。

3 结果对比分析

平均风速为额定风速的湍流风况是风力机极为恶劣的工况，在这种工况下风力机叶片的叶根载荷和叶尖变形量最大，叶片也容易产生气弹稳定性[7]～[9]。因此本文选取平均风速为9.5 m/s的湍流风对叶片进行动态载荷分析，湍流风时间为600 s，其湍流强度纵向为20.61%、横向为17.56 %、垂直方向为11.04 %。SIMPACK的仿真时间设为100 s，选择SODASRT_2求解器对风力机模型进行离线积分求解，求解步长设为0.01 s。选取50～100 s时间段进行数据对比分析，该时间段内轮毂高度处的风速变化如图4所示。

图4 轮毂高度处的风速变化Fig.4Wind speed variation athub height

风力机开始运行时处于不稳定阶段，所以选择图4时间段内的分析结果与行业普遍认可的GH-Bladed软件的分析结果进行对比分析，以验证本文方法的可靠性。

在挥舞、摆振及扭转方向，2MW叶片叶尖截面（距叶根53m处）变形的分析结果对比见图5。

图5 湍流风下叶尖动态气弹变形Fig.5 Dynamic aerodynamic deformation of blade tip under turbulentwind

由图5（a）可知：SIMPACK计算的叶尖平均挥舞位移为5.1 m，最大挥舞位移为7.4 m，在该时间段内，叶尖挥舞位移的变化趋势与图4中X方向的风速变化趋势一致，说明叶尖挥舞变形主要受X方向的风速影响；GH-Bladed计算的叶片平均挥舞位移为5.3 m，最大挥舞位移为5.8 m。两种软件计算的挥舞位移变化趋势一致，且相差不大。

由图5（b）可知，叶尖摆振变形范围在-1～0.8 m，两种软件计算的摆振位移很接近。

图5（c）为叶尖在扭转方向的变形量，其中叶片扭转方向的变形量是叶片在工作中各截面变形后相对于叶根坐标系欧拉变换角φ，SIMPACK的风机模块可以直接输出这个计算结果，从而可以非常直接地评估叶片各截面的气弹扭转变形。通过分析发现，在风力机运行稳定后，叶片截面扭转变形通常维持在某一角度，并随着风速、风向的变化上下波动。由图5（c）可知：叶尖截面的平均扭转角度为3.5 °，且在50～100 s时间段内扭转角度波动幅度为0～7°；而GH-Bladed计算的平均扭转角度为3.3 °，在该时间段内扭转角度的波动幅度也为0～7°。

4 气弹稳定性分析

在实际工作中，大型风力机叶片截面会产生较大的弯扭弹性变形，特别是叶尖部位的截面。叶片截面的弹性变形会引起截面翼型攻角的变化，翼型攻角的变化使得叶片截面的气弹载荷发生改变，进而对其弹性变形产生影响。在极端工况下，这种气弹耦合响应会使叶尖部位产生失速或气弹失稳现象，因此对叶尖翼型的气弹稳定性分析对大型风电叶片的设计具有重要的意义。

叶尖失速是叶片失稳的重要特征之一。对于叶尖翼型，首先通过空气动力学软件AeroDyn分析叶片在旋转过程中的攻角变化情况，并计算攻角对应的升力系数；最后，定义1 s为一个周期来研究叶尖翼型的攻角-升力系数变化情况，在一个周期内，同一攻角可对应多个升力系数值，取这些升力系数中的最大值和最小值，这两个值之间的差值越大，说明叶片在这个周期内经历了较大的气动载荷变化。

对于本文研究的2MW风力机叶片，在前文所述的湍流风况下，叶尖截面翼型（距叶根53m处）在控制系统正常和控制系统失效两种工况下的攻角-升力系数关系如图6所示。

图6 叶尖截面翼型的攻角—升力系数曲线图Fig.6 Angle of attack—lift coefficient curve of blade tip section

由图6可知：在控制系统正常的情况下，叶尖截面在同一攻角下对应的升力系数比较接近，升力系数之差的最大值为0.046；而在控制系统失效的情况下，最大升力系数与最小升力系数相差0.2，此时叶尖截面呈现较典型的失速状态[10]，[11]。

为了进一步分析叶片的气弹稳定性，在两种工况下，对叶尖在摆振和扭转方向的变形进行频域分析，分析结果如图7所示。

图7 叶尖变形的频域分析Fig.7 Frequency domain analysis of blade tip deformation of the pre-bend blades

叶片一阶摆振的固有频率为0.962 Hz，在控制系统失效的情况下，图7中叶片在0.962 Hz频率处产生较大的扭转和摆振变形，表明叶片在该工况下产生了明显的气弹失稳现象[8]，[12]。

5 结论

本文以某2MW风力机为研究对象，基于SIMPACK软件了构建了整机刚柔耦合的气-弹-控联合仿真模型和方法，通过与行业内公认的风力机动力学分析软件GH-Blaed相比，本文提出的风力机气弹性能分析方法和程序可准确分析风轮叶片的载荷和气弹变形量。相比于GH-Bladed软件，本文采用与AeroDyn联合气动仿真的方法，可直观地分析叶片在一个周期内的升力系数变化情况，进而实现对叶片叶尖截面失速失稳现象的预测分析。