謝 浩，李德源，曲珍壯，黃 維

基于幾何精確梁的風力機葉片動力學建模及動態(tài)特性分析*

謝 浩，李德源?，曲珍壯，黃 維

（廣東工業(yè)大學 機電工程學院，廣州 510006）

基于幾何精確梁理論，結(jié)合廣義-時間預(yù)測法的迭代算法，考慮葉片復(fù)合鋪層材料的各向異性特性，建立了大型風力機葉片的幾何非線性動力學模型并導(dǎo)出了相應(yīng)的特征方程，編制了數(shù)值仿真程序。通過對幾何非線性梁標準算例和某10 MW柔性風力機葉片動力特性的模擬分析，驗證了動力學模型的正確性，以及幾何精確梁模型對分析葉片幾何非線性大變形及其所導(dǎo)致的非線性動力學效應(yīng)的有效性。葉片在靜止和轉(zhuǎn)動工況下的模態(tài)分析結(jié)果表明，在動力剛化效應(yīng)作用下，葉片的固有頻率會隨著轉(zhuǎn)速的增加而增大，動力剛化效應(yīng)在揮舞方向比在擺振方向更明顯，在低階模態(tài)比在高階模態(tài)更明顯。

風力機葉片；幾何精確梁；動力剛化；模態(tài)分析

0 引 言

隨著全球能源危機與環(huán)境污染等問題的日益加劇，以風能為代表的各種清潔可再生能源逐漸被世界各地所重視。為了提高風力機的風能利用率、降低發(fā)電成本，風力機葉片尺寸的大型化已經(jīng)成為風力機的重要發(fā)展方向[1]。

風力機葉片具有細長結(jié)構(gòu)的特點，在研究分析葉片結(jié)構(gòu)時，葉片通常被近似假設(shè)為梁模型來進行研究。在工程上，葉片梁模型的分析常使用基于小變形假設(shè)的軟件，如Bladed、HAWC2等[2-3]。對于短葉片，小變形假設(shè)能獲得較為準確的結(jié)果和較高的計算效率。然而隨著風力機的大型化，小變形假設(shè)已經(jīng)難以表達更為細長的柔性葉片大變形所帶來的幾何非線性現(xiàn)象[4]。對大變形葉片的幾何非線性問題研究已經(jīng)成為風力機領(lǐng)域的重點研究方向之一[5]。

現(xiàn)在常用的非線性梁模型主要有多體動力學模型[6]、絕對節(jié)點坐標模型[7]以及幾何精確梁模型[8-10]。幾何精確梁模型能以較少的節(jié)點自由度獲得較高的計算精度，能準確地用于葉片的非線性動力學分析[11]，采用線性化系統(tǒng)得到的解析形式的模態(tài)便于獲得系統(tǒng)的低自由度的非線性模態(tài)振動方程，為大型柔性葉片的非線性振動研究提供了有效途徑。

本文基于幾何精確梁理論，使用四元數(shù)表達截面轉(zhuǎn)動，采用Gauss-Lobatto積分法則，結(jié)合具有二階精度的廣義-時間預(yù)測法[12]迭代算法，考慮葉片復(fù)合鋪層材料的各向異性特性，建立大型風力機葉片的非線性動力學模型，并對該模型進行線性化求解。通過與丹麥技術(shù)大學（Technical University of Denmark, DTU）所提供的10 MW風力機葉片數(shù)據(jù)[13]作對比，對葉片動力學模型的正確性進行驗證，研究幾何精確梁模型對分析非線性大變形和動力學特性的可行性，為大型風力機葉片的結(jié)構(gòu)分析及后續(xù)的氣彈耦合分析提供一個計算效率高且準確的非線性動力學模型。

1　風力機葉片非線性單元模型

1.1　葉片構(gòu)型及變形

風力機葉片是由一定幾何形狀的翼型與多種材料鋪層所形成的，由于截面形狀較為特殊，且復(fù)合材料鋪層的各向異性同一截面的鋪層厚度不一，因此具有重心、剪心、扭心不共點的特性和材料的各向異性特性。圖1為10 MW風力機葉片的主要截面形狀。

圖1 風力機葉片截面

在葉片根部建立全局坐標系（G，G，G），G與風輪平面垂直，G與葉根平面垂直。以各截面的剪心作為原點建立各截面的局部彈性坐標系（E，E，E），三個坐標軸的指向與全局坐標系一致，定義此時的葉片狀態(tài)為直葉片的參考構(gòu)型。

將各截面的剪心連成一線，該連線則為葉片的軸線?；趲缀尉_梁的平截面假設(shè)，葉片的變形和葉片的預(yù)彎預(yù)扭初始構(gòu)型，可視作是由各截面EE平面上的平移與繞E軸的轉(zhuǎn)動所形成的。

1.2　截面轉(zhuǎn)動描述

幾何精確梁理論基于由Timoshenko梁理論發(fā)展而成的Reissner梁理論[14]，同樣采用平截面假設(shè)，即梁的位移形變由梁軸線位移與截面轉(zhuǎn)動所形成，還使用三維有限轉(zhuǎn)動理論，精確得到梁在大位移大變形條件下的應(yīng)變?位形關(guān)系，如圖2所示。

圖2 截面轉(zhuǎn)動在梁變形前后的變化關(guān)系

圖2中，（1,2,3）為全局坐標系的基向量，當梁發(fā)生變形后，同一截面的局部坐標系基向量(0)變成了(1)，其中的轉(zhuǎn)動變化可通過三階方向余弦矩陣來表示

方向余弦矩陣有9個分量，計算效率低，截面轉(zhuǎn)動的更新也比較復(fù)雜，歐拉四元數(shù)是空間三維轉(zhuǎn)動非奇異的最少參數(shù)表達，因此本文采用歐拉四元數(shù)描述截面的轉(zhuǎn)動[15]。四元數(shù)與方向余弦矩陣的關(guān)系為

由歐拉四元數(shù)性質(zhì)可知，0可通過計算所得，則可以用節(jié)點線位移、節(jié)點歐拉參數(shù)，定義節(jié)點的位移向量為

通過四元數(shù)法與增量轉(zhuǎn)動法得到增量四元數(shù)可實現(xiàn)旋轉(zhuǎn)的可加性，用于確定截面轉(zhuǎn)動。四元數(shù)增量可以使用指數(shù)映射的方式

節(jié)點的角速度和角加速度與四元數(shù)的關(guān)系為

對于10 MW葉片截面節(jié)點的選擇，本文采用Legendre-Gauss-Lobatto積分法則中的積分點。

1.3　動力學平衡方程

幾何精確梁理論采用一階假設(shè)的Reissner應(yīng)變，定義Reissner應(yīng)變?yōu)槠絼討?yīng)變向量，應(yīng)變?yōu)檗D(zhuǎn)動應(yīng)變向量，則截面等效力、與應(yīng)變、的關(guān)系，即本構(gòu)關(guān)系為

根據(jù)DTU所提供的10 MW風力機葉片數(shù)據(jù)，建立葉片三維模型，再根據(jù)Lobatto積分點選擇相應(yīng)葉片截面，提取對應(yīng)葉片截面的有限元單元幾何數(shù)據(jù)，加上葉片鋪層材料的數(shù)據(jù)，導(dǎo)入專業(yè)分析截面軟件VABS計算得到葉片截面剛度矩陣[16]。

結(jié)合虛功原理，得到葉片單元的內(nèi)力、外力和慣性力的虛功平衡方程為

1.4　平衡方程的線性化

當葉片從靜止加速到所設(shè)定的轉(zhuǎn)速且方程達到平衡狀態(tài)時，在平衡位置處線性化質(zhì)量矩陣及剛度矩陣，由特征值分解計算轉(zhuǎn)動葉片的模態(tài)的公式如下

1.5 穩(wěn)態(tài)位置求解

葉片非線性系統(tǒng)動力學平衡方程的求解需要使用具有預(yù)測?校正特性的方法，此處采用廣義-時間預(yù)測法[12]。

廣義-時間預(yù)測?校正的方法對增量平衡方程的求解為

由位移增量得到速度和加速度增量為：

2 算例驗證

本節(jié)根據(jù)以上幾何精確梁理論公式，編寫對應(yīng)的MATLAB程序，構(gòu)建幾何精確梁模型，并用以下三個算例驗證幾何精確梁理論對靜力學及動力學特性描述的準確性。

2.1　純彎曲梁

為了體現(xiàn)幾何精確梁理論描述非線性大變形的準確性，以純彎曲懸臂梁的靜力學變形作為算例。根據(jù)純彎理論可得，梁僅在自由端受到一定彎矩后，梁軸線會變形成一段圓弧，其對應(yīng)曲率為()，為彎矩，為梁截面剛度?，F(xiàn)對梁長1 m的懸臂梁自由端施加的彎矩，則得到變形后梁軸線的曲率解析值為。

圖3為幾何精確梁理論分析不同的所對應(yīng)的彎矩對懸臂梁變形的結(jié)果，與相應(yīng)的解析值作為對比，兩者結(jié)果基本符合。

圖3 純彎曲梁變形

2.2　葉片靜撓度

為了體現(xiàn)幾何精確梁理論能對結(jié)構(gòu)不規(guī)則且具有復(fù)合材料鋪層的10 MW風力機葉片進行準確的變形描述，對葉片末端分別在揮舞（垂直于風輪平面的方向）、擺振方向（風輪平面所在的方向）施加?100 000 N的力，利用幾何精確梁理論計算其靜撓度，并與使用軟件Abaqus的計算結(jié)果做對比，分別如圖4所示，兩者結(jié)果基本符合。

2.3　旋轉(zhuǎn)梁的動力特性

為了體現(xiàn)幾何精確梁理論對旋轉(zhuǎn)物體動力學特性的準確性，以定軸旋轉(zhuǎn)梁為研究對象，此矩形截面梁的幾何尺寸比例與葉片的尺寸比例相近，長度= 7 m，橫截面寬= 0.5 m、高= 0.15 m，材料為葉片結(jié)構(gòu)常用的環(huán)氧樹脂，彈性模量= 4.0 GPa，剪切模量= 1.481 4 GPa，密度= 1 140 kg/m3，泊松比= 0.35。

表1為使用幾何精確梁模型，結(jié)合21個Lobatto積分點，計算該定軸旋轉(zhuǎn)梁分別在0 r/min、10 r/min、20 r/min轉(zhuǎn)速下的前五階固有頻率結(jié)果，以及使用軟件ANSYS Workbench對該梁做有限元分析的結(jié)果，該有限元模型采用的是Beam188梁單元，該梁單元考慮了剪切效應(yīng)，每個節(jié)點有3個平動自由度，3個轉(zhuǎn)動自由度，共6個自由度。算例單元數(shù)為234個。

表1 旋轉(zhuǎn)梁的固有頻率對比

表1的對比結(jié)果表明了幾何精確梁模型可以用較少的自由度得出旋轉(zhuǎn)梁的高精度動力學特性。

3 葉片的模態(tài)分析

上一節(jié)已證明幾何精確梁理論對描述復(fù)合材料非線性大變形及動力學特性表達的準確性，現(xiàn)使用幾何精確梁理論，同樣基于21個Lobatto積分點，對靜止的風力機葉片進行模態(tài)分析，研究對象采用DTU所提供的10 MW風力機葉片，得到各階模態(tài)相應(yīng)固有頻率如表2所示，相應(yīng)振型如圖5所示。

表2 靜止葉片的模態(tài)分析結(jié)果

圖5 靜止葉片的前六階振型

進一步通過幾何精確梁理論對風力機葉片分別在轉(zhuǎn)速為3.0 r/min、6.0 r/min、9.6 r/min（額定工況轉(zhuǎn)速）、15.0 r/min、19.2 r/min時進行模態(tài)分析，得到各轉(zhuǎn)速對應(yīng)的各階模態(tài)固有頻率如表3所示，對應(yīng)變化趨勢如圖6所示。

表3 旋轉(zhuǎn)葉片的固有頻率對比

圖6 旋轉(zhuǎn)葉片的固有頻率對比

動力剛化效應(yīng)是指物體在轉(zhuǎn)動過程中自身的固有頻率隨著轉(zhuǎn)速的增加而升高。由表3和圖6可得，隨著葉片轉(zhuǎn)速的增大，葉片的固有頻率越來越高，動力剛化效應(yīng)變得更加明顯。當葉片從靜止加速到額定工況轉(zhuǎn)速9.6 r/min時，前六階固有頻率的增量分別為6.14%、2.12%、2.21%、0.76%、1.04%、0.44%，而加速到兩倍額定工況轉(zhuǎn)速19.2 r/min時，增量分別為22.16%、8.09%、8.49%、2.98%、4.04%、1.70%。由此可見，葉片的揮舞方向的動力剛化效應(yīng)比擺振方向更加明顯，且動力剛化效應(yīng)對一階揮舞及一階擺振的影響比對高階振型的影響更加大。

4 結(jié) 論

基于幾何精確梁理論，結(jié)合廣義-時間預(yù)測法，建立梁的幾何非線性動力學模型，并通過不同的算例驗證了該模型對描述非線性大變形及動力學特性的準確性。然后在此模型的基礎(chǔ)上，加入10 MW風力機的數(shù)據(jù)，考慮葉片復(fù)合鋪層材料的各向異性特性，建立大型風力機葉片的非線性動力學模型，將葉片非線性動力學模型線性化，對各轉(zhuǎn)速的葉片進行模態(tài)分析。通過基于幾何精確梁理論對大型風力機葉片的建模及動力學特性分析可得出以下結(jié)論：

（1）所建的幾何精確梁模型能對柔性梁進行準確的動力學分析。

（2）模型考慮葉片復(fù)合鋪層材料的各向異性特性，能對大型風力機葉片進行準確的模態(tài)分析。

（3）大型風力機葉片隨著轉(zhuǎn)速的增大，其動力剛化效應(yīng)越明顯，揮舞方向的動力剛化效應(yīng)比擺振方向更明顯，低階振型動力剛化效應(yīng)比高階振型更明顯。

（4）幾何精確梁模型能夠?qū)Υ笮惋L力機葉片做非線性大變形及動力學特性分析，且能以較少的自由度對葉片進行計算效率高且準確的動力學特性分析，為葉片非線性振動和氣彈耦合響應(yīng)的研究提供了可行方案。

[1] 樂威. 新能源背景下我國風力發(fā)電現(xiàn)狀和未來發(fā)展方向探索[J]. 綠色環(huán)保建材, 2020(11): 165-166. DOI: 10.16767/j.cnki.10-1213/tu.2020.011.080.

[2] 趙俊杰. 2MW風力發(fā)電機葉片設(shè)計與分析[D]. 邯鄲: 河北工程大學, 2018.

[3] PASSON P, KüHN M, BUTTERFIELD S, et al. OC3—benchmark exercise of aero-elastic offshore wind turbine codes[J]. Journal of physics: conference series, 2007, 75: 012071. DOI: 10.1088/1742-6596/75/1/012071.

[4] 黃俊東, 夏鴻建, 李德源, 等. 大型風力機柔性葉片非線性氣彈模態(tài)分析[J]. 機械工程學報, 2020, 56(14): 180-187. DOI: 10.3901/JME.2020.14.180.

[5] HANSEN M O L, S?RENSEN J N, VOUTSINAS S, et al. State of the art in wind turbine aerodynamics and aeroelasticity[J]. Progress in aerospace sciences, 2006, 42(4): 285-330. DOI: 10.1016/j.paerosci.2006.10.002.

[6] 莫文威. 基于多體模型的水平軸風力機氣彈耦合分析[D]. 廣州: 廣東工業(yè)大學, 2013.

[7] 張海波. 絕對節(jié)點坐標法在大型風力機葉片結(jié)構(gòu)分析中的應(yīng)用研究[D]. 廣州: 廣東工業(yè)大學, 2021. DOI: 10.27029/d.cnki.ggdgu.2021.000773.

[8] 肖乃佳. 基于幾何精確梁理論的框架的弱形式求積元分析[D]. 北京: 清華大學, 2011.

[9] 呂品, 廖明夫, 尹堯杰. 考慮幾何非線性的風力機葉片氣彈模型[J]. 機械科學與技術(shù), 2015, 34(12): 1805-1812. DOI: 10.13433/j.cnki.1003-8728.2015.1201.

[10] REISSNER E. On one-dimensional large-displacement finite-strain beam theory[J]. Studies in applied mathematics, 1973, 52(2): 87-95. DOI: 10.1002/sapm197352287.

[11] ROMERO I. A comparison of finite elements for nonlinear beams: the absolute nodal coordinate and geometrically exact formulations[J]. Multibody system dynamics, 2008, 20(1): 51-68. DOI: 10.1007/s11044-008-9105-7.

[12] 黃正. 基于等幾何配點法的幾何精確Euler-Bernoulli梁幾何非線性分析[D]. 武漢: 華中科技大學, 2017. DOI: 10.7666/d.D01313662.

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

[14] SHABANA A A, YAKOUB R Y. 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.

[15] ZUPAN E, SAJE M, ZUPAN D. The quaternion-based three-dimensional beam theory[J]. Computer methods in applied mechanics and engineering, 2009, 198(49/52): 3944-3956. DOI: 10.1016/j.cma.2009.09.002.

[16] YU W B, HODGES D H, HO J C. Variational asymptotic beam sectional analysis-An updated version[J]. International journal of engineering science, 2012, 59: 40-64. DOI: 10.1016/j.ijengsci.2012.03.006.

Wind Turbine Blade Dynamics Modeling and Dynamic Characteristics Analysis Based on Geometrically Exact Beam Theory

XIE Hao, LI De-yuan, QU Zhen-zhuang, HUANG Wei

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

Based on the geometrically exact beam theory, combined with the iterative algorithm of generalized-time prediction method, the geometric nonlinear dynamics model of large wind turbine blade was established considering the anisotropic characteristics of composite laminates. The corresponding characteristic equations were derived and a numerical simulation program was developed. The correctness of the dynamic model and the effectiveness of the geometrically exact beam model in analyzing the large deformation of the geometrically nonlinear blade and its nonlinear dynamic effects were verified by the simulation analysis of the standard geometric nonlinear beam model and the dynamic characteristics of a 10 MW flexible wind turbine blade. The results of the modal analysis of the blade at rest and rotation conditions showed that the natural frequency of the blade increased with the increase of the rotational speed under the dynamic stiffening effect, and the dynamic stiffening effect was more obvious in flapwise direction than in edgewise direction, and more obvious in the low order modal than in the high order modal.

wind turbine blade; geometrically exact beam; dynamic stiffening; modal analysis

2095-560X（2022）03-0203-06

TK83；O313.7

A

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

2022-03-11

2022-04-06

李德源，E-mail：lidey@gdut.edu.cn

謝 浩（1997-），男，碩士研究生，主要從事基于幾何精確梁模型的風力機葉片研究。

李德源（1965-），男，教授，碩士生導(dǎo)師，主要從事大型風力機氣動與結(jié)構(gòu)分析、風力機系統(tǒng)測試和計算機軟件的開發(fā)與應(yīng)用等方面的研究。