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

基于三维壳有限元?梁模型的风力机叶片位移计算方法
【专利摘要】本发明公开了一种基于三维壳有限元?梁模型的风力机叶片位移计算方法，在该计算方法中将叶片三维壳有限元模型等效为梁模型，梁模型的节点与截面参考中心重合，对叶片叶根进行完全约束，在壳有限元模型的自由端分别施加6个独立的单位载荷，基于梁模型柔度矩阵的节点位移计算公式，只需要求解一次各梁模型柔度矩阵，在载荷已知情况下，通过做简单的矩阵乘法，可以直接计算出叶片各气动中心的位移，极大提高在叶片气弹耦合、弯扭耦合特性分析计算中叶片弹性变形的计算效率。
【专利说明】
基于三维壳有限元-梁模型的风力机叶片位移计算方法
技术领域
[0001] 本发明涉及一种基于三维壳有限元-梁模型的风力机叶片位移计算方法，用于计 算叶片等效梁模型的节点位移，同时也用于气弹分析中的位移计算及弯扭耦合特性计算。
【背景技术】
[0002] 叶片变形是风力机叶片进行结构设计和气弹耦合分析的重要参数。叶片刚度设计 需要准确地计算变形;气弹耦合分析中，需要反复迭代计算叶片变形及气动性能。因此，准 确、快速地计算出叶片变形是叶片设计中的重要问题。
[0003] 风力机叶片具有长厚比、长宽比大的特点，常被简化为梁模型，叶片变形计算采用 梁模型时自由度可以很好适应基于叶素动量理论得到的叶片气动载荷，同时梁模型由于具 有单元数量少、计算时间短的特点，是目前应用最广泛的描述叶片的模型之一。在风力机的 设计过程中，为了分析叶片的动力学行为，需要进行风力机多体系统的气弹分析，采用梁模 型使得具有大自由度的多体系统分析成为可能。
[0004] 现代大型风力机叶片具有复杂的截面形状、多种类的铺层材料、丰富的主梁结构 等特点。采用三维壳单元建立叶片模型比采用传统的一维梁单元建立模型能更准确更真实 地涵盖大型叶片的这些特点。此外，三维分析结果能够提供叶片更丰富的变形、应变等信 息，但是叶片三维模型建立的效率低。
[0005] 综上所述，在叶片三维壳有限元模型的建模过程中，由于未引入任何假设，因此具 有很高的建模精度。但是，叶片壳模型建模过程较为繁琐，特别是需要花费大量的时间设置 壳单元的铺层信息，这导致壳有限元模型的建模效率很低。此外，叶片壳有限元模型的自由 度总数很大，无法高效地用于风力机多体系统的气弹分析中。因此，如何从三维壳有限元模 型数据直接得到多体系统分析所需的梁模型数据非常值得研究。

【发明内容】

[0006] 本发明要解决的技术问题是提供一种基于三维壳有限元-梁模型的风力机叶片 位移计算方法，该方法基于梁模型柔度矩阵的节点位移计算公式，在已知叶片各叶素上的 载荷后，可以直接计算出叶片各参考中心的位移，极大提高在叶片气弹耦合、弯扭耦合特性 计算中叶片弹性变形计算效率。
[0007] 为解决上述技术问题，本发明所采取的技术方案是:一种基于三维壳有限元-梁模 型的风力机叶片位移计算方法，其特征在于包括如下计算步骤：
[0008] 1)依据叶片的铺层、翼型和材料数据，建立叶片三维壳有限元模型，壳有限元模型 分为η个横向的截面，将叶片三维壳有限元模型等效为梁模型，梁模型的节点与截面参考中 心重合，对叶片叶根进行完全约束，在壳有限元模型的自由端分别施加6个独立的单位载 荷，其中包括3个相互垂直的力载荷和3个绕力载荷逆时针方向的力矩载荷，形成6种工况；
[0009] 2)计算每种载荷工况下，叶片各个截面等效节点的位移，
[0011]式中，（u，v，w)为叶片壳有限元模型上任意节点(x，y，z)的位移；（U，V，W)为截面参 考中心的线位移，(Φχ，ι1ν，Φζ)为截面参考中心的角位移；
[0012]将各个截面参考中心的位移代入下式，得到每个截面的位移矩阵cU;
[0014] 3)根据壳有限元模型的自由端施加的单位载荷，计算每个横向截面的载荷矩阵fi
[0016]式中，L1为第i个截面距离叶片约束端的距离，i = 1~η;
[0017] 4)根据位移矩阵Cl1和载荷矩阵乜，计算出叶片各悬臂梁模型的6X6柔度矩阵匕
[0018] Ci = Clifr1
[0019] 5)根据每种工况下作用于叶片叶素的载荷、叶片梁模型的柔度矩阵C1和载荷矩阵 h，计算得到梁模型各节点的刚体位移及变形rff
[0022]式中，Mlj为转换矩阵，其计算公式如下
[0020]
[0021]
[0023]
[0024] 式中，Lj为第j个截面距离叶片约束端的距离，j = l~η;
[0025] 6)根据刚体位移及变形 <，计算梁模型各节点的位移Cl1
[0026]
[0027] 7)重复步骤5)和6)，可计算不同载荷下的叶片梁模型节点位移。
[0028]对上述方法作进一步补充，所述的叶片三维壳有限元模型包括叶片坐标信息，其 中坐标原点位于叶根截面的参考中心，沿叶片长度方向上、各截面参考中心连线为ζ轴，且 沿叶片展向为正，xoz平面位于旋转平面内，y为挥舞方向，X为摆振方向。
[0029] 对上述方法作进一步补充，所述的步骤4)中，柔度矩阵匕可表示为
[0030]
[0031] 式中：各元素为梁模型在叶片坐标系下的柔度，Cn和C22为剪切柔度，C33为拉伸柔 度，C44为X向弯曲柔度，C 55为y向弯曲柔度，C66为扭转柔度;非对角线元素表示耦合柔度。
[0032] 采用上述技术方案所产生的有益效果在于：本发明用于风力机叶片位移计算方 法，从叶片三维壳有限元模型出发，采用梁模型柔度定义，已知叶片各叶素上的载荷后，可 以直接计算出叶片各参考中心的位移，极大提高在叶片气弹耦合、动态分析计算中叶片弹 性变形计算效率;在计算过程中，本发明提出方法中，只需要求解一次各悬臂梁模型的柔度 矩阵，即可计算不同载荷作用下的梁模型节点位移，计算得到的叶片气动中心变形与三维 壳有限元计算得到的叶片气动中心的变形具有相同精度，为叶片的气弹耦合分析中位移的 计算或叶片弯扭耦合特性的计算提供便利。
【附图说明】
[0033]下面结合附图和【具体实施方式】对本发明作进一步详细的说明。
[0034]图1是本发明中风力机叶片三维壳有限元模型；
[0035]图2是本发明中风力机叶片梁模型节点载荷示意图；
[0036]图3是本发明中风力机叶片梁模型节点位移示意图；
[0037]图4是风力机叶片叶素载荷示意图；
[0038]图5是各向同性箱型叶片主梁的三维有限元模型；
[0039]图6是图5中叶片主梁的自由端局部放大图；
[0040]图7是主梁梁模型刚度结果对比示意图，其中图7a为扭转刚度结果，图7b为拉伸刚 度结果，图7c为弯曲刚度结果；
[0041 ]图8是叶片主梁的载荷简图；
[0042]图9是叶片主梁变形对比，其中图9a为挠度对比，图9b为角位移对比；
[0043]图10是5MW风力机叶片壳有限元模型；
[0044]图11是5MW风力机叶片壳有限元模型在额定风速工况下的变形；
[0045]图12是5MW风力机叶片梁模型在额定风速工况下变形，其中图12a为角位移，图12b 为线位移。
【具体实施方式】
[0046]本发明具体为一种基于三维壳有限元-梁模型的风力机叶片位移计算方法，用于 计算叶片等效梁模型的节点位移，包括如下计算步骤：
[0047] 1)依据叶片的铺层、翼型和材料数据，建立叶片三维壳有限元模型，壳有限元模型 分为η个横向的截面，如附图1所示。xoz平面位于旋转平面内，y为挥舞方向，X为摆振方向， 叶片各截面参考中心的连线为z轴，且沿叶片展向为正。叶片壳有限元模型的叶根固定，壳 模型上任意节点(x，y，z)的位移为(u, V，w)。将叶片三维壳有限元模型等效为梁模型，如图1 中黑色粗实线所示，即梁模型的节点与截面参考中心重合，梁模型位于z轴上。
[0048]为了得到叶片等效梁模型第i个节点（i截面参考中心）的位移，i节点到叶根的距 离为U，将从叶根至截面i的叶片壳有限元模型段截取出来，如图2和3所示。那么，该叶片段 可以简化为一个悬臂梁模型，见图中黑色粗实线，该悬臂梁模型由一个梁单元组成，且梁单 元的左节点固定(位移为0)。
[0049] 对壳有限元模型的自由端分别施加6个独立的单位载荷，其中包括3个相互垂直的 力载荷和3个绕力载荷逆时针方向的力矩载荷，形成6种工况；
[0050] 假设悬臂梁模型i的右节点（i截面参考中心）的载荷f和位移向量d如图2和3所示， 为
[0051]
[0052]
[0053] 式中，（U，V，W)为截面参考中心的线位移，（ΦΧ，ΦΥ，ΦΖ)为参考中心的角位移；（f x，fy， fz)为截面力，(mx，my，mz)为截面力矩。
[0054] 平衡方程为
[0055] di = Cifi
[0056] 式中，Ci为6X6的梁模型柔度矩阵，
[0057]
[0058] 式中：各元素为梁模型在叶片坐标系下的柔度，如Cn和C22为剪切柔度，C33为拉伸 柔度，C44为X向弯曲柔度，C 55为y向弯曲柔度，C66为扭转柔度;非对角线元素表示耦合柔度， 如C46为弯曲扭转耦合柔度。
[0059] 2)假设叶片在变形前后，截面均为平面，且形状保持不变，即忽略面内、面外翘曲。 则得到由截面参考中心线位移(U，V，W)和角位移(φ χ，!?ν，φζ)表示的叶片三维壳有限元模型 各截面节点的位移。计算每种载荷工况下，叶片各个截面等效节点的位移，其表达式如下
[0060]
[0061]式中，（u，v，w)为叶片壳有限元模型上任意节点(x，y，z)的位移；（U，V，W)为截面参 考中心的线位移，(Φχ，ι1ν，Φζ)为截面参考中心的角位移；
[0062]在上式中，给出了叶片壳有限元模型某截面节点位移和该截面参考中心位移之间 的关系。该方程有6个未知变量(υ，ν，￥，φχ，ΦΥΒ，ΦΒζ)，而方程数目为3η个，为超静定方程。采 用最小二乘法，可求得6个未知变量。
[0063]分别计算在这六个独立工况下叶片壳有限元模型的变形，并将各工况下的壳有限 元模型节点位移代入上式，计算得到这六个工况下悬臂梁模型的i截面参考中心的位移列 向量，组装位移列向量，得到每个截面的位移矩阵di;
[0064]

[0065] 3)根据叶片的三维壳有限元模型，在叶尖施加沿叶片坐标系X、Y和Z轴方向的单位 力或单位力矩，则在这六个独立工况下，由静力平衡方程计算得到第i个梁模型自由端的载 荷矩阵f，
[0066]
[0067]式中，L为叶片长度;L1为第i个悬臂梁模型的长度，即第i个截面到叶片根部的距 离。
[0068] 4)根据位移矩阵Cl1和载荷矩阵fi，计算出叶片各悬臂梁模型的6X6柔度矩阵匕
[0069] Ci = Clifr1
[0070] 5)由叶素动量理论，可以得到叶片在相应载荷工况下的各叶素的气动载荷，如附 图4所示，
[0071] 根据每种工况下作用于叶片叶素的载荷、叶片梁模型的柔度矩阵C1和载荷矩阵A， 计算得到作用在第」0 = 1，2，一，卜1)个截面的气动载荷在1截面引起的刚体位移^/，和作 用在第」0 = 1，1+1，一，111)个截面的气动载荷作用在1截面引起的变形<，
[0072]
[0073]
[0074]
[0075]
[0076]式中，Lj为第j个截面距离叶片约束端的距离，j = l~η;
[0077] 6)由叠加原理，根据刚体位移rff及变形<，计算梁模型各节点的位移Cl1
[0078]
[0079] 7)重复步骤5)和6)，可计算不同载荷下的梁模型节点位移。
[0080] 由以上求解梁模型节点位移的计算步骤可以看出，基于三维壳有限元-梁模型节 点位移计算方法，只需要求解一次各梁模型柔度矩阵，在载荷已知情况下，通过做简单的矩 阵乘法，即可求得叶片变形。该方法可以快速计算不同载荷下的变形，可用于气弹耦合等分 析中。
[0081] 在气弹分析过程中，需要反复迭代计算叶片在不同载荷作用下的变形，采用给出 的方法，只需要求解一次各悬臂梁模型的柔度矩阵，然后重复步骤5)和6)，即可计算不同载 荷作用下的梁模型节点位移。因此，本发明的方法具有较高的计算效率。此外，从本发明方 法的推导过程可以看出，该方法建立了由叶片三维壳模型等效为叶片梁模型的方法，没有 引入假设直接建立梁模型，因此，计算精度较高。此方法定义的叶片变形计算模型与叶片动 量叶素理论中的自由度相适应，没有涉及到载荷从一维到三维的转换，非常方便。
[0082] 下面给出两个计算实例，来说明本发明中计算方法的准确性和计算精度。
[0083]实例一为各向同性材料箱型叶片主梁。
[0084]各向同性材料箱型叶片主梁的宽度和高度分别为1.5m和3m，厚度为0.05m，长度为 30m。箱型叶片主梁材料为钢材，弹性模量为200GPa，泊松比为0.3,密度为7800kg/m3。
[0085]采用有限元ANSYS软件建立的根部固定端约束的箱型叶片主梁三维壳有限元模 型，如图5所示。其中，为了减小边界效应的影响，在自由端处耦合了一段梁单元，用于施加 载荷，如图6所示。
[0086] 按照本发明的方法，计算得到各梁模型的柔度矩阵C，代入式
[0087] Si ~ Ci1
[0088]中得到梁模型的刚度矩阵S。
[0089]为了验证基于三维壳有限元-梁模型风力机叶片位移计算方法的正确性，将刚度 矩阵S的各系数与经典梁的刚度矩阵公式进行比较。
[0090] 经典梁的刚度矩阵公式为
[0091]
[0092] 式中，EA为截面拉伸刚度;EIx和EIy为弯曲刚度，GJ为扭转刚度，各刚度的计算公式 如下，
[0093]
[0094] 式中，Ω为箱型梁中面所围面积，δ为壁厚。
[0095] 将本发明的方法计算得到的刚度矩阵与经典梁的刚度矩阵计算得到的结果相比 较，得到
[0096]
[0097] 将本发明的方法计算得到的刚度矩阵与经典梁的刚度矩阵计算得到Sn、S22、S33和 S66进行比较，如图7所示。
[0098] 从图7可以看出，采用本发明计算方法与经典梁的刚度矩阵计算得到的拉伸、扭转 刚度几乎相同，误差均小于5%。其中由于边界效应的影响，弯曲刚度在约束部位附近误差 大于其他截面。
[0099] 箱型叶片主梁受到载荷如图8所示，采用本发明的计算步骤，计算得到叶片各截面 变形。将计算结果与梁模型、壳有限元模型、材料力学解析解进行比较如图9所示。
[0100]从图9中可以看出，各种计算方法得到的挠度和弯曲角的变化趋势是相同的。采用 该方法计算得到的叶片气动中心变形与三维壳有限元模型计算得到的叶片气动中心的变 形具有相同精度，验证了计算方法的正确性。
[0101]实例二为5MW风力机叶片变形计算。
[0102]建立的5丽风力机叶片模型，5丽叶片由蒙皮、梁帽和腹板组成，叶片长度为61.5m， 由8种不同的翼型组成，在20%附近的位置达到最大弦长，为4.7m。
[0103] 表1给出了5MW叶片所用复合材料的力学特性。采用ANSYS软件的shell99单元建立 的三维叶片壳有限元模型如图10所示，计算了5MW叶片在额定风速为11.4m/s工况下的叶片 壳有限元模型的位移如图11所示。
[0104]
[0105] 按照位移计算步骤，计算了5MW叶片在额定风速为11.4m/s工况下的叶片梁模型的 位移如图12所示。从图12可以看出，在11.4m/s风速下，叶片拉伸变形及扭转变形均很小，通 常忽略不计;叶尖挥舞方向变形为4.6m，达到了合理安全距离，避免了与塔架发生碰撞。
[0106] 由上述两个实例可知，本发明的计算方法定义的叶片变形计算模型与叶片动量叶 素理论中的自由度相适应，没有涉及到载荷从一维到三维转换，非常方便。而且该方法基于 叶片三维模型，未引入各种假设直接建立梁模型，因此，计算结果比传统方法精度更高。
[0107] 本发明的计算方法兼具了三维壳有限元模型的高精度和叶片梁模型的高效率。在 已知叶片各叶素上的载荷的情况下，利用本计算方法可以直接计算出叶片梁模型节点的位 移。为叶片的气弹耦合分析中气动载荷的计算或叶片弯扭耦合特性的计算提供便利。
【主权项】
1. 一种基于三维壳有限元-梁模型的风力机叶片位移计算方法，其特征在于包括如下 计算步骤： 1) 依据叶片的铺层、翼型和材料数据，建立叶片三维壳有限元模型，壳有限元模型包含 n个横向的截面，将叶片三维壳有限元模型等效为梁模型，梁模型的节点与截面参考中心重 合，对叶片叶根进行完全约束，在壳有限元模型的自由端分别施加6个独立的单位载荷，其 中包括3个相互垂直的力载荷和3个绕力载荷逆时针方向的力矩载荷，形成6种工况； 2) 计算每种载荷工况下，叶片各个截面参考中心的位移，式中，（u，v，w)为叶片壳有限元模型上任意节点(x，y，z)的位移；（U，V，W)为截面参考中 心的线位移，(t，如，也)为截面参考中心的角位移； 将各个截面参考中心的位移代入下式，得到每个截面的位移矩阵d1;3) 根据壳有限元模型的自由端施加的单位载荷，计算每个横向截面的载荷矩阵fi式中，U为第i个截面距离叶片约束端的距离，i = 1~n; 4) 根据位移矩阵cU和载荷矩阵h，计算出叶片各悬臂梁模型的6X6柔度矩阵匕 Ci = difi-1 5) 根据每种工况下作用于叶片叶素的载荷、叶片梁模型的柔度矩阵Q和载荷矩阵h，计 算得到梁模型各节点的刚体位移及变形<式中，为转换矩阵，其计算公式如下式中，L」为第j个截面距离叶片约束端的距离，j = l~n; 6) 根据刚体位移及变形rff，计算梁模型各节点的位移cU dt = df + rff 4 7) 重复步骤5)和6)，可计算不同载荷下的叶片梁模型节点位移。2. 根据权利要求1所述的基于三维壳有限元-梁模型的风力机叶片位移计算方法，其特 征在于所述的叶片三维壳有限元模型包括叶片坐标信息，其中坐标原点位于叶根截面的参 考中心，沿叶片长度方向上、各截面参考中心连线为z轴，且沿叶片展向为正，xoz平面位于 旋转平面内，y为挥舞方向，x为摆振方向。3. 根据权利要求2所述的基于三维壳有限元-梁模型的风力机叶片位移计算方法，其特 征在于所述的步骤4)中，柔度矩阵匕可表示为式中：各元素为梁模型在叶片坐标系下的柔度，Cll和C22为剪切柔度，C33为拉伸柔度，C44 为x向弯曲柔度，C55为y向弯曲柔度，C66为扭转柔度;非对角线元素表示耦合柔度。
【文档编号】G06F17/50GK106055764SQ201610356352
【公开日】2016年10月26日
【申请日】2016年5月26日
【发明人】周邢银, 安利强
【申请人】华北电力大学（保定）