基于参数建模的风力机预弯叶片结构优化设计*

陈进 郭小锋 李松林 孙振业

(重庆大学 机械传动国家重点实验室，重庆400044)

叶片是风力机最主要的部件之一，为了减小叶片的重量，当前多对大型叶片进行预弯设计，以避免因叶片柔性变形后与塔筒发生碰撞.对叶片进行预弯设计的目的是在满足材料强度的前提下，增大叶片的柔性变形许可，减轻叶片的重量，因此，弯曲型线与叶片的铺层结构存在耦合关系，如何在叶片的结构设计中实现预弯型线与结构铺层的协同是一个十分有意义的课题.

在预弯叶片研究方面，LM 公司最早提出了预弯设计的概念，并申请了发明专利［1］，Riziotis 等［2-3］对预弯叶片的气动性能和气弹性能进行了研究，提出了风力机预弯叶片的气弹分析方法.但在预弯叶片的弯曲型线表达及其与叶片铺层结构的协同设计方面，未见相关研究文献;在对非预弯叶片的复合材料铺层结构优化设计方面，Ashuri 等［4-6］将叶片等效为悬臂梁，运用材料力学方法近似计算叶片的强度，提供了一种简单实用的叶片结构设计方法.秦志文等［7-9］对风力机叶片的有限元建模和加载方法进行了研究，表明有限元分析方法较工程计算方法有更高的计算精度.还有一些学者［10-11］在叶片的结构设计方面也做出了贡献，值得借鉴.Chen等［12］直接采用有限元建模分析方法对叶片的铺层结构进行优化设计，在对叶片各部位铺层顺序和层数进行初步设计之后，以单层玻璃纤维布厚度的改变来等效实现叶片各部位铺层厚度的变化，取得了较好的优化结果.然而在实际生产中，单层玻璃纤维布的厚度是固定的，以单层厚度的改变来等效实现铺层层数的变化会影响有限元分析结果的准确性，并且该优化设计方法需要预先确定好叶片各部位的铺层层数，优化结果的好坏依赖于前期的铺层结构设计.因此，对叶片的参数化建模分析方法还有待于进一步研究.

为解决以上问题，文中提出了一种预弯叶片弯曲型线的表达方法;研究叶片铺层结构设计模型，提出以叶片展向关键点处的铺层层数为变量来表征叶片的整体铺层结构;对某850 kW 叶片进行参数化建模和有限元分析，验证叶片参数化有限元建模分析方法的正确性;构建叶片预弯型线和铺层结构的一体化优化设计模型，对850 kW 叶片进行预弯和铺层结构优化设计，并给出了优化结果和评价.

1 叶片的铺层结构设计

风力机叶片由玻璃纤维布和夹层泡沫材料铺设而成，为了实现对叶片铺层层数的优化设计，首先需对叶片的铺层结构设计模型进行研究，以实现用若干个变量来表征叶片的整体铺层结构.典型的叶片结构模型如图1所示，叶片截面由前缘加强面板(LRP)、前缘面板(LEP)、梁冒面板(CAP)、后缘面板(TEP)、后缘加强面板(TRP)和大小腹板(WEB)7 个区域组成，前缘面板、后缘面板的铺层结构如图1(c)，梁冒的铺层结构如图1(d)，大小腹板的铺层结构如图1(e)，前缘加强和后缘加强的铺层结构如图1(f).

图1 叶片结构模型Fig.1 Structural modal of the blade

叶片截面各区域由不同类型的层合板组成，从结构形态上，可将其划分为6 种类型，分别为壳体层、翼面加强层、单向布加强层、翼面泡沫和腹板泡沫，其功能和形态可描述如下:

(1)壳体层 由两层交叉铺设的双轴布材料组成(叶片铺层的0 度方向为叶片展长方向)，主要功能是为了保证叶片外形的光滑连贯性，在叶片的中后部位也起到保证叶片强度的作用;

(2)翼面加强层 由若干个子层合板组成，每个子层合板又由一层双轴布材料和一层单轴布材料构成，每相邻子层合板中的双轴布材料铺层方向成交叉布置，翼面加强层的主要功能是增强前缘和后缘的强度;

(3)单向布加强层 由若干层单向布材料组成，铺设在梁冒的主要功能是抵抗叶片的挥舞弯矩，铺设在叶片前后缘的功能是抵抗叶片的摆振弯矩、增加前后缘强度和防止胶接脱裂的作用;

(4)翼面泡沫 主要功能是增加叶片的厚度，防止叶片表面的屈曲变形;

(5)腹板泡沫 一方面用于防止腹板发生屈曲变形，另一方面也用于抵抗叶片截面受到的剪切力;

(6)叶根增强层 由若干层三轴向玻璃布材料组成，用于增强叶片根部的强度.

叶片设计中所用到的铺层材料如表1所示.叶片根部的铺层结构由叶片外表面至内依次为:外壳体层、由32 层双轴布和32 层单轴布材料交替铺设而成的翼面加强层、由32 层单向布铺设而成的单向布加强层和内壳体层.距叶根0.55 ～26.75 m 之间的截面铺层结构如图1所示，沿叶片展向，叶片翼面加强层和梁冒单向布加强层的铺层层数如图2所示.根据其铺层规律，设置若干个关键点，将关键点处的铺层层数作为设计变量，叶片其他位置处的铺层层数可通过线性插值的方式得到.

为了增加叶片前后缘的强度和防止翼面胶接处脱裂，需要对叶片的前后缘进行加强.在叶片的过渡段(距叶根1.55 ～7.25 m 之间的截段)，前后缘处翼面加强层和单向布加强层的铺层层数为4层，在其它截段，前后缘处翼面加强层和单向布加强层的层数为2 层.文中在优化设计时，前后缘加强层的层数仍采用原叶片的铺层层数设计，不对其进行优化.

图2 850 kW 叶片铺层结构及变量设置Fig.2 Lay-up of 850 kW blade and setting of design variables

2 预弯叶片的参数化建模

2.1 预弯叶片的参数化表达

叶片外形是由若干个截面扫掠而成的，对叶片各个截面的空间坐标进行参数化转换和表达就能实现对整个叶片的参数化表达，叶片各个截面坐标的转换过程如图3所示.

图3 叶片截面坐标变换Fig.3 Coordinate transformation of blade section

在图3中，翼型1 为标准翼型，翼型1 经过气动中心变换后得到翼型2，翼型2 经过弦长变换后得到翼型3，翼型3 经过扭角变换后得到翼型4，翼型4经过预弯变换后得到翼型5，其变换过程为［13］

式中，(X0，Y0，Z0)为基本翼型坐标，(Xi，Yi，Zi)为变换后的叶片截面坐标，T1、T2、T3、T4、T5为变换矩阵.

2.2 预弯型线设计

文中提出了一种简单实用的预弯型线函数表达方法，该型线函数在预弯段为幂函数，由幂函数的特性可知，预弯型线函数在弯曲起始点处与x 轴相切，因此能保证叶片在预弯过渡段处的平滑过渡.此外，采用这种型线函数还可以减少优化设计中设计变量的个数.函数表达式为

式中，A 为幂函数预弯型线的指数，R 为叶片的长度，B 为叶尖的预弯量，b(x)为叶片各个截面外形在图3中由翼型4 到翼型5 的纵向偏移量，x 为叶片展长.预弯变形起始于距离叶片根部R/3 处，通过改变式(2)中A 和B 的值就可以得到不同的预弯型线.对于文中作为实例研究的850 kW 叶片，当A =3，B =1时，叶片的预弯型线和有限元模型如图4所示.

图4 预弯型线及预弯后的850 kW 叶片Fig.4 Pre-bending profile and 850 kW blade model

2.3 有限元建模

为了满足叶片结构优化设计的需要，文中采用Matlab 与Ansys APDL 语言相结合的方法编写程序建立叶片参数化有限元模型，构建的850 kW 叶片的有限元模型如图5所示.

图5 850 kW 叶片模型Fig.5 850 kW wind blade model

叶片参数化建模和分析的方法可概括如下:

(1)基于叶片外形的参数化表达方法，准备叶片外形参数的文本文件;

(2)根据叶片展向关键点处的铺层层数，插值得到叶片各个构成面上的铺层层数，用Matlab 程序自动生成叶片铺层结构的APDL 程序文本文件;

(3)编写叶片建模的主程序，包含叶片外形建模、复合材料铺层和结构分析3 个子模块.在叶片外形建模子模块中，调用步骤(1)中生成的叶片外形参数文本文件生成叶片模型的线和面;在复合材料铺层模块中，读取并执行步骤(2)中生成的铺层结构APDL 程序文本文件，对叶片各区域面赋予材料并对其进行网格划分;结构分析子模块实现对叶片强度、变形、重量和自振频率等结构参数的计算和输出.

2.4 载荷计算与加载

风力机叶片在实际运行中，要经历各种复杂的工况，为了保证叶片各个部位在这些工况下都不会发生结构破坏，在叶片的强度校核时需要让叶片每个截面都承受所有计算工况下的最大载荷.为此，文中依据GL 2010 标准［14］将风力机运行工况细分为240 种典型工况，运用风力机载荷计算软件GHBladed 计算得到叶片各个截面在240 种典型工况下的极限载荷，将极限载荷施加在有限元模型的节点上，以分析叶片的最大应力和最大变形.图6为计算得到的850 kW 叶片极限弯矩载荷分布.在对叶片的结构强度计算时，分别考虑挥舞方向载荷最大和摆振方向载荷最大两种情况.计算得到的极限载荷为集中载荷，在有限元分析时需要选用适当的加载方式以避免应力集中.文献［15］中在对叶片进行静载分析时，将相邻两个截面上的弯矩增量平均加载在两个截面之间的所有节点上，但作者在叶片建模分析时发现，这种加载方式会产生局部应力集中，导致在叶片根部过渡段处应力明显偏大.文中在综合文献［15］和现行叶片试验测试方法的基础上，采用如下加载方式:将叶片每两个相邻截面上的挥舞弯矩增量平均加载在两个截面之间的梁冒节点上，以模拟分析叶片在挥舞方向上的极限载荷;将50 kN的集中力平均加载在距离叶根19.75 ～20.75 m 截段的梁冒节点上，以模拟叶片根部及过渡段在摆振方向承受的极限载荷，该载荷的施加方向为使叶片朝向尾缘摆振变形的方向.

图6 叶片极限弯矩分分布Fig.6 Distribution of ultimate bending moment of the blade

2.5 有限元模型的验证

在叶片的有限元模型中，采用shell181 单元对叶片的铺层结构进行模拟，为了验证有限元模型的正确性，对850 kW 叶片进行有限元建模和分析，并与试验报告数据进行对比.原叶片在申请GL 认证时参照GL 2010 标准中的实验规范进行了静力实验［14］，其中在校核最大挥舞弯矩载荷下叶片的强度时，在距叶根20.25 m 处加载121 kN 的集中力，载荷施加方向为使叶片产生挥舞变形、叶片压力面受拉伸的方向.文中在验证有限元模型时，为避免应力集中，将121 kN 的集中力平均加载在距离叶根19.75 ～20.75 m 叶片截段的梁冒上，计算叶片叶尖的最大变形，此外，也对模型叶片进行了质量求解和模态分析，计算和对比结果如表2所示.从中可以看出，叶片有限元模型的分析结果与试验结果比较接近，各项指标相差最大不超过2%，误差的产生可能是有限元建模时对叶片的结构铺层进行了一定的简化，与叶片实际铺层不完全相同所致.误差控制在比较小的范围内，验证了采用Matlab 编程结合Ansys 二次开发建立叶片有限元模型对叶片进行结构分析的可行性.此外，文中也对有限元模型的网格无相关性进行了分析，最终确定叶片有限单元的个数为24 498.

表2 有限元模型分析与测试值的对比Table2 Comparison between the results obtained by FEM modal analysis and test

3 优化设计

3.1 设计变量

按照图2中叶片铺层关键点的设置，对叶片关键点处梁冒单向层的铺层层数和翼面加强层的层数进行优化，设置设计变量如下:距离叶根1.75、2.75、3.75 和4.75 m 处翼面加强层的层数分别为变量x1、x2、x3、x4;距离叶根1.75、2.75、4.75 和26.75 m处梁冒单向层的层数分别为变量x5、x6、x7、x8.变量x1-x8的变化范围为3 ～45.

选取式(2)中的指数A 和叶尖预弯量B 为叶片预弯型线变量参数.A 的变化范围为1.5 ～4.0，B 的变化范围为0.8 ～1.8.

3.2 设计约束

(1)设计叶片的强度校核以叶片最外层双轴布材料不失效为设计准则.在求得叶片表面所有单元的拉应力、压应力和剪切应力后，运用Tsai-Wu 失效准则式(3)对叶片每个单元的强度进行校核，并统计不满足约束条件的单元个数，由式(4)计算不满足强度设计约束的惩罚因子P1:

式(3)、(4)中，F1、F2、F11、F12、F22、F66为复合材料的强度因子;σ1、σ2、σ6为叶片有限单元应力，N 为不满足Tsai-Wu 失效准则的单元个数.

式(3)中各强度因子的计算依赖于铺层材料的许用强度和各单元处的应力值.对于铺设在叶片最外层的双轴布，其许用强度的计算方法为

式中，UTS、UCS 和USS 分别为复合材料的极限拉伸、压缩和剪切强度，ATS、ACS 和ASS 分别为复合材料的许用应力，γm为安全因子，其计算公式为

式中，C1a为材料的安全系数，C2a为所有材料的蠕变强度因子，C3a为所有材料的温度影响因子，C4a为预浸或手敷铺层，C5a为控制固化层压板的热处理系数.

对本优化设计叶片，计算得到设计安全因子γm为2.45，双轴布材料的许用强度如表3所示.

表3 双轴布材料的许用强度Table3 Allowable strength of the biaxial material

(2)为了避免共振，考虑原叶片的结构动力学特性，在优化时，要求设计叶片的一阶挥舞和一阶摆振频率不能在以下范围内:

若叶片一阶挥舞频率ω1f和1 阶摆振频率ω1e不在式(7)范围内，由式(8)计算不满足结构动力学要求的惩罚因子P2.

(3)叶片通过轮毂固定在塔架上，为了避免叶片在极限工况载荷下与塔筒发生碰撞，在优化设计时要求叶尖最大变形量不超过叶片未变形时叶尖与塔架距离的70%［14］.叶片在静止状态下，叶尖与塔筒壁之间的距离xtip，max为

式中，xhub为轮毂中心到塔筒中心的距离，Rrotor为叶片的长度，αtitle为风轮安装倾角，αcone为叶片锥角，dtower，tip为叶尖相应高度处塔筒的直径，dtip，max为叶尖的最大变形量.对于原850 kW 叶片，叶尖与塔筒之间的距离为4.731 m，计算得到叶尖最大变形量为3.312 m.在考虑叶片预弯后，假设叶片的叶尖预弯量为B，优化设计中叶片的最大变形量约束dtip，max为

若叶片的最大变形量超过式(11)的限制，由式(10)计算不满足最大变形要求的惩罚因子P3.

(4)为避免优化进程的局部收敛，假定优化设计的理想目标为2 350 kg，对大于2 350 kg 的计算适应度值，计算未达到理性质量要求的惩罚因子P4:

式中，mblade为叶片的实际质量.

3.3 优化目标

在叶片外形一定的情况下，结构铺层设计中所使用的复合材料质量越小越有利于减轻叶片的制造成本和载荷.因此，在叶片的优化设计中，以复合材料叶片的质量最小为优化目标.

式中，P1-P4为计算适应度值的惩罚因子，F(x)为计算适应度值.

3.4 优化算法

图7为基于改进粒子群算法的叶片优化设计流程图，对于每一个个体，在解码得到叶片设计变量参数后，生成叶片有限元建模和结构分析的APDL 程序文件;主优化程序调用Ansys 软件执行生成的叶片APDL 程序文件，对叶片进行结构特性分析;读取计算结果文件并计算个体的计算适应度值.优化设计中，在每一代个体计算完成后，对粒子群算法的计算参数进行自适应调整，以避免优化程序的过早收敛.文中在优化设计中，对粒子群算法中的主要参数作如下设置:惯性权重w 取0.9，学习因子C1、C2均取0.5，变量维数为10，种群大小为20.

图7 优化算法流程图Fig.7 Flow chat of the design method for airfoil

4 优化结果

850 kW 叶片的优化设计进程如图8所示.在优化初期，因惩罚因子的引入，计算适应度值远大于叶片初始质量，经过80 代的计算，计算适应度值逐渐逼近于叶片实际质量，但因为仍小于目标预期值2350 kg而引入惩罚因子P4，计算适应度值的收敛值小于叶片实际质量的收敛值.整个计算过程耗时20.3 h.

优化叶片预弯型线的幂函数指数A 为2.22，叶尖预弯量为1.474m.优化叶片与原叶片梁冒加强层层数和翼面加强层层数的对比如图9所示.

从图9中可以看出，沿叶片展向梁冒的铺层层数和叶根处翼面加强层的层数均有一定的减少，表明对原叶片进行预弯设计后原有的铺层结构过于保守，必须对其进行优化设计.

图8 进化历程图Fig.8 Iterative course for the optimization

图9 优化前后铺层层数的对比Fig.9 Layers comparison of the original and the optimized blades

叶片的质量、叶尖最大变形、自振频率的对比情况如表4所示.从表4中可以看出，优化后叶片的质量为2539.2 kg，比原叶片减轻了15.26%.原叶片的最大蔡武失效因子为0.890，发生在叶根过渡处，优化叶片的最大蔡武失效因子为0.997，发生在叶根过渡段和叶片中后段梁冒位置.优化叶片的一阶挥舞频率和一阶摆振频率满足设计要求，表明优化叶片具有较好的结构动力学性能.

表4 优化结果Table4 Optimization results

图10为优化叶片沿展长方向的应力分布图，对比分析原叶片和优化叶片的应力分布情况，沿叶片展向的压应力为两只叶片的主要应力.原叶片的最大应力值为72.4MPa，发生在叶根过渡段处，叶片中后段(距叶根14 ～26.75m)的最大应力值只有53.6 MPa，说明原叶片的铺层设计在叶片中后段有较高的安全冗余度，在中后段铺设冗余的梁帽铺层是为了防止叶片变形过大.优化叶片的最大应力为81.6MPa，同时发生在叶根过渡段和叶片中厚段梁帽部位，表明优化设计在减轻叶根过渡段质量的同时，对叶片中后段的铺层结构也进行了减重优化.原叶片为非预弯叶片，优化叶片因为进行了预弯设计，所以在满足强度的条件下可以容许更大的变形，其叶尖最大变形量为4.786 m，满足设计要求.

图10 优化叶片的应力分布图Fig.10 Stress distribution of the optimized blade

5 结语

文中提出了一种预弯叶片弯曲型线的参数化表达方法;构建了叶片的铺层结构设计模型，采用Matlab 和Ansys APDL 语言编程实现叶片的参数化建模;实例验证了复杂叶片的参数化有限元建模方法;提出了一种叶片预弯型线与铺层结构的一体化优化设计方法，以叶片展向关键点处翼面加强层层数和梁帽单向布加强层层数为变量参数，以叶片的质量最小为优化目标，对某850 kW 叶片进行了优化设计，有效地减轻了叶片的重量.用经过验证的有限元建模分析方法对叶片进行计算和分析，所得到优化结果具有较高的可信度.对于大型风力机叶片，在完成气动外形优化设计之后，可运用文中提出的设计方法对叶片的预弯型线和铺层结构进行一体化设计和优化.本研究为后续预弯叶片的气动、结构和气弹多学科优化设计奠定了坚实的基础.

［1］LM 玻璃纤维有限公司.风车转子和用于该风车转子的机翼型叶片:中国，98808849.5［P］.2000-10-11.

［2］Riziotis，Voutsinas，Manolas，et al.Aeroelastic analysis of precurved rotor blades ［C］∥Proceedings of EWEC 2010.Warsaw:EWEC，2010:3-10.

［3］郭婷婷，吴殿文，张强，等.风力机叶片预弯技术的数值模拟［J］.太阳能学报，2011，32(7):1011-1025.Guo Ting-ting，Wu Dian-wen，Zhang Qiang，et al.Numerical researches of the wind turbine blade’s preflex tech-nology［J］.Acta Energiae Solaris Sinica，2011，32(7):1011-1025.

［4］Ashuri T，Zaaijer M B，Martins J R R A，et al.Multidisplinary design optimization of offshore wind turbines for minimum levelized cost of energy［J］.Renewable Energy，2014，68(2014):893-905.

［5］王珑，王同光，吴江海，等.基于改进NSGA-Ⅱ算法的风力机叶片多目标优化设计［J］.南京航空航天大学学报，2011，43(5):672-676.Wang Long，Wang Tong-guang，Wu Jiang-hai，et al.Multiobjective optimization design of wind turbine blade based on improved NSGA-Ⅱ［J］.Journal of Nanjing University of Aeronautics ＆ Astronautics，2011，43(5):672-676.

［6］蔡新，朱杰，潘盼.水平轴风力机叶片最优体型设计［J］.工程力学，2013，30(2):477-480.Cai Xin，Zhu Jie，Pan Pan.The best shape design of horizontial axis wind turbine blade［J］.Engineering Mechanics，2013，30(2):477-480.

［7］秦志文，徐宇，杨科，等.风力机叶片FRP 铺层结构强度与失效预测［J］.工程热物理学报，2012，33(4):591-594.Qin Zhi-wen，Xu Yu，Yang Ke，at al.Structural strength and failure prediction on FRP laminate of wind turbine blade［J］.Journalof Engineering Thermophysics，2012，33(4):591-594.

［8］李慧，王同光.风力机叶片分布力函数加载方式有限元分析［J］.太阳能学报，2013，34(3):485-489.Li Hui，Wang Tong-guang.FEM simulation of distributed force function loading applied on wind turbine blade［J］.Acta Energiae Solaris Sinica，2013，34(3):485-489.

［9］宗楠楠，董湘怀.小型风力机叶片强度的有限元分析［J］.太阳能学报，2010，31(6):764-768.Zong Nan-nan，Dong Xiang-hui.Finite element analysis of blade strength for small wind turbine［J］.Acta Energiae Solaris Sinica，2010，31(6):764-768.

［10］Todd Griffith D，Ashwill Thomas D.The sandia 100-m all-glass baseline wind turbine blade［R］.New Mexico:Sandia National Laboratories，2004.

［11］Kong C，Bang J，Sugiyama Y.Structural investigation of composite wind turbine blade considering various load cases and fatigue life［J］.Energy，2005，30(11/12):2101-2114.

［12］Chen Jin，Wang Quan，Shen Wen-Zhong，et al.Structural optimization study of composite wind turbine blade［J］.Materials and Design，2013，46(2013):247-255.

［13］包飞.风力机叶片几何设计与空气动力学仿真［D］.大连:大连理工大学机械工程学院，2009.

［14］GL 2010，Guideline for the certification of wind turbines［S］.

［15］王锋，陈坤.CGI56-850 kW 聚酯叶片设计与强度分析［R］.绵阳:中国空气动力研究与发展中心计算空气动力研究所，2011:1-35.