切出风速下风力机变桨故障叶片气动特性及准静态结构响应分析

刘一鸣，李 春，王渊博，李 根，闫阳天

（上海理工大学 能源与动力工程学院/上海市动力工程多相流与传热重点实验室，上海 200093）

化石燃料的过度使用及能源结构的不合理导致温室效应和环境污染等问题，且在不久的将来还可能导致资源枯竭之境况。为此，国家“十四五”规划纲要对发展风电等清洁能源提出了明确要求。不仅我国，世界各国均遭遇同样问题并都积极开发风能资源。此外，现代风力机通常均安装叶片变桨系统，然而变桨系统为风力机所有部件中故障率最高，且一旦发生故障将导致风力机无法正常运行时间最长的关键子系统之一，但鲜见与此相关气动及结构方面的研究。因此，在同时兼顾流体域和固体域基础上，对风力机特定状态下变桨故障叶片进行数值计算以探究其气动特性及结构响应显得尤为重要。

对于风力机变桨故障，国内外学者进行了一些研究。Ronsten通过实测相关数据提出风力机停机后叶片状态将影响其力矩大小及压力分布。据此不难推断，停机后变桨故障与变桨成功叶片由于状态差别极大，受载亦将存在较大差异。Tian等基于数值模拟方法并采用独立变桨技术控制大型三叶片水平轴风机运行，结果证明该策略可有效减小叶片载荷变化范围及风轮转矩波动幅度。因此，该独立变桨技术能保证叶片时刻处于最佳攻角，旋转过程中叶片状态始终较为稳定，而风力机发生变桨故障时则与此情况相反。Chen等通过研究台风“天兔”袭击汕尾红海湾风电场事件，指出变桨系统故障与风力机大面积严重损毁情况存在很大联系。闫学勤等采取将叶片方位角权系数分配与叶根气动载荷反馈相结合的变桨新技术进行叶根挥舞力矩及轮毂倾覆力矩控制，结果显示新型独立变桨技术可保证叶片具有低变桨故障率和小变桨误差。任年鑫等研究了台风达到兆瓦级水平轴风力机切出风速时，风力机顺桨停机后变桨故障与变桨成功叶片的受载情况，结果证明，风力机叶轮方位角及台风切入方向会对叶片受力产生影响，且变桨故障叶片载荷总大于变桨成功叶片。

由上述研究可知，风力机变桨故障对叶片受载影响极大，故此时叶片结构亦最危险。然而，在绝大多数研究中仅考虑了气动侧载荷却未涉及结构侧响应。此外，切出风速下变桨故障叶片受载及入流风速极大，但该风况下的风力机变桨故障研究鲜见相关报道。因此，本文针对NREL 5 MW风力机在切出风速下变桨故障叶片展开研究，对比分析变桨成功及各典型方位角下变桨故障叶片气动特性，然后开展变桨故障叶片准静态结构响应分析，以期为叶片结构设计提供更全面的参考数据。

1 研究对象及数值计算模型

1.1 叶片参数化建模及铺层

叶片是风力机最为关键的部件。本文利用Unigraphics NX软件建立NREL 5 MW风力机叶片三维实体模型，如图1所示，其中对关键部位进行了放大显示。

图1 NREL 5 MW风力机叶片模型Fig. 1 Blade model of NREL 5 MW wind turbine

由于本文结构侧计算涉及有限元分析，因此在Unigraphics NX软件中完成NREL 5 MW风力机三维实体建模后，还需对其赋予相关材料属性。共选取9种单一材料，通过对材料种类、纤维方向以及铺层顺序进行组合便可得到适用于风力机叶片各部位的复合材料。例如，叶片主梁帽需要承担弯曲载荷，可以选择抗拉性强的碳纤维和玻璃纤维并使纤维方向顺着叶片展向进行铺层，便可满足梁帽的结构力学性能要求，其他部位铺层方法类似。最终叶片铺层方案如图2所示。通过Abaqus软件Property模块的Composite Layup功能进行NREL 5 MW风力机叶片复合材料铺层。这与真实风力机叶片结构设计和制造工艺相同，从而可保证计算结果能够直接为相关设计制造及实际工程问题提供参考。

图2 叶片铺层方案Fig. 2 Blade layout

1.2 气动侧模型

气动载荷分析需通过与叶片同步旋转的局部坐标系定义各方向力矩。叶片坐标系及各方向力矩如图3所示。

图3 叶片坐标系及各方向力矩Fig. 3 Blade coordinate system and torque in all directions

为降低计算域边界对数值模拟的影响并减小风洞阻塞效应，需尽可能扩大计算域尺寸以增加风力机叶片与各边界间的距离。本文中计算域以及边界条件如图4所示，图中右下角坐标系为数值模拟全局坐标系，为风轮直径126 m。入口边界条件为指数率风剪切速度入口，参高度处是NREL 5 MW风力机25 m·s的切出风速，地面粗糙度指数选为0.16，出口为压力出口，其值为1个准大气压.地面为无滑移壁面，顶面是对称平面。两侧面采用周期性边界条件，即计算域一侧面网格节点与另一侧面网格节点一一对应，某侧周期边界计算域外“镜像单元”信息由紧邻另一侧周期边界计算域内的单元提供。CCM+前处理创建高质量多面体网格，且壁面边界层采用低处理（＜1）。近壁面网格高度为5×10m，边界层网格共24层，总厚度为0.1 m。同时为使风力机数值模拟准确度更高，对其尾迹

图4 计算域以及边界条件Fig. 4 Calculation domain and boundary conditions

计算域最终网格分布如图5所示。由图中可看出，相较于采用四面体网格，全部采用多面体网格类型能以更少的网格数量达到相同的计算精度。本文中气动侧均利用计算流体力学软件STAR-区域作了加密处理。

图5 网格分布Fig. 5 Grid distribution

数值计算采用隐式非定常算法（implicit unsteady），时间步长为0.001 s，对流项与时间项的离散均使用二阶格式，压力-速度解耦基于SIMPLE算法，湍流处理采用SST模型，湍流强度取0.05，湍流长度尺度为1.0 m。当数值模拟运行至载荷处于相对稳定或呈周期性振荡时，则判定为已收敛。

1.3 结构侧模型

图6为流固耦合结构侧计算模型。风力机叶片几何模型由片体（无厚度曲面）构建，这是由于有限元分析中若三维结构某一方向尺寸远小于另外两个方向尺寸即可采用片体模型，如此既能保证精度又可节省相当多的计算资源。同样为节省计算资源，将风力机叶片划分为结构化网格。与非结构化网格相比，结构化网格生成速度快且质量高，基于结构化网格数值计算的数据存储结构更简单且占用计算内存较少。有限元分析中片体结构化网格适合使用S4R单元，其为一种通用壳单元类型，既可应用于厚壳问题求解，亦能应用于薄壳模型计算。此外，将叶根一圈通过运动耦合约束在叶根圆形的中心点（即约束点）上，然后在该约束点施加ENCASTRE边界条件对六个自由度同时进行约束。

图6 结构侧计算模型Fig. 6 Calculation model for blade structure

2 准静态结构分析方法

2.1 双向弱流固耦合

双向弱流固耦合流程如图7所示。该流程共有8个步骤，对各步骤的解释分别为：①对固体域进行网格划分及物理模型设置，并参考气动计算结果施加一个近似载荷条件，然后单独运行一段时间以验证结构侧相关设置的准确性及可行性；②流体域保留气动计算数据结果，并对多面体网格应用Morphing功能以控制其变形；③判断流固两侧柔性结构是否完全匹配，若一切正常便可通过Import mode程序导出流体域柔性结构所受压力和黏性剪力，并通过abaqus关键字加载至固体域柔性结构上；④对结构侧进行计算，柔性结构变形等数据保存至ODB文件；⑤首先将结构侧的ODB文件加载至Import mode程序，其次由Import mode程序提取其中变形数据传递给STAR-CCM + ，最后STAR-CCM +对柔性结构形状进行更新，同时利用Morphing功能控制计算域网格变形以适应当前柔性结构形状；⑥对流体域进行计算并得到柔性结构新的压力以及黏性剪力；⑦重复前面第③～⑥步直至流固两侧相关特征不再有明显变化为止（本文中流固两侧数据共交换4次）；⑧双向弱流固耦合数值模拟完成。

图7 双向弱流固耦合流程Fig. 7 Flow chart of two-way weak fluid-structure coupling

2.2 屈曲分析

细长结构体由于长度方向尺寸远大于另外两个方向尺寸而表现出柔性特征，当受到微小载荷或扰动后便易发生失稳。当这些外部微小影响消失后，细长弹性体能够恢复至稳态，但若载荷或扰动较大时，细长弹性体则很可能发生屈曲。目前风力机大型化趋势非常明显，风力机越来越长导致其柔性增加，NREL 5 MW风力机叶片长达61.5 m，属于标准的细长弹性体，因此有必要对其进行屈曲分析。屈曲分析一般分为线性与非线性，基于线性屈曲所研究的结构并不会发生塑性变形，而是在结构弹性阶段产生弹性失稳。线性屈曲分析最终所得临界载荷可能会稍大于结构发生屈曲时的实际值，但线性屈曲分析计算效率非常高，且可作为结构分析时的前期预估值，故本文采用线性屈曲分析。

3 结果分析

3.1 单叶片模型可行性验证

首先对比全叶轮模型（full rotor model，FRM）与单叶片模型（single blade model，SBM）的计算结果，以便为减少计算量的策略提供理论依据。图8为SBM与FRM挥舞力矩和摆振力矩时域图。对比图8（a）、（b）中的力矩可知：SBM与FRM挥舞力矩和摆振力矩时域曲线基本重合，初步证明了基于SBM计算的可行性。

图8 SBM与FRM挥舞力矩和摆振力矩时域图Fig. 8 Time domain spectra of flapwise moment and edgewise moment for SBM and FRM

将SBM和FRM挥舞力矩和摆振力矩脉动时域序列进行快速傅里叶变换，可得到表1中的SBM与FRM载荷脉动主要刺激频率。由表中可以看出，SBM与FRM挥舞力矩和摆振力矩主要刺激频率相同或极其接近，这更进一步证明了水平轴风力机SBM与FRM载荷脉动差异非常小。

表1 SBM与FRM载荷脉动主要刺激频率Tab. 1 Main excitation frequency of load pulsation for SBM and FRM (单位：Hz)

经过叶片顶点且与叶片坐标系轴垂直平面上SBM与FRM在轴方向的叶尖涡量场如图9所示。从图中可观察到，SBM与FRM叶尖涡大小及位置几乎完全相同，说明SBM与FRM流场流动状况极为接近，此亦为SBM与FRM载荷频域主要刺激频率以及力矩时域曲线差别甚微的根本原因。综上可知，使用SBM进行数值模拟具有较高的准确性及可行性，故后文均基于SBM开展计算。

图9 SBM与FRM在Z轴方向的叶尖涡量场Fig. 9 Blade tip vortex of SBM and FRM in Z-axis direction

3.2 变桨故障叶片气动特性

图10为NREL 5 MW风力机在切出风速下变桨故障与变桨成功叶片挥舞力矩和摆振力矩。对比图10（a）、（b）可看出，挥舞力矩比摆振力矩大1个数量级，但都表现出变桨故障叶片力矩大于变桨成功叶片的特征，其中变桨故障叶片挥舞力矩平均值为变桨成功叶片的13.8倍。

为进一步分析变桨故障叶片与变桨成功叶片的差异，对周围流场速度云图进行分析，结果如图11所示。由图中可知，变桨故障叶片的尾迹较变桨成功叶片的更为明显。这是因为变桨故障叶片从来流风中吸收了更多能量，因此变桨故障叶片的受载比变桨成功叶片的大，这一结果与图10中的数据非常吻合。综上，后文仅分析变桨故障叶片特征。

图10 变桨故障与变桨成功叶片挥舞力矩和摆振力矩Fig. 10 Flapwise moment and edgewise moment of a blade under fault and normal status of pitch system

图11 变桨故障叶片与变桨成功叶片周围流场速度云图Fig. 11 Velocity contour around a blade under fault and normal status of pitch system

NREL 5 MW风力机达到切出风速停机后叶轮方位角是随机的。图12为0°、60°、120°、180°、240°和300°共6个方位角下变桨故障叶片挥舞力矩和摆振力矩时域图。从图12（a）中可看出，0°方位角下挥舞力矩最大，其次是60°和300°方位角下的，然后是120°、240°、180°方位角下挥舞力矩最小。挥舞力矩在各方位角下的差异主要是由风剪切造成风速从地面沿竖直方向越来越大所致，故0°方位角叶片平均来流风速最大，180°方位角叶片平均来流风速最小。此外，60°与300°方位角以及120°与240°方位角叶片平均来流风速两两相同，因此它们的力矩大小也一样，故后续研究仅分析60°和120°方位角，对300°和240°方位角的结果不再赘述。

图12 各方位角下变桨故障叶片挥舞力矩及摆振力矩时域图Fig. 12 Time domain spectra of flapwise moment and edgewise moment of a blade with pitch fault at different azimuths

3.3 准静态结构响应分析

图13为双向弱流固耦合叶片应力及叶尖位移。由图中可知，切出风速下4个典型方位角的变桨故障叶片均出现了应力集中现象，其中0°方位角下变桨故障叶片应力集中最强，为9.886 MPa，而180°方位角下应力集中最弱，为6.932 MPa，后者较前者降低了29.8%。叶尖位移变化趋势与应力集中变化趋势相同，从0°、60°、120°到180°方位角下叶尖位移逐渐减小，其中180°方位角下变桨故障叶片叶尖位移较0°方位角下的减少了32.7%。此外，还发现双向弱流固耦合结构侧结果与图12中气动力矩情况吻合较好，证明了本文计算结果的准确性。

图13 双向弱流固耦合叶片应力及叶尖位移Fig. 13 Stress and tip displacement of a blade under two-way weak fluid-structure coupling

图14为0°方位角下变桨故障叶片前六阶屈曲因子和屈曲模态。屈曲因子的物理意义是载荷与其之乘积即为叶片发生屈曲的临界载荷。从图中可看出，最小的屈曲因子是一阶屈曲因子1.234 5，说明叶片并未发生屈曲。此外，仔细观察图14能够发现，随着屈曲阶数的增加，叶片发生屈曲部位的面积越来越大，且屈曲因子逐渐提升，但实际工程中最关注的是一阶屈曲因子，因为一阶屈曲发生之后其他阶结构屈曲将无意义。此外，其他方位角结构屈曲情况与0°方位角变化趋势一样，且随着方位角增加叶片屈曲越弱，其中180°方位角下变桨故障叶片的一阶屈曲因子较0°方位角下的增大20.2%，结构安全性更高。需要说明的是，本文结构侧屈曲分析载荷来自双向弱流固耦合最后一次流固两次数据交换并计算结束后气动侧的载荷数据。

图14 0°方位角下变桨故障叶片前六阶屈曲Fig. 14 The first six order buckling of a blade with pitch fault at 0°azimuth

4 结 论

本文基于计算流体力学方法对NREL 5 MW风力机变桨故障/成功叶片气动侧状态进行对比分析，并利用双向弱流固耦合及曲屈分析对典型方位角下变桨故障叶片特征展开研究。主要结论有：

（1）SBM与FRM挥舞力矩和摆振力矩时域序列、频域特征及涡量场强度和分布差别均非常小或相同，故采用节省大量计算资源的SBM模型可行。

（2）切出风速下变桨故障叶片较变桨成功叶片受载更大，且前者流场尾迹比后者亦更明显，说明变桨故障叶片从来流风中吸收了更多能量。

（3）180°方位角下变桨故障叶片较0°方位角下变桨故障叶片应力及叶尖位移分别减小29.8%和32.7%，一阶屈曲因子增加20.2%。