基于气动弹性剪裁的大型风力机弯扭耦合叶片力学性能分析

2021-08-23缪维跑张万福闫阳天

动力工程学报 2021年8期

张立，缪维跑，李春,2，张万福,2，闫阳天，王博

(1．上海理工大学能源与动力工程学院，上海 200093； 2．上海市动力工程多相流动与传热重点实验室，上海 200093)

由于传统化石能源不可再生及其引起的环境污染问题日益加剧，开发利用零污染可再生能源逐渐受到各国重视[1]。风能作为最具前景的替代能源之一，据世界气象组织估计，地球上可利用风能总量约为2×1010MW，为水能总量的20倍，其作为未来清洁能源的重要组成部分，对人类发展及国家能源战略具有深远意义[2]。风力机为吸收利用风能的主要设备，现阶段研发的兆瓦级风力机风轮直径超过200 m，成为人类建造的外形尺寸最大的动力机械。

叶片作为风力机汲取风能的关键部件及风力机源动力输入的首要载体，直接影响风力机气动性能及结构特性，对风力机经济性和安全运行均有重要影响[1]。为实现叶片轻量化设计的同时提高其承载能力，现代大型风力机叶片多采用高比强度、高比刚度且低密度的复合材料制造[3]。气动弹性剪裁源自飞机机翼设计，目的在于使机翼在飞行中产生理想的耦合形变[4]。20世纪90年代末引入风力机叶片设计，应用于风力机的减载和降振控制。气动弹性剪裁利用复合材料显著的各向异性，通过对叶片上、下表面纤维角度进行偏轴镜像对称铺设，使叶片气弹中心与扭转轴不平衡而产生弯扭耦合特性。

目前，国内外学者针对弯扭耦合叶片开展了诸多研究，主要集中于弯扭耦合叶片减载性能、材料特性对叶片弯扭耦合性能影响以及弯扭耦合对叶片结构性能影响等方面的研究。刘旺玉等[5]采用有限元方法研究了偏置角度及材料体积比对叶片弯扭耦合特性的影响，并将主梁及蒙皮镜像偏置不同角度以获取更优的弯扭耦合特性。刘宇航等[6]通过镜像对称铺设蒙皮纤维实现叶片气动弹性剪裁，采用节点位移法研究了蒙皮偏置角度对叶片弯扭耦合系数的影响。Gözcü等[7]基于多体动力学模型研究了弯扭耦合叶片降低叶根载荷的效果以及其对减轻齿轮箱和变桨控制器等系统的损伤等效疲劳载荷的影响。Hayat等[8]采用ABAQUS及ADAMS研究了铺层材料、铺层厚度及纤维角度3种铺层不平衡对弯扭耦合叶片颤振性能的影响。Sener等[9]研究了玻璃钢-碳纤维混合弯扭耦合叶片对降低风力机系统部件(如叶片、轴承、齿轮箱及偏航转轴等)内部疲劳损伤的影响。Miao等[10]基于弯扭耦合叶片被动自适应特性，建立了不同结构特性弯扭耦合叶片模型并研究其在极端台风环境下减载、降振及抗台风性能。Shakya等[11]通过参数化方法探究蒙皮及主梁铺层材料、铺层厚度以及铺设先后顺序对弯扭耦合叶片颤振不稳定性的影响。

由于叶片为展向长、弦向短的细长弹性体，对其进行结构分析及优化设计时常采用线性梁模型[12]，但难以获取具体叶片结构信息，如层间应力应变分布和结构屈曲失稳现象等[13]。此外，一些文献中常忽略重力载荷作用[14]，而风力机趋于大型化发展，叶片重力对整体结构影响较大，结构响应求解时不应忽略。基于此，笔者以NREL 5 MW风力机叶片为研究对象，通过对CAD软件NX二次开发建立叶片三维模型，结合复合材料铺层并采用CFD方法求解叶片表面压力，研究主梁偏置角度对叶片力学性能的影响。

1 研究对象与模型

1.1 研究对象

研究对象选取美国可再生能源实验室公开的NREL 5 MW风力机，其主要结构及运行参数见表1。

叶片由叶根至叶尖可分为叶根区、过渡区及叶尖区，叶根区主要用来保证叶片结构刚度，叶尖区主要用以提供更大气动力。叶片主梁帽为抵抗挥舞力矩的主要承载结构，压力面及吸力面主梁帽与腹板共同构成箱型结构，主要承受剪切及扭转力矩。为使叶片具有足够强度，叶片内布置2个剪切腹板。为使水平轴风力机叶片在设计工况下达到最佳攻角，叶片展向不同位置处翼型以变桨轴为扭转中心，预设不同扭角。NREL 5 MW风力机叶片截面翼型扭角及弦长沿展向的变化趋势如图1所示。由图1可见，叶片截面翼型扭角由叶根至叶尖逐渐减小，最大扭角为13.308°，叶片截面翼型弦长由叶根至叶尖先增大后减小。

表1 NREL 5 MW风力机主要设计参数

图1 风力机叶片扭角及弦长沿展向的分布

1.2 模型建立

精确的叶片几何模型对气动载荷及结构响应计算尤为重要，直接影响计算结果精度、收敛及过程稳定性。风力机叶片外形为一系列翼型曲线构成的空间曲面，为满足叶片气动及结构性能要求，除叶根处为圆或近似圆外，其他部位均由不同翼型组成，且扭角及弦长沿展向不断变化，一方面增大了叶片实体建模的难度，另一方面也更突显基于CAD平台二次开发参数化建模的重要性。为此，笔者采用NX/Open Grip语言对三维建模软件NX进行二次开发，建立包含NREL-S、DU、NACA、Riso及FFA-W等风力机叶片常用翼型的翼型数据库，自动调用NX曲线组命令来快速构建叶片几何模型。风力机叶片二次开发程序流程图如图2所示。

图2 风力机叶片二次开发程序流程图Fig.2 Flow chart of secondary development program

具体开发步骤为：选择具体翼型并对其进行平移、缩放及旋转操作，以确定翼型展向位置、弦长及扭角等参数；获得叶片各截面翼型后，自动调用NX曲线组命令得到叶片几何模型。该程序可生成主流有限元软件所支持的众多文件格式，能准确快速建模的同时又具有一定通用性，可为后续流场及结构计算打下良好基础。图3为采用NX二次开发程序创建的风力机叶片各截面翼型及几何模型。

(a) 各截面翼型曲线

(b) 叶片实体模型图3 叶片各截面翼型及实体模型Fig.3 Airfoil of blade cross section and model of the blade

2 叶片气动载荷求解

2.1 计算方法

叶片气动载荷求解主要采用叶素-动量(BEM)模型、涡尾迹模型、致动模型及CFD模型4种。1935年，Glauert首次将叶素理论与动量理论相结合提出BEM模型，之后学者提出诸多修正方法对其进行完善与改进，如叶尖损失修正、推力系数修正及径向流动修正等[15]。相比其余模型，BEM模型在可靠翼型气动数据下可快速提供较精确结果，但无法实现叶片周围流场的可视化。涡尾迹模型常用于叶片复杂非定常尾迹研究中[16]，相比BEM模型，需占用更多计算资源且存在发散的可能。此外，由于涡尾迹模型忽略黏性效应，也限制了其在风力机某些方面的应用[17]。致动模型分为致动盘模型、致动线模型和致动面模型，叶片分布载荷用盘、线和面表示，可较好地反映叶片三维流动和尾迹动态发展，但其仍基于叶素理论，需可靠翼型气动数据，相比BEM模型，耗时更久且计算结果精度更低[17]。随着计算机技术的不断发展，CFD方法得到广泛重视，其优点在于无需预先给定翼型气动数据，可实现叶片周围流场的可视化,且计算结果与实验值较为吻合[17]，因此选用CFD方法求解叶片气动载荷。

CFD方法中，基于时均统计的雷诺时均法(Reynolds-Averaged Navier-Stokes, RANS)因计算量较小而在工程实际中应用较广，其通过采用湍流模型使Navier-Stokes方程封闭。目前，基于RANS的湍流模型已广泛应用于风力机数值模拟中，如S-A湍流模型、k-ε湍流模型和SSTk-ω湍流模型等。其中，SSTk-ω湍流模型被誉为二维翼型和三维叶片数值模拟的最佳湍流模型[18]。故选用SSTk-ω湍流模型求解叶片气动载荷，压力速度耦合使用SIMPLE算法，对流项差分格式采用二阶迎风格式。

笔者采用有限元软件ANSYS中CFX模块计算流体侧叶片表面压力，并将其加载至结构侧复合材料有限元模型，运用有限元方法对叶片进行结构响应计算，具体分析流程见图4。

2.2 数值计算模型

叶片计算域包括外流域及旋转域，由于风轮流场具有周期性特点，取1/3风轮流场作为计算域以降低计算量，计算域边界条件设置如图5所示，图中R指叶片长度。进口设为速度入口，进口长度取2倍叶片长度，速度大小为11.4 m/s；出口为压力出口，出口长度取5倍叶片长度，压力大小为标准大气压(即101.325 kPa)；两侧矩形区域为周期性边界条件，宽取5倍叶片长度；叶片表面为无滑移壁面；旋转域施加转速模拟叶片旋转状态，转速大小为12.1 r/min。

图4 分析流程图Fig.4 Analysis flow chart

图5 计算域边界条件Fig.5 Boundary conditions of the calculation domain

为保证整体网格质量，采用非结构网格进行划分，其对叶片等复杂几何形状有良好适应性。由于外流域对最终计算结果影响较小，故适当放大网格尺寸以减少计算量。叶片表面覆盖边界层网格以确保叶片壁面边界层流动模拟的准确性。通过网格镜像对称设置，使周期性边界两侧网格节点一一对应，周期性边界一侧计算域“镜像单元”信息可由另一侧周期性边界计算域单元获得。划分网格后，整体网格畸变度(Skewness)为0.215 95，网格数约1 900万，流场网格划分见图6。叶片气动载荷分布如图7所示。

图6 流场网格划分Fig.6 Meshing of fluid domain

3 叶片结构响应计算

3.1 铺层设计原则及材料

风力机叶片铺层设计时应遵循以下原则[19]：(1)铺层定向原则，叶片所采用铺层方向数应尽可能少，简化铺层设计工作量且利于生产制造，常用铺层方向主要有-45°、0°、+45°和90°等。(2)按内力方向铺层取向原则，对于承受压缩载荷或单轴拉伸载荷的部件，其铺层方向应与载荷方向相同。叶根增强结构和叶片主梁主要承载压缩载荷和轴向拉伸载荷，铺层方向应取0°或90°。叶片及剪切腹板蒙皮主要承受剪切载荷，应采用-45°及+45°对称铺层。(3)变厚度设计原则，铺层厚度与载荷大小呈正相关，所受载荷越大，铺层越厚。叶片所受弯曲载荷由叶根至叶尖逐渐减小，故主梁层合板厚度也趋于减小，过渡段变厚度设计时各层台阶宽度应不小于2.5 mm。为避免各层台阶处产生剥离破坏现象，叶片表面最终应由连续铺层覆盖。

(a) 压力面

(b) 吸力面图7 叶片气动载荷分布Fig.7 Blade aerodynamic load distribution

选用胶衣、泡沫、玻璃钢、碳纤维、双轴向蒙皮和三轴向蒙皮6种层合板结构作为NREL 5 MW叶片铺层材料。其中，胶衣、泡沫、玻璃钢和碳纤维均为[0]2单轴向层合板，双轴向蒙皮为[±45]4双轴向层合板，三轴向蒙皮为[±45]2[0]2三轴向层合板，各层合板材料力学性能参数见表2[20]。表2中，Ex、Ey分别为展向和弦向弹性模量；Gxy为剪切模量；μxy为泊松比；ρ为材料密度;TL、CL分别为极限纵向拉伸和抗压强度。

表2 叶片层合板材料力学性能Tab.2 The mechanical properties of blade laminate

3.2 叶片铺层策略

由于叶片各部位受力不同，将叶片分割为多区域以便铺设不同材料及厚度的层合板，使叶片同时满足强度、刚度、稳定性及疲劳等要求[21]。为此，将叶片吸力面和压力面分别沿展向分割为25段，沿弦向依次分割为前缘、前缘镶板、主梁、尾缘镶板、尾缘加强及尾缘6个区域。叶片内部结构及材料如图8所示,图中各区域左侧为叶片外层铺层材料。叶片压力面和吸力面对称铺层，最终叶片共分为302块区域。

由图8可见，叶片表面及腹板均采用三明治夹芯结构，铺层方向为叶片外侧指向内侧。叶片上下表面主梁区采用镜像铺层方式，在两侧纤维共同引导下叶片将向预设方向扭转，产生弯扭耦合特性，其他区域采用传统叶片常采用的对称铺设方式，对称及镜像铺层方式如图9所示,图中θ为偏置角度。叶片前缘镶板区、尾缘镶板区及腹板夹芯层均为PVC泡沫，主梁区夹芯层为玻璃钢或碳纤维，尾缘加强区夹芯层为玻璃钢和PVC泡沫，叶根铺设多层三轴向蒙皮以加强叶根性能。叶片表面铺设蒙皮以保证其几何气动外形。外侧蒙皮由辅助层和强度层构成:辅助层为最外层的胶衣层，可避免叶片受海洋环境腐蚀以及紫外线长期照射老化，并提供光滑气动表面；强度层为三轴向蒙皮，以防止纤维布开裂和出现局部凹陷并承担叶片扭转载荷。内层蒙皮采用三轴向蒙皮。剪切腹板两侧均铺设双轴向蒙皮。铺层结果模型如图10所示。

图8 叶片内部结构及材料Fig.8 Internal structure and material of blade

(a) 对称铺层

(b) 镜像铺层图9 叶片铺层方式Fig.9 Blade laying method

图10 复合材料叶片模型Fig.10 Composite blade model

3.3 结构网格划分

借助Ansys中的复合材料铺层模块(Ansys Composite PrePost, ACP)实现叶片铺层设计。在叶片根部施加固定约束，采用SHELL 181层合单元对叶片进行剖分，其为四节点六自由度单元，可用于复合材料壳结构线性、非线性及扭转分析等，最多可允许铺设250层材料，设置后可进行真实叶片厚度距离偏置和层压板铺层方向调整，全局网格尺寸取0.1 m，网格单元质量为0.946 31，共计划分约9万个网格节点及9万个网格单元，叶片结构网格划分如图11所示。

图11 结构网格划分Fig.11 Meshing of the structure

4 计算结果与分析

4.1 模态分析

模态是指与外部条件无关的结构固有属性，通过模态分析可确定结构固有频率及该频率下的结构响应，常用于预测结构是否产生共振及共振下的位移形变等[22]。

将CFD方法计算所得叶片气动载荷加载至复合材料有限元模型，以y轴正方向施加重力加速度来模拟重力载荷作用[23]，并在叶片根部添加固定约束后对叶片进行模态分析。研究表明偏置角度在-45°～45°时，叶片弯扭耦合特性最佳[5]，故选取-45°、-30°、-15°、0°、15°、30°及45° 7个角度，以探究主梁偏置角度变化对弯扭耦合叶片力学性能的影响，其中0°铺层主梁叶片为未耦合的传统叶片，其他角度均为弯扭耦合叶片，偏置角度θ定义如图12所示。

图12 叶片主梁偏置角度Fig.12 Offset angle of main beam of the blades

为验证所建复合材料叶片模型的准确性，将传统叶片固有频率与文献[24]和文献[21]中所获固有频率进行对比，前六阶固有频率如图13所示。Resor等[24]采用NuMAD设计叶片铺层，因其参考表面选择中层而存在一定误差，而缪维跑等[21]采用ABAQUS设计叶片铺层，其参考表面则选择顶层，较为精确。对比计算值与参考值，可见三者上升趋势一致且数值较为吻合，证明所建复合材料叶片模型具有一定准确性。由于该风力机额定转速为12.1 r/min，其对应叶片旋转频率1P、通过频率3P分别为0.202 Hz和0.605 Hz，而传统叶片一阶固有频率为0.881 64 Hz，与1P、3P相差大于10%，可见其不会产生共振破坏。

图13 叶片前六阶固有频率Fig.13 First six natural frequencies of the blades

图14给出了弯扭耦合叶片前六阶固有频率的相对变化量，该相对变化量均以传统叶片固有频率为基准。由图14可见，随偏置角度逐渐增大，各阶固有频率均先增大后下降，即将主梁材料偏轴镜像铺设后，叶片各阶固有频率均有不同程度下降，且正向及反向偏置相同角度后各阶固有频率下降量较接近。对比各阶固有频率相对变化量可知，叶片一、三和五阶固有频率下降最多，二、四和六阶固有频率下降较少，而一、三和五阶振型为挥舞振动，即主梁偏轴镜像铺设对叶片挥舞振动影响较大。主梁偏轴镜像铺设使得叶片弹性模量发生变化，而结构刚度矩阵K受弹性模量的影响，因此影响了叶片模态参数。分析各叶片一阶固有频率可知，θ=±45°时固有频率下降最多(约50%)，而θ=+15°和θ=-15°时固有频率下降最少(约20%)，对应一阶固有频率分别为0.693 69 Hz和0.683 40 Hz，与叶片通过频率3P较为接近，易发生共振现象。因此，通过偏轴镜像铺设实现叶片气动弹性剪裁的同时，应对模态固有频率予以重点关注，避免叶片发生共振破坏。

图14 弯扭耦合叶片固有频率的相对变化量Fig.14 Relative variation in natural frequency of the bend-twist coupling blades

图15给出了θ=0°及θ=45°时主梁叶片前六阶模态振型。由图15可见，两叶片低阶振型类似而高阶振型有所不同，而叶片振型主要有挥舞振动(垂直于旋转平面弯曲运动)、摆振振动(旋转平面内弯曲运动)及扭转振动(绕变桨轴扭转运动)3种。其中，一、二阶振型为一阶挥舞及摆振振动，三、四阶振型为二阶挥舞及摆振振动，五阶振型为三阶挥舞振动，六阶振型为一阶扭转振动，可见叶片抗扭转能力均较强。

4.2 静力学分析

结构线性静力学分析为基本力学分析过程，在一定理论假设下实际问题均可简化为线性问题，其主要用以分析特定载荷作用下的结构变形情况，经静力学分析后可获得结构基本信息(应力、应变和应变能等)。

对叶片进行静力学分析时可整体分析叶片结构应力、应变和应变能等分布规律，并为进一步优化弯扭耦合叶片铺层设计等提供参考。

图16给出了各偏置角度下叶片表面及内部腹板的最大应力及应变值。由图16可知，叶片主梁偏置角度变化会对叶片表面应力和应变产生一定影响。θ=-15°弯扭耦合叶片的表面最大应力最低，约为63.52 MPa，而θ=-45°弯扭耦合叶片的表面最大应力最高，约为90.41 MPa，与传统叶片分别相差约14.78%和21.31%。传统叶片表面最大应变最低，约为0.002 995 1，而θ=-45°弯扭耦合叶片的表面最大应变最高，约为0.005 381 7，较传统叶片增大约79.68%。

(a) θ=0°一阶

(b) θ=45°一阶

(d) θ=45°二阶

(e) θ=0°三阶

(f) θ=45°三阶

(g) θ=0°四阶

(h) θ=45°四阶

(i) θ=0°五阶

(j) θ=45°五阶

(k) θ=0°六阶

(l) θ=45°六阶图15 叶片前六阶模态振型Fig.15 First six vibration modes of the blades

图16 不同偏置角度下叶片最大应力和应变

叶片主梁偏置角度变化对叶片内部腹板应力和应变也会产生一定影响。传统叶片腹板最大应力最低，约为71.76 MPa，而θ=30°弯扭耦合叶片的腹板最大应力最高，约为103.73 MPa，较传统叶片增大约44.55%。传统叶片腹板最大应变最低，约为0.004 818 2，而θ=-45°弯扭耦合叶片的腹板最大应变最高，约为0.006 660 5，较传统叶片增大约38.27%。综合分析可知，传统叶片整体结构最大应力及应变均最低，将主梁材料镜像对称铺设后，叶片整体结构最大应力及应变均不同程度增大。

表3和图17给出了各偏置角度下叶片应变能数值及分布。由表3可知，叶片最大应变能及总应变能均随着偏置角度的增大呈先减小后增大的趋势，即将主梁镜像铺设后叶片最大应变能及总应变能均有不同程度提升，弯扭耦合叶片的最大应变能和总应变能较传统叶片最大增幅分别约为185.64%和288.10%。由图17可见，由于主梁为叶片主要承载部件，承担叶片大部分弯曲载荷，叶片应变能主要集中于此处，与叶根相距7～50 m内为应变能最大区域。由叶根至叶尖，叶片总应变能逐渐减小。叶片未产生弯扭耦合特性时(θ=0°)，最大应变能处与叶根相距10.4 m，而叶片产生弯扭耦合特性后，最大应变能处与叶尖更接近。

表3 不同偏置角度下叶片应变能

4.3 屈曲分析

当受载细长结构存在微小横向干扰时，其可能因不稳定而发生大尺度形变失效，即屈曲失稳。通过屈曲分析可确定结构开始不稳定时对应的临界屈曲载荷以及发生屈曲时的模态形状[25]。

屈曲失稳是致使叶片结构强度不足的重要原因之一，研究表明较之叶片强度不足导致的疲劳失效及断裂失效，叶片结构屈曲失稳失效显得更为重要。因此，将CFD方法计算所得叶片表面气动载荷映射至有限元模型，对叶片结构进行屈曲稳定性分析。

表4给出了未耦合的传统叶片前六阶屈曲因子。屈曲因子与初始外载荷相乘即为促使叶片发生屈曲失稳的临界屈曲载荷。由表4可知，屈曲因子随阶次升高逐渐增大，工程中最关注的1阶屈曲因子为2.530 8。根据风力发电机组风轮叶片标准[26]所述，若进行线性稳定性分析时，应对所得结果附加安全系数1.25，而对1阶屈曲因子附加1.25安全系数后其值变为2.024 64，仍大于1，可见传统叶片具有较高的抗屈曲能力，不会发生整体屈曲失稳现象。

(a) θ=-45°

(b) θ=-30°

(d) θ=0°

(e) θ=15°

(f) θ=30°

(g) θ=45°图17 不同偏置角度下叶片应变能分布

表4 传统叶片屈曲因子Tab.4 Buckling factor of the conventional blade

图18给出了弯扭耦合叶片前六阶屈曲因子相对变化量，均以传统叶片屈曲因子为基准。由图18可知，随着偏置角度逐渐增大，各阶屈曲因子均先增大后减小，即将主梁偏轴镜像铺设后，叶片各阶屈曲因子均有不同程度减小，且正向和反向偏置同角度后各阶屈曲因子减小量较为接近。主梁偏置角度变化会对叶片抗屈曲能力产生影响。对比各阶屈曲因子相对变化量可知，同偏置角度下弯扭耦合叶片各阶屈曲因子减小量较为接近，θ=±15°时屈曲因子减小最少(约39%)，而θ=±45°时屈曲因子减小最多(约78%)，此时一阶屈曲因子分别为0.574 87和0.575 54，可见其已低于1，即叶片初始外载荷F0均已达到临界屈曲载荷Fcr。因此，为使叶片产生弯扭耦合特性，叶片屈曲因子必然会有所减小，通过气动弹性剪裁实现弯扭耦合的同时,应关注该叶片是否满足屈曲稳定性设计要求，必要时应优化铺层结构以提高弯扭耦合叶片整体结构抗屈曲能力。