极限载荷下的风力机叶片复合材料损伤过程研究

时间：2024-07-28

马铁强，林耀坤，孙传宗，苏龙，单光坤

(沈阳工业大学机械工程学院，辽宁沈阳 110870)

0 前言

叶片成本约占风电机组总成本的20%。高性价比叶片是风电机组性能和效益的重要保障[1]。为了降低叶片成本、提高叶片的可靠性，玻璃纤维、碳纤维等复合材料被广泛用于叶片的制造领域。复合材料对于实现叶片的高可靠性、低成本和轻量化，一直发挥着至关重要的作用。然而，也有部分资料显示[2]：每年约有3 800余起复合材料叶片失效和故障发生，诱因存在于叶片的加工、运输、安装等环节。此外，长期服役条件下的叶片也会出现疲劳损伤和极限强度破坏，较为典型的损伤包括脱层、裂纹等[3-4]。为了有效降低复合材料叶片损伤和破坏发生的概率，深入掌握风力机叶片材料损伤的机理和损伤发生及演化过程显得十分必要。这将为风力机叶片的形貌、结构和铺层材料的优化设计提供依据、明确方向[5]。

文献[6]采用纤维失效和纤维间失效分析方法，研究了风力机叶片复合材料的分层和基体失效问题，利用FF和IFF分析、估计结构的破坏模式及发生位置。文献[7]分析了风力机叶片的结构失效特性，在复合加载条件下开展了52.5 m叶片结构的倒塌试验，观测了极限载荷作用下叶片损伤全过程，分析了叶片失稳后的失效区域和试验过程中的失效模式。文献[8]分析了临界载荷作用下风力机叶片的损伤行为。以上研究对复合材料风力机叶片损伤问题进行了初步探索。

然而，风力机叶片失效形式复杂，多种失效形式同时发生、彼此复合，从单一尺度难以清楚得出损伤发生机理，有必要从细观到宏观的多个不同尺度联合审视和分析风力机叶片损伤发生的渐进演化过程，才能真正揭示风力机叶片损伤从孕育、发生到发展的过程和机理。文献[9-10]研究表明：高密度的复合材料基体裂纹是分层等严重损伤模式发生的前兆。文献[11]阐述了风力机叶片损伤演化过程，用多尺度渐进损伤模型评估了风力机叶片结构的次临界损伤演化和刚度退化。文献[12]提出了一种多尺度损伤模型，描述了复合材料风力机叶片在准静态和循环载荷作用下的渐进层裂和接头脱粘，可预测叶片的渐进失效和相应的承载能力损失。上述文献在风力机叶片多尺度分析方法方面做了大量探索性研究，但对叶片损伤机理和演化历程缺少更深入的分析和研究。在文献[13-15]研究中，某些复合材料叶片模型仅考虑叶片结构中的损伤起始和演化，而其他模型仅考虑叶片在变形和应力水平方面的整体响应[16-17]。此外，虽然在一些叶片模型中模拟了损伤演化，但它们通常基于失效准则，如蔡吴或蔡希尔，这些准则并未考虑任何物理损伤或相关材料刚度降低[18]。为了准确地预测叶片结构分析中所用材料的非线性响应，必须考虑材料的行为。

本文的风力机叶片复合材料渐进损伤分析新方法，首先定义了材料的本构关系并得到材料损伤后积分点的本构方程，然后基于连续介质损伤力学对HASHIN准则进行修正并定义材料的损伤演化规律，得到渐进损伤模型并在有限元程序中作为用户定义的材料子程序(UMAT)实现。本文以某型号复合材料风力机叶片为例以此方法研究极限风荷载作用下风力机叶片复合材料的力学行为，对叶片材料损伤的演化历程进行数值模拟。

1 渐进损伤分析模型

1.1 风力机叶片复合材料渐进损伤分析流程

首先，目标极限载荷分若干单步按一定增量逐渐增加载荷，每步载荷增量对应的材料状态保持不变。由于风力机叶片由不同材料和铺设角度按一定规律铺设制成故将风力机叶片复合材料视为层合板结构，建立对应的有限元方程并求解。计算得到的有限元模型各单元“应力-应变”状态代入材料失效准则，判断单元是否失效。

如果单元失效，则根据损伤演化规律计算损伤状态，并依据材料退化方案对材料刚度进行退化处理。此时，前一方程组的解已不再满足当前结构的平衡方程，因此根据退化后的材料属性，在相同载荷条件下建立新的有限元平衡方程。

重复应力求解、失效准则判断及材料刚度退化等步骤，直至风力机叶片复合材料结构中不再出现新的损伤。持续增加载荷，重复迭代过程，直至风力机叶片复合材料结构最终失效。在复合材料渐进损伤分析过程中，如果载荷增量足够小，那么损伤处的材料刚度在退化后无需建立新的有限元方程组。图1为具体的叶片复合材料渐进损伤分析方法。

图1 复合材料叶片渐进损伤分析方法

1.2 材料本构关系

复合材料作为正交各向异性材料,其包含三个相互正交的弹性对称面，所以材料有9个独立的弹性常数，材料的应力应变本构方程可以表示如式(1)所示。

[σ]=[C][ε]

(1)

(2)

根据文献[19]提出的渐进损伤模型，在材料发生损伤后，引入损伤矩阵M[D]，则损伤后等效应力可以表示成式(3)。

(3)

(4)

其中，

(5)

(6)

(7)

相应的损伤刚度矩阵为

[Cd]=[M-1][C][MT-1]

(8)

所以可得材料损伤后积分点的本构方程为

[σ]=[Cd][ε]

(9)

1.3 复合材料失效判断准则

研究复合材料风力机叶片失效机理和损伤演化常用的失效判定准则包括最大应力应变、蔡-希尔(Tsai-Hill)、蔡-吴(Tsai-Wu)等失效准则。上述失效准则考虑了不同应力及相互作用的影响，但却缺少对不同失效模式的描述[20]。为了评估风力机叶片复合材料中纤维和基体的损伤情况，本文用HASHIN准则作为损伤起始判据和失效准则。该判据和准则可预测复合材料的层内损伤响应。HASHIN准则涵盖了纤维拉伸失效、纤维压缩失效、基体拉伸失效、基体压缩失效、法向拉伸失效等5种失效模式。

(1)纤维拉伸失效(σ11≥0)。

(10)

(2)纤维压缩失效(σ11<0)。

(11)

(3)基体拉伸失效(σ22>0)。

(12)

(4)基体压缩失效(σ22<0)。

(13)

(5)法向拉伸失效(分层)。

(14)

式(10)～(14)中，σii，σij是单层复合材料各个主方向以及相应面内的剪切应力；Xa，Ya，Za,Sij分别是单层复合材料各主方向的强度，拉伸情况下a为T，压缩情况下a为C；α为材料非线性因子。

1.4 材料损伤演化与退化模型

在满足HASHIN失效准则条件下，需要对材料属性按退化方案进行刚度折减。多数复合材料结构渐进损伤分析模型的材料退化方案是根据失效模式及特点，按比例对材料弹性常数直接进行折减[21-22]。文献[23-24]采用连续介质损伤力学中内部损伤的状态变量表征材料损伤状态，刚度矩阵的退化是基于损伤状态变量对材料单元节点计算得到的。不同状态变量表示不同损伤模式，材料的临界断裂应变能释放率决定其损伤状态变量。

与连续介质损伤力学方法相比，直接刚度折减模型中的材料性能折减系数由试验或经验确定，而有限元模型只能用于特定材料体系和特定结构。此外，材料性能折减系数选用不当，会导致材料退化刚度矩阵奇异，最终将导致数值模拟过程不收敛。为此，本文用连续介质损伤力学方法表征材料退化。

损伤演化过程可描述为应变能释放过程。若材料点应变能释放密度等于断裂能密度，则该材料点失效。该过程中的材料以线性或指数形式软化，宏观上表现为弹性模量退化和承载能力下降。本文假定损伤起始发生后材料满足线性等效“应力—应变”本构关系，如图2所示。图中A、B、C点分别为材料损伤起始点、材料完全失效点、损伤过程中的当前时刻。

图2 材料应力-应变本构关系曲线

本文模型中对材料的本构关系的定义和对材料的损伤演化通过自编用户子程序UMAT实现。UMAT子程序能够实现对单元积分点应力分析、失效判定、损伤状态变量计算、损伤本构建立，以及损伤状态信息的反馈。

近年来，Linde[25]等对损伤后的单元刚度进行非线性退化，考虑了损伤累积对刚度的影响，同时网格划分对计算的收敛性有较大影响，为了减小对网格划分的依赖性引入等效位移，定义为

δeq,iq=εijLC

(15)

定义损伤状态变量变化规律为

(16)

(17)

(18)

2 复合材料风力机叶片模型

本文研究所用叶片为美国国家可再生能源实验室(NREL)公开发表，用于研究比较5MW风机模型所属的叶片如图3所示。该风力机叶片采用DU翼型和双腹板截面支撑结构。如图3b所示，该叶片包括叶根、叶尖、前缘和后缘，内部有主梁、腹板等支撑结构，主梁与腹板采用C型连接，构成支撑结构。市场上梁常用D型、O型、矩形和C型等型式。D型、O型、矩形腹板其特点是重量轻，对叶片运输要求较高。由于叶片前缘强度和刚度较低，在运输过程中局部易于损坏。同时这种叶片整体刚度较低，运行过程中叶片变形较大，必须选择高性能的结构胶，否则极易造成后缘开裂。C型梁其优点是主梁和腹板在前缘粘接部位常重叠，以便增加粘接面积。在后缘粘接缝，由于粘结角的产生而变坚固了。在有扭曲变形时，粘接部分不会产生剪切损坏。叶片整体强度和刚度较大，在运输、使用中安全性好。因此本文应用的双腹板支撑结构具有一定的典型性和代表性。

图3 叶片结构图

2.1 风力机叶片有限元模型

本文建立了风力机叶片翼型截面草图，通过放样建立叶片气动外形的三维几何模型。考虑到叶尖对有限元分析结果影响较小，并且会造成建模过程更复杂，几何建模过程中应去掉叶尖结构。文献[26]指出了该叶片内部组件的形状和位置，其中腹板位于横截面翼弦的15%和50%处；前后腹板之间为主梁。主梁是叶片的主要支撑结构，承受叶片所受到的主要弯曲荷载。为了正确模拟复合材料风力机叶片的力学性能，用壳单元对叶片的几何模型进行有限元分网，同时假定壳单元呈平面应力状态。图4所示为风力机叶片在某截面处的有限元网格模型。该网格采用了壳单元，因此上壳和下壳连接边缘被视为永久粘合的界面。

图4 宏观尺度叶片截面分析模型

2.2 材料参数与层压策略

风力机叶片的蒙皮和翼梁通常用多个不同的层压材料和夹层板制造。为了简化风力机叶片模型，假设蒙皮、前腹板和后腹板均采用由玻璃纤维与芯材按照[45°/0°/-45°/C16]s铺层增强策略组成的夹芯结构。根据表1[13]推荐的材料，定义风力机叶片有限元模型中的材料属性。如表1所示，复合材料有10种材料属性，例如弹性模量、泊松比、剪切模量、密度等。风力机叶片共计用到了混合纤维材料(Hybrid Fiber Mat)、玻璃纤维(GFRP)、聚乙烯(Polyethylene)、梁帽混合物(Spar Cap Mixture)等4种材料，分别用于叶片外侧铺层材料、增强纤维材料、铺层夹芯材料和腹板夹芯材料。

表1 材料及其属性

3 极限载荷与加载

3.1 极限载荷

标准工况条件下，空气动力、重力和离心力是风力机叶片的主要载荷来源。为了计算和分析简便，此处忽略离心力及其惯性力影响，空气动力成为引起风力机叶片发生变形、损伤和破坏的主要外力来源。沿叶片展向各叶素上所受的空气动力可由叶素-动量理论(BEM)[27-28]评估和计算。

叶素周围气体流动速度W可分解为叶片挥舞方向分量vrx和叶片摆振方向分量vry。即风速可w表示为

(19)

式中，v为风轮远前方气流速度,m/s；vr为风轮处气流速度,m/s；a为轴向诱导因子；a′为周向诱导因子；Ω为风力机旋转角速度,rad/s；r为叶素旋转形成的圆环半径,m。

综合动量定理和叶素理论，可求出周向诱导因子a和周向诱导因子a′，进而在叶片展向方向上，按照公式(20)～(23)对各叶素所受的空气动力和动力矩进行积分，求出叶片所受到的空气动力载荷。

(20)

(21)

(22)

(23)

式中，R0为轮毂半径，m；ρ为空气密度，kg/m3；t为叶素的弦长，m；ct为切向力系数；cn为法向力系数；FN为叶片挥舞方向作用力，N；FT为叶片摆振方向作用力，N；MN为叶片挥舞方向力矩，N·m；MT为叶片摆振方向力矩，N·m。

3.2 极限载荷加载

该风力机叶片由DU30、DU40、DU93等3种翼型构成，并按照翼型的展向分布可将叶片分为不同分段。各分段处的弦长不同，所受载荷也有较大差异。本文所述风力机叶片沿展向共分7段，各段分别施加载荷。DU40翼型、DU30翼型分别占据风力机叶片上一段，DU93翼型占据5段，如图5a所示。风力机叶片各段单元节点上施加的风载荷及其分布情况如图5b所示。为了研究极限风载荷作用下的风力机叶片复合材料损伤演化过程，按时间历程改变逐渐增加风载荷。假定时间间隔为1 s，作用于风力机叶片上的风速从1 m/s到22 m/s逐渐加大。

图5 叶片加载部位示意图

4 分析与讨论

仍以典型的5 MW复合材料风力机叶片为例，并按本文方法进行非线性有限元建模、风况加载和应力计算，得出不同风况条件下复合材料风力机叶片的应力分布云图。分析可知复合材料风力机叶片中纤维和基体的损伤与风速及极限风载荷之间的关系，并分析了复合材料风力机叶片中纤维和基体的损伤形成机理。

分析可知，空气动力使叶片弯曲和扭转；离心力使叶片承受拉伸、弯曲和扭转；重力使叶片承受拉压、弯曲和扭转。叶片的前缘采用了0°、﹢45°、-45°三种角度铺层。0°铺层主要支撑离心力和气动力产生的叶片轴向拉压应力。﹢45°和-45°铺层承受扭矩产生的剪应力。

4.1 纤维损伤演化过程

图6为叶片不同位置的纤维损伤与风速之间关系。本文定义了损伤变量的含义，损伤变量值为 0 代表材料未受损伤，值为 1 时候代表材料发生失效，依据HASHIN准则，图6中数值达到1时，即表明纤维发生损伤。结合图7所示有限元模型，分析可知复合材料叶片各部位纤维损伤在细观层面的演化和发展过程。

图6 不同位置纤维损伤单元随风速变化图

图7a～图7d是叶片纤维材料损伤变量云图。如图7a所示，当风速达到14 m/s时，从叶跟向叶片过渡的后缘区域纤维最先出现损伤。如图7b所示，当风速达到16 m/s时，叶根区域纤维出现轻微损伤，在叶根和过渡区域的主梁也出现损伤。随风速增大，叶根与后缘的损伤开始延翼展方向延伸，并逐渐连接。图7c中当风速达到19 m/s时，后缘表面出现新损伤。

图7 叶片纤维损伤与风速关系图

如图8所示，主梁与后缘的过渡区域在主应力方向出现应力集中，由于主梁较后缘更厚、强度更高。因此，主梁与后缘过渡区域的纤维易出现损伤如图7c所示。此处纤维损伤将随风速增大，逐渐沿着翼展方向扩展到叶尖如图7d 所示。

图8 风速为19 m/s时，叶片过渡区域应力图

图9为纤维损伤处单元应力随风速变化关系图。由图9可知，不同位置纤维应力变化趋势大致相同。后缘过渡区域的前半程应力值始终保持最大，并最先达到最大；但当风速达到14 m/s时，却出现骤降。此时，后缘过渡区域纤维发生损伤、断裂，并失去承载能力，单元应力值迅速下降。叶根和后缘中段区域的纤维也将发生损伤、断裂，失去承载能力，应力曲线下降。值得注意的是，某一单元应力消失即纤维损伤失效，势必造成周围单元应力增加，从而引起相应纤维的加速损伤和失效。

图9 不同位置纤维损伤单元应力与风速关系图

4.2 基体损伤演化过程

图10为不同位置基体损伤与风速关系图。分析图11可知，基体损伤首先发生在叶根处，并沿着翼展方向逐渐扩展，与文献[12]的研究结论基本一致。

图10 不同位置基体损伤单元随风速变化图

如图11a 所示，当风速达到12 m/s时，叶根处基体出现轻微损伤。这意味着叶片根部应力大于其他区域应力。如图11b所示，当风速达到14 m/s时，后缘过渡和前缘区域处发现轻微的基体损伤。如图11c所示，当风速增大到15 m/s时，叶片后缘中段发现轻微的基体损伤。如图11d所示，当风速增大到19 m/s时，叶片中段的基体损伤开始分别向叶尖、叶根逐步扩展。

图11 叶片基体损伤与风速关系图

在基体损伤演化过程中，当风速达到14 m/s时，基体损伤开始逐渐加重。这一现象是由纤维断裂后失去承载能力所引起。纤维的断裂失效会导致周围基体材料承受更大的载荷而出现应力集中，从而加速基体损伤失效过程。

图12为基体损伤处单元应力随风速的变化关系图。

图12 不同位置基体损伤单元应力与风速关系图

图12中，后缘和叶根的应力曲线与纤维损伤的曲线特性一致，纤维损伤后应力曲线也相应下降。但前缘和后缘中段两条曲线与叶根和后缘过渡区域两条曲线不同，在风速增大过程中出现一定程度的应力下降。分析可知，该现象是由于基体损伤并快速扩散所导致，形成曲线上第一次下降，随着风速增大单元应力再次增大，而当风速达到发生纤维损伤值时应力曲线骤降形成第二次下降。