基于緊致差分格式的風(fēng)電葉片抗彎剛度分布辨識

摘 要：針對葉片抗彎剛度分布難以辨識的實際問題，提出一種基于單點靜載撓度擬合和緊致差分格式的葉片抗彎剛度辨識方法，建立葉片靜態(tài)標定工況彎曲變形數(shù)學(xué)模型，進一步推導(dǎo)出各截面撓度和抗彎剛度表達式。通過對多支不同型號葉片的分析結(jié)果表明，葉中部分抗彎剛度辨識誤差均小于5%，驗證了該方法的有效性。

關(guān)鍵詞：風(fēng)電葉片；抗彎剛度；曲線擬合；標定；緊致差分格式

中圖分類號：TM315" " " " " " " " " " " " 文獻標志碼：A

DOI：10.19912/j.0254-0096.tynxb.2023-0545

文章編號：0254-0096（2024）08-0523-06

同濟大學(xué)機械與能源工程學(xué)院，上海 201804

0 引 言

風(fēng)電葉片抗彎剛度是指葉片在受到外力作用時，抵抗彎曲變形的能力。風(fēng)力機運行時葉片承受風(fēng)力、重力和慣性力等作用，需要具備足夠的抗彎剛度以保證不發(fā)生過度變形、分離、損壞和掃塔等故障。葉片采用變截面、變剛度設(shè)計，通過調(diào)整葉片截面形狀、尺寸、材料來優(yōu)化不同位置的抗彎剛度以減輕質(zhì)量、提高性能。在葉片設(shè)計和實際制造過程中，葉片材料、截面形狀和鋪層結(jié)構(gòu)的復(fù)雜性、制造過程中可能存在的誤差等使得葉片實際剛度與設(shè)計剛度值不一致，進而對葉片性能和使用壽命造成不利影響。

為了在設(shè)計開發(fā)階段評估葉片模型的可信度，國際標準IEC 61400-5：2020［1］提出通過與全尺寸葉片測試結(jié)果比較來驗證全局模型，其模型驗證要求在IEC 61400-23：2014［2］規(guī)定的試驗荷載水平下，允許葉片最靠近葉尖加載位置的整體彎曲撓度偏差小于等于±7%。基于該誤差要求，文獻［3-5］采用有限元法和有限差分法計算葉片設(shè)計模型撓度，與實際葉片測試數(shù)據(jù)進行對比。風(fēng)電葉片具有變截面、變剛度特性，盡管可通過靜態(tài)測試確定被測葉片的撓度間接驗證其綜合抗彎剛度與設(shè)計模型的一致性，但難以確定葉片各截面抗彎剛度的分布。

對于實際葉片的抗彎剛度分布辨識，IEC 61400-23：2001［6］提出剛度試驗思路，在葉尖施加載荷，測量葉片位移或應(yīng)變，根據(jù)各截面的彎矩及其形成的曲率計算得到詳細的剛度分布，其中位移法更簡單快捷。文獻［7］推導(dǎo)了葉尖施加集中載荷工況下截面剛度的計算式。文獻［8-9］應(yīng)用位移法辨識了65 kW葉片兩個截面的抗彎剛度，監(jiān)測在役葉片剛度退化情況。文獻［10］根據(jù)鏡面法［11］和三點彎曲技術(shù)［12］法推導(dǎo)了葉片單點載荷彎曲剛度表達式，采用千分表測量撓度，用于辨識1.873 m葉片抗彎剛度分布，葉片均勻區(qū)域（葉中至葉尖）揮舞剛度誤差為2.37%，擺振剛度誤差為14.3%；針對該葉片，文獻［13］使用數(shù)字圖像相關(guān)技術(shù)測量變形，葉片均勻區(qū)域的剛度辨識誤差小于7%。葉片抗彎剛度分布還可通過計算每個截面的慣性矩和材料的彈性模量得到，文獻［14-15］將8.325 m葉片切片計算各截面慣性矩，基于靜力和模態(tài)測試數(shù)據(jù)校準葉片模型梁單元楊氏模量，進而得到抗彎剛度分布。

基于此，本文提出一種基于緊致差分格式的大型風(fēng)電葉片抗彎剛度辨識方法。根據(jù)靜態(tài)標定撓度數(shù)據(jù)進行高階多項式擬合，采用緊致差分格式求解葉片彎曲變形方程，從而得到葉片抗彎剛度分布的辨識值?？箯潉偠缺孀R值可用于評估葉片結(jié)構(gòu)的實際剛度性能、校準葉片設(shè)計模型，并可為葉片后續(xù)疲勞試驗的載荷計算提供實際剛度參數(shù)。同時，抗彎剛度辨識值還可應(yīng)用于葉片的運維和損傷監(jiān)測等需要剛度指標的場景。

1 葉片抗彎剛度分布辨識原理

1.1 撓度測量數(shù)據(jù)擬合

在葉片上施加單點靜態(tài)載荷進行葉片靜態(tài)標定［16-17］，通過測量葉片各截面撓度對被測葉片的抗彎剛度分布進行辨識。測試系統(tǒng)由試驗臺、被測葉片、夾具、加載裝置、測量儀器等組成，如圖1所示，采用試驗臺坐標系，其中原點位于葉片軸線根部。

針對可能產(chǎn)生的測量誤差，需對標定中測量的撓度值進行擬合處理，以獲得更準確數(shù)值。根據(jù)葉片懸臂梁理論模型，撓度測量值擬合函數(shù)需滿足：

1）邊界條件，葉根處撓度和轉(zhuǎn)角為零；

2）加載點至葉尖部分無彎曲變形，曲率為零；

3）葉根至加載點部分有彎曲變形，截面的轉(zhuǎn)角、曲率不為零；

4）數(shù)據(jù)擬合平滑處理，無需過所有測量點。

1.2 四階緊致差分格式彎曲變形方程

將被測葉片離散化，葉根截面編號記為0，加載力截面為[n]，截面間距為[h]，如圖2所示，由材料力學(xué)撓曲軸近似微分方程得：

[d2xidz2=±MiEIi] （1）

式中：[xi]——第[i]（[0≤i≤n]）截面撓度，m；[z]——沿葉片展向方向，m；[Mi]——截面彎矩，N·m；[E]——楊氏模量，Pa；[Ii]——截面慣性矩，[m4]。

在單點載荷作用下，葉片將產(chǎn)生相應(yīng)的形變，可通過施加的載荷和力臂計算得到葉片各截面的彎矩分布，其表達式為：

[Mi=F×Li] （2）

式中：[F]——與葉片軸線垂直的載荷，N；[Li]——第[i]截面的力臂，m。

有限差分法的主要原理是對微分方程中的微分項進行直接差分近似，從而將微分方程轉(zhuǎn)化為代數(shù)方程組求解［18］。文獻［19］發(fā)現(xiàn)同樣的節(jié)點數(shù)量下四階緊致差分比中心差分精度高，將其應(yīng)用至葉片靜載彎曲變形：

[xi-1+10xi+xxi+1=12h2xi-1-2xi+xi+1] （3）

式中：[h]——各截面間距，m。

根據(jù)上述擬合函數(shù)葉根邊界條件，即[x0=0]，[x′0=0]。由式（3）可得葉片第[i]截面（[1≤i≤n-1]）的四階緊致差分彎曲變形方程組，寫為矩陣形式：

[Ax=BX+C] （4）

式中，[A=1010…01101??0???0??11010…0110]，[x=x″1x″2?x″n-1]，[B=12h2-210…01-21??0???0??1-210…01-2]，[C=00?12h2xn-x″n]，[X=x1x2?xn-1]。

式（3）不包含葉根第0截面，為確定葉根抗彎剛度[EI0]，在懸臂梁截面0向左延伸一個“虛擬”截面，編號記為-1，由一階中心差分格式，葉根第0截面的轉(zhuǎn)角為：

[x′0=x1-x-12h=0] （5）

對第0截面采用二階中心差分格式，其方程為：

[x0=1h2x-1-2x0+x1] （6）

聯(lián)立式（4）～式（6）求解線性方程組，可得葉片抗彎剛度分布[Si=S0，S1，S2，…，Sn-1T]為：

[Si=M0h22x1" " " " " " " " " " " ，" i=0Mih212xi-1-2xi+xi+1" " " ，" 1≤i≤n-1] （7）

2 仿真驗證與分析

2.1 撓度驗證

葉片靜態(tài)標定彎曲變形數(shù)學(xué)模型中，撓度計算是抗彎剛度辨識的逆向問題，由四階緊致差分格式彎曲變形方程組式（4），可得葉片撓度[xk=x0，x1，x2，…，xnT]為：

[xk=0，" k=0k2D0+i=1k-1（k-i）Ei+1" " " ，1≤k≤n] （8）

式中：[Ei=112Di-2+10Di-1+Di]； [Di=h2MiEIi]。

以某型73 m葉片靜態(tài)標定工況為例，取不同截面間距1.5、3.0、6.0 m，應(yīng)用式（8）計算葉尖撓度，與有限元仿真結(jié)果進行對比，以驗證葉片靜態(tài)標定彎曲變形數(shù)學(xué)模型的準確性。如表1所示，基于四階緊致差分格式的計算結(jié)果與有限元計算相比誤差小于2%，且截面間距越小，計算精度越高。

2.2 高階多項式擬合

以4支不同型號的葉片為例，葉片參數(shù)及加載信息如表2所示，理論撓度曲線如圖3所示。考慮工程實際測量儀器的精度，將靜載工況仿真得到的理論撓度取至毫米級，用于后續(xù)剛度辨識仿真計算。對4號葉片的撓度進行擬合，以驗證葉片抗彎剛度辨識方法中高階多項式擬合撓度的必要性。

為了得到最優(yōu)撓度擬合，提出撓度擬合均方根誤差指標[Ed]，其表達式為：

[Ed=1ni=1nxi-xi2] （9）

式中：[n]——截面數(shù)量；[xi]——第[i]截面撓度擬合值，m。

為了評估剛度辨識誤差，提出均方根誤差[Ee]指標，其表達式為：

[Ee=1ni=1nSi-Si2] （10）

式中：[ Si]——第[i]截面辨識剛度，N/m；[Si]——第[i]截面理論剛度，N/m。

4號葉片撓度測量值未經(jīng)擬合及三階擬合的剛度辨識誤差[Ee]較大；高階多項式撓度擬合誤差[Ed]及剛度辨識誤差[Ee]較小，其中四階擬合的誤差最小，如表3所示。高階多項式函數(shù)的優(yōu)點是更易擬合，當撓度擬合最優(yōu)時，剛度辨識誤差最小。因此，對撓度測量值進行高階多項式擬合可有效提高葉片抗彎剛度辨識精度。

2.3 抗彎剛度辨識

為了驗證前述方法的通用性，對表2中4支葉片進行抗彎剛度辨識。在葉根部分［0，0.1］區(qū)間，其理論抗彎剛度數(shù)值和變化較大，而撓度測量值為0.000 m，導(dǎo)致剛度辨識效果不佳，如圖4所示；葉根部分［0.1，0.3］區(qū)間存在葉片弦長最大區(qū)域，大部分截面誤差在10%以內(nèi)；在葉中部分［0.3，0.7］區(qū)間剛度辨識效果較好，誤差在5%以內(nèi)；在葉尖部分［0.8，0.95］區(qū)間，葉片理論抗彎剛度數(shù)值較小、變化較大，大部分截面誤差在10%以內(nèi)；因無法辨識加載點處的剛度值，導(dǎo)致［0.95，1］區(qū)間內(nèi)的剛度辨識效果不佳。

總體而言，當葉片靜態(tài)標定撓度測量值精確到毫米時，基于撓度高階多項式擬合和四階緊致差分格式的剛度辨識方法，在葉片長度的［0.1，0.9］區(qū)間抗彎剛度辨識誤差為10%，在葉片長度的［0.3，0.7］區(qū)間抗彎剛度辨識誤差為5%。

3 實測數(shù)據(jù)分析

以某70 m級大型風(fēng)電葉片為例，其撓度測試歸一化結(jié)果如表4所示。測試值與理論值的最大偏差為4.5%，測試葉片的變形與理論偏差滿足規(guī)范±7%要求［1］。為檢驗軟件工具包輸出參數(shù)和施加載荷的正確性，需進一步辨識被測葉片詳細的抗彎剛度分布。被測葉片靜態(tài)標定現(xiàn)場如圖5所示，加載點為59 m。

被測葉片葉中部分撓度實測數(shù)據(jù)歸一化結(jié)果如表5所示。對被測葉片4個工況的撓度測量值進行擬合，最優(yōu)擬合曲線如圖6所示。

剛度辨識值曲線較為光滑連續(xù)，而理論值曲線在23 m和35 m處有明顯突變（圖7a），這是為了后續(xù)的靜力和疲勞測試，葉片在該位置進行了局部結(jié)構(gòu)補強，其剛度值會顯著增加?？拷~根部位弦長最大的截面形狀的突變也會造成剛度值的突變［20］。針對剛度突變截面，可在剛度突變處增加測量點，并采取分段擬合的方法，提高突變處剛度辨識的精度。葉片揮舞方向23～53 m、擺振方向25～49 m的剛度辨識值和理論值誤差在10%以內(nèi)，如圖7b所示，說明理論葉片和實際葉片抗彎剛度在葉中部分具有一致性?？筛鶕?jù)辨識的抗彎剛度分布，進一步修正設(shè)計模型參數(shù)，獲得更精準的葉片模型和測試數(shù)據(jù)。

4 結(jié) 論

本文研究了全尺寸葉片揮舞和擺振方向抗彎剛度辨識和測試方法，建立了基于四階緊致差分格式的葉片靜載彎曲變形數(shù)學(xué)模型，提出撓度測量數(shù)據(jù)采用多項式函數(shù)擬合方法，推導(dǎo)理論撓度和抗彎剛度分布計算表達式，通過仿真與試驗分析，得出如下主要結(jié)論：

1）通過在靜態(tài)標定過程中增加撓度測量截面數(shù)量，可達到辨識實際葉片抗彎剛度分布目的。

2）單點靜載工況下，基于四階緊致差分格式的葉片撓度計算速度快，其結(jié)果與有限元模型仿真誤差小于2%，截面間距越小，撓度計算越精確。

3）對于精確到毫米的撓度測量值，采用高階多項式擬合撓度，能有效提高葉片抗彎剛度辨識精度。

4）通過對4支不同型號葉片計算與分析，葉中部分（［0.3，0.7］）剛度誤差小于5%，驗證了上述抗彎剛度辨識方法的有效性。對于葉根和補強等剛度突變位置的精確辨識，還需進一步探究。

［參考文獻］

［1］ IEC 61400-5：2020， Wind energy generation systems-part 5： wind turbine blades［S］.

［2］ IEC 61400-23：2014， Wind turbines-part 23：Full-scale structural testing of rotor blades［S］.

［3］ 張顏明， 李國勇， 王艷， 等. 全尺寸葉片變形一致性測試技術(shù)研究［J］. 風(fēng)能， 2022（1）： 68-71.

ZHANG Y M， LI G Y， WANG Y， et al. Research on deformation consistency testing technology of full-scale blades［J］. Wind energy， 2022（1）： 68-71.

［4］ 靳交通， 劉力根， 鄧航， 等. 基于有限差分法和有限元法的大型風(fēng)電葉片剛度計算［J］. 機械設(shè)計與研究， 2012， 28（6）： 83-85.

JIN J T， LIU L G， DENG H， et al. Stiffness calculation of large-scale wind turbine blade based on the finite difference method and the finite element method［J］. Machine design amp; research， 2012， 28（6）： 83-85.

［5］ 靳交通， 梁鵬程， 曾竟成， 等. 復(fù)合材料風(fēng)電葉片有限元剛度分析［J］. 武漢理工大學(xué)學(xué)報， 2009， 31（21）： 133-136.

JIN J T， LIANG P C， ZENG J C， et al. Finite element analysis for composite material wind rotor blade stiffness［J］. Journal of Wuhan University of Technology， 2009， 31（21）： 133-136.

［6］ IEC 61400-23： 2001， wind turbines-part 23： full-scale structural testing of rotor blades［S］.

［7］ 傅程， 張方慧. 風(fēng)輪葉片疲勞試驗加載設(shè)計［J］. 中國機械工程， 2018， 29（6）： 743-747.

FU C， ZHANG F H. Loading design of wind turbine blade fatigue tests［J］. China mechanical engineering， 2018， 29（6）： 743-747.

［8］ 白學(xué)宗， 安宗文， 侯運豐， 等. 低速沖擊載荷下的風(fēng)電葉片剛度退化規(guī)律［J］. 太陽能學(xué)報， 2022， 43（1）： 132-139.

BAI X Z， AN Z W， HOU Y F， et al. Stiffness degradation rule of wind turbine blade under low-velocity impact load［J］. Acta energiae solaris sinica， 2022， 43（1）： 132-139.

［9］ 白學(xué)宗， 安宗文， 侯運豐， 等. 風(fēng)電葉片剛度檢測值和退化速率的溫度效應(yīng)［J］. 太陽能學(xué)報， 2021， 42（12）： 245-252.

BAI X Z， AN Z W， HOU Y F， et al. Temperature effect of stiffness detection value and degradation rate of wind turbine blades［J］. Acta energiae solaris sinica， 2021， 42（12）： 245-252.

［10］ JUNG S N， LAU B H. Determination of HART Ⅰ blade structural properties by laboratory testing［R］. NASA/CR-2012-216039， 2012.

［11］ CHANDRA R， STEMPLE A D， CHOPRA I. Thin-walled composite beams under bending， torsional， and extensional loads［J］. Journal of aircraft， 1990， 27（7）： 619-626.

［12］ OGORKIEWICZ R M， MUCCI P E R. Testing of fibre-plastics" " " composites" " " in" " " three-point" " "bending［J］. Composites， 1971， 2（3）： 139-145.

［13］ PFLUMM T， GAUGELHOFER L， HEUSCHNEIDER V， et al. Hybrid experimental measurement of sectional stiffness properties of the merit rotor blade with digital image correlation［C］// 47th European Rotorcraft Forum， United Kingdom，7-9th September， 2021.

［14］ GRIFFITH D T， PAQUETTE J A， CARNE T G. Development of validated blade structural models［C］//46th AIAA Aerospace Sciences Meeting and Exhibit. Reno， Nevada， USA， 2008.

［15］ GRIFFITH D T， CARNE T G. Experimental modal analysis of 9-meter research-sized wind turbine blades［M］. New： Springer New York， 2011： 1-14.

［16］ IEC 61400-13： 2015， Wind turbines- Part 13： Measurement of mechanical loads［S］.

［17］ 烏建中. 風(fēng)電葉片全尺寸結(jié)構(gòu)測試理論基礎(chǔ)、系統(tǒng)設(shè)計和試驗方法［M］. 上海： 同濟大學(xué)出版社， 2021.

WU J Z. Full-scale structural testing theoretical basis， system design and experimental method of wind turbine blades［M］. Shanghai： Tongji University Press， 2021.

［18］ 王漢權(quán)， 成蓉華. 微分方程數(shù)值方法： 有限差分法［M］. 北京： 科學(xué)出版社， 2020.

WANG H Q， CHENG R H. Numerical method of differential" "equation：" "finite" "difference" "method［M］. Beijing： Science Press， 2020.

［19］ LELE S K. Compact finite difference schemes with spectral-like" "resolution［J］." "Journal" "of" "computational physics， 1992， 103（1）： 16-42.

［20］ POST N. Fatigue test design： scenarios for biaxial fatigue testing of a 60-meter wind turbine blade［R］. NREL/TP-5000-65227， 1271941， 2016： NREL/TP-5000-65227， 1271941.

IDENTIFICATION OF WIND TURBINE BLADE ANTI-BENDING

STIFFNESS BASED ON COMPACT DIFFERENCE SCHEME

Ma Yi，Zhou Aiguo，Shi Jinlei，Zhao Shiwen，Zhu Yutian

（College of Mechanical Engineering， Tongji University， Shanghai 201804， China）

Abstract：Identifying the distribution of wind turbin blade flexural stiffness in practical scenarios poses a significant challenge. To address this issue， this study proposes a novel blade flexural stiffness identification method based on two key components： single point static load deflection fitting and a compact finite difference scheme. By establishing a comprehensive mathematical model that captures blade bending deformation under static load calibration conditions， this study derives precise expressions for deflection and anti-bending stiffness for each section of the blade. Extensive analysis on various blade types reveals that the anti-bending stiffness identification error in the midsection of the blade remains below 5%， thus validating the effectiveness of the proposed method.

Keywords：wind turbine blades； anti-bending stiffness； curve fitting； calibration； compact difference scheme