绝对节点坐标法大型风力机复合材料叶片的动态特性*

胡锡彪，夏鸿建，李德源

绝对节点坐标法大型风力机复合材料叶片的动态特性*

胡锡彪，夏鸿建†，李德源

（广东工业大学 机电工程学院，广州 510006）

为研究大型风力机旋转复合材料叶片的非线性结构动态特性，采用绝对节点坐标法，并基于一般连续介质力学理论和复合材料本构关系，推导广义质量与刚度矩阵，建立大型风力机旋转复合材料叶片的非线性动力学模型。为描述叶片的转动，引入大范围转动连体坐标系，结合摄动原理及线性化方法，构建旋转叶片的状态方程。以非线性梁标准算例验证了模型的准确性，并以DTU-10MW风力机复合材料叶片为对象，分析了静止与旋转条件下的动态特性。结果表明，由于叶片的大型化和非线性变形特点，叶片静止时的挥舞、摆振和扭转振型存在显著耦合；由于动力刚化效应，叶片的各阶频率会随着转速的提高而改变，并对振型的挥舞和摆振产生较大影响。

风力机叶片；复合材料；绝对节点坐标法；动态特性分析

0 引 言

叶片是风力机发电系统中的两个主要柔性部件之一，其结构性能的优劣将直接影响风力发电机的功率和风能利用率，由于复合材料具有质量轻、耐腐蚀、比强度高及减震性能好等特点[1]，采用新型复合材料及合理的铺层优化可有效提高刚度并降低叶片重量，进而降低成本[2]，现有大型风力机叶片多采用复合材料铺层制造。国内外学者对风力机复合材料叶片的动态特性开展了许多研究。张立等[3]结合复合材料铺层设计，采用有限元法对复合材料叶片模型进行模态和屈曲分析，得到叶片自重对其结构特性影响。吕计男等[4]建立了复合材料风力机叶片的实体模型，研究了影响叶片频域的因素。刘雄等[5]考虑叶片气动阻尼及旋转产生的离心刚化作用，研究叶片弯曲振动的动力特性，但其将叶片简化为二节点梁单元。李德源等[6]将多体动力学理论与有限元结合起来，研究风力机叶片的振动特性，揭示了动力刚化效应对叶片固有频率的影响规律。周鹏展等[7]采用ANSYS软件对复合材料风力机叶片进行结构分析，结果表明，该叶片的振型以一阶挥舞和一阶摆振为主。BAYOUMY等[8]利用绝对节点坐标公式建立一个三维弹性风力机叶片连续网单元模型，并用六节点板单元对叶片结构进行建模。吕品等[9]采用几何精确梁理论，建立一种用于叶片的非线性单元，计算风力机叶片和整机的频域特性。刘宇航等[10]采用节点位移法分析铺层结构对弯扭耦合叶片整体性能的影响，发现镜像对称铺层结构可实现叶片弯扭耦合。

综上，目前对风力机叶片动态特性的分析方法主要为有限元法、几何精确梁法以及绝对节点坐标法。有限元法在分析柔长叶片的大旋转、大变形问题时，为建立准确模型而增加其自由度，难以将其用于多场耦合分析。几何精确梁法考虑叶片的弯曲和挠曲变形，能够较好地描述叶片的非线性变形特性，但其对插值策略和数值积分方法有较高的要求。绝对节点坐标法[11-12]基于一般连续介质力学理论，质量矩阵为常数，不存在离心力和科氏力，可用于模拟叶片的弯曲和扭转以及非线性材料的应力应变行为，从而更准确地预测叶片的动力响应。

为此，利用丹麦技术大学（Technical University of Denmark, DTU）所提供的10 MW风力机叶片模型[13]，采用绝对节点坐标法，基于一般连续介质力学理论和复合材料本构关系，建立大型风力机旋转复合材料叶片的非线性动力学模型。为了描述叶片的转动，引入大范围转动连体坐标系，结合摄动原理及线性化方法，构建旋转叶片的状态方程，研究叶片在静止与不同转速下的模态振型变化，提供算例验证该模型的正确性，为改进柔性叶片的结构设计提供参考。

1 绝对节点坐标法下的动力学方程

绝对节点坐标法采用连续介质力学理论，通过单元形状函数和节点坐标向量定义梁上任意点的全局坐标向量。

图1 三维梁单元模型

如图1所示，柔性梁单元的所有坐标都建立在全局坐标系下，三维二节点梁单元上的任意一点对应的绝对位置矢量可以表示为

1.1 单元质量矩阵

通过节点的速度表达式，可计算得到梁单元的动能：

1.2 单元刚度矩阵

在一般连续介质力学方法中，任意节点的变形梯度可由非线性应变−位移关系计算，表达式为

通过变形梯度可计算出梁的格林−拉格朗日应变张量：

可简写为

根据应变表达式和应力应变关系式可得到单元应变能为

式中：

1.3 系统动力学方程

绝对节点坐标下三维旋转梁的动力学方程为

2 风力机叶片非线性动力学建模

2.1 叶片建模及材料铺层设计

以DTU的10 MW风力机叶片为研究对象[13]，叶片总长为86.366 m，最大弦长为6.026 m，其叶片翼型主要使用的是FFA-W3-xxx翼型系列。采用玻璃纤维复合材料和巴沙木（basal wood）作为夹芯材料，其叶片的复合层由多向层的堆积顺序来定义，而多向层板铺层结构主要由两个E型玻璃纤维组成的单向层板和环氧树脂基体组成。每个层都有多个方向角度的纤维，利用微观力学方程和经典的层压理论[16]，根据组成材料的力学性能，可推导出多向层板的力学性能。

叶片的各部分分别由不同材料铺层所组成，根据文献[13]中叶片的铺层和材料参数，多向层板的纤维方向与力学性能如表1所示。

表1 多向层板的纤维方向和力学性能

2.2 叶片翼型截面信息处理

由于叶片翼型截面形状不规则且复合材料铺层的分区不同，在计算质量矩阵和刚度矩阵的积分时，往往需要计算复合函数的定积分，而这些复杂函数的原函数很难求得，而且计算会消耗大量的计算资源，因此引入数值积分，利用网格划分软件对叶片模型进行三角形单元的划分，将叶片分成10个梁单元，每个单元沿着风力机叶片展向方向选取4个点，此时可将风力机叶片截取出40个截面，如图2a所示。图2b为第10截面的形状，而表2、表3为获取的第10截面三角形单元节点部分数据，将各截面的节点参数导入数值软件中，通过高斯积分和Hammer积分公式相结合求得质量矩阵和刚度矩阵。

图2 叶片截面位置及形状

表2 第10截面部分单元信息

表3 第10截面部分节点信息

在叶片截面中，使用Hammer积分进行求解时，需要在同一个三角形单元进行积分后再累加计算所有的三角形单元，得到总体积分。通过对单元信息找到三个节点的序号，如表2第1号单元，由节点2、14、8组成，根据这些节点信息找到对应坐标值，表3中、方向表示该节点的绝对坐标，夹角表示该节点复合材料铺层方向与变桨中心坐标系的夹角。其Hammer积分形式为

在叶片单元展向方向使用高斯积分求解，其积分形式一般如下

2.3 旋转动力学方程的建立

为了描述叶片的转动，引入大范围浮动坐标系，建立关于浮动坐标系下的动力学方程，绝对坐标需要用新的坐标向量来表示。

将式（19）代入式（16）中，将其转化为浮动坐标系下旋转梁的动力学方程：

2.4 动力学方程线性化

可将式（24）改为矩阵表达式，如下式所示：

求解得到状态方程：

3 仿真算例

根据上述绝对节点坐标法构建的动力学模型，编写对应的MATLAB程序，通过两个算例验证本模型的准确性，以此来分析风力机旋转复合材料叶片的动力学特性。

3.1 算例一：定转速旋转梁

由于风力机叶片旋转与旋转梁的运动较为相似，梁的一端固定，绕着定轴按照设定角度旋转，所以选用该算例来验证模型准确性。在该算例中，使用文献[18]相同参数，梁绕轴旋转的转动角度随时间变化的函数如下：

3.2 算例二：复合材料梁的固有频率

通过文献[19]中的复合材料梁参数，结合复合材料本构关系，计算复合材料梁的固有频率。表4列出了计算结果与文献结果的对比，可看出在梁的低阶频率误差较小，以此验证复合材料本构关系与绝对节点坐标法建立的动力学模型是准确的。

表4 复合材料的固有频率

4 风力机叶片动力特性分析

线性化方程（31），获得其状态方程，并求解出叶片静止时以及不同转速下的固有频率和振型。

4.1 叶片静止时的模态

得到各阶模态对应的固有频率如表5，并与文献[13]的结果进行对比可以看出，结果基本一致，说明基于绝对节点坐标法的复合材料叶片单元对叶片的结构处理是合理的。图4为叶片前八阶模态振型。

表5 静止叶片的模态分析结果

图4 静止时叶片的前八阶振型

由图4可见，叶片振型随着阶数的增大而存在显著耦合。低阶振型主要为挥舞和摆振，高阶时叶片在第六阶开始出现扭转，第七阶为挥舞、摆振、扭转三者之间的耦合运动。

4.2 不同转速下叶片的固有频率

为研究叶片柔性变形与大范围刚性转动对叶片固有频率的影响，进一步通过绝对节点坐标法分别对3.0、6.0、9.6（额定转速）、12.0、15.0、19.2 r/min时进行模态特性分析，前八阶固有频率随转速的变化趋势如图5所示。

由图5可以看出，叶片旋转时，由于动力刚化效应，叶片的各阶频率会随着转速的提高而改变，并对振型的挥舞和摆振产生较大影响。其中当转速达19.2 r/min时，相对于静止时的频率，一阶、二阶挥舞频率分别增加了0.131 Hz和0.145 Hz，增加幅度为21.96%和8.32%；摆振方向分别增加了0.024 Hz和0.063 Hz，增加幅度为2.57%和2.18%。

图5 叶片在不同转速下的固有频率

5 结 论

采用绝对节点坐标法，基于一般连续介质力学理论和复合材料本构关系，建立柔性梁的非线性动力学模型，并通过不同算例验证该模型动态特性的正确性。为了描述叶片的转动，结合浮点坐标法，建立风力机旋转复合材料叶片的非线性动力学模型，并使用线性化方法得到其状态方程，将该模型应用于10 MW风力机叶片的动态特性分析，计算了该叶片在不同转速下的固有频率和振型。通过绝对节点坐标法对风力机旋转复合材料叶片的建模及其动态特性进行分析，结论如下：

（1）通过两个算例验证了所建立的绝对节点坐标法模型能够分析几何非线性变形问题，并应用到叶片的模态分析中，所得结果与参考文献结果一致。

（2）在绝对节点坐标法基础上结合高斯积分与Hammer积分，通过对叶片截面信息处理可进行复杂积分的运算，实现复杂截面叶片模型动力学方程的建立。

（3）在静止时，叶片振型随着阶数的增大而出现耦合现象。低阶振型主要为挥舞振动和摆振振动，叶片在第六阶开始出现扭转，第七阶为挥舞、摆振、扭转三者之间的耦合运动。

（4）在不同转速下，由于动力刚化效应，在挥舞方向比摆振方向更容易受叶片转速的影响。

[1] ZHENG Y Q, MA H D, WEI J F, et al. Robust optimization for composite blade of wind turbine based on kriging model[J]. Advanced composites letters, 2020, 29: 2633366X2091463. DOI: 10.1177/2633366X20914631.

[2] SIEROS G, CHAVIAROPOULOS P, SØRENSEN J D, et al. Upscaling wind turbines: theoretical and practical aspects and their impact on the cost of energy[J]. Wind energy, 2012, 15(1): 3-17. DOI:10.1002/we.527.

[3] 张立, 缪维跑, 闫阳天, 等. 考虑自重影响的大型风力机复合材料叶片结构力学特性分析[J]. 中国电机工程学报, 2020, 40(19): 6272-6283.

[4] 吕计男, 刘子强. NREL S809风力机叶片动力学特征分析[C]//第十一届全国空气弹性学术交流会会议论文集. 昆明: 中国力学学会, 2009: 579-584.

[5] 刘雄, 李钢强, 陈严, 等. 水平轴风力机叶片动态响应分析[J]. 机械工程学报, 2010, 46(12): 128-134, 141. DOI:10.3901/JME.2010.12.128.

[6] 李德源, 池志强, 汪显能, 等. 10MW风力机叶片设计与动力特性[J]. 沈阳工业大学学报, 2016, 38(3): 241-246. DOI:10.7688/j.issn.1000-1646.2016.03.01.

[7] 周鹏展, 肖加余, 曾竟成, 等. 基于ANSYS的大型复合材料风力机叶片结构分析[J]. 国防科技大学学报, 2010, 32(2): 46-50.

[8] BAYOUMY A H, NADA A A, MEGAHED S M. A continuum based three-dimensional modeling of wind turbine blades[J]. Journal of computational and nonlinear dynamics, 2013, 8(3):031004.

[9] 吕品, 廖明夫, 尹尧杰. 考虑几何非线性的风力机叶片气弹模型[J]. 机械科学与技术, 2015, 34(12): 1805-1812. DOI:10.13433/j.cnki.1003-8728.2015.1201.

[10] 刘宇航, 王渊博, 李春, 等. 应用于气动弹性剪裁的大型风力机叶片弯扭耦合性能研究[J]. 动力工程学报, 2018, 38(12): 1016-1021. DOI:10.3969/j.issn.1674-7607. 2018.12.011.

[11] SHABANA A A. Computational continuum mechanics[M]. New York: Cambridge University Press, 2008.

[12] SHABANA AA, YAKOUB RY. Three dimensional absolute nodal coordinate formulation for beam elements: theory[J]. Journal of mechanical design, 2001, 123(4): 606-613. DOI:10.1115/1.1410100.

[13] BAK C, ZAHLE F, BITSCHE R, et al. The DTU 10-MW reference wind turbine[R]. Fredericia: DTU, 2013.

[14] 范纪华, 章定国, 谌宏. 基于绝对节点坐标法的弹性线方法研究[J]. 力学学报, 2019, 51(5): 1455-1465. DOI:10.6052/0459-1879-19-076.

[15] KOLLÁR L P, SPRINGER G S. Mechanics of composite structures[M]. New York: Cambridge University Press, 2003.

[16] 步鹏飞, 任辉启, 阮文俊. 铺层角度对碳纤维/环氧树脂基复合材料板等效刚度的影响[J]. 南京理工大学学报, 2021, 45(5): 537-544. DOI:10.14177/j.cnki.32-1397n. 2021.45.05.003.

[17] MAQUEDA L G, BAUCHAU O A, SHABANA A A. Effect of the centrifugal forces on the finite element eigenvalue solution of a rotating blade: a comparative study[J]. Multibody system dynamics,2008, 19(3): 281-302.DOI:10.1007/s11044-007-9070-6.

[18] FAN W, ZHU W D. An accurate singularity-free formulation of a three-dimensional curved Euler– Bernoulli beam for flexible multibody dynamic analysis[J]. Journal of vibration and acoustics, 2016, 138(5): 051001. DOI:10.1115/1.4033269.

[19] HODGES D H, ATILGAN A R, FULTON M V, et al. Free-vibration analysis of composite beams[J]. Journal of the American helicopter society, 1991, 36(3): 36-47. DOI:10.4050/JAHS.36.36.

Dynamic Characteristics of Large Wind Turbine Composite Blades Using Absolute Nodal Coordinate Formulation

HU Xi-Biao, XIA Hong-jian, LI De-yuan

(School of Electromechanical Engineering, Guangdong University of Technology, Guangzhou 510006)

To investigate the nonlinear structural dynamic characteristics of rotating composite blades in large-scale wind turbines, the absolute nodal coordinate method was used to derive the generalized mass and stiffness matrices based on the theory of general continuous media mechanics and the constitutive relationship of composite materials, and a nonlinear dynamic model for rotating composite blades of large-scale wind turbines was established. To describe blade rotation, a large-scale rotating body coordinate system was introduced, and the state equation of the rotating blade was constructed by combining the perturbation principle and linearization method. The model’s accuracy was verified using a nonlinear beam standard example, and the dynamic characteristics under stationary and rotating conditions of the DTU-10MW wind turbine composite blades were analyzed. The results showed that due to the blades’ large scale and nonlinear deformation characteristics, there was significant coupling between the flapwise, edgewise, and torsional modes when the blades were stationary. Due to the dynamic stiffening effect, the frequencies of each order of the blades changed with the increase of rotation speed, significantly impacting the flapwise and edgewise modes.

wind turbine blade; composite material; absolute nodal coordinate method; dynamic characteristic analysis

2095-560X（2023）03-0247-08

TK83；O313.7

A

10.3969/j.issn.2095-560X.2023.03.007

2023-03-08

2023-03-28

国家自然科学基金项目（51776044，51105079）；广东省自然科学基金项目（2020A1515010844）

夏鸿建，E-mail：hjxia@gdut.edu.cn

胡锡彪, 夏鸿建, 李德源.绝对节点坐标法大型风力机复合材料叶片的动态特性[J]. 新能源进展, 2023, 11(3): 247-254.

HU Xi-Biao, XIA Hong-jian, LI De-yuan. Dynamic characteristics of large wind turbine composite blades using absolute nodal coordinate formulation[J]. Advances in new and renewable energy, 2023, 11(3): 247-254.

胡锡彪（1996-），男，硕士研究生，主要从事绝对节点坐标法在风力机叶片中的应用研究。

夏鸿建（1978-），男，副教授，主要从事风力机气动与结构和计算机软件的开发与应用等方面的研究。