基于三维壳有限元‑梁模型的风力机叶片位移计算方法

摘要

本发明公开了一种基于三维壳有限元‑梁模型的风力机叶片位移计算方法，在该计算方法中将叶片三维壳有限元模型等效为梁模型，梁模型的节点与截面参考中心重合，对叶片叶根进行完全约束，在壳有限元模型的自由端分别施加6个独立的单位载荷，基于梁模型柔度矩阵的节点位移计算公式，只需要求解一次各梁模型柔度矩阵，在载荷已知情况下，通过做简单的矩阵乘法，可以直接计算出叶片各气动中心的位移，极大提高在叶片气弹耦合、弯扭耦合特性分析计算中叶片弹性变形的计算效率。

说明书

技术领域

本发明涉及一种基于三维壳有限元-梁模型的风力机叶片位移计算方法，用于计算叶片等效梁模型的节点位移，同时也用于气弹分析中的位移计算及弯扭耦合特性计算。

背景技术

叶片变形是风力机叶片进行结构设计和气弹耦合分析的重要参数。叶片刚度设计需要准确地计算变形；气弹耦合分析中，需要反复迭代计算叶片变形及气动性能。因此，准确、快速地计算出叶片变形是叶片设计中的重要问题。

风力机叶片具有长厚比、长宽比大的特点，常被简化为梁模型，叶片变形计算采用梁模型时自由度可以很好适应基于叶素动量理论得到的叶片气动载荷，同时梁模型由于具有单元数量少、计算时间短的特点，是目前应用最广泛的描述叶片的模型之一。在风力机的设计过程中，为了分析叶片的动力学行为，需要进行风力机多体系统的气弹分析，采用梁模型使得具有大自由度的多体系统分析成为可能。

现代大型风力机叶片具有复杂的截面形状、多种类的铺层材料、丰富的主梁结构等特点。采用三维壳单元建立叶片模型比采用传统的一维梁单元建立模型能更准确更真实地涵盖大型叶片的这些特点。此外，三维分析结果能够提供叶片更丰富的变形、应变等信息，但是叶片三维模型建立的效率低。

综上所述，在叶片三维壳有限元模型的建模过程中，由于未引入任何假设，因此具有很高的建模精度。但是，叶片壳模型建模过程较为繁琐，特别是需要花费大量的时间设置壳单元的铺层信息，这导致壳有限元模型的建模效率很低。此外，叶片壳有限元模型的自由度总数很大，无法高效地用于风力机多体系统的气弹分析中。因此，如何从三维壳有限元模型数据直接得到多体系统分析所需的梁模型数据非常值得研究。

发明内容

本发明要解决的技术问题是提供一种基于三维壳有限元-梁模型的风力机 叶片位移计算方法，该方法基于梁模型柔度矩阵的节点位移计算公式，在已知叶片各叶素上的载荷后，可以直接计算出叶片各参考中心的位移，极大提高在叶片气弹耦合、弯扭耦合特性计算中叶片弹性变形计算效率。

为解决上述技术问题，本发明所采取的技术方案是：一种基于三维壳有限元-梁模型的风力机叶片位移计算方法，其特征在于包括如下计算步骤：

1)依据叶片的铺层、翼型和材料数据，建立叶片三维壳有限元模型，壳有限元模型分为n个横向的截面，将叶片三维壳有限元模型等效为梁模型，梁模型的节点与截面参考中心重合，对叶片叶根进行完全约束，在壳有限元模型的自由端分别施加6个独立的单位载荷，其中包括3个相互垂直的力载荷和3个绕力载荷逆时针方向的力矩载荷，形成6种工况；

2)计算每种载荷工况下，叶片各个截面等效节点的位移，

式中，(u,v,w)为叶片壳有限元模型上任意节点(x,y,z)的位移；(U,V,W)为截面参考中心的线位移，(ψ_x,ψ_y,ψ_z)为截面参考中心的角位移；

将各个截面参考中心的位移代入下式，得到每个截面的位移矩阵d_i；

3)根据壳有限元模型的自由端施加的单位载荷，计算每个横向截面的载荷矩阵f_i

式中，L_i为第i个截面距离叶片约束端的距离，i＝1～n；

4)根据位移矩阵d_i和载荷矩阵f_i，计算出叶片各悬臂梁模型的6×6柔度矩阵C_i

C_i＝d_if_i^-1

5)根据每种工况下作用于叶片叶素的载荷、叶片梁模型的柔度矩阵C_i和载荷矩阵f_i，计算得到梁模型各节点的刚体位移及变形

式中，M_ij为转换矩阵，其计算公式如下

式中，L_j为第j个截面距离叶片约束端的距离，j＝1～n；

6)根据刚体位移及变形，计算梁模型各节点的位移d_i

7)重复步骤5)和6)，可计算不同载荷下的叶片梁模型节点位移。

对上述方法作进一步补充，所述的叶片三维壳有限元模型包括叶片坐标信息，其中坐标原点位于叶根截面的参考中心，沿叶片长度方向上、各截面参考中心连线为z轴，且沿叶片展向为正，xoz平面位于旋转平面内，y为挥舞方向， x为摆振方向。

对上述方法作进一步补充，所述的步骤4)中，柔度矩阵C_i可表示为

式中：各元素为梁模型在叶片坐标系下的柔度，C₁₁和C₂₂为剪切柔度，C₃₃为拉伸柔度，C₄₄为x向弯曲柔度，C₅₅为y向弯曲柔度，C₆₆为扭转柔度；非对角线元素表示耦合柔度。

采用上述技术方案所产生的有益效果在于：本发明用于风力机叶片位移计算方法，从叶片三维壳有限元模型出发，采用梁模型柔度定义，已知叶片各叶素上的载荷后，可以直接计算出叶片各参考中心的位移，极大提高在叶片气弹耦合、动态分析计算中叶片弹性变形计算效率；在计算过程中，本发明提出方法中，只需要求解一次各悬臂梁模型的柔度矩阵，即可计算不同载荷作用下的梁模型节点位移，计算得到的叶片气动中心变形与三维壳有限元计算得到的叶片气动中心的变形具有相同精度，为叶片的气弹耦合分析中位移的计算或叶片弯扭耦合特性的计算提供便利。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明中风力机叶片三维壳有限元模型；

图2是本发明中风力机叶片梁模型节点载荷示意图；

图3是本发明中风力机叶片梁模型节点位移示意图；

图4是风力机叶片叶素载荷示意图；

图5是各向同性箱型叶片主梁的三维有限元模型；

图6是图5中叶片主梁的自由端局部放大图；

图7是主梁梁模型刚度结果对比示意图，其中图7a为扭转刚度结果，图7b为拉伸刚度结果，图7c为弯曲刚度结果；

图8是叶片主梁的载荷简图；

图9是叶片主梁变形对比，其中图9a为挠度对比，图9b为角位移对比；

图10是5MW风力机叶片壳有限元模型；

图11是5MW风力机叶片壳有限元模型在额定风速工况下的变形；

图12是5MW风力机叶片梁模型在额定风速工况下变形，其中图12a为角位移，图12b为线位移。

具体实施方式

本发明具体为一种基于三维壳有限元-梁模型的风力机叶片位移计算方法，用于计算叶片等效梁模型的节点位移,包括如下计算步骤：

1)依据叶片的铺层、翼型和材料数据，建立叶片三维壳有限元模型，壳有限元模型分为n个横向的截面，如附图1所示。xoz平面位于旋转平面内，y为挥舞方向，x为摆振方向，叶片各截面参考中心的连线为z轴，且沿叶片展向为正。叶片壳有限元模型的叶根固定，壳模型上任意节点(x,y,z)的位移为(u,v,w)。将叶片三维壳有限元模型等效为梁模型，如图1中黑色粗实线所示，即梁模型的节点与截面参考中心重合，梁模型位于z轴上。

为了得到叶片等效梁模型第i个节点(i截面参考中心)的位移，i节点到叶根的距离为L_i，将从叶根至截面i的叶片壳有限元模型段截取出来，如图2和3所示。那么，该叶片段可以简化为一个悬臂梁模型，见图中黑色粗实线，该悬臂梁模型由一个梁单元组成，且梁单元的左节点固定(位移为0)。

对壳有限元模型的自由端分别施加6个独立的单位载荷，其中包括3个相互垂直的力载荷和3个绕力载荷逆时针方向的力矩载荷，形成6种工况；

假设悬臂梁模型i的右节点(i截面参考中心)的载荷f和位移向量d如图2和3所示，为

式中，(U,V,W)为截面参考中心的线位移，(ψ_x,ψ_y,ψ_z)为参考中心的角位移；(f_x,f_y,f_z)为截面力，(m_x,m_y,m_z)为截面力矩。

平衡方程为

d_i＝C_if_i

式中，C_i为6×6的梁模型柔度矩阵，

式中：各元素为梁模型在叶片坐标系下的柔度，如C₁₁和C₂₂为剪切柔度，C₃₃为拉伸柔度，C₄₄为x向弯曲柔度，C₅₅为y向弯曲柔度，C₆₆为扭转柔度；非对角线元素表示耦合柔度，如C₄₆为弯曲扭转耦合柔度。

2)假设叶片在变形前后，截面均为平面，且形状保持不变，即忽略面内、面外翘曲。则得到由截面参考中心线位移(U,V,W)和角位移(ψ_x,ψ_y,ψ_z)表示的叶片三维壳有限元模型各截面节点的位移。计算每种载荷工况下，叶片各个截面等效节点的位移，其表达式如下

式中，(u,v,w)为叶片壳有限元模型上任意节点(x,y,z)的位移；(U,V,W)为截面参考中心的线位移，(ψ_x,ψ_y,ψ_z)为截面参考中心的角位移；

在上式中，给出了叶片壳有限元模型某截面节点位移和该截面参考中心位移之间的关系。该方程有6个未知变量(U,V,W,ψ_x,ψ_yB,ψB_z)，而方程数目为3n个，为超静定方程。采用最小二乘法，可求得6个未知变量。

分别计算在这六个独立工况下叶片壳有限元模型的变形，并将各工况下的壳有限元模型节点位移代入上式，计算得到这六个工况下悬臂梁模型的i截面参考中心的位移列向量，组装位移列向量，得到每个截面的位移矩阵d_i；

3)根据叶片的三维壳有限元模型，在叶尖施加沿叶片坐标系X、Y和Z轴方向的单位力或单位力矩，则在这六个独立工况下，由静力平衡方程计算得到第i个梁模型自由端的载荷矩阵f_i

式中，L为叶片长度；L_i为第i个悬臂梁模型的长度，即第i个截面到叶片根部的距离。

4)根据位移矩阵d_i和载荷矩阵fi，计算出叶片各悬臂梁模型的6×6柔度矩阵C_i

C_i＝d_if_i^-1

5)由叶素动量理论，可以得到叶片在相应载荷工况下的各叶素的气动载荷，如附图4所示，

根据每种工况下作用于叶片叶素的载荷、叶片梁模型的柔度矩阵C_i和载荷矩阵f_i，计算得到作用在第j(j＝1,2,…,i-1)个截面的气动载荷在i截面引起的刚体位移，和作用在第j(j＝i,i+1,…,m)个截面的气动载荷作用在i截面引起的变形，

式中，M_ij为转换矩阵，其计算公式如下

式中，L_j为第j个截面距离叶片约束端的距离，j＝1～n；

6)由叠加原理，根据刚体位移及变形，计算梁模型各节点的位移d_i

7)重复步骤5)和6)，可计算不同载荷下的梁模型节点位移。

由以上求解梁模型节点位移的计算步骤可以看出，基于三维壳有限元-梁模型节点位移计算方法，只需要求解一次各梁模型柔度矩阵，在载荷已知情况下，通过做简单的矩阵乘法，即可求得叶片变形。该方法可以快速计算不同载荷下的变形，可用于气弹耦合等分析中。

在气弹分析过程中，需要反复迭代计算叶片在不同载荷作用下的变形，采用给出的方法，只需要求解一次各悬臂梁模型的柔度矩阵，然后重复步骤5)和6)，即可计算不同载荷作用下的梁模型节点位移。因此，本发明的方法具有较高的计算效率。此外，从本发明方法的推导过程可以看出，该方法建立了由叶片三维壳模型等效为叶片梁模型的方法，没有引入假设直接建立梁模型，因此，计算精度较高。此方法定义的叶片变形计算模型与叶片动量叶素理论中的自由度相适应，没有涉及到载荷从一维到三维的转换，非常方便。

下面给出两个计算实例，来说明本发明中计算方法的准确性和计算精度。

实例一为各向同性材料箱型叶片主梁。

各向同性材料箱型叶片主梁的宽度和高度分别为1.5m和3m，厚度为0.05m，长度为30m。箱型叶片主梁材料为钢材，弹性模量为200GPa，泊松比为0.3，密度为7800kg/m³。

采用有限元ANSYS软件建立的根部固定端约束的箱型叶片主梁三维壳有限元模型，如图5所示。其中，为了减小边界效应的影响，在自由端处耦合了一段梁单元，用于施加载荷，如图6所示。

按照本发明的方法，计算得到各梁模型的柔度矩阵C，代入式

中得到梁模型的刚度矩阵S。

为了验证基于三维壳有限元-梁模型风力机叶片位移计算方法的正确性，将刚度矩阵S的各系数与经典梁的刚度矩阵公式进行比较。

经典梁的刚度矩阵公式为

式中，EA为截面拉伸刚度；EI_x和EI_y为弯曲刚度，GJ为扭转刚度，各刚度的计算公式如下，

式中，Ω为箱型梁中面所围面积，δ为壁厚。

将本发明的方法计算得到的刚度矩阵与经典梁的刚度矩阵计算得到的结果相比较，得到

将本发明的方法计算得到的刚度矩阵与经典梁的刚度矩阵计算得到S₁₁、S₂₂、S₃₃和S₆₆进行比较，如图7所示。

从图7可以看出，采用本发明计算方法与经典梁的刚度矩阵计算得到的拉伸、扭转刚度几乎相同，误差均小于5％。其中由于边界效应的影响，弯曲刚度在约束部位附近误差大于其他截面。

箱型叶片主梁受到载荷如图8所示，采用本发明的计算步骤，计算得到叶片各截面变形。将计算结果与梁模型、壳有限元模型、材料力学解析解进行比较如图9所示。

从图9中可以看出，各种计算方法得到的挠度和弯曲角的变化趋势是相同的。采用该方法计算得到的叶片气动中心变形与三维壳有限元模型计算得到的叶片气动中心的变形具有相同精度，验证了计算方法的正确性。

实例二为5MW风力机叶片变形计算。

建立的5MW风力机叶片模型,5MW叶片由蒙皮、梁帽和腹板组成，叶片长度为61.5m，由8种不同的翼型组成，在20％附近的位置达到最大弦长，为4.7m。

表1给出了5MW叶片所用复合材料的力学特性。采用ANSYS软件的shell99单元建立的三维叶片壳有限元模型如图10所示，计算了5MW叶片在额定风速为11.4m/s工况下的叶片壳有限元模型的位移如图11所示。

按照位移计算步骤，计算了5MW叶片在额定风速为11.4m/s工况下的叶片梁模型的位移如图12所示。从图12可以看出，在11.4m/s风速下，叶片拉伸变形及扭转变形均很小，通常忽略不计；叶尖挥舞方向变形为4.6m，达到了合理安全距离，避免了与塔架发生碰撞。

由上述两个实例可知，本发明的计算方法定义的叶片变形计算模型与叶片动量叶素理论中的自由度相适应，没有涉及到载荷从一维到三维转换，非常方便。而且该方法基于叶片三维模型，未引入各种假设直接建立梁模型，因此，计算结果比传统方法精度更高。

本发明的计算方法兼具了三维壳有限元模型的高精度和叶片梁模型的高效率。在已知叶片各叶素上的载荷的情况下，利用本计算方法可以直接计算出叶片梁模型节点的位移。为叶片的气弹耦合分析中气动载荷的计算或叶片弯扭耦合特性的计算提供便利。