彎扭耦合對(duì)風(fēng)力機(jī)葉片氣彈模態(tài)特性的影響*

曲珍壯，夏鴻建，李德源，林偉豪

彎扭耦合對(duì)風(fēng)力機(jī)葉片氣彈模態(tài)特性的影響*

曲珍壯，夏鴻建?，李德源，林偉豪

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

彎扭耦合對(duì)風(fēng)力機(jī)葉片氣動(dòng)彈性響應(yīng)有重要影響。為研究彎扭耦合對(duì)風(fēng)力機(jī)葉片氣彈模態(tài)特性的影響，采用幾何精確梁理論和葉素動(dòng)量定理，建立風(fēng)力機(jī)葉片的氣彈耦合方程，在Beddoes-Leishman動(dòng)態(tài)失速模型的基礎(chǔ)上對(duì)氣彈耦合方程線性化，構(gòu)建特征平衡方程。由復(fù)特征值分析法計(jì)算DTU 10 MW葉片的氣彈頻率和氣彈阻尼比。通過在截面剛度矩陣中引入耦合系數(shù)，考慮彎扭耦合對(duì)氣彈模態(tài)的影響。揮舞順槳耦合會(huì)增加揮舞模態(tài)氣彈頻率，降低揮舞模態(tài)氣彈阻尼比。揮舞失速耦合會(huì)降低揮舞模態(tài)氣彈頻率，增大揮舞模態(tài)氣彈阻尼比。揮舞順槳耦合會(huì)降低風(fēng)力機(jī)葉片的顫振風(fēng)速。

風(fēng)力機(jī)葉片；幾何精確梁；彎扭耦合；顫振分析

0 引 言

彎扭耦合對(duì)風(fēng)力機(jī)葉片氣動(dòng)彈性響應(yīng)有重要影響。葉片在氣動(dòng)力作用下彎曲后，會(huì)改變攻角的大小，進(jìn)而改變?nèi)~片的氣動(dòng)力，形成反饋循環(huán)。彎扭耦合可以由葉片的幾何形狀（幾何耦合）或各向異性葉片材料（材料耦合）引起。幾何耦合是葉片幾何形狀彎曲的結(jié)果（預(yù)彎曲、后掠或負(fù)載偏移），從而引起額外的扭轉(zhuǎn)。材料耦合是葉片表皮或梁帽的纖維方向偏離主方向，從而引起彎曲和扭轉(zhuǎn)的耦合。

彎扭耦合可以用來減輕風(fēng)力機(jī)葉片的負(fù)載。根據(jù)LOBITZ[1]、BOTTASSO[2]等的研究，揮舞扭轉(zhuǎn)順槳耦合會(huì)使風(fēng)力機(jī)葉片的疲勞載荷降低10% ～ 20%。然而彎扭耦合除了降低葉片負(fù)載外，還會(huì)改變?nèi)~片的氣彈模態(tài)特性。LOBITZ[3]研究了彎扭耦合葉片與非彎扭耦合葉片發(fā)生顫振失穩(wěn)臨界風(fēng)速，發(fā)現(xiàn)彎扭耦合葉片發(fā)生顫振失穩(wěn)的臨界風(fēng)速會(huì)比非彎扭耦合葉片降低12%左右。LARWOOD[4]研究后掠葉片的顫振特性，發(fā)現(xiàn)后掠葉片顫振失穩(wěn)的模態(tài)是二階揮舞與一階扭轉(zhuǎn)的耦合模態(tài)。

隨著風(fēng)力機(jī)葉片長(zhǎng)度和柔性的增加，葉片的非線性效應(yīng)越加顯著。KALLES?E[5]采用非線性二階 Bernoulli-Euler梁理論，研究了5 MW葉片彎曲變形對(duì)氣彈穩(wěn)定性的影響，表明葉片的幾何非線性對(duì)其氣彈穩(wěn)定性有明顯的影響。傳統(tǒng)采用歐拉梁結(jié)構(gòu)模型[6]或模態(tài)法[7]構(gòu)建葉片結(jié)構(gòu)模型，難以表達(dá)柔性葉片的幾何非線性變形。由于風(fēng)力機(jī)葉片氣彈耦合狀態(tài)方程的建立比較復(fù)雜，基于頻域的顫振分析研究相對(duì)較少，黃俊東等[8]采用多體動(dòng)力學(xué)方法和葉素動(dòng)量理論的方法建立葉片氣彈耦合方程，研究了美國(guó)可再生能源實(shí)驗(yàn)室（National Renewable Energy Laboratory, NREL）發(fā)布的5 MW風(fēng)力機(jī)氣彈穩(wěn)定性問題。和其他非線性模型相比，幾何精確梁模型可以更為準(zhǔn)確方便地計(jì)算葉片的變形位置以及葉片扭轉(zhuǎn)的角度。通過在截面剛度矩陣中加入耦合項(xiàng)，研究由葉片材料引起的彎扭耦合問題。在結(jié)構(gòu)動(dòng)力學(xué)求解時(shí)，對(duì)結(jié)構(gòu)動(dòng)力學(xué)方程進(jìn)行線性化處理，可以方便地與氣動(dòng)力線性化結(jié)合起來，進(jìn)行特征值分析。

本文首先根據(jù)幾何精確梁理論，結(jié)合葉素動(dòng)量定理，建立風(fēng)力機(jī)葉片氣彈耦合方程，找出葉片在不同工況下的平衡位置。然后對(duì)氣動(dòng)力進(jìn)行線性化，并引入Beddoes-Leishman（B-L）動(dòng)態(tài)失速模型，構(gòu)建特征平衡方程，采用特征值分析方法，研究風(fēng)力機(jī)葉片的氣彈模態(tài)特性。通過在截面剛度矩陣中引入耦合系數(shù)[9]，以丹麥技術(shù)大學(xué)（Technical University of Denmark, DTU）10 MW RWT[10]葉片為參考葉片，研究了揮舞扭轉(zhuǎn)耦合對(duì)氣彈模態(tài)特性的影響。

1 幾何精確梁變形的描述

幾何精確梁采用平截面假設(shè)，梁的位形由梁軸線的位移和相應(yīng)截面的轉(zhuǎn)動(dòng)決定。如圖1有全局坐標(biāo)系（不考慮葉片預(yù)彎預(yù)扭以及截面偏心等因素）、不變形位置坐標(biāo)系（考慮葉片真實(shí)幾何形狀）和變形位置坐標(biāo)系（考慮氣彈耦合后葉片的變形）。

圖1 梁的轉(zhuǎn)動(dòng)坐標(biāo)系

在變形位置，梁上任意一點(diǎn)的位置向量可以寫為

2 結(jié)構(gòu)動(dòng)力學(xué)方程的建立

2.1 幾何精確梁本構(gòu)關(guān)系

幾何精確梁理論基于Reinssner應(yīng)變，根據(jù)參考文獻(xiàn)[12-13]，定義應(yīng)變?yōu)槠絼?dòng)應(yīng)變向量，應(yīng)變?yōu)檗D(zhuǎn)動(dòng)應(yīng)變向量，截面等效應(yīng)力、和應(yīng)變、之間的關(guān)系為

應(yīng)變與變形的關(guān)系為

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

首先建立風(fēng)力機(jī)葉片弱形式動(dòng)力學(xué)平衡方程[14]，然后通過數(shù)值積分和微分方法離散求解。

對(duì)內(nèi)力沿單元軸線積分，得到內(nèi)力虛功為[15]

對(duì)氣動(dòng)力沿軸線積分，重力進(jìn)行體積積分，得到單元的外力虛功

各點(diǎn)的速度進(jìn)行體積分，得到單元的慣性力虛功

根據(jù)虛功原理，單元內(nèi)力、外力和慣性力虛功滿足

通過Gauss-Lobatto積分和微分求積法[14]離散平衡方程，得到基于節(jié)點(diǎn)虛位移表達(dá)的內(nèi)力、外力和慣性力虛功[15]。

單元節(jié)點(diǎn)虛位移表示為

把式（12）代入式（11），并在全局坐標(biāo)系下進(jìn)行組裝，得到等效節(jié)點(diǎn)內(nèi)力、外力和慣性力的離散平衡方程

3 動(dòng)力學(xué)方程線性化

為了構(gòu)建葉片氣彈線性平衡方程，形成特征值求解，需要對(duì)動(dòng)力學(xué)平衡方程進(jìn)行線性化處理。分為兩個(gè)部分，第一部分是結(jié)構(gòu)動(dòng)力學(xué)方程線性化，第二部分是氣動(dòng)力線性化。

3.1 結(jié)構(gòu)動(dòng)力學(xué)線性化

分別對(duì)節(jié)點(diǎn)內(nèi)力、外力和慣性力矩陣線性化[15]得到

將式（15）～ 式（17）代入平衡方程（14）可以得到

3.2 氣動(dòng)力線性化

根據(jù)葉素動(dòng)量定理，作用在葉片第個(gè)截面的氣動(dòng)力可以寫為

圖2 葉片翼型截面

氣動(dòng)載荷的線性化表示為[8]

對(duì)狀態(tài)變量線性化[6]得到

聯(lián)立式（18）、式（22）和式（23），氣彈線性平衡方程可以寫為

4 特征值復(fù)模態(tài)計(jì)算

根據(jù)式（24）建立風(fēng)力機(jī)顫振分析的特征值平衡方程

利用Matlab軟件提供的函數(shù)計(jì)算特征值和特征向量，得到葉片的氣彈頻率、氣彈阻尼比和振型，其中特征值的虛部為氣彈頻率，特征值的實(shí)部為氣彈阻尼比。

式中：Im為復(fù)特征值的虛部；Re為復(fù)特征值的實(shí)部。

5 計(jì)算結(jié)果

5.1 旋轉(zhuǎn)梁算例

旋轉(zhuǎn)梁作為能否處理梁大變形標(biāo)準(zhǔn)算例，受到廣泛的關(guān)注。通過與參考文獻(xiàn)[18]對(duì)比來驗(yàn)證動(dòng)力學(xué)方程的正確性，梁的旋轉(zhuǎn)角度隨時(shí)間變化的函數(shù)為

與只有轉(zhuǎn)動(dòng)沒有變形的剛性梁相比，柔性梁的自由端偏轉(zhuǎn)位移如圖3所示。圖中可見，計(jì)算結(jié)果與文獻(xiàn)[18]吻合，說明本文的動(dòng)力學(xué)方程是正確的。

5.2 運(yùn)行工況下氣彈阻尼和氣彈頻率

由于材料的各向異性，風(fēng)力機(jī)葉片截面耦合特性比較復(fù)雜，為了單獨(dú)研究揮舞扭轉(zhuǎn)耦合對(duì)葉片顫振特性的影響，將丹麥技術(shù)大學(xué)[10]提供的10 MW風(fēng)力機(jī)葉片（截面剛度不包含耦合項(xiàng)）設(shè)為參考葉片。在葉片中間部分某截面的剛度參數(shù)如式（30）。

圖4顯示了運(yùn)行工況下，揮舞扭轉(zhuǎn)耦合對(duì)一階揮舞模態(tài)氣彈頻率和氣彈阻尼的影響。在風(fēng)力機(jī)運(yùn)行工況下，三種葉片的一階揮舞氣彈頻率和氣彈阻尼比的變化趨勢(shì)基本相同。具體為在風(fēng)速4 ～ 7 m/s，一階揮舞氣彈阻尼比保持不變，在風(fēng)速7 ～ 12 m/s，為上升趨勢(shì)，風(fēng)速12 m/s后保持平穩(wěn)，此時(shí)風(fēng)力機(jī)達(dá)到額定轉(zhuǎn)速運(yùn)行工況下，一階揮舞氣彈頻率變化趨勢(shì)與氣彈阻尼比相似。相比于參考葉片，揮舞順槳耦合增大了一階揮舞氣彈頻率，降低了一階揮舞氣彈阻尼比，揮舞失速耦合降低了一階揮舞氣彈頻率，增大了一階揮舞氣彈阻尼比。

圖4 運(yùn)行工況下一階揮舞氣彈頻率和氣彈阻尼比

揮舞扭轉(zhuǎn)耦合對(duì)一階擺振氣彈頻率和氣彈阻尼的影響如圖5所示。在運(yùn)行工況下，三種葉片的一階擺振氣彈頻率基本重合，一階氣彈阻尼也比較接近，表明揮舞扭轉(zhuǎn)耦合對(duì)一階擺振模態(tài)氣彈特性影響較小。在運(yùn)行工況下，三種葉片的一階擺振氣彈頻率變化較小，變化區(qū)間為0.945 ～ 0.979 Hz。一階擺振氣彈阻尼比變化趨勢(shì)為，在風(fēng)速為4 ～ 11 m/s時(shí)緩慢上升，然后短暫下降后又繼續(xù)保持上升，一直到風(fēng)速25 m/s時(shí)達(dá)到最大值，此時(shí)氣彈阻尼比大約為4%，與風(fēng)速11 m/s時(shí)的氣彈阻尼接近。

圖5 運(yùn)行工況下一階擺振氣彈頻率和氣彈阻尼比

5.3 定葉尖速比工況下顫振風(fēng)速

為了研究彎扭耦合對(duì)風(fēng)力機(jī)葉片顫振風(fēng)速的影響，將風(fēng)力機(jī)葉片運(yùn)行工況設(shè)為定葉尖速比10，即風(fēng)力機(jī)葉尖速度為來流風(fēng)速的10倍，并將葉片的槳距角設(shè)為0。結(jié)果如圖6。

圖6 定葉尖速比工況下一階扭轉(zhuǎn)氣彈頻率和氣彈阻尼比

揮舞順槳耦合葉片一階扭轉(zhuǎn)氣彈阻尼比在轉(zhuǎn)速1.38 rad/s時(shí)出現(xiàn)負(fù)值，此時(shí)一階扭轉(zhuǎn)氣彈頻率為7.51 Hz，參考葉片大約在轉(zhuǎn)速1.97 rad/s時(shí)出現(xiàn)負(fù)值，此時(shí)一階扭轉(zhuǎn)氣彈頻率為8.47 Hz，揮舞扭轉(zhuǎn)失速耦合葉片沒有出現(xiàn)負(fù)值。表明揮舞順槳耦合會(huì)降低風(fēng)力機(jī)葉片的顫振風(fēng)速。這種現(xiàn)象的原因可能是，揮舞順槳耦合葉片變形小，造成扭轉(zhuǎn)剛度小于參考葉片，根據(jù)文獻(xiàn)[19]，扭轉(zhuǎn)剛度降低會(huì)使葉片的顫振風(fēng)速降低。

6 結(jié) 論

采用非線性幾何精確梁模型和葉素動(dòng)量定理建立風(fēng)力機(jī)葉片的氣彈模型，在Beddoes-Leishman動(dòng)態(tài)失速模型的基礎(chǔ)上對(duì)氣彈模型線性化，構(gòu)建特征平衡方程，通過特征值分析，計(jì)算風(fēng)力機(jī)葉片的氣彈頻率和氣彈阻尼比。通過在截面剛度矩陣中引入耦合系數(shù)，研究風(fēng)力機(jī)葉片揮舞扭轉(zhuǎn)耦合對(duì)氣彈模態(tài)特性的影響。主要結(jié)論如下：

（1）在運(yùn)行工況下，揮舞扭轉(zhuǎn)耦合主要影響揮舞模態(tài)，對(duì)擺振模態(tài)影響較小。揮舞順槳耦合增大了一階揮舞模態(tài)氣彈頻率，降低了一階揮舞模態(tài)氣彈阻尼比；揮舞失速耦合降低了一階揮舞模態(tài)氣彈頻率，增大了一階揮舞模態(tài)氣彈阻尼比。

本文建立的氣彈模型還可以用于分析風(fēng)力機(jī)葉片其他的氣彈問題。

[1] LOBITZ D W, VEERS P S, LAINO D J. Performance of twist-coupled blades on variable speed rotors[R]. Albuquerque: Sandia National Labs, 1999.

[2] BOTTASSO C L, CAMPAGNOLO F, CROCE A, et al. Optimization-based study of bend-twist coupled rotor blades for passive and integrated passive/active load alleviation[J]. Wind energy, 2013, 16(8): 1149-1166. DOI: 10.1002/we.1543.

[3] LOBITZ D W. Aeroelastic stability predictions for a MW-sized blade[J]. Wind energy, 2004, 7(3): 211-224. DOI: 10.1002/we.120.

[4] LARWOOD S. Flutter of variations on a 5 MW swept wind turbine blade[J]. Journal of solar energy engineering, 2016, 138(2): 024504. DOI: 10.1115/1.4032475.

[5] KALLES?E B S. Effect of steady deflections on the aeroelastic stability of a turbine blade[J]. Wind energy, 2011, 14(2): 209-224. DOI: 10.1002/we.413.

[6] HANSEN M H. Aeroelastic properties of backward swept blades[C]//49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition. Orlando: AIAA, 2011. DOI: 10.2514/6.2011-260.

[7] 劉雄, 李鋼強(qiáng), 陳嚴(yán), 等. 水平軸風(fēng)力機(jī)葉片動(dòng)態(tài)響應(yīng)分析[J]. 機(jī)械工程學(xué)報(bào), 2010, 46(12): 128-134, 141. DOI: 10.3901/JME.2010.12.128.

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

[9] LOBITZ D W, VEERS P S. Aeroelastic behavior of twist-coupled HAWT blades[C]//1998 ASME Energy Symposium. Reno: ASME, 1998: 75-83. DOI: 10.2514/6.1998-29.

[10] BAK C, ZAHLE F, BITSCHE R, et al. Description of the DTU 10 MW reference wind turbine[R]. Roskilde: DTU, 2013.

[11] 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.

[12] REISSNER E. On fifinite deformation of space-curved beams[J]. Z. angew. Math Phys., 1981, 32: 734-744. https://doi.org/10.1007/BF00946983.

[13] SIMO J C. A finite strain beam formulation. The three-dimensional dynamic problem. Part I[J]. Computer methods in applied mechanics and engineering, 1985, 49(1): 55-70. DOI: 10.1016/0045-7825(85)90050-7.

[14] ZHANG R, ZHONG H Z. Weak form quadrature element analysis of planar slender beams based on geometrically exact beam theory[J]. Archive of applied mechanics, 2013, 83(9): 1309-1325. DOI: 10.1007/s00419-013-0748-3.

[15] 呂品, 廖明夫, 徐陽, 等. 基于幾何精確梁理論的風(fēng)力機(jī)葉片單元[J]. 太陽能學(xué)報(bào), 2015, 36(10): 2422-2428. DOI: 10.3969/j.issn.0254-0096.2015.10.017.

[16] ZHONG H Z, YU T. Flexural vibration analysis of an eccentric annular Mindlin plate[J]. Archive of applied mechanics, 2007, 77(4): 185-195. DOI: 10.1007/s00419-006-0083-z.

[17] HANSEN M H, GAUNAA M, MADSEN H A. A Beddoes-Leishman type dynamic stall model in state-space and indicial formulations[R]. Denmark: Risφ National Laboratory, 2004.

[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] ZHOU X D, HUANG K F, LI Z. Effects of bend-twist coupling on flutter limits of composite wind turbine blades[J]. Composite structures, 2018, 192: 317-326. DOI: 10.1016/j.compstruct.2018.02.071.

Effect of Bending-Torsion Coupling on Aeroelastic Modal Characteristics of Wind Turbine Blades

QU Zhenzhuang, XIA Hongjian?, LI Deyuan, LIN Weihao

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

The bending-torsion coupling has an important influence on the aeroelastic response of wind turbine blades. To study the effect of bending-torsion coupling on the aeroelastic modal characteristics of wind turbine blades, the geometrically exact beam theory and blade element momentum theory are used to establish the aeroelastic coupling equation for wind turbine blades. Based on the Beddoes-Leishman dynamic stall model, the aeroelastic coupling equation is linearized and the characteristic equilibrium equation is constructed. The DTU 10 MW blade’s aeroelastic frequencies and damping ratio are calculated using complex eigenvalue analysis. The influence of bending-torsion coupling on the aeroelastic modes is considered by introducing coupling coefficients into the section stiffness matrix. The flap-twist to feather coupling increases the flapwise aeroelastic frequency and decreases the flapwise aeroelastic damping ratio, while the flap-twist to stall coupling decreases the flapwise aeroelastic frequency and increases the flapwise aeroelastic damping ratio. The flap-twist to feather can reduce the flutter speed of wind turbine blades.

wind turbine blade; geometrically exact beam; bending-torsion coupling; flutter analysis

2095-560X（2023）04-0333-07

TK83；O313.7

A

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

2023-03-08

2023-04-12

國(guó)家自然科學(xué)基金項(xiàng)目（51776044，51105079）；廣東省自然科學(xué)基金項(xiàng)目（2020A1515010844）

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

曲珍壯, 夏鴻建, 李德源, 等. 彎扭耦合對(duì)風(fēng)力機(jī)葉片氣彈模態(tài)特性的影響[J]. 新能源進(jìn)展, 2023, 11(4): 333-339.

：QU Zhenzhuang, XIA Hongjian, LI Deyuan, et al. Effect of bending-torsion coupling on aeroelastic modal characteristics of wind turbine blades[J]. Advances in new and renewable energy, 2023, 11(4): 333-339.

曲珍壯（1996-），男，碩士研究生，主要從事基于幾何精確梁模型的風(fēng)力機(jī)葉片研究。

夏鴻建（1978-），男，博士，副教授，碩士生導(dǎo)師，主要從事風(fēng)力機(jī)氣動(dòng)與結(jié)構(gòu)等研究。