收稿日期：2022-11-02

基金項目：國家自然科學(xué)基金（51976113）

通信作者：陳 嚴(yán)（1964—），男，博士、教授、博士生導(dǎo)師，主要從事風(fēng)能利用及風(fēng)力機械方面的研究。ychen@stu.edu.cn

DOI：10.19912/j.0254-0096.tynxb.2022-1667 文章編號：0254-0096（2024）02-0143-09

摘 要：針對超長柔性風(fēng)電葉片在傳統(tǒng)線性（小變形）條件下發(fā)生動力學(xué)分析失效的問題，提出考慮大變形結(jié)構(gòu)的柔性葉片模型并進行氣動彈性分析。首先將葉片簡化成非均勻懸臂梁，利用大變形Euler-Bernoulli梁模型進行離散化建模；后結(jié)合葉素動量理論，建立葉片的非線性結(jié)構(gòu)動力學(xué)方程，并運用Newton-Raphson和逐步積分相結(jié)合的方式進行時域仿真。整個過程重點考慮大變形結(jié)構(gòu)的剛度變化、以及柔性葉片截面扭轉(zhuǎn)運動和外形變化。以DTU 10 MW風(fēng)力機為算例，比較葉片在大變形和小變形結(jié)構(gòu)下的動態(tài)氣動特性，結(jié)果表明：動態(tài)響應(yīng)中，大變形結(jié)構(gòu)的剛度變化導(dǎo)致?lián)]舞方向的葉尖位移明顯增加，略微減小擺振方向的位移，同時葉尖扭角也會顯著減小。此外，葉根面外載荷在大變形條件下會顯著下降，但葉根面內(nèi)載荷則會略微增加。這與傳統(tǒng)線性條件下的氣彈現(xiàn)象具有明顯差異。該文的葉片模型能有效預(yù)測長柔葉片的動態(tài)氣動特性。

關(guān)鍵詞：風(fēng)力機；柔性葉片；大變形結(jié)構(gòu)；氣動彈性；截面扭轉(zhuǎn)

中圖分類號：TK83""""""""nbsp; """"""""" """"""文獻標(biāo)志碼：A

0 引 言

風(fēng)能作為一種取之不盡、用之不竭的清潔能源，有助于減少人類對化石能源的依賴［1］。特別是海上風(fēng)能，由于具有資源豐富、風(fēng)況穩(wěn)定等優(yōu)點，有望在未來能源領(lǐng)域發(fā)揮重要作用。當(dāng)前，風(fēng)力機是在海上提取風(fēng)能的最有效手段之一［2］。為提高有限土地資源和海洋資源的空間風(fēng)能利用效率、降低度電成本，風(fēng)電單機容量大型化是必然發(fā)展趨勢。由此帶來的重量、制造成本極具增長問題，使得葉片被設(shè)計得更加輕柔，導(dǎo)致葉片更易產(chǎn)生變形。

在風(fēng)力機運行過程中，風(fēng)電葉片將同時承受氣動力、慣性力和彈性力的耦合作用，使得葉片產(chǎn)生多自由度耦合振動，加劇葉片的非線性變形。特別是在極端風(fēng)況條件下，葉片將面臨更加嚴(yán)峻的幾何大變形問題以及非線性結(jié)構(gòu)和復(fù)雜流場相互耦合下的氣動彈性問題［3］。因此，對風(fēng)電葉片大變形條件下的氣彈效應(yīng)了解不足已成為限制大型風(fēng)力機進一步發(fā)展的“瓶頸”，研究建立合適的葉片動力學(xué)模型并分析柔性大變形下的氣動彈性問題勢在必行。

目前用于建立葉片動力學(xué)模型的方法中，常見的有一維有限元模型（finite element method，F(xiàn)EM）。一維FEM是基于梁模型來建立葉片結(jié)構(gòu)，但目前葉片采用的一維FEM大部分是基于小變形梁結(jié)構(gòu)，對于大尺寸柔性細(xì)長葉片，存在大變形失效問題?？紤]一維FEM具有描述葉片非線性變形的潛力，諸多學(xué)者都在其基礎(chǔ)上進行了研究。Kalles?e［4］基于Hodegs［5］推導(dǎo)的偏微分運動方程，采用線性化方法建立大變形結(jié)構(gòu)葉片的氣動模型，研究表明穩(wěn)態(tài)條件下葉片大變形對擺陣方向氣動阻尼具有顯著影響；Manolas等［6］在其基礎(chǔ)上繼續(xù)研究結(jié)構(gòu)非線性對風(fēng)電葉片載荷的影響，研究表明彎扭耦合是引起載荷變化的主要非線性因素；Hodegs推導(dǎo)的模型最早用于研究直升機轉(zhuǎn)子葉片大彎曲和扭轉(zhuǎn)耦合的非線性問題上，該模型性能良好但空間結(jié)構(gòu)復(fù)雜。陳嚴(yán)等［7］考慮了葉片截面的扭轉(zhuǎn)運動，研究不同柔性變形條件對葉片動態(tài)氣動特性的影響，結(jié)果表明揮舞傾角和截面扭角是影響氣動特性的主要因素；張新榃等［8］考慮歐拉梁的幾何轉(zhuǎn)角建立大變形結(jié)構(gòu)葉片，研究某些極端工況下葉根的彎曲載荷情況；曹九發(fā)等［9］在歐拉梁基礎(chǔ)上加入應(yīng)變二階項，建立葉片的大變形結(jié)構(gòu)模型并進行穩(wěn)態(tài)響應(yīng)分析，結(jié)果表明結(jié)構(gòu)大變形對風(fēng)機穩(wěn)態(tài)條件下的氣動特性和葉片載荷具有較大影響；付志超等［10］采用本征梁模型研究大展弦比機翼幾何非線性條件下的穩(wěn)態(tài)氣彈特性，研究表明幾何非線性結(jié)構(gòu)對機翼氣彈特性具有顯著影響。但是無論國內(nèi)外，基于一維FEM的非線性方法受限于結(jié)構(gòu)方程的復(fù)雜程度，大部分只停留在穩(wěn)態(tài)氣動特性研究上。因此，本文將重點研究大變形條件下，柔性葉片模型的動態(tài)氣動特性。

本文將使用考慮橫縱向非線性位移場和平面 Green應(yīng)變的大變形Euler-Bernoulli梁模型，并考慮隨動坐標(biāo)系，通過將葉片離散成非均勻懸臂梁，結(jié)合葉片截面扭轉(zhuǎn)和外形變化，建立柔性葉片的非線性結(jié)構(gòu)動力學(xué)方程。后通過將Newton-Raphson法和Newmark法應(yīng)用在非線性結(jié)構(gòu)方程求解中，實現(xiàn)葉片氣彈耦合下的時域仿真。

1 結(jié)構(gòu)模型

1.1 大變形Euler-Bernoulli模型

區(qū)別線性梁模型［11］，在大變形條件下，等截面梁單元的變形如圖1所示。單元形心軸與x軸重合，單元局部坐標(biāo)系的原點[O]與單元左端點重合。形心軸上任意點的位移可分解為縱向位移、橫向位移和轉(zhuǎn)角，其中轉(zhuǎn)角可由橫向位移的導(dǎo)數(shù)進行描述；而梁單元的變形可分為縱向拉壓變形和橫向彎曲變形兩部分?？v向拉壓變形主要影響單元左右節(jié)點的縱向位移；而橫向彎曲變形除了影響單元左右節(jié)點的橫向位移，還顯著影響節(jié)點的縱向位移，因此必須考慮［12］。

首先對梁單元形心軸的橫向位移場進行描述，為了同時考慮橫向位移和轉(zhuǎn)角，采用3次Hermite插值進行建立，有橫向位移場[uy]：

[uy=Nyue]""""" （1）

式中：[ue]——單元節(jié)點的位移列陣，表達式為[ue=u1u2u3u4u5u6T]；[Ny]——梁單元橫向位移場對應(yīng)的形函數(shù)矩陣，[Ny=Ny1Ny2Ny3Ny4Ny5Ny6]，表達為：

[Ny1=Ny4=0Ny2=1-3ξ2+2ξ3Ny3=lξ-2ξ2+ξ3Ny5=3ξ2-2ξ3Ny6=l-ξ2+ξ3]"""" （2）

接著對梁形心軸的縱向位移場進行描述，包括縱向伸縮變形和橫向彎曲變形引起的附加縱向位移?？v向伸縮變形使用一次Lagrange函數(shù)建立，滿足：

[ux=Nxue]""""" （3）

式中：[Nx]——梁單元橫向伸縮位移場對應(yīng)的形函數(shù)矩陣，表達式為[Nx=1-ξ00ξ00]。

考慮因橫向彎曲變形產(chǎn)生的附加縱向位移，如圖1所示，取單元軸線上某微段[M1M2=dδ0≤δ≤x]，其移動后的位置為[M1′M2′]，則該微段的縱向附加位移為：

[dΔ=N1N2=M1M2-M1′M2′2-M2′N22]"""" （4）

由于純彎曲變形后微段長度不變，因此有[M1′M2′=dδ， M2′N2=duy，]由幾何關(guān)系推導(dǎo)，式（4）可表示為：

[dΔ=dδ1-1-?uy?δ2]"""""" （5）

進一步簡化，對式（5）的根號項進行展開，并忽略二階以上的高次項，有：

[dΔ=12?uy?δ2dδ]" （6）

對上式積分可得單元上任意一點的附加縱向位移[Δ=0x12?uy?δ2dδ]。

綜上，有梁單元形心軸的位移場為：

[u=ux-120x?uy?δ2dδuy=Nue]""""" （7）

式中：[N]——梁單元完整的形函數(shù)矩陣。

由二階位移場可知，梁單元的應(yīng)變可分為縱向拉壓應(yīng)變[εa]和橫向彎曲應(yīng)變[εb]兩部分。當(dāng)考慮結(jié)構(gòu)的幾何大變形時，基于平面Green應(yīng)變，工程應(yīng)變[ε]為：

[ε=ε0+εL]"""""" （8）

式中：[ε0、][εL]——工程應(yīng)變的線性和非線性部分，[ε0=ε0aε0b=duxdx-yd2uydx2，εL=12duydx20]。

該梁單元的勢能為：

[Πe=Ue-We=VεTσdV-ueTp]""" （9）

式中：[Ue]——應(yīng)變能；[We]——外力功；[σ]——應(yīng)力列陣；[ue]——單元節(jié)點位移；[p]——梁單元上作用的外力。

基于本構(gòu)關(guān)系，有應(yīng)力列陣[σ=Dε]。[D]為單元的彈性矩陣。

由最小勢能原理，即δ=0，滿足：

[VεTσdV-δueTp=0]" （10）

將式（7）代入式（8），可將應(yīng)變表達式轉(zhuǎn)變?yōu)椋?/p>

[ε=Bue] （11）

式中：[B]——單元的應(yīng)變矩陣。

式（10）可表示為：

[VBTσdV=p]"""""" （12）

根據(jù)式（8），應(yīng)變列陣[ε]含有非線性項[εL]，由此可知應(yīng)變列陣[ε]是單元節(jié)點位移列陣[ue]的函數(shù)。因此，可將單元應(yīng)變矩陣[B]進行分解，滿足：

[B=B0+BL]""" （13）

式中：[B0]——單元的線性應(yīng)變矩陣，矩陣元素與節(jié)點位移無關(guān)，屬于小變形條件下的應(yīng)變矩陣；[BL]——單元的非線性應(yīng)變矩陣，矩陣元素與節(jié)點位移相關(guān)。對式（12）取變分，滿足：

[VδBLTσdV+VBTδσdV=δp]"""""" （14）

將式（13）代入式（14），則式（14）等號左邊第一項可表示為：

[VδBLTσdV=kσδue]""""" （15）

式中：[kσ]——單元的幾何剛度矩陣，即“離心剛化”項，它表示單元縱向應(yīng)力對單元剛度矩陣的影響。

式（14）等號左邊第二項可表示為：

[VBTδσdV=VBT0DB0dVδue+VBT0DBLdV+VBTLDB0dV+VBTLDBLdVδue"""""" =k0δue+kLδue]

（16）

式中：[k0]——單元的線性剛度矩陣，屬于小變形條件下的剛度矩陣；[kL]——大變形矩陣，表示變形后的單元剛度矩陣相較于對變形前所產(chǎn)生的附加剛度矩陣。

綜上，有大變形條件下，單元增量平衡方程為：

[keδue=δp] （17）

式中：[ke]——梁單元的剛度矩陣，表達式為["ke=k0+kL+kσ]。

在非線性結(jié)構(gòu)分析中，本文采用修正拉格朗日法，即在增量載荷條件下，采用Newton-Raphson法進行迭代，并在迭代過程中，以結(jié)構(gòu)上一時刻的構(gòu)型為參考構(gòu)型，最終求解結(jié)構(gòu)的非線性響應(yīng)。同時，非線性響應(yīng)是針對完整結(jié)構(gòu)，需在整體坐標(biāo)系下對所有單元剛度矩陣進行組裝，考慮上述的單元切線剛度矩陣是在局部坐標(biāo)系下推導(dǎo)的，因此需轉(zhuǎn)換到整體坐標(biāo)系中。

另外，上文的大變形Euler-Bernoulli梁理論是在小轉(zhuǎn)角假設(shè)基礎(chǔ)上建立的，為減少小轉(zhuǎn)角假設(shè)對幾何非線性的影響，更好地描述結(jié)構(gòu)的大位移情況，引入隨動坐標(biāo)系，即在每次迭代計算過程中，都會對結(jié)構(gòu)的坐標(biāo)進行修正，通過隨動的局部坐標(biāo)對單元的位移情況進行更詳盡地描述，此種方法稱為共旋坐標(biāo)法［13］。

1.2 葉片非線性結(jié)構(gòu)動力學(xué)方程

研究風(fēng)電葉片的動態(tài)響應(yīng)時，本文固定葉片根部，將其近似看作懸臂梁并進行單元離散。隨后基于虛功原理，對離散化單元建立結(jié)構(gòu)動力學(xué)方程［14］：

[Mv+Cv+Kv=F]"""" （18）

式中：[M]——質(zhì)量矩陣；[C]——阻尼矩陣；[K]——剛度矩陣；[F]——外載荷向量；[v]——葉片變形量。

葉片結(jié)構(gòu)采用2節(jié)點Euler-Bernoulli梁單元，每個節(jié)點含有3個自由度［[u1" u2" u3]］，包括展向位移、揮舞（或擺振）方向位移和轉(zhuǎn)角。將1.1節(jié)推導(dǎo)的非線性剛度矩陣[ke]代入式（18），當(dāng)考慮大變形Euler-Bernoulli梁單元時，有增量條件下的單元非線性動力學(xué)方程為：

[meqe+ceqe+ktedqe=dfe]""""" （19）

式中：[dqe]——增量條件下的單元不平衡位移；[dfe]——單元增量載荷；[me、][ce、][kte]——局部坐標(biāo)系下的單元質(zhì)量矩陣、單元阻尼矩陣和單元切線剛度矩陣。

上述的切線剛度矩陣，是基于梁截面第一主軸的彈性剛度分布推導(dǎo)得到，而實際葉片的翼型截面沿展向呈螺旋狀分布，因此必須考慮截面扭角的影響，將梁的彈性剛度轉(zhuǎn)到葉片中。如圖2所示，[O]為翼型截面形心，由幾何關(guān)系，截面在兩個主軸方向x、y上的抗彎剛度分別滿足：

[EIx=EIFcos2θb+EIcsin2θbEIy=EIFsin2θb+EIccos2θb]""" （20）

式中：[IF]——葉片截面對軸[ξ]的慣性矩，[IF=Aη2dA，]m4；[Ic]——葉片截面對軸[η]的慣性矩，[Ic=Aξ2dA，]m4；[θb]——截面扭角，（ °）。

通過坐標(biāo)變換將單元動力學(xué)方程轉(zhuǎn)到整體坐標(biāo)系下并組裝，即可得到完整結(jié)構(gòu)在整體坐標(biāo)系下的增量動力學(xué)方程：

[MQ+CQ+KtdQ=dF]"""" （21）

式中：[M]——組裝后的葉片質(zhì)量矩陣；[C]——組裝后的葉片阻尼矩陣；[Kt]——組裝后的葉片切線剛度矩陣；[Q]——組裝后的葉片不平衡位移；[dF]——組裝后的外載荷向量。

2 葉片載荷模型

風(fēng)電葉片的氣彈動力學(xué)分析基礎(chǔ)是計算作用在風(fēng)輪上的氣動載荷；此外，如離心力、重力等非結(jié)構(gòu)化載荷也具有重要影響。同時，為了更準(zhǔn)確地描述風(fēng)電葉片的結(jié)構(gòu)特性，基于IEC 61400-1標(biāo)準(zhǔn)［15］，對風(fēng)力機建立統(tǒng)一的參考坐標(biāo)系。如圖3所示，有X-Y-Z是塔架坐標(biāo)系，[Xn-Yn-Zn]是機艙坐標(biāo)系，[Xh-Yh-Zh]是輪轂坐標(biāo)系，[Xb-Yb-Zb]是葉片坐標(biāo)系。輪轂主軸與機艙主軸間的夾角[η]，稱為俯仰角；葉片轉(zhuǎn)動到各個位置與豎直向上的夾角[ψ]，稱為方位角；當(dāng)葉片在擺振方向振動時，對應(yīng)方向會有1個振動傾角[k]；同理在揮舞方向有1個振動傾角[β]；當(dāng)葉片發(fā)生扭轉(zhuǎn)振動時，展向會有1個附加的扭轉(zhuǎn)振動。

結(jié)合風(fēng)機坐標(biāo)系變換，有葉片坐標(biāo)系下，入流風(fēng)速[U]沿[x、][y]方向的分量為：

[U=UxUy=110TbladeTrotorThub0U∞0]"""" （22）

式中：[Tblade、][Trotor、][Thub]——輪轂到葉片、機艙到輪轂、塔架到機艙的坐標(biāo)變換矩陣。

2.1 葉片的氣動載荷

在氣動載荷方面，葉素動量理論（blade element moment，BEM）是應(yīng)用最廣泛的理論，其效率高、穩(wěn)定性好，適合非線性結(jié)構(gòu)復(fù)雜的迭代計算流程，因此本文基于BEM理論來建立氣動模型。在穩(wěn)態(tài)條件下，葉素動量理論包含“動量理論”和“葉素理論”。根據(jù)動量理論，有作用在致動盤上的單元軸向推力[FT]和扭矩[FQ]分別表示為：

[dTdr=4πρU2ra1-a]""" （23）

[dQdr=4πρΩUr3a′1-a]"" （24）

式中：[ρ]——空氣密度，kg/m3；r——到旋轉(zhuǎn)軸的展向距離，m；[a]、[a′]——軸向和切向誘導(dǎo)因子；[Ω]——風(fēng)輪角速度，rad/s。

在葉素理論中，作用在葉片單元上的氣動載荷由法向分量和切向分量以及氣動扭矩組成，分別表示為：

[dFxdr=12ρW2cCLcosφ+CDsinφ]"" （25）

[dFydr=-12ρW2cCLsinφ-CDcosφ]"" （26）

[dMaerodr=12ρW2ceCLcosα+eCDsinα+CM]" （27）

式中：[W]——相對入流風(fēng)速，m/s；[CL]、[CD]——截面的升力和阻力系數(shù)；[φ]——入流角，（ °）；[c]——翼型弦長，m。

葉片截面參數(shù)如圖3b所示，[E]為截面扭轉(zhuǎn)剛心（距離前沿1/4處），[o]為氣動中心，[e]為剛心和氣動中心間的距離，[α]為攻角，[θ]為扭角，[Δθ]為考慮截面扭轉(zhuǎn)的附加扭角［7］。對于柔性葉片，有考慮葉片扭轉(zhuǎn)運動和外形變化下的氣動參數(shù)為：

[φ=tan-1cosβcosη-sinβcosψsinη1-a-sinσsinβcosη--cosσsinψ+sinσcosβcosψsinη+ΩrUcosβ1+a′]"""""" （28）

[α=φ-θ-Δθ]""" （29）

[a=18πrcosβcosσ?F?sin2φBcCLcosφ+CDsinφ+1]" （30）

[a′=18πcosβcosσ?F?sinφcosφBcCLsinφ-CDcosφ-1]""" （31）

式中：[F]——損失因子；[B]——葉片數(shù)量。

由式（28）～式（31）可知，當(dāng)柔性葉片因揮舞、擺振、扭轉(zhuǎn)方向振動而導(dǎo)致外形變化時，氣動參數(shù)會隨之變化，從而改變?nèi)~片氣動載荷。

2.2 葉片的非結(jié)構(gòu)化載荷

在非結(jié)構(gòu)化載荷方面，有葉片單元相較于旋轉(zhuǎn)軸的離心力為：

[Fc=miΩ2r]"" （32）

式中：[mi]——葉片單元質(zhì)量，kg。

考慮葉片在揮舞、擺振方向上彎曲變形的影響，有分段葉片的有效半徑為[re=rcosβcosΔσ]；同時，由于葉片的離心力計算是在輪轂坐標(biāo)系下進行的，必須轉(zhuǎn)換到葉片坐標(biāo)系下，因此有：

[FcedgeFcflapFcaxial=Tblade00Fc=FcsinσcosβFcsinβFccosσcosβ]"" （33）

式中：[Fcedge]、[Fcflap]、[Fcaxial]——變換后在擺振、揮舞和展向3個方向上的離心力，N。

由于葉片截面上存在著質(zhì)心、剛心，而葉片離心力是基于翼型的質(zhì)心來計算的，因此考慮在揮舞方向上，離心載荷會對剛心產(chǎn)生一個力矩[Mc]，滿足：

[Mc=Fcflapecgcosθ]""" （34）

式中：[ecg]——截面上質(zhì)心和剛心間的距離，m。

由于葉片重力載荷隨著葉片旋轉(zhuǎn)方位角的變化而改變，考慮塔架-機艙-輪轂-葉片坐標(biāo)系的變化，滿足：

[FgedgeFgflapFgaxial=TbladeTrotorThub00-mig]"""""" （35）

[Mg=Fgflapecgcosθ]""" （36）

式中：[g]——重力加速度，m/s2。

在綜合考慮氣動載荷、重力和離心力載荷之后，即可得到柔性葉片各方向載荷的準(zhǔn)確表達。

3 算例分析

3.1 程序的非線性迭代流程

結(jié)合第1節(jié)和第2節(jié)，即可得到葉片完整的非線動力學(xué)方程，隨后進行葉片動態(tài)響應(yīng)分析。動態(tài)響應(yīng)的基本求解思路：在一個時間步內(nèi)，首先計算葉片所受的完整載荷，然后計算葉片的剛度矩陣、質(zhì)量矩陣和阻尼矩陣；基于上述參數(shù)建立動力學(xué)方程并求解動態(tài)響應(yīng)，響應(yīng)結(jié)果會反饋到葉片氣動載荷中，對葉片載荷更新并進入下個時間步的計算，由此循環(huán)［11］。

區(qū)別線性動力學(xué)方程，在非線性結(jié)構(gòu)中，剛度矩陣是位移響應(yīng)的函數(shù)，計算過程中會不斷更新，因此只能迭代求解。本文采用Newmark-β和增量Newton-Raphson方法［15］相結(jié)合的方式，在一個時間步內(nèi)，通過將外載荷劃分多個加載步形式，每個加載步再進行多次迭代，每次迭代對剛度矩陣更新并基于Newmark-β求解響應(yīng)，直至響應(yīng)結(jié)果滿足收斂條件。圖4為葉片非線性動態(tài)響應(yīng)的計算流程。

3.2 DTU10 MW氣動模型

基于上述推導(dǎo)的非線性氣彈響應(yīng)理論，以DTU 10 MW水平軸風(fēng)力機為對象，研究大變形結(jié)構(gòu)下葉片模型的動態(tài)響應(yīng)。葉片各截面翼型和詳細(xì)參數(shù)同文獻［17］。

圖5為各穩(wěn)定風(fēng)況下指定翼型截面的升力系數(shù)[CL]和阻力系數(shù)[CD]。

為了確保風(fēng)力機氣動模型的準(zhǔn)確性，首先通過Matlab對風(fēng)力機葉片進行參數(shù)建模，并將計算結(jié)果與專業(yè)軟件Bladed 4.3計算結(jié)果進行對比分析。

不難看出，本文計算結(jié)果與Bladed結(jié)果具有較好的契合度，說明本文建立的風(fēng)力機氣動模型具有一定的準(zhǔn)確性和可靠性。

3.3 DTU 10 MW線性模型驗證

首先在線性條件下驗證葉片結(jié)構(gòu)模型的準(zhǔn)確性。通過將某一時長50 s的湍流風(fēng)場仿真數(shù)據(jù)，輸入本人編寫的Matlab程序，并將結(jié)果與Bladed結(jié)果比較。如圖6a和圖6b所示，葉尖在揮舞和擺陣方向上的線性響應(yīng)變形量均與Bladed吻合，說明基礎(chǔ)線性模型的有效性；稍微的數(shù)據(jù)差異主要考慮Bladed內(nèi)置使用模態(tài)疊加法，從而在結(jié)果上受到模態(tài)截斷的影響。

線性葉片結(jié)構(gòu)模型的準(zhǔn)確性，也從側(cè)面說明本文基于線性Euler-Bermoulli梁推導(dǎo)的大變形結(jié)構(gòu)葉片模型可靠。

3.4 計算結(jié)果及結(jié)論

3.4.1 穩(wěn)態(tài)響應(yīng)結(jié)果

穩(wěn)態(tài)響應(yīng)工況取常數(shù)風(fēng)速11.4 m/s，槳距角1.9°，風(fēng)輪轉(zhuǎn)速0.99 rad/s進行計算。

圖7a～圖7c分別為穩(wěn)態(tài)條件下，葉尖位移和扭轉(zhuǎn)隨方位角的變化情況。

在揮舞方向上，對比圖7a中的計算結(jié)果曲線，可看出當(dāng)風(fēng)機轉(zhuǎn)動一周時，葉尖揮舞位移均為正值；在整個周期內(nèi)（特別是在50°～300°方位角范圍內(nèi)），有考慮結(jié)構(gòu)大變形的葉尖揮舞位移明顯增加，并在209°時達到峰值7.48 m，最大差值約為0.47 m。

由圖7b可看出，在擺振方向上，葉尖位移結(jié)果在-1.0～0.4 m之間波動。對比兩條計算結(jié)果曲線可知，有考慮結(jié)構(gòu)大變形的葉尖位移絕對值明顯減小，葉尖位移振幅降低，在256°達到最大振幅-0.75 m，最大差值約為0.26 m。對比圖7c中的計算結(jié)果曲線可知，在扭轉(zhuǎn)方向上，當(dāng)考慮結(jié)構(gòu)大變形時，葉尖扭角變化量明顯減小。

對比穩(wěn)態(tài)條件下葉尖位移的計算結(jié)果，有無考慮結(jié)構(gòu)大變形對柔性葉片葉尖的位移和扭角具有顯著影響。發(fā)生顯著差異的重要原因是大變形條件下，每個時刻的葉片結(jié)構(gòu)剛度矩陣受上一時刻響應(yīng)位移的影響，且這種影響是逐漸累積并影響到后續(xù)時刻，即葉片當(dāng)前構(gòu)型與初始構(gòu)型直接掛鉤，理論上更能準(zhǔn)確反映葉片的結(jié)構(gòu)特性，進而影響到葉片外形姿態(tài)（揮舞傾角、擺振傾角和扭角），最終影響到葉尖響應(yīng)結(jié)果。同時，葉片每個單元均通過隨動坐標(biāo)系更準(zhǔn)確地描述了局部坐標(biāo)系下的大位移情況，即變形前后的位姿狀態(tài)，使計算結(jié)果進一步精確。

3.4.2 動態(tài)響應(yīng)結(jié)果

動態(tài)響應(yīng)工況取平均風(fēng)速11.4 m/s，湍流強度15.21，仿真時間150 s進行計算。

圖8a～圖8c分別為湍流工況下，葉尖位移和扭角隨時間的動態(tài)響應(yīng)情況。

在揮舞方向上，對比圖8a兩種計算結(jié)果曲線可知，是否考慮大變形對葉尖位移振型并無明顯影響，整體葉尖位移均為正值。但在位移變化量方面，結(jié)合圖9葉尖揮舞位移的差值曲線放大圖可知，當(dāng)考慮結(jié)構(gòu)大變形時，葉尖位移的振幅明顯增加（峰值更大、谷值更?。?，有大變形條件下葉尖揮舞方向的平均位移量為5.79 m，而小變形條件下的葉尖平均位移量為5.63 m。動態(tài)響應(yīng)結(jié)果表明葉片的柔性大變形會導(dǎo)致葉尖揮舞位移量增大，這對動態(tài)響應(yīng)中葉片和塔架的間隙提出了更大的要求。

在擺振方向上，由圖8b可知，是否考慮大變形對葉尖位移振型同樣影響不大，葉尖整體位移量偏小，在-1.2～0.6 m間正負(fù)波動，說明對擺振位移起決定作用的仍是重力載荷。在位移變化量方面，當(dāng)考慮結(jié)構(gòu)大變形時，葉尖位移的振幅相對減小，有大變形條件下葉尖的平均位移絕對值為0.24 m，小變形條件下的葉尖平均位移絕對值為0.33 m。動態(tài)響應(yīng)結(jié)果表明葉片的柔性大變形會導(dǎo)致葉尖擺振位移量減小。

扭轉(zhuǎn)方向，由圖8c可知，當(dāng)考慮結(jié)構(gòu)大變形時，葉尖扭轉(zhuǎn)變化量顯著減小。

圖10為兩種模型狀態(tài)下，葉根面外和面內(nèi)方向的載荷差值變化。由圖10b可知，大變形條件下葉根面外載荷結(jié)果明顯小于小變形結(jié)果，下降幅度集中在5～15 kN，這主要是葉根面外載荷主要來源于氣動載荷，由于葉片在大變形條件下整體姿態(tài)發(fā)生顯著變化，從而影響了氣動載荷并導(dǎo)致葉根面外載荷下降。由圖10b可知，大部分時間內(nèi)，大變形條件下的葉根面內(nèi)載荷相對增加，主要集中在5 kN以內(nèi)，其他時間變化不明顯，這考慮是葉片姿態(tài)變化對重力載荷產(chǎn)生了一定影響。對比兩個計算結(jié)果表明，考慮柔性大變形能顯著降低動態(tài)響應(yīng)過程中葉根面外載荷，但會一定程度地增加葉根面內(nèi)載荷。

綜上，對于早期風(fēng)電葉片，由于葉片較短，彎曲變形較小，故一般分析時不考慮結(jié)構(gòu)柔性大變形對葉片動態(tài)響應(yīng)的影響。但隨著風(fēng)力機單機容量大型化，特別是海上風(fēng)力機快速發(fā)展，結(jié)構(gòu)柔性和大變形對葉片揮舞、擺振和扭轉(zhuǎn)方向的穩(wěn)定性分析、動態(tài)性能、載荷計算的影響將不容忽視，這對未來風(fēng)力機的指導(dǎo)設(shè)計和減重降本具有重要意義。

4 結(jié) 論

本文通過對超長柔性風(fēng)電葉片進行氣動和大變形結(jié)構(gòu)建模，研究大變形結(jié)構(gòu)下葉片在動態(tài)響應(yīng)中的氣彈特性，得到如下結(jié)論：

1）本文通過湍流風(fēng)場下的線性響應(yīng)結(jié)果與Bladed結(jié)果進行比較，說明本文基于線性結(jié)構(gòu)的葉片模型的具有較高精度；而大變形結(jié)構(gòu)模型是在小變形結(jié)構(gòu)基礎(chǔ)上進行建立，表明本文基于大變形結(jié)構(gòu)的葉片模型能較好地預(yù)測長柔葉片的氣動特性。

2）動態(tài)響應(yīng)過程中，是否考慮結(jié)構(gòu)大變形對葉尖位移振型影響不大，但會顯著影響葉尖位移量。采用大變形結(jié)構(gòu)會導(dǎo)致?lián)]舞方向的葉尖位移明顯增加，略微減小擺振方向的葉尖位移，同時葉尖扭角變化量也會明顯減小，這主要認(rèn)為是大變形條件下葉片結(jié)構(gòu)剛度受前面時刻的影響，且這種影響持續(xù)累計，直接影響葉片結(jié)構(gòu)特性并導(dǎo)致葉片姿態(tài)發(fā)生改變，最終影響葉尖位移響應(yīng)。這對動態(tài)響應(yīng)中葉片和塔架的間隙提出了更大的要求。

3）動態(tài)響應(yīng)過程中，采用大變形結(jié)構(gòu)會顯著降低葉根面外載荷，下降幅度集中在5～15 kN；相反，在大部分時間內(nèi)葉根面內(nèi)載荷會略微增加，增幅在5 kN以內(nèi)。其原因是葉根面外載荷主要來源于氣動載荷，由于葉片在大變形條件下整體姿態(tài)發(fā)生顯著變化，從而影響了氣動載荷并導(dǎo)致葉根面外載荷下降；而面內(nèi)方向考慮是葉片姿態(tài)變化對重力載荷產(chǎn)生了一定影響導(dǎo)致。這對未來風(fēng)力機減重降本具有重要意義。

［參考文獻］

［1］"""" JOSELIN HERBERT G M， INIYAN S， SREEVALSAN E， et al. A review of wind energy technologies［J］. Renewable and sustainable energy reviews， 2007， 11（6）： 1117-1145.

［2］"""" LEUNG D Y C， YANG Y. Wind energy development and its environmental impact： a review［J］. Renewable and sustainable energy reviews， 2012， 16（1）： 1031-1039.

［3］"""" VEERS P， DYKES K， LANTZ E， et al. Grand challenges in the science of wind energy［J］. Science， 2019， 366（6464）： eaau2027.

［4］"""" KALLES?E B S. Equations of motion for a rotor blade， including gravity， pitch action and rotor speed variations［J］. Wind energy， 2007， 10（3）： 209-230.

［5］"""" HODGES D H. Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades［R］. NASA TN D-7818， 1974.

［6］"""" MANOLAS D I， RIZIOTIS V A， VOUTSINAS S G. Assessing the importance of geometric nonlinear effects in the prediction of wind turbine blade loads［J］. Journal of computational and nonlinear dynamics， 2015， 10（4）： 041008.

［7］"""" 陳嚴(yán)， 方郁鋒， 劉雄， 等. 考慮風(fēng)力機葉片截面扭轉(zhuǎn)運動及外形變化的氣動彈性分析［J］. 太陽能學(xué)報， 2014， 35（4）： 553-561.

CHEN Y， FANG Y F， LIU X， et al. Aeroelastic analysis on section twist and shape change of wind turbine blade［J］. Acta energiae solaris sinica， 2014， 35（4）： 553-561.

［8］"""" 張新榃， 楊超， 吳志剛. 基于大變形梁模型的風(fēng)力機葉片氣動彈性分析［J］. 北京航空航天大學(xué)學(xué)報， 2014， 40（3）： 327-332.

ZHANG X T， YANG C， WU Z G. Aeroelastic analysis method for wind turbine blade based on large deformation beam model［J］. Journal of Beijing University of Aeronautics and Astronautics， 2014， 40（3）： 327-332.

［9］"""" 曹九發(fā)， 王同光， 王梟. 大型風(fēng)力機的幾何非線性動態(tài)響應(yīng)研究［J］. 振動與沖擊， 2016， 35（14）： 182-187， 214.

CAO J F， WANG T G， WANG X. Geometric nonlinear dynamics response of large-scale wind turbine［J］. Journal of vibration and shock， 2016， 35（14）： 182-187， 214.

［10］""" 付志超， 陳占軍， 劉子強. 大展弦比機翼氣動彈性的幾何非線性效應(yīng)［J］. 工程力學(xué)， 2017， 34（4）： 231-240.

FU Z C， CHEN Z J， LIU Z Q. Geometric nonlinear aeroelastic behavior of high aspect ratio wings［J］. Engineering mechanics， 2017， 34（4）： 231-240.

［11］""" 陳嚴(yán)， 張林偉， 劉雄， 等. 水平軸風(fēng)力機柔性葉片氣動彈性響應(yīng)分析［J］. 太陽能學(xué)報， 2014， 35（1）： 74-82.

CHEN Y， ZHANG L W， LIU X， et al. Response analysis of aeroelasticity for hawt flexible blade［J］. Acta energiae solaris sinica， 2014， 35（1）： 74-82.

［12］""" 夏擁軍， 繆謙. Euler-Bernoulli梁單元的完整二階位移場［J］. 中國工程機械學(xué)報， 2011， 9（4）： 416-420.

XIA Y J， MIAO Q. Complete second-order displacement field of Euler-Bernoulli beam element［J］. Chinese journal of construction machinery， 2011， 9（4）： 416-420.

［13］""" GALVANETTO U， CRISFIELD M A. An energy-conserving" co-rotational" procedure" for" the" dynamics" of planar"" beam"" structures［J］." International"" journal"" for numerical methods in engineering， 1996， 39（13）： 2265-2282.

［14］""" CLOUGH R W， PENZIEN J. Dynamics of structures［M］. 2nd ed. New York： McGraw-Hill， 1993.

［15］""" IEC 61400-1： 2019， Wind energy generation systems—part 1： design requirements.

［16］""" 陳至達. 桿、板、殼大變形理論［M］. 北京： 科學(xué)出版社， 1994.

CHEN Z D. Large deformation theory of rod， plate and shell［M］. Beijing： Science Press， 1994.

［17］""" BAK C， ZAHLE F， BITSCHE R D， et al. The DTU 10-MW" reference" wind" turbine［C］//Danish" Wind" Power Research 2013. Trinity， Fredericia， Denmark， 2013.

DYNAMIC RESPONSE OF FLEXIBLE WIND TURBINE BLADES CONSIDERING LARGE-DEFLECTION

Wang Zedong1，Wang Dian1，Qi Liangwen2，Chen Yan3，Huang Zhihong3

（1. CRRC Zhuzhou Institute Co.，Ltd.， Zhuzhou 412000， China；

2. Institute of Engineering Thermophysics， Chinese Academy of Sciences， Beijing 100049， China；

3. Key Laboratory of Intelligent Manufacturing Technology， Ministry of Education， Shantou University， Shantou 515063， China）

Abstract：Aiming at the problem of dynamic analysis failure of ultra-long flexible wind turbine blades under traditional linear （small deformation） conditions， a flexible blade model with the large-deflection structure is developed to study the dynamic response. The blade is simplified as a cantilever and discretely modeled by applying the large-deflection model which is derived from the Euler-Bernoulli beam. Then， the nonlinear aeroelastic dynamics equations of the blade model are established by using blade element momentum （BEM） theory， and time-domain simulations are performed by using the nonlinear Newmark method and Newton-Raphson method. The changes of the structural stiffness and the pose of the blade are considered throughout the process. Finally， the nonlinear model is applied to the DTU 10 MW wind turbine. The results show that the large-deflection model significantly increases the flapwise deflection and torsion angle of the blade tip in dynamic response， while slightly reduce the edgewise deflection. Meawhile， the large-deflection model significantly reduces the flapwise load of the blade root in dynamic response， while slightly increase the edgewise load of the blade root， which is different from the result of the linear model.

Keywords： wind turbines； flexible blade； large-deflection； aeroelasticity； cross-section twist