風電葉片低階頻率識別的近似解法及試驗

黃雪梅，張磊安

(山東理工大學 機械工程學院，山東 淄博 255049)

風電葉片是風力發(fā)電機的關(guān)鍵核心部件之一，隨著發(fā)電機功率的不斷增大，要求葉片長度不斷增加，對葉片動態(tài)響應的研究也越來越重要。模態(tài)分析作為研究葉片動態(tài)響應最常用的方法，主要是為確定葉片的第1階、第2階頻率，可分為模態(tài)計算方法與模態(tài)試驗方法［1－3］。模態(tài)計算方面，由于風電葉片形狀為不規(guī)則的多自由度構(gòu)件，對葉片結(jié)構(gòu)的處理也比較粗糙，所以其模態(tài)參數(shù)解不易準確求得。模態(tài)試驗由于該方法比較成熟、可靠，在工程上已經(jīng)得到了廣泛的應用。由于MW級風電葉片具有尺寸大、價格高等特點，因而現(xiàn)有的試驗模態(tài)方法大多局限于長度不超過1 m的小葉片上［4］。吳春梅等［5］采用錘擊法和正弦激勵法對長度為0.6 m，額定功率為300 W的葉片進行了試驗模態(tài)分析。劉旺玉等［6］對一個翼型為SG6050，半徑為1 m的小型風力機葉片采用仿生學原理進行了模態(tài)試驗。歐陽濤等［7－8］通過多傳感器辨識和葉尖定時方法對葉片振動進行測試。

基于上述原因，本文將傳遞矩陣法引入到葉片的模態(tài)參數(shù)識別計算中。將葉片整體結(jié)構(gòu)的力學分析問題轉(zhuǎn)化為若干單元對接與傳遞的力學分析問題，利用矩陣相乘的計算方法實現(xiàn)對葉片的模態(tài)分析。引入邊界條件后，確定該葉片的低階頻率計算值。最后通過構(gòu)建葉片試驗模態(tài)試驗系統(tǒng)，對構(gòu)件施加正弦激勵，使之產(chǎn)生自由振動，通過安裝在某些確定位置的加速度記錄相應振動信號，再進行快速傅里葉變換，得到試驗值，并將兩者結(jié)果進行比較。

1 低頻識別算法

1.1 Bernoulli-Euler Beam模型

風電葉片根部通過螺栓法蘭盤被固定在筒型加載支座上，通常將其視為一個沿展向變質(zhì)量的懸臂梁，主要變形為自由彎曲。目前，識別自由懸臂梁模態(tài)參數(shù)的主要方法有解析法、有限元法［9］、直接動剛度法［10－11］等。

建立風電葉片歐拉－伯努利(Bernoulli-Euler Beam)模型［12］:

式中:y為相對于根部固定端的位移;E為葉片材料彈性模量;J為葉片截面慣性矩;ρ為葉片密度;A為橫截面面積;P(x，t)為單位長度受到外力;m(x，t)為單位長度受到外力矩。

風電葉片做自由振動時，有:

此時:

可知葉片的自由彎曲振動為四階微分方程，由于葉片并不是真實等截面，式(3)的解析解不易準確求得，本文將傳遞矩陣算法應用于風電葉片的模態(tài)參數(shù)識別。

1.2 矩陣傳遞算法

將葉片等效成一個根部固定的變截面懸臂梁，將其分割成若干段，利用二節(jié)點梁做離散化建模，離散后模型如圖1所示。圖中，mi為該段的質(zhì)量，li為該段長度，EiIi為該段的抗彎剛度。

圖1 二節(jié)點梁模型Fig.1 Two-node beam model

葉片做彎曲振動時，離散質(zhì)量塊mi除了有撓度yi和轉(zhuǎn)角 θi，還有與之相對應的剪力Qi和彎矩Mi。這4個變量構(gòu)成狀態(tài)矢量:

從葉片中截取第i個離散質(zhì)量塊，每個質(zhì)量塊受力分別為彎矩Mi和剪力Qi。建立其受力模型，如圖2所示。

圖2 二節(jié)點梁受力模型圖Fig.2 Two-node beam force model picture

對葉片第i個離散質(zhì)量塊mi進行受力分析，左邊梁li對其作用力表示為彎矩MLi和剪力QLi，右段梁 li對其作用力表示為彎矩MRi和剪力QRi。當質(zhì)量塊mi做簡諧運動的位移為yi時，產(chǎn)生的慣性力可表示為miω2yi。

由于質(zhì)量塊兩邊的位移yi和轉(zhuǎn)角θi相等，寫成點矩陣形式，表示為:

為了建立質(zhì)量塊mi兩邊的撓度 yi和轉(zhuǎn)角 θi的關(guān)系，引入材料力學關(guān)系式［13］:

聯(lián)立以上各式，建立葉片第i個質(zhì)量塊mi的狀態(tài)矢量與第i－1個質(zhì)量塊mi－1的狀態(tài)矢量之間關(guān)系，表示為:

根據(jù)式(7)，以 aeroblade3.6－56.4風電葉片為研究對象，沿葉片展向離散為25個元件，26個聯(lián)結(jié)點。定義葉片根部為輸入端，尖部為輸出端，假設(shè)Pij表示第i個元件和第j個元件的聯(lián)結(jié)點，每個聯(lián)結(jié)點Pij對應一個狀態(tài)矢量Zij。根據(jù)矩陣傳遞特性，每個元件的傳遞方程為:

式(8)進一步表示為:

將 H25，H24…H0代入式(9)，得到總傳遞矩陣方程:

根據(jù)實際情況，邊界條件分別取:

葉片根部:

葉片尖部:

根據(jù)懸臂梁特性，將式(10)的第一行與第二行展開，得到風電葉片低階頻率識別方程，表示為:

2 低頻識別計算

風電葉片由截面翼型組成，葉片根部截面為圓形且較厚，中間為氣動翼型，往葉尖方向逐漸變薄，主要為了滿足氣動性能要求。各截面翼型材料主要分為:① 葉片根部為玻璃纖維/環(huán)氧樹脂;② 葉片中部為(玻璃纖維/環(huán)氧樹脂)/輕木/(環(huán)氧樹脂/玻璃纖維);③葉片尖部為(玻璃纖維/環(huán)氧樹脂)/PVC塑料/(環(huán)氧樹脂/玻璃纖維)。質(zhì)量作為風電葉片的重要特征，為了保證葉片盡可能輕，內(nèi)部通常采用中空構(gòu)型。風電葉片的4個加載面分別定義為，最大面向(Maximum flapwise)、最小面向(Minimum flapwise)、最大弦向(Maximun edgewise)和最小弦向(Minimum edgewise)。按照結(jié)構(gòu)分為:筋板－承載結(jié)構(gòu)、殼體－氣動結(jié)構(gòu)和葉根－連接結(jié)構(gòu)，總體結(jié)構(gòu)如下。

圖3 風電葉片截面示意圖Fig.3 Wind turbine blade cross－ section schematic

本文研究對象為aeroblade3.6－56.4風電葉片，其額定功率為3.6 MW，長度為56.4 m，質(zhì)量為17 277 kg，重心為距離葉片根部19.05 m處。該款葉片各截面翼型主要設(shè)計參數(shù)(預彎、扭角、抗彎剛度和彈性模量等)如表1所示。

表1 葉片截面翼型設(shè)計參數(shù)Tab.1 The design parameters of each section airfoil in the blade model

根據(jù)表1中葉片各截面的設(shè)計參數(shù)，編制沿aeroblade3.6－56.4風電葉片展向的矩陣傳遞程序，計算出Max flapwise方向的低階固有頻率曲線，如圖4所示。

圖4 Aeroblade3.6－56.4風電葉片計算頻率曲線Fig.4 Aeroblade3.6－56.4 wind turbine blade calculate frequency curve

從圖中得出，aeroblade3.6－56.4風電葉片 Max flapwise方向的低階固有頻率計算結(jié)果分別為:一階頻率:f1st=0.68 Hz;二階頻率:f2st=1.76 Hz;三階頻率:f3st=6.46 Hz。

3 低頻識別試驗

為了與上述計算結(jié)果進行比較，構(gòu)建一套葉片模態(tài)參數(shù)識別試驗系統(tǒng)。

3.1 FIR濾波器設(shè)計

模態(tài)試驗系統(tǒng)采集的原始數(shù)據(jù)中可能混有噪聲，進行采樣之前，需要對連續(xù)信號用低通過濾器進行濾波。基于運算精度與穩(wěn)定性方面考慮，本文設(shè)計一個有限脈沖濾波器(FIR)。FIR濾波器的基本結(jié)構(gòu)為一個分節(jié)的延時線，將每一節(jié)的輸出值加權(quán)累加，得到濾波器的輸出值。數(shù)學表示為［14－15］:

式中:N為FIR濾波器的抽頭數(shù);x(i)為第i時刻的輸入樣本;h(i)為FIR濾波器的第i級抽頭系數(shù)。

其傳遞函數(shù)為:

為了將傅里葉反變換導出的無限長序列變?yōu)殚L度為N的序列，采用矩形窗函數(shù)法進行截取。矩形窗函數(shù)表達式為:

3.2 模態(tài)試驗平臺構(gòu)建

被試驗葉片目前屬于國內(nèi)最大的風力機發(fā)電葉片，為了對該葉片的低階固有頻率進行試驗識別，構(gòu)建了一個模態(tài)參數(shù)識別試驗平臺。筒型加載支座固定于地面，葉片根部通過88個高強度螺桿與加載支座相連，葉片處于自由懸臂狀態(tài)，如圖5所示。為了消除了外界因素(譬如，風載荷)影響，該試驗平臺布置在室內(nèi)，系統(tǒng)的主要硬件包括加速度計、電荷放大器、數(shù)據(jù)采集儀和計算機等，主要技術(shù)指標如表2。

圖5 模態(tài)參數(shù)識別平臺Fig.5 Modal parameter identify test platform

表2 硬件設(shè)備技術(shù)指標Tab.2 The hardware device technical property

3.3 模態(tài)試驗過程

將aeroblade3.6－56.4風電葉片固定在筒型加載支座上，加速度傳感器置于葉片尖部，試驗現(xiàn)場如圖6所示。采用單點激勵中的“不測力法”，對葉片施加一個正弦激勵，引起葉片自由振動，電荷放大器將加速度信號放大并傳送給信號采集裝置，如圖7所示。

3.4 試驗結(jié)果及分析

限于篇幅，僅給出 aeroblade3.6－56.4風電葉片Max flapwise與Max edgewise方向的低階固有頻率試驗結(jié)果。將數(shù)據(jù)采集儀得到的加速度信號進行FIR濾波及快速FFT變換，頻譜圖如圖8所示。

根據(jù)圖8得出，aeroblade3.6－56.4風電葉片的低階固有頻率分別為，Max flapwise1st=0.64 Hz，Max flapwise 2nd=1.68 Hz;Max edgewise1st=1.12 Hz，Max edgewise 2nd=3.64 Hz。將揮舞方向的計算結(jié)果與試驗結(jié)果進行比較，兩者誤差不超過6%。

圖8 aeroblade3.6－56.4風電葉片試驗頻率曲線Fig.8 Aeroblade3.6－56.4wind turbine blade test frequency curve

4 結(jié)論

本文將模態(tài)計算方法與試驗方法相結(jié)合對風電葉片的低階固有頻率進行識別。通過將葉片沿展向離散化，建立葉片的離散化模型。編制了矩陣傳遞算法，將其應用于aeroblade3.6－56.4風電葉片的低階固有頻率識別計算。構(gòu)建了一套大葉片模態(tài)參數(shù)識別試驗平臺，將模態(tài)試驗結(jié)果與計算結(jié)果進行比較，兩者數(shù)值基本一致，說明采用矩陣傳遞算法能較好地計算出葉片的低階固有頻率值，也驗證了該算法的有效性。將葉片兩個方向的試驗結(jié)果對比，得出擺振方向的同階固有頻率遠高于揮舞方向。另外，本文通過模態(tài)試驗方法求出的葉片低階固有頻率，也為后續(xù)的風電葉片疲勞加載系統(tǒng)研發(fā)及試驗奠定了基礎(chǔ)。

［1］ 毛火軍，石可重，李宏利，等.大型風電葉片的模態(tài)測試與數(shù)值模擬［J］.工程熱物理學報，2009，30(4):601－604.MAO Huo-jun，SHI Ke-chong，LI Hong-li，et al.Modal testing and numerical simulation of large wind turbine blade［J］.Journal of Engineering Thermophysics，2009，30(4):601－604.

［2］ Cha P D，de Pillis L G.Model updating by adding known masses［J］.International Journal for Numerical Method in Engineering，2001，50:2547－3571.

［3］ Ryu J，Kim Sang，Kim Sung Soo.Ageneral approach to stress stiffening effents on flexible multibody dynamic systems［J］.Mech Struct＆Mach，1994，22(2):157－180.

［4］ 葉枝全，馬昊旻，丁 康，等.水平軸風力機槳葉的實驗模態(tài)分析［J］.太陽能學報，2001，22(4):473－476.YE Zhi-quan，MA Hao-min，DING Kang，et al.Experimental model ananlysis of the rotor blade of the horizontal axis wind turbine［J］.Acta Energiae Solaris Sinica，2001，22(4):473－476.

［5］ 吳春梅，田 瑞，劉 博，等.小型水平軸風力機葉片的振動研究［J］.能源技術(shù)，2006，27(5):205－207.WU Chun-mei，TIAN Rui，LIU Bo，et al.Study on small size blade vibration characteristic of the horizontal axis wind turbine［J］.Energy Technology，2006，27(5):205－207.

［6］ 劉旺玉，張海全，李明剛.基于仿生設(shè)計的風力發(fā)電機葉片力學性能的實驗研究［J］.實驗力學，2009，24(2):121－125.LIU Wang-yu，ZHANG Hai-qian，LI Ming-gang.Experimental study of mechanical performance of wind turbine blade based on bionic design［J］.Journal of Experimental Mechanics，2009，24(2):121－125.

［7］ 李孟麟，段發(fā)階，歐陽濤，等.基于葉尖定時的旋轉(zhuǎn)機械葉片振動頻率識別 ESPRIT方法［J］.振動與沖擊，2010，29(12):18－21.LI Meng-lin，DUAN Fa-jie，OUYANG Tao，et al.ESPRIT Method for blade vibration frequency identification in rotating machinery based on blade tip-timing measurement［J］.Journal of Vibration and Shock，2010，29(12):18－21.

［8］ 歐陽濤，郭文力，段發(fā)階，等.基于葉尖定時的旋轉(zhuǎn)葉片同步振動辨識新方法［J］.振動與沖擊，2011，30(8):247－252.OUYANG Tao，GUO Wen-li，DUAN Fa-jie，et al.New method for identifying rotating blades synchronous vibration based on tip-timing［J］.Journal of Vibration and Shock，2011，30(8):247－252.

［9］ Martinerz C A E，Museros P，Castillo L A.Semranalytic solution in the time domain for non－uniform multispan Bernoulli－ Euler beams traversed by moving loads［J］.Sound Vib，2006，294(1):277－278.

［10］ Henchi K，F(xiàn)afard M.Dynamic behaviour of mulitispan beams under moving loads［J］.Sound Vib，1997，199(1):33－35.

［11］ Dugush Y A，Eisenberger M.Vibrations of nonuniform continuous beams under moving loads［J］.Sound Vib，2002，254(5):9－11.

［12］ 王光遠.建筑結(jié)構(gòu)的振動［M］.北京:科學出版社，1978.

［13］ 賀興書.機械振動學［M］.上海:上海交通大學出版社，1985.

［14］ 彭紅平，楊福寶.基于Matlab的FIR數(shù)字濾波器設(shè)計［J］.武漢理工大學學報(信息與管理工程版)，2005，27(5):275－278.PENG Hong-ping，YANG Fu-bao.Design of FIR digital filter based on matlab［J］.Journal of WUT(Information ＆ Managementent Engin-eering)，2005，27(5):275－278.

［15］ 李 輝，張 安，趙 敏，等.粒子群優(yōu)化算法在FIR數(shù)字濾波器設(shè)計中的應用［J］.電子學報，2005，33(7):1338－1341.LI Hui，ZHANG An，ZHAO Min，et al.Particle swarm optimization algorithm for FIR digital filters design［J］.Acta Electronica Sinca，2005，33(7):1338－1341.