基于幾何精確梁的風(fēng)力機(jī)葉片動(dòng)力學(xué)建模及動(dòng)態(tài)特性分析*

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

基于幾何精確梁的風(fēng)力機(jī)葉片動(dòng)力學(xué)建模及動(dòng)態(tài)特性分析*

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

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

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

風(fēng)力機(jī)葉片；幾何精確梁；動(dòng)力剛化；模態(tài)分析

0 引 言

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

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

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

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

1　風(fēng)力機(jī)葉片非線(xiàn)性單元模型

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

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

圖1 風(fēng)力機(jī)葉片截面

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

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

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

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

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

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

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

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

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

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

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

1.3　動(dòng)力學(xué)平衡方程

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

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

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

1.4　平衡方程的線(xiàn)性化

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

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

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

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

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

2 算例驗(yàn)證

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

2.1　純彎曲梁

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

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

圖3 純彎曲梁變形

2.2　葉片靜撓度

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

（1）所建的幾何精確梁模型能對(duì)柔性梁進(jìn)行準(zhǔn)確的動(dòng)力學(xué)分析。

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

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

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

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

[2] 趙俊杰. 2MW風(fēng)力發(fā)電機(jī)葉片設(shè)計(jì)與分析[D]. 邯鄲: 河北工程大學(xué), 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] 黃俊東, 夏鴻建, 李德源, 等. 大型風(fēng)力機(jī)柔性葉片非線(xiàn)性氣彈模態(tài)分析[J]. 機(jī)械工程學(xué)報(bào), 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] 莫文威. 基于多體模型的水平軸風(fēng)力機(jī)氣彈耦合分析[D]. 廣州: 廣東工業(yè)大學(xué), 2013.

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

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

[9] 呂品, 廖明夫, 尹堯杰. 考慮幾何非線(xiàn)性的風(fēng)力機(jī)葉片氣彈模型[J]. 機(jī)械科學(xué)與技術(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] 黃正. 基于等幾何配點(diǎn)法的幾何精確Euler-Bernoulli梁幾何非線(xiàn)性分析[D]. 武漢: 華中科技大學(xué), 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-），男，碩士研究生，主要從事基于幾何精確梁模型的風(fēng)力機(jī)葉片研究。

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