曹九發(fā)， 王同光， 王　梟

(1. 揚州大學 水利與能源動力工程學院，江蘇 揚州　225127; 2.南京航空航天大學 江蘇省風力機設計高技術(shù)研究重點實驗室,南京　210016)

?

大型風力機的幾何非線性動態(tài)響應研究

曹九發(fā)1， 王同光2， 王梟2

(1. 揚州大學 水利與能源動力工程學院，江蘇 揚州225127; 2.南京航空航天大學 江蘇省風力機設計高技術(shù)研究重點實驗室,南京210016)

摘要：隨著風力機大型化、柔性化發(fā)展，風力機葉片的變形越來越大，幾何非線性效應也越來越明顯。采用自由渦尾跡方法提高風力機氣動載荷計算精度，應用有限元方法建立了風力機葉片結(jié)構(gòu)梁模型。應用伽遼金方法建立葉片的幾何非線性剛度矩陣，同時構(gòu)建了葉片的非線性動力學方程。采用Newmark直接積分法進行時間步長推進，并采用修正的Newton-Raphson方法進行增量迭代計算，實現(xiàn)了大型風力機葉片氣動與非線性結(jié)構(gòu)耦合計算模擬。最后，通過美國可再生能源實驗室的 Phase VI風力機實驗數(shù)據(jù)驗證了計算模型的準確性和可靠性。同時以大型風力機NH1500為算例，計算了葉片的線性與非線性動態(tài)響應，分析了葉片幾何非線性對動態(tài)響應和氣動性能的影響。這對于提高大型風力機設計水平和各個部件載荷計算準確度具有重要意義。

關(guān)鍵詞：風力機；自由渦尾跡方法；幾何非線性；Newton-Raphson方法；動態(tài)響應

隨著風力機大型化、柔性化發(fā)展，其葉片在氣動力的作用下，揮舞和擺振方向變形進一步增大，這不僅會影響風輪氣動性能，而且會改變?nèi)~片表面載荷分布，影響結(jié)構(gòu)安全性，在一定情況下，甚至會引起結(jié)構(gòu)破壞。因此，能夠準確地計算出風力機的動態(tài)響應，對風力機的發(fā)展起著至關(guān)重要的作用。在葉片大尺度變形的情況下，將其簡化為線性梁模型的方法不再適用，而采用幾何非線性[1-2]假設能更準確地計算風力機動態(tài)響應。

葉素動量理論簡單快捷，但準確度不高。CFD方法雖然能夠較為準確的模擬葉片上的復雜流動，但計算效率低，消耗資源巨大。渦尾跡方法具有旋渦特性，計算成本較CFD更小，準確度較葉素動量理論更高，而且能夠較好的模擬風力機尾流場。預定渦尾跡方法是根據(jù)大量的尾流場實驗數(shù)據(jù)，建立尾跡形狀的經(jīng)驗描述函數(shù)。相比之下，自由渦尾跡方法[3-5]不需要渦元位置的先驗數(shù)據(jù)，而且尾跡渦元允許在當?shù)厮俣葓龅挠绊懴伦杂勺冃??？紤]到葉片變形會導致風場及葉片氣動外形變化，自由渦尾跡方法更適用于風力機氣動結(jié)構(gòu)耦合計算。

將風力機葉片簡化為歐拉梁模型，其響應計算方法可分為工程梁方法，有限元方法[6-7]和多體動力學方法[8-10]。有限元方法求解動力學方程時，可采用振型疊加法[11]和直接積分法。振型疊加法對動力學方程進行模態(tài)分析，能夠有效地減少結(jié)構(gòu)模型的自由度數(shù)和矩陣規(guī)模，提高計算速度，但只適用于線性計算。直接積分法如Newmark方法，計算量比振型疊加法較大，但在節(jié)點數(shù)較少時，二者差別不明顯，而且直接積分法適用于非線性計算。本文在Newmark直接積分法的基礎上采用修正的Newton-Raphson增量計算法[12-13]，用來計算非線性動態(tài)響應。

本文采用自由渦尾跡方法計算了風力機的氣動載荷，并與實驗值進行對比，驗證了氣動計算方法的可靠性與準確性。然后，通過自由渦尾跡方法與有限元方法耦合計算，實現(xiàn)了大型風力機NH1500的動態(tài)響應求解。

1自由渦尾跡方法

1.1尾流場描述

自由渦尾跡方法中，對風力機流場作不可壓和位流假設，氣動模型可以簡化為來流、葉片附著渦線和自由渦面的總和。如圖1所示，葉片附著渦線位于1/4弦線處，并采用“arc-cosine”法離散，每段附著渦線代替每段葉素，葉素控制點位于3/4弦線處，從而葉片被模擬成一個Weissinger-L升力面模型。自由渦面是由葉片尾緣拖出渦線形成，可分為尾隨渦線和脫體渦線，分別模擬附著環(huán)量在空間和時間上的變化。尾隨渦強度定義為相鄰葉素的附著環(huán)量之差：

(1)

式中：i=1,2,…,NE，體現(xiàn)不同葉素，Γt為尾隨渦環(huán)量強度，Γb為葉素附著渦環(huán)量強度。

脫體渦強度為相鄰方位角上葉素附著環(huán)量之差，則第j個方位角的脫體渦強度為：

(2)

式中：j=1,2,…,NT，體現(xiàn)不同方位角，Γs為脫體渦環(huán)量強度。

圖1　風力機尾跡離散描述示意圖Fig.1 Wake discrete description of wind turbine

1.2渦線控制方程

本文模型中每根渦線均在遠場截斷，流場中渦線隨當?shù)亓魉僖苿訒杂删砥稹u線的偏微分控制方程可寫為：

(3)

推導尾跡控制方程中的對時間步微分方程的差分，令Δψ=Δζ，可得到控制方程的離散格式：

預估步：

(4)

校正步：

2幾何非線性動力學方程

風力機向大型化、柔性化方向發(fā)展，葉片變形進一步增大，盡管應變很小，沒有超過彈性極限，但是位移較大，為了更準確的計算風力機動態(tài)響應，動力學方程應該建立在變形后的位形上。因此，本文考慮的幾何非線性問題源于葉片的大位移、小應變問題。

風力機為旋轉(zhuǎn)機械，動力學方程建立在隨葉片轉(zhuǎn)動的葉根坐標系上[14]。計算葉片變形時，將其根部固定，建立懸臂梁的動力學方程；進行非線性求解時，則在每個方位角時，采用完全拉格朗日格式(T.L.)進行增量分析。

2.1結(jié)構(gòu)動力學方程

求解葉片動態(tài)響應，需要構(gòu)造葉片的動力學方程：

(6)

式中：M為質(zhì)量矩陣，C為阻尼矩陣，K為剛度矩陣，x為葉片變形量。

將風力機葉片簡化為歐拉梁模型[15]，采用有限元方法進行線性與非線性分析[16]，構(gòu)造動力學方程所需的質(zhì)量矩陣，剛度矩陣和阻尼矩陣。采用伽遼金方法計算單元矩陣，其在力學范圍內(nèi)體現(xiàn)為虛功原理。質(zhì)量矩陣在線性與非線性分析時，是沒有區(qū)別的。線性與非線性分析的主要區(qū)別在于有沒有考慮應變的二次項，即非線性項，其體現(xiàn)在剛度矩陣的求解過程中。阻尼矩陣表示為質(zhì)量矩陣與剛度矩陣的線性組合，即C=c1M+c2K，c1和c2為比例常數(shù)，可通過實驗確定。

忽略葉片扭轉(zhuǎn)，梁模型節(jié)點有五個自由度：展向位移u，揮舞方向位移v和轉(zhuǎn)角α，擺振方向位移w和轉(zhuǎn)角β。由此，單元任一點位移表示為[u,v,w]T=N(x)δ。其中，N(x)為形函數(shù)矩陣，δ=[ui,vi,αi,wi,βi,uj,vj,αj,wj,βj]T為單元節(jié)點位移列陣。

單元動能為：

(7)

2.1.1線性分析

線性分析時，單元應變能為：

(8)式中：u′)為u對x方向一次導數(shù)，表示x方向應變；v″,w″為v,w對x方向二次導數(shù)，表示曲率；E為楊氏彈性模量。

外力做功為：

W=Fδ

(9)

由哈密頓原理得：

(10)

式中：K0為線性單元剛度矩陣。

得到單元質(zhì)量矩陣和單元剛度矩陣之后，通過有限元方法，將各單元矩陣組合為一致質(zhì)量矩陣和一致剛度矩陣，并求得一致阻尼矩陣，建立動力學方程。然后，采用Newmark方法[16]進行時間步長推進。

2.1.2非線性分析

非線性分析時，單元應變能為：

(11)

由哈密頓原理得：

(12)

式中：KN為非線性剛度矩陣，其中包含了與當前位形相關(guān)的元素，K0+KN為全量剛度矩陣。

由于每次位移和轉(zhuǎn)角變化，都會導致非線性剛度矩陣變化，因此在時域分析時，宜采用增量計算的方法，需求解切線剛度矩陣：

(13)

式中：(K0+KT)為切線剛度矩陣。

時間步長推進時，在Newmark方法的基礎上采用Newton-Raphson方法，從t時刻向t+Δt時刻遞推公式如下：

(14)

式中變量可參考文獻[13,16]。由于方程左端的切線剛度矩陣是非線性的，所以每一個時間步內(nèi)都需要迭代求解。

2.2氣動結(jié)構(gòu)耦合計算實現(xiàn)

在一個時間步內(nèi)，首先應用自由渦尾跡方法進行氣動計算，然后將氣動載荷轉(zhuǎn)移到結(jié)構(gòu)模型中，再考慮葉片所受重力與離心力，作為葉片受到的激勵。然后計算葉片的剛度矩陣，質(zhì)量矩陣和阻尼矩陣，建立動力學方程，求解葉片響應，并反饋到氣動計算中，重新構(gòu)建氣動外形并進行下一個時間步的計算。

線性計算時，剛度矩陣在整個時間流程中是不變的。非線性計算時，每次迭代都會更新葉片節(jié)點位移與轉(zhuǎn)角，全量剛度矩陣和切線剛度矩陣也會不斷更新，直到計算結(jié)果達到收斂條件。由于每次迭代都需計算切線剛度矩陣，計算量較大，因此采用修正的Newton-Raphson方法，即每個時間步中，全量剛度矩陣每次迭代都進行更新，而切線剛度矩陣只在第一個迭代步更新。這樣雖然可能需要更多的迭代次數(shù)才能收斂，但是總的計算量減小了，而且此方法更穩(wěn)定，不易發(fā)散。非線性動態(tài)響應計算流程圖如圖2所示。

圖2　耦合計算流程圖Fig.2 Coupling calculation flow chart

3計算結(jié)果及分析

3.1計算方法驗證

通過美國可再生能源實驗室進行的NREL Phase VI風力機實驗數(shù)據(jù)[17]，驗證本文計算模型的可靠性與準確性。NREL Phase VI風力機為2葉片風力機，風輪直徑為10.058 m，輪轂中心高度為12.192 m，額定轉(zhuǎn)速為72 r/min，翼型采用S809翼型，額定功率為19.8 kW，屬于失速型功率控制。

圖3是各個穩(wěn)定風速工況下葉片上翼型的法向系數(shù)Cn和切向系數(shù)Ct，同時給出它們隨著展向方向的變化曲線圖。從圖3可以看出，葉片展向位置各個截面的計算值與實驗值具有較好的吻合度，特別是在低風速時?？拷~尖處95%截面，整個風速區(qū)域Cn和Ct與實驗值的吻合度比葉片中部和葉根要好，這是由于靠近葉尖處的誘導速度主要是葉尖渦在起作用，而且靠近葉尖處槳距角小，是整個葉片展向入流迎角最小的位置，失速程度較弱。但是，在高風速段，由于葉根處總是容易處于深失速的狀態(tài)，失速程度較嚴重，因此，相對于葉片其他部位和低風速，靠近葉根位置預測計算精確度較弱。從整體上看，本文計算模型具有一定可靠性和準確性。

3.2大型風力機葉片非線性動態(tài)響應計算

首先計算了大型風力機葉片NH1500不同葉尖速比的響應，并與實驗結(jié)果進行對比分析，如圖4所示。NH1500風力機是變槳變速風力機，具體的風力機性能和幾何參數(shù)見表1。

圖3　穩(wěn)態(tài)風速下葉片展向的Cn和CtFig.3 Variation of normal and tangential coefficient at different span with steady wind speed

本文計算結(jié)果與實驗值比較吻合，但比實驗值偏大。主要是因為本文氣動計算方法只在渦核處考慮了黏性，而沒有考慮全流場的黏性，而且實驗采用縮比模型，雷諾數(shù)偏小，導致功率系數(shù)偏小。

然后，分別采用線性分析與非線性分析計算了風速10 m/s時NH1500動態(tài)響應，并將葉尖揮舞、擺振位移以及風力機氣動性能進行對比分析。

表1　 NH1500風力機設計參數(shù)

圖4　NH1500計算驗證Fig.4 Validation of NH1500

對于葉尖揮舞位移，非線性與線性計算結(jié)果有明顯差別。通過圖5和表2可以看出，無論是平衡位置還是振幅，非線性計算結(jié)果都比線性計算結(jié)果要小。這是由于非線性分析過程中，會考慮到葉片由于揮舞和擺振位移而導致的展向長度減小和變形導致的彎曲剛度增強，其中后者為主要影響因素，因此葉尖的變形會比線性計算結(jié)果較小。

表2　 葉尖揮舞特性對比

圖5　葉尖揮舞位移對比Fig.5 Comparison of the blade tip wave displacement

對于葉片擺振位移，非線性與線性計算結(jié)果差別較小,見圖6。這是由于擺振方向葉片剛度較大，變形較小，而非線性剛度矩陣是與形變量相關(guān)的，在葉片擺振變形不明顯，即形變量很小的情況下，非線性全量剛度矩陣與線性剛度矩陣差別非常小。因此，對于剛度較大，變形較小的擺振方向，幾何非線性效應不明顯。

圖6　葉尖擺振位移對比Fig.6 Comparison of the blade tip shimmy displacement

從表3可以看出，對于風力機的氣動性能，非線性與線性計算結(jié)果有一定差別。風力機氣動性能主要與葉片控制點入流速度和展向位置有關(guān)。從圖7(a)、圖7(b)和表4可以看出，線性與非線性計算所得葉片揮舞速度有明顯差別，擺振速度差別較小。并且非線性分析會考慮葉片由于彎曲變形導致的展向位置變化，顯然葉片彎曲會導致葉片控制點距離葉根的展向長度減小，即氣動力力臂減小，從而導致功率的減小，而風輪平面的減小則會導致風輪推力減小。對于兆瓦級以上的大型風力機，氣動性能的稍許變化，就會對風能的捕獲產(chǎn)生重要影響。因此，能夠更準確的計算風力機氣動性能，對于提高風力機效率及相關(guān)的控制策略，有著重要作用。

表3　氣動性能對比

圖7　葉尖速度對比Fig.7 Comparison of the blade tip velocities

然后對線性與非線性分析得到的葉根剪力進行對比，如圖8，表5所示。揮舞方向葉根剪力的非線性計算結(jié)果明顯小于線性計算結(jié)果，擺振方向則差別較小。這是由于擺振方向葉根剪力的來源主要是重力載荷，而揮舞方向葉根剪力主要為氣動力，其受葉片變形影響較大，因此線性與非線性計算結(jié)果差別明顯。非線性分析能夠更準確的計算葉片變形及其對風力機氣動載荷產(chǎn)生的影響，從而更準確的計算各部件載荷。因此，非線性分析對于風力機載荷設計與結(jié)構(gòu)強度設計有重要意義。

表4　葉尖速度振幅對比

圖8　葉根剪力對比Fig.8 Comparison of the blade root shear forces

線性非線性差別平衡位置/N71000705360.65%振幅/N28172821-0.14%

4結(jié)論

本文計算了Phase VI葉片與NH1500葉片的穩(wěn)態(tài)性能，并與實驗結(jié)果進行對比，驗證了本文氣動計算方法的可靠性與準確性。然后計算了NH1500在線性與非線性分析下的動態(tài)響應，并進行對比分析，得出了以下結(jié)論：

(1) 非線性分析所得葉片揮舞位移明顯小于線性分析結(jié)果；非線性與線性分析所得葉片擺振位移差別較小，這是由于擺振方向剛度較大，葉片變形較小，幾何非線性效應不明顯造成的。

(2) 非線性與線性分析所得葉片揮舞速度差別明顯，并且非線性分析能更準確地計算控制點位置變化，而風力機氣動性能主要與葉片控制點位置及其入流速度有關(guān)，因此兩種分析方法得到的氣動性能差別明顯。

(3) 非線性分析能夠更準確的計算風力機部件載荷，對于風力機載荷設計與結(jié)構(gòu)強度設計有重要意義，能夠更好的保障風力機安全運行。

(4) 非線性分析能夠更好的考慮由于葉片變形帶來的幾何非線性效應，更準確地計算風力機動態(tài)響應，更適用于大型風力機模擬計算。

風力機大型化、柔性化發(fā)展，隨之葉片變形越來越大，線性分析勢必會造成資源不合理分配和計算準確度低，采用非線性分析能夠更準確地模擬葉片變形，對于更合理的分配資源，提高結(jié)構(gòu)強度與壽命，以及準確模擬風力機氣動性能，有著重要作用。

在考慮了幾何非線性后，葉片揮舞位移計算結(jié)果明顯減小，因此可以在保證安全系數(shù)的同時，進一步降低葉片的柔性或質(zhì)量，從而降低葉片生產(chǎn)成本；功率系數(shù)的改變則需要考慮新的功率控制方案來改善風力機性能。對于更大型的風力機，或者柔性更強的葉片，非線性與線性響應的差別則會更大。因此，在風力機大型化、柔性化的趨勢下，非線性分析對模擬風力機響應是必不可少的。

參 考 文 獻

[1] Yuan W B. Geometric nonlinear bending response and buckling of angle-section beams under pure bending [J]. Journal of Constructional Steel Research, 2013, 85:73-80.

[2] Almeida Carlos A, Albino Juan C R, Menezes Ivan F M, et al. Geometric nonlinear analysis of functionally graded beams using a tailored Lagrangian formulation [J].Mechanics Research Communications, 2011, 38(8):553-559.

[3] Sebastian T, Lackner M A. Development of a free vortex wake method code for offshore floating wind turbines [J]. Renewable Energy, 2012, 46: 269-275.

[4] Jeong M S, Yoo S J, Lee I. Wind turbine aerodynamics prediction using free-wake method in axial flow[C]//4th International Symposium on Physics of Fluids. Lijiang, China, 2012.

[5] Zhou W P, Tang S L, Lu H. Computation on aerodynamic performance of horizontal axis wind turbine based on time-marching free vortex method [J]. Proceedings of the CSEE, 2011,31(29):124-130.

[6] Zheng Y Q, Zhao R Z, Liu H. Finite element modal analysis of large-scale composite wind turbine blade [J]. Advanced Materials Research, 2013,694/695/696/697:453-457.

[7] Zhu S F, Rustamov I. Structural design and finite element analysis of composite wind turbine blade [J]. Key Engineering Materials, 2013,525/526: 225-228.

[8] Wu A, Sun W L, Wang P P. Large wind turbine dynamic modeling and simulation [J].Advanced Materials Research, 2010, 34/35:1098-1103.

[9] Riziotis V A, Politis E S, Voutsinas S G, et al. Stability analysis of pitch-regulated, variable-speed wind turbines in closed loop operation using a linear eigenvalue approach [J]. Wind Energy, 2008, 11(5):517-535.

[10] 莫文威，李德源，夏鴻建，等. 水平軸風力機柔性葉片多體動力學建模與動力特性分析[J]. 振動與沖擊,2013,32(22):99-105.

MO Wen-wei, LI De-yuan, XIA Hong-jian, et al. Multibody dynamic modeling and dynamic characteristics analysis of flexible blades for a horizontal axis wind turbine [J]. Journal of Vibration and Shock, 2013, 32(22):99-105.

[11] 馬劍龍，汪建文，董波，等. 風力機風輪低頻振動特性的實驗模態(tài)研究[J]. 振動與沖擊,2013,32(16):164-170.

MA Jian-long, WANG Jian-wen, DONG Bo, et al.Experimental modal analysis on low-frequency vibration characteristics of wind turbine [J]. Journal of Vibration and Shock, 2013, 32(16):164-170.

[12] 巴斯 K J. 工程分析中的有限元法[M]. 北京：機械工業(yè)出版社,1991.

[13] 王勖成. 有限單元法[M]. 北京：清華大學出版社,2003.

[14] Martin O L H. Aerodynamics of wind turbines [M]. Oxford:Taylor & Francis Group,2008.

[15] 程危危，曹登慶，初世明. 風力機葉片動力學建模研究[J]. 動力學與控制學報, 2011, 9(4):342-347.

CHENG Wei-wei, CAO Deng-qing, CHU Shi-ming. Dynamics modeling for the blade of a wind turbine [J]. Journal of Dynamics and Control,2011, 9(4):342-347.

[16] 彭細榮，楊慶生，孫卓. 有限單元法及其應用[M].北京:清華大學出版社,2012:369-373.

[17] Hand M M, Simms D A, Fingersh L J, et al.Unsteady aerodynamics experiment phases Ⅵ: Wind tunnel test configurations and available data campaigns[R]. Golden: National Renewable Energy Laboratory, NREL/TP-500-29955, 2001.

基金項目：國家973計劃項目(2014CB046200)大型風力機的關(guān)鍵力學問題研究及設計實現(xiàn)；江蘇省自然科學基金(BK20140059)；江蘇高校優(yōu)勢學科建設工程資助項目

收稿日期：2015-04-28修改稿收到日期：2015-08-06

通信作者王同光 男，博士，研究員，1962年生

中圖分類號:O355;TK89

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.14.030

Geometric nonlinear dynamics response of large-scale wind turbine

CAO Jiu-fa1, WANG Tong-guang2, WANG Xiao2

(1. School of Hydraulic, Energy and Power Engineering, Yangzhou University, Yangzhou 225127, China;2. Jiangsu Key Laboratory of Hi-tech Research for Wind Turbine Design, Nanjing University of Aeronautics & Astronautics, Nanjing 210016, China)

Abstract:With the development of large-scale flexible wind turbine, the deflection of the blades gets larger, which makes the geometric nonlinear effect more obvious. A free vortex wake model was presented to improve the accuracy of the aerodynamics load. The blade was assumed to be a beam model to analyze the structural response with the finite element method. Then a nonlinear kinetic equation was built, in which the geometric nonlinear stiffness matrix was derived by the Galerkin method. The equation was solved with the Newmark immediate integration method and the corrected Newton-Raphson method for time marching. Thus, the coupling simulation between the aerodynamics performance and the geometric nonlinear structural response of the large-scale wind turbine blade was realized. The computation results of the NREL Phase VI rotor agree well with the experimental data, validating the prediction accuracy of the calculation model. The linear and nonlinear dynamics responses of the large-scale wind turbine NH1500 were then calculated. The results of the geometric nonlinear analysis indicate the influence on the structural response and aerodynamics performance. The model is helpful to improve the design standard of the large-scale wind turbine and the accuracy of the load calculation.

Key words:wind turbine; free vortex wake method; geometric nonlinearity; Newton-Raphson method; dynamics response

第一作者 曹九發(fā) 男，博士生，1986年生