随着风力发电行业的迅速发展，叶片长度不断增加，结构设计安全余量逐步降低。叶片作为核心的捕风装置，其结构安全稳定是风机正常运行的前提。根据叶片的结构特点，屈曲已成为叶片极限工况下的主要的失效形式[1,2]。由于全尺寸试验检验叶片的承载能力与结构可靠性成本过高，故基于有限元的仿真分析一直是验证叶片结构稳定性的主要手段。

 

风力叶片研发伊始，研究者基于有限元方法针对叶片屈曲失效进行了大量的工作。薛彩虹等[3]基于 Abaqus 采用线性特征值屈曲方法分析了 40m 级叶片的稳定性；李成友等[4]应用 Nastran软件中特有的特征值抽取算法判断 34m 级叶片的临界失稳点，并计算出屈曲因子；张立[5]等基于 ANSYS 软件研究了 NREL 5MW 风力机叶片承受弯扭耦合载荷作用下的线性屈曲特征。不难发现，大部分研究都基于特征值屈曲的分析方法，特征值屈曲以小位移、小应变的线弹性理论为基础，忽略了结构在受载变形过程中结构构形的变化[6,7]，而叶片结构复杂、铺层繁多、变形很大，基于线性方法难免存在误差，故 IEC61400-5 和 DNV-0376 均提出了基于非线性的叶片屈曲分析方法，以准确地预测叶片的承载能力，进而指导叶片结构铺层设计[8]。在大兆瓦、小节圆、长叶片设计成为当前研发主流的情况下，非线性屈曲非常具有工程应用价值。

 

基于上述讨论，本文应用 Abaqus 有限元软件对时代新材 90m 级叶片进行屈曲分析,采用线性、非线性多种方法对大挥工况的临界载荷进行求解，并对两种方法比较分析。

 

2. 有限元模型

时代新材某 90m 级叶片为双腹板结构，如图 2-1 所示的叶片截面，蒙皮、主梁是叶片最主要的部分，主梁承担叶片绝大部分的弯曲载荷由，而蒙皮起气动作用且承担部分弯曲载荷及大部分的剪切载荷，其中弯曲载荷大约为总载荷的 20% ～ 30%。

 

 

图 2-1 90m 级叶片截面图

 

2.1 叶片有限元模型

基于外部气动外形的输入，应用 Foucs 进行模型建立，依据铺层设计划分结构区域，并划分网格，网格尺寸为 50mm*50mm，如图 2-2 所示。将模型导入 Abaqus 有限元软件中，应用python 与 abaqus 的接口进行编程并对叶片模型进行铺层，其中应用的为 S4r 壳体单元以及C3D8R实体单元。

 

 

图 2-2 90m 级叶片有限元模型

 

2.2 边界条件

在实际运行中叶片与轮毂通过螺栓相连，结构分析需将其简化为悬臂梁模型，故对叶根部进行固定全约束。对于叶片的载荷形式，根据整机气动模拟给到的弯矩载荷等效转换为集中力，施加在叶片主梁位置，如图 2-3 所示。

 

 

图 2-3 叶片有限元模型边界条件

 

3. 屈曲分析方法

屈曲分析主要用于研究结构在特定载荷下的稳定性，以及确定结构失稳的临界载荷，它包括线性屈曲和非线性屈曲分析。屈曲稳定性分析的意义就是为进一步优化结构提供依据。

 

3.1 特征值屈曲分析法

其中线性特征值屈曲分析法，它适用于对理想弹性结构的理想屈曲强度进行预测，主要是使用特征值公式计算造成结构负刚度的应力刚度阵的比例因子。结构在达到屈曲载荷之前其位移—变形曲线表现出线性关系，达到屈曲以后曲线将跟随另外的路径。发生转折的这一点称为分支点，分支点的载荷称为屈曲载荷，结构发生屈曲的不同形态成为屈曲模态。

 

 

 

入即上述方程的最小特征值解，代表临界载荷比例因子，最小特征值乘以初始载荷得到临界载荷，即抗屈曲值。

 

3.1 非线性 Riks 分析法

Riks 是求解非线性屈曲应用较为广泛的方法，其基本思想是假定载荷沿弧长方向增加或减小，Abaqus 中提供了修正的 Riks 算法，通过沿切线移动到给定距离到当前解点限制增量大小，进而结合获得的点且与同一切线正交的平面中进行搜索取得平衡[7]，具体如下：

 

以Pⁿ表示加载模式，N表示模型的自由度，入载荷放大因子，故模型实时运算的载荷为

入Pⁿ为对应的位移，正则化为uⁿ，实时速度uo，同理，以矢量 (uⁿ；入)表达当前弧长步计算结果

 

 

图 3-1 修正的 Riks 法位移-载荷示意曲线

 

 

 

其中， Abaqus/Standard 自动负载增量算法可根据收敛率针对静态问题进行调整，之后，同一弧长步内，取得当前弧长步最后的解并进行迭代，进而得到最后结果。

 

4. 屈曲仿真结果与分析

以大挥方向的工况作为研究对象，首先基于 Abaqus 的 Buckle 求解整机厂输入的载荷进行特征值屈曲分析，求解结果如图 4-1 和图 4-2所示，在这里列举了一阶模态及二阶模态的屈曲线性计算结果，一阶模态屈曲特征值为 1.978，二阶模态屈曲特征值为 1.997，均在 1.965以上，满足 GL2015 及 IEC61400-5 标准，也就证明叶片铺层结构对该载荷结构相对稳定，但值得注意的线性分析结果表面相对容易出现屈曲的位置靠近叶尖。

 

 

图 4-1 叶片一阶模态特征值屈曲

 

 

图 4-2 叶片二阶模态特征值屈曲

 

之后，应用非线线性弧长法进行计算，考虑到叶片由于生产制造过程中可能存在的误差，引入几何缺陷，应用 abaqus关键字语句*IMPERFECTION,FILE=Job-name,STEP=1 将特征值计算得到前两阶模态的节点扰动引入到修正的 Riks后屈曲计算中，依据叶片设计 IEC61400-5标准对缺陷进行缩放。非线性屈曲求解结果如下，当载荷放大因子为时间载荷的1.653倍时停止计算，图 4-3 和 4-4 为计算结束后的应力及应变云图。

 

 

图 4-3 修正的 Riks 求解屈曲应力云图

 

 

图 4-4 修正的 Riks 求解屈曲应变云图

 

屈曲位置为 PS 面靠近叶根整体叶片长度 5.2%处，由压应力主导致使该位置出现局部屈曲进而导致整体屈曲。

 

我们选取叶片处展向长度 15.2%处及展向 30.4%处绘制其载荷位移曲线，如图 4-5 所示，在OA 及 OC 段载荷系数较低，叶片各截面挥舞方向位移与加载力成线性关系，当载荷加载到达 A点及 C 点之后，位移对于加载力的响应更加敏感且为非线性关系，这表明在 A 点的载荷加载历程中叶片出现了局部屈曲，当载荷继续加载到 B、D 点后，在加载力不变的情况下挥舞方向位移也会增加，Abaqus 停止继续计算，表明此时已达到叶片的整体屈曲。

 

 

图 4-5 叶片展向长度 15.2%及 30.4%位置加载力-位移关系图

 

为更直观的讨论叶片局部屈曲与整体屈曲，比较局部屈曲位置单元结点与结构稳定区域单元节点(随机选取展向长度 54.3%位置主梁区域的节点)的应变及应力演化历程，如图 4-6 至图 4-9。

 

 

图 4-6 屈曲位置应变-弧长曲线

 

 

图 4-7 叶片展向长度 54.3%处应变-弧长曲线

 

 

图 4-8 屈曲位置应力-弧长曲线

 

 

图 4-7 叶片展向长度 54.3%处应力-弧长曲线

 

通过上述对比可以发现，屈曲位置节点在弧长为1.2左右应力应变均有明显的突变，而稳定部位节点虽有斜率的增加，但变化微小。此外，从图 4-6 及 4-8 可以看出，叶根屈曲位置受两个方向的压应力，以及较高的剪切力，足以将单元压溃，致使发生屈曲。

 

对比特征值屈曲以及 Riks 非线性方法，计算预测到的屈曲位置有显著差异，且 Riks 法得到的载荷系数也较特征值法小，这主要是因为改进的 Riks 方法因考虑结构的几何缺陷以及结构的非线性关系 ，得到的结果更加准确可靠，同时也表明特征值只提供可能发生屈曲的位置，对于大变形的计算，该方法适用性不理想。。

 

结论

基于 Abaqus有限元软件采用线性特征值以及非线性弧长法对时代新材 90m级叶片挥舞工况进行了屈曲分析，得到的结论如下：

 

1. 通过特征值屈曲方法得到的屈曲因子为 1.978，满足了 GL2015 及 IEC61400-5 标准的要求，验证了本司叶片结构设计的可靠性。

 

2.通过非线性 Riks 方法考虑了肯能存在的制造缺陷以及结构的非线性，得到的屈曲因子为

1.653，同样满足叶片设计标准。

 

3. 通过线性及非线性两种方法计算结果的对比，表明特征值只提供可能发生屈曲的位置，对于大变形、大挠度的计算，非线性屈曲更加准确。

 

4.应用 Riks 方法可较准确的捕捉到叶片出现局部屈曲的位置，对叶片结构设计具有重要指导意义，可在叶片结构校核工作中进行推广。

 

资料来源：达索官方