大型风力机叶片净空计算方法研究

董笑宇 曲春明 葛文澎 夏 冰 苗得胜

明阳智慧能源集团股份公司

0 引言

随着风能技术的不断发展，风电机组越来越大型化，叶轮直径也越来越大，在现有加工工艺和材料技术未完成革命性突破前，叶片越来越注重轻量化设计，从而使叶片更加柔性化［1］。轻型化、柔性化的叶片在实际运行中发生几何非线性大变形的几率增大，会出现叶片与塔筒相撞的扫塔现象，严重时叶片会断裂，造成巨大的经济损失。为保障风电机组大型化对叶片结构强度的要求及风电机组的安全运行，对风机叶片的净空研究十分重要。

叶片净空计算方法很多，风力机载荷计算方面主要分为叶素动量理论BEM、计算流体力学CFD、和升力线方法等；叶片变形方面，按几何结构可分为梁单元法、壳单元法和实体单元法，按方法原理可分为多体动力学MED 和有限元结构力学FEM。Liu 等人［2］利用计算流体力学和多体运动系统耦合方法计算了海上风力机组运行状态下的动态结构响应，Yang 等人［3］研究了在旋转轮毂上欧拉-伯努利梁的柔性运动，推导了一组描述梁轴向、横向和旋转运动全耦合的非线性积分微分方程。

本文比较了几种叶片净空计算方法，首先从风力机空气动力计算方法进行比较，再从叶片结构响应方法比较，最后从工程实用角度出发给出推荐的计算方法。

1 叶片模型

本文以某三叶片海上机组为研究对象，选取方位角为180°的叶片（叶轮呈“Y”形时叶尖向下的叶片）为计算模型，机组及叶片参数见表1，叶片模型见图1。

图1 机组叶片模型

表1 机组及叶片参数

有研究表明，当机组处于以下两种状态时，叶片净空距离最小：①风轮达到额定转速、环境风速接近额定风速、叶片尚未变桨；②机组运行存在极端负切变，叶片载荷过大。鉴于海上风切变较小，本文选取第一种状态作为分析工况。

本研究叶片采用纤维增强型复合材料，根据叶片弦长、扭角、预弯、厚度、主梁宽度及铺层情况，FOCUS软件可输出叶片各截面的挥舞刚度、摆振刚度、轴向刚度和扭转刚度，同时输出用于结构计算的截面质量属性。本文采用已知的气动外形参数和FOCUS输出结构参数作为叶片净空计算的数据来源。

2 载荷计算方法

风力机气动性能及载荷计算的方法有计算流体力学CFD、叶素动量理论BEM、升力线、升力面、面元法及区域求解法等［4-5］，其中BEM 理论和CFD 计算方法较普遍，以下为两种方法的对比和结果分析。

2.1 叶素动量理论BEM

叶素动量理论BEM 由英国科学家Glauert 提出，将动量理论与叶素理论相结合求解轴向诱导因子a和周向诱导因子a′，进而计算风轮叶片上的力和力矩。叶素动量理论建立后，一直用于风力机气动特性的研究，正确性已得到充分验证［6］，该理论计算步骤如下［7-8］：

（1）a、a′赋初值，作为迭代第一步；

（2）根据公式（1）计算入流角：

其中：Ω 为叶轮转速，r为流体微团所在的叶轮展向位置，u为流体沿流线方向的运动速度，下标∞为流体位于上游无穷远处；

（3）利用公式（2）计算攻角：

θ为局部扭角和局部变桨角之和；

（4）根据翼型气动特性曲线得到叶素的升力系数CL和阻力系数CD；

（5）通过轴向和周向诱导因子公式（3）、（4）计算a和a′：

其中，B为叶片数，F为考虑叶根和叶尖流动损失的总损失系数，H为推力系数的修正；

（6）根据更新后的a和a′判断迭代过程的收敛性；

（7）由此得到叶轮的推力系数和功率系数，见公式（5）（6）：

ρ为微团密度，T为叶轮上的推力，M为叶轮上的转矩。

2.2 计算流体力学CFD

计算流体力学CFD 就是通过计算机和数值方法求解流体力学控制方程，获取各条件下流动的数据和作用在绕流物体上的力、力矩、流动图像和热量等。

流动的数值计算一般有三个方法：直接数值模拟法DNS，大涡模拟法LES 和雷诺平均方程RANS，但解决流体动力学及空气动力学的流动和工程实际问题的主要是雷诺平均法。雷诺平均法的基本思想是将Navier-Stokes 方程时均化，用时均值来代表瞬时值。该方法可节省大量时间与资源，因此在工程中被广泛应用［9］。

2.3 方法验证

BEM 方法在无黏假设下的一维动量理论与实际三维黏性流动之间存在差异，为此，CFD 设置了整机、风轮、单叶片三种计算模型，图2、图3 为叶片方位角为180°时BEM和CFD载荷计算结果。

图2 叶片推力

图3 叶片力矩

由图2 可知，两种计算方法得到的沿叶展方向的推力整体变化趋势一致，均先增大后急剧减小，但单叶片模型推力分布区别较大，BEM 在叶片中段位置较大，在叶尖位置较小，且推力沿叶展方向更光滑，即推力值沿叶展方向增长且增长率较平稳均匀，CFD 单叶片整体推力较BEM 存在-2.09%差值；由图3可见，两种计算方法得到的沿叶展方向的力矩整体变化趋势一致，但CFD 叶根与叶尖位置比BEM 大，CFD 单叶片模型力矩较BEM 转矩存在8.03%差值。

CFD 整机、三叶片风轮和单叶片模型沿叶片展向推力分布一致，单叶片推力最大，整机最小。与风轮模型相比，整机模型存在一定的塔影效应，使位于塔筒前的叶片所受推力减小；与单叶片模型相比，三叶片风轮模型叶片推力更小，这是因为当叶片以较高线速度旋转时，每个叶片都会对轨迹后方叶片所处流场造成一定的影响，因此三个叶片气动性能均存在一定程度的损失。

CFD 整机、三叶片风轮和单叶片模型力矩结果中，塔筒与叶片对气动性能干扰作用更加明显，尤其是叶片间流场的相互干扰，很大程度降低了叶片力矩。

图4 为整机模型流场图。由图可见，由于塔筒的存在使叶片与塔筒之间流场方向变得紊乱且流速变慢。图5、图6 分别为叶片展向50%和80%位置翼型周向压力分布图，取单叶片模型与整机模型为对象，弦长方向位置0 处为前缘，1 处为后缘。由图可见，单叶片与整机模型翼型周向压力分布规律一致，压力面前缘压力最大，后缘稍低，中段为负压；吸力面均为负压，前缘压力最小，由前向后逐渐升高。

图5 叶片展向50%位置翼型周向压力分布

图6 叶片展向80%位置翼型周向压力分布

此外，从图5中还可看到，单叶片压力面前缘压力较整机大，吸力面则小，因此单叶片模型前缘压差更大。该规律在图6 展向80%位置处依然可得，即在更靠近叶尖处，单叶片模型前缘压差更大，两种模型在压力面压差最大达155 Pa，吸力面达-970 Pa，说明靠近叶尖位置受塔筒和其他叶片的干扰更明显。单叶片模型前缘位置压差更大说明单叶片比整机所受推力更大。

3 结构变形计算方法

现阶段对风力机叶片振动及变形的研究大多采用梁模型或多体动力系统。本研究采用传统有限元方法FEM 进行叶片静力学分析：将叶片简化为梁模型，通过调整叶片截面参数使梁模型获得更准确的挥舞、摆振及拉伸方向刚度等参数，同时采用多体动力学方法开展叶片结构分析，即将叶片离散为多个弹簧与质量的组成。

3.1 FEM方法

FEM为将所求解连续区域的弹性体离散为有限个单元，各个单元通过相邻节点连接，使复杂边界连续体无限自由度问题简化为单元节点有限自由度问题，从而得到整个区域的近似解，如果采用的单元满足收敛要求，那么单元数越多越趋近于精确解。

风力机叶片具有细长结构的特点，在分析叶片结构时，叶片通常被假设为梁模型(1D-FEM)，在有限元中，梁单元模型的基本假设是原截面受弯曲力后依然是平面。

本文采用的线性梁单元B31 是考虑剪切变形的Timoshenko 梁单元［10］，受力前截面和中性面是垂直的，弯曲后不再垂直。该单元适用于模拟剪切变形起重要作用的深梁和模拟剪切变形不太重要的细长梁。

3.2 MBD方法

MBD 方法为将柔性叶片离散为若干个有一定质量与形状的体，体与体之间用旋转铰、弹簧及阻尼器连接，组成弹簧质量多体系统，通过系统的位移表征叶片的运动和变形。

为准确模拟柔性叶片弹性变形及力学特性，根据最小势能原理和伯努利梁假设确定弹簧刚度系数。多体系统动力学拉格朗日数学模型可通过虚功形式的动力学普遍方程建立，即多体系统在虚加惯性力、真实主动力与理想约束力的作用下处于静平衡，联立运动约束方程和受力平衡方程对静力学进行计算［11-12］。

3.3 方法验证

首先对梁模型与多体模型进行模态分析，验证两种模型质量刚度分布的准确性，模态分析结果见表2。

根据模态分析结果可看出，两种模型的挥舞频率较接近，梁模型模态二阶挥舞频率比多体动力系统模态结果误差更大，相差0.9%。模态结果误差的原因是模型没有采用叶片的实际外形（以近似叶片模型的方法），故叶片整体质量与刚度等效存在一定的误差，但与叶片实际测量值误差均在2%以内，因此认为模型质量刚度分布符合精度要求。

为对比多体动力学方法与有限元方法计算叶片变形情况，采用BEM方法得到的载荷值对梁模型和多体系统进行静力分析，叶尖变形结果见表3。

表3 叶尖变形结果

此机组叶片最大净空距离为34 m，由表3 可见，MBD 和FEM 两种静力分析方法结果均未超过34 m，且十分接近，仅相差0.48，与最大净空距离相比不超过20.59%，采用FEM 梁单元分析叶片变形稍大于MBD结果。

图7 为1D-FEM 方法计算得到的叶片变形图。由图可见，取叶片展向50%与80%两个位置a、b，原模型a 处叶片切向夹角θa为1.5°，加载变形后切向夹角θa′为-8.5°，角度变化10°，位移2.92 m，原模型b 处叶片切向夹角θb为7.0°，加载后b 处切向夹角θb′为-20.3°，角度变化27.3°，位移14.12 m。

图7 叶片变形图

4 叶片变形敏感因素分析

为分析叶片变形影响因素，本文从材料刚度与载荷角度出发，采用BEM 耦合FEM 的方法计算叶片空气动力特性和结构响应，并将BEM载荷计算结果按比例作为叶片静力分析的输入。载荷对叶尖变形分析结果见表4，变形曲线见图8。

图8 叶尖变形随载荷变化曲线

表4 叶尖变形随载荷变化结果

由图8 可见，叶尖变形随载荷增大呈单调递增正比关系，当载荷比值为1.24 时，叶片将发生扫塔现象。将材料刚度（包括弹性模量和剪切模量）按比例取值，分析刚度对叶尖变形影响，结果见表5，变形曲线见图9。

图9 叶尖变形随刚度变化曲线

表5 叶尖变形随刚度变化结果

由图9可见，叶尖变形随刚度增大而减小，当刚度比值为0.83时，机组将发生扫塔现象。

5 结论

本文以某108 m 长叶片为模型，对不同叶片净空计算方法进行了研究和对比，计算方法包括气动性能计算方法与结构响应计算方法，结论如下：

1）不同计算方法得到沿叶展方向推力与力矩均呈先增大后急剧减小的趋势。相比CFD，BEM 载荷结果在叶片中段更大，沿叶展方向均匀增长，在叶尖位置更小。由于叶片气动外形，导致CFD 推力在根部及中段较小，在叶尖有明显增加。

2）塔筒效应和叶片间尾流相互影响均会导致叶片气动性能的损失。

3）MBD 二阶模态分析误差较FEM 低0.9%，相同载荷下，MBD 叶尖变形27.00 m，相较FEM 变形少0.48 m。

4）载荷增大会引起叶尖变形正比增大，净空减小，而刚度增大则会使叶尖变形减小，净空增大，且变形随刚度变化逐渐降低。

5）CFD 方法由于考虑了叶片实际外形而精确度更高，但计算耗时更长，BEM 方法与CFD 方法差值较小，但计算效率更高；结构分析方面，MBD与FEM 在计算效率方面差距不大，但模态分析MBD 比FEM 准确性更高一些，因此采用BEMMBD 计算叶片净空可快速得到准确结果，较适合于实际工程。