錢曉航，郜志騰，王同光，*，王 瓏，柯世堂，2

(1. 南京航空航天大學(xué) 江蘇省風(fēng)力機(jī)設(shè)計(jì)高技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，南京 210016；2. 南京航空航天大學(xué) 民航學(xué)院，南京 210016)

0 引言

隨著風(fēng)力機(jī)復(fù)合材料葉片尺寸的增加，幾何非線性效應(yīng)、截面面內(nèi)和面外翹曲等非經(jīng)典的效應(yīng)，對(duì)葉片結(jié)構(gòu)動(dòng)力學(xué)動(dòng)態(tài)響應(yīng)產(chǎn)生顯著影響。然而由于彈性耦合影響，復(fù)合材料結(jié)構(gòu)的分析要比各向同性結(jié)構(gòu)的更難[1]。傳統(tǒng)的模態(tài)疊加法仍然適用于計(jì)算小功率風(fēng)力機(jī)的彈性變形。然而對(duì)于大型風(fēng)力機(jī)而言，在復(fù)雜工況下可能會(huì)發(fā)生大的變形，在這時(shí)，傳統(tǒng)的線性方法無(wú)法再準(zhǔn)確地預(yù)測(cè)葉片的動(dòng)態(tài)響應(yīng)。

為帶有預(yù)扭和彎曲的復(fù)合材料葉片開(kāi)發(fā)一種考慮幾何非線性的梁模型是風(fēng)力機(jī)領(lǐng)域的一個(gè)新的焦點(diǎn)。葉片截面上的應(yīng)變通常是小應(yīng)變，因此，幾何非線性[2-3]主要是由葉片截面的有限旋轉(zhuǎn)造成的。傳統(tǒng)的線性分析方法是將葉片簡(jiǎn)化為歐拉伯努利梁模型[4]，并采用模態(tài)疊加法進(jìn)行求解，但此方法沒(méi)有考慮扭轉(zhuǎn)自由度，精度不夠。文獻(xiàn)[5]中采用微分求積單元求解幾何非線性，而本文應(yīng)用了基于Lengendre譜有限元的幾何精確梁理論[6-8]，這個(gè)理論是基于經(jīng)典鐵木辛柯梁理論[9-10]且考慮了截面旋轉(zhuǎn)發(fā)展而來(lái)。這個(gè)結(jié)構(gòu)模型包含了耦合的揮舞、擺振和扭轉(zhuǎn)自由度。目前，氣動(dòng)計(jì)算中常用的方法是AL-LES方法和葉素動(dòng)量理論等。由于AL-LES方法[11-12]計(jì)算的復(fù)雜性，于是在本文中用幾何精確梁理論與葉素動(dòng)量理論[13-15]耦合來(lái)建立風(fēng)力機(jī)葉片氣動(dòng)彈性模型。這個(gè)模型能準(zhǔn)確計(jì)算在氣動(dòng)載荷下的葉片變形并且充分考慮變形對(duì)氣動(dòng)彈性穩(wěn)定性造成的幾何非線性影響。

本文主要將幾何非線性分析方法應(yīng)用在兩種不同功率的大型風(fēng)力機(jī)中的大柔性葉片，首先采用幾何精確梁理論計(jì)算懸臂梁變形，并與理論值進(jìn)行對(duì)比，驗(yàn)證結(jié)構(gòu)計(jì)算方法的可靠性與準(zhǔn)確性。然后通過(guò)葉素動(dòng)量理論與幾何精確梁方法耦合計(jì)算，實(shí)現(xiàn)大型風(fēng)力機(jī)5 MW和15 MW的動(dòng)態(tài)響應(yīng)計(jì)算。

1 理論描述

1.1 葉素動(dòng)量理論

葉素動(dòng)量方法實(shí)際上是葉素理論和動(dòng)量理論的結(jié)合。氣流在流管內(nèi)流動(dòng)滿足動(dòng)量定理，但是氣流在流經(jīng)葉片時(shí)會(huì)受到擾動(dòng)，從而導(dǎo)致了氣流的切向和軸向速度發(fā)生變化，通常引入切向和軸向誘導(dǎo)因子來(lái)反應(yīng)氣流在通過(guò)風(fēng)輪平面時(shí)速度的損失。葉素理論將葉片沿展向離散為有限數(shù)量的葉素，葉素隨著風(fēng)輪旋轉(zhuǎn)形成了一個(gè)圓環(huán)，沿展向?qū)ι?、阻力積分便可求得氣動(dòng)推力和扭矩。葉片一個(gè)葉素微元的受力如圖1所示。

圖1 葉素上的來(lái)流速度與氣動(dòng)力示意圖Fig. 1 Inflow velocity and aerodynamic forces on the blade

W為入流合速度，速度合成關(guān)系為：

式中 ?為入流角，V0為入流速度， ?為風(fēng)輪轉(zhuǎn)速，a為軸向誘導(dǎo)因子，a′為切向誘導(dǎo)因子。

葉素升力、阻力表達(dá)式為：

式中：ρ為空氣密度，c為二維翼型弦長(zhǎng)，CL、CD為二維翼型升力、阻力系數(shù)。

通過(guò)將葉段升、阻力分解到風(fēng)輪旋轉(zhuǎn)平面及垂直于風(fēng)輪旋轉(zhuǎn)平面的方向，得到了切向力Ftq和法向力Fth：

葉素的剖面是一系列基本翼型，而基本翼型的氣動(dòng)數(shù)據(jù)一般作為輸入條件來(lái)確定不同入流條件下的氣動(dòng)力，輸入的氣動(dòng)數(shù)據(jù)一般為不同攻角α下的升力系數(shù)CL和阻力系數(shù)CD，翼型的升、阻力系數(shù)根據(jù)局部風(fēng)速和攻角，通過(guò)翼型數(shù)據(jù)表線性插值獲得，翼型軸向力系數(shù)Cth、切向力系數(shù)Ctq與升、阻力系數(shù)的關(guān)系為：

考慮葉片數(shù)量B，微元上受到的推力和轉(zhuǎn)矩通過(guò)葉素理論描述為：

根據(jù)動(dòng)量理論和葉素理論下的推力和轉(zhuǎn)矩公式相等，可得軸向誘導(dǎo)因子a和 切向誘導(dǎo)因子a′的表達(dá)式：

式中，σ為葉片局部長(zhǎng)度， σ=Bc/(2πr)。

1.2 幾何精確梁理論

幾何精確梁理論以受初始彎扭的梁能承受大的位移和旋轉(zhuǎn)的能力為特點(diǎn)。通過(guò)一個(gè)三維橫截面分析，六個(gè)自由度的所有耦合效應(yīng)，包括拉伸、剪切、扭轉(zhuǎn)、彎曲、扭轉(zhuǎn)翹曲、面內(nèi)翹曲等，都能被幾何精確梁理論涵蓋?！皫缀尉_”指的是在公式中沒(méi)有對(duì)初始幾何形狀和變形后的幾何形狀進(jìn)行近似[16]，如圖2所示。

圖2 梁變形狀態(tài)示意圖Fig. 2 Schematic of the beam in deformed states

幾何精確梁理論的非線性運(yùn)動(dòng)控制方程如下：

式中h和g為 線動(dòng)量和角動(dòng)量；t為時(shí)間；F和M為截面上的力和力矩；u為參考軸上的一維位移；x0為沿參考線點(diǎn)的初始位置向量；f和m為施加在梁上的分布力和分布力矩； (·)′表示對(duì)s求導(dǎo)。

基于小應(yīng)變假設(shè)，動(dòng)量與速度、應(yīng)變與截面力之間的關(guān)系為：

式中M1和K為截面質(zhì)量和剛度矩陣；ε和 κ為一維應(yīng)變和一維曲率；v和 ω為線速度和角速度。

一維應(yīng)變與曲率定義為：

式中R表 示旋轉(zhuǎn)張量；R0表 示初始旋轉(zhuǎn)張量；k表示截面的曲率向量；l1表示沿s方向的單位向量。

梁的非線性控制運(yùn)動(dòng)方程通過(guò)Newton-Raphson方法來(lái)求解，并用Lengendre譜有限元方法對(duì)增量方程離散化。線性分析中假設(shè)位移無(wú)窮小，因而物體的位形保持不變，但在大變形下必須建立參考位形為已知位形的方程。想推導(dǎo)由線性化得出的近似解的控制方程，可將應(yīng)力應(yīng)變參照于已知位形之上，首先需要線性化處理。非線性運(yùn)動(dòng)控制方程的線性化形式如下所示：

幾何精確梁理論的時(shí)間積分采用廣義α?xí)r間積分器來(lái)計(jì)算。由非線性控制運(yùn)動(dòng)方程定義的廣義α?xí)r間積分系統(tǒng)需要非線性系統(tǒng)在每個(gè)時(shí)間步長(zhǎng)進(jìn)行求解，相鄰兩個(gè)步數(shù)之間的差別在一定范圍內(nèi)時(shí)判定為收斂。幾何精確梁理論應(yīng)用的是能量停止準(zhǔn)則，用每一次迭代時(shí)的內(nèi)能增量與初始內(nèi)能增量相比來(lái)判斷是否收斂，該準(zhǔn)則提供了當(dāng)位移和力接近其平衡值時(shí)的度量：

式中： |·|表 示絕對(duì)值； ?U為位移向量的增量；R為外部施加的節(jié)點(diǎn)載荷向量；F為對(duì)應(yīng)于內(nèi)部單元應(yīng)力的節(jié)點(diǎn)力向量；εE為預(yù)設(shè)的能量容差。變量左側(cè)的上標(biāo)表示時(shí)間值，說(shuō)明處于動(dòng)態(tài)分析中，右側(cè)上標(biāo)表示Newton-Raphson迭代次數(shù)。

幾何精確梁理論是基于由鐵木辛柯梁發(fā)展而來(lái)的梁理論，也采用了平截面假設(shè)，但二者對(duì)于截面轉(zhuǎn)動(dòng)的處理方式不同[17]。一般的線性梁模型通過(guò)小變形假設(shè)忽略了方程中與截面轉(zhuǎn)動(dòng)相關(guān)的正、余弦項(xiàng)和高次冪項(xiàng)來(lái)簡(jiǎn)化為線性方程，而幾何精確梁理論考慮了截面轉(zhuǎn)動(dòng)和扭轉(zhuǎn)自由度，其中的三維旋轉(zhuǎn)可表示為Wiener-Milenkovic參數(shù)，用以下方程來(lái)表示：

式中，φ為旋轉(zhuǎn)角，n為旋轉(zhuǎn)軸的單位向量。

將風(fēng)力機(jī)葉片簡(jiǎn)化為梁?jiǎn)卧?，輸入相?yīng)結(jié)構(gòu)參數(shù)后，對(duì)節(jié)點(diǎn)自由度通過(guò)Legendre譜有限元進(jìn)行數(shù)值實(shí)現(xiàn)，采用梯形求積法 ，用單個(gè)單元對(duì)風(fēng)力機(jī)葉片進(jìn)行建模，用Wiener-Milenkovic參數(shù)表示三維旋轉(zhuǎn)，得到葉片各個(gè)節(jié)點(diǎn)的線位移、角位移，其次對(duì)非線性運(yùn)動(dòng)控制方程采用Newton-Raphson求解，線性化處理后，采用廣義α?xí)r間積分器判斷方程是否收斂。

1.3 氣動(dòng)彈性響應(yīng)計(jì)算

輸入風(fēng)模型后，在一個(gè)時(shí)間步長(zhǎng)內(nèi)，BEM計(jì)算葉片的氣動(dòng)載荷，然后計(jì)算葉片氣動(dòng)力、重力、離心力合力。通過(guò)輸入的結(jié)構(gòu)數(shù)據(jù)構(gòu)建質(zhì)量、剛度矩陣，建立葉片動(dòng)力學(xué)方程，通過(guò)Newton-Raphson方法求解葉片響應(yīng)，將葉片當(dāng)前狀態(tài)反饋到氣動(dòng)模型中，根據(jù)葉片變形后當(dāng)?shù)氐膩?lái)流條件和攻角確定升、阻力系數(shù)，葉片氣動(dòng)外形更新后進(jìn)行下一個(gè)時(shí)間步長(zhǎng)的計(jì)算。氣動(dòng)彈性響應(yīng)計(jì)算流程圖見(jiàn)圖3。

圖3 氣彈響應(yīng)計(jì)算流程圖Fig. 3 Flow chart of the aeroelastic response calculation

2 計(jì)算結(jié)果與分析

2.1 懸臂梁驗(yàn)證

通過(guò)與文獻(xiàn)[18]中理論值進(jìn)行對(duì)比，來(lái)驗(yàn)證本文幾何精確梁計(jì)算模型的可靠性與準(zhǔn)確性。文獻(xiàn)[18]中采用的是一個(gè)帶有預(yù)彎的懸臂梁。梁的橫截面是正方形，彈性模量為68.9 GPa。懸臂梁變形時(shí)的狀態(tài)如圖4所示，預(yù)彎梁與x軸夾角為45°，向葉尖分別施加1335 N和2670 N的力。

圖4 預(yù)彎梁的變形示意圖Fig. 4 Schematic of the curved beam in deformed states

表1和表2分別是施加1335 N和2670 N的力時(shí)與文獻(xiàn)中的葉尖位移結(jié)果對(duì)比。從表中可以看出，本文采用的幾何精確梁算法與文獻(xiàn)[18]中的解析解吻合良好，說(shuō)明幾何精確梁理論對(duì)于帶有預(yù)彎的梁位移的預(yù)測(cè)有較高的精度。對(duì)于風(fēng)力機(jī)葉片來(lái)說(shuō)，文獻(xiàn)[19]中已經(jīng)證明了幾何精確梁模型的結(jié)果與2.3 MW風(fēng)力機(jī)葉片的實(shí)驗(yàn)值更加吻合。

表1 施加1 335 N時(shí)葉尖位移對(duì)比Table 1 Comparison of the blade tip displacements under a force of 1 335 N

表2 施加2 670 N時(shí)葉尖位移對(duì)比Table 2 Comparison of the blade tip displacements under a force of 2 670 N

2.2 穩(wěn)態(tài)風(fēng)工況

本文研究大型水平軸風(fēng)力機(jī)柔性葉片，選用NREL 5 MW[20]和IEA 15 MW[21]風(fēng)力機(jī)，風(fēng)力機(jī)的葉片主要參數(shù)見(jiàn)表3和表4。

表3 NREL 5 MW葉片主要參數(shù)Table 3 Key parameters of the NREL 5 MW blade

表4 IEA 15 MW葉片主要參數(shù)Table 4 Key parameters of the IEA 15 MW blade

分別采用線性模態(tài)疊加法和幾何精確梁方法，對(duì)均勻來(lái)流條件下的5 MW和15 MW風(fēng)力機(jī)組的功率和推力模擬進(jìn)行對(duì)比分析，在單個(gè)風(fēng)速下的功率值和推力值均采用風(fēng)力機(jī)穩(wěn)定運(yùn)行周期的平均值。圖5～圖12中，Linear代表模態(tài)疊加法，Nonlinear代表幾何精確梁方法。

圖5給出了5 MW和15 MW機(jī)組在不同風(fēng)速下的功率曲線。從圖5（a）中可以看出，對(duì)于5 MW風(fēng)力機(jī)，額定風(fēng)速在11.4 m/s附近。對(duì)于功率來(lái)說(shuō)，在低于額定風(fēng)速下，非線性結(jié)果均略小于線性結(jié)果，但總體差距不明顯。從圖5（b）中可以看出，對(duì)于15 MW風(fēng)力機(jī)組，額定風(fēng)速在10.56 m/s附近；非線性方法下的額定風(fēng)速在11.5 m/s附近，并且對(duì)于功率來(lái)說(shuō)，低于額定風(fēng)速下的結(jié)果也均小于線性結(jié)果，但其差值比5 MW機(jī)組的要大，分析認(rèn)為這是由于葉片的扭轉(zhuǎn)變形可能會(huì)造成輸出功率的降低，因?yàn)閹в信まD(zhuǎn)變形的模型具有較低的載荷并且發(fā)電量也較低。越接近11.5 m/s風(fēng)速，兩種方法的功率差值越大，這說(shuō)明，模態(tài)疊加法和幾何精確梁理論兩種方法對(duì)低風(fēng)速下小變形葉片的風(fēng)力機(jī)發(fā)電功率的預(yù)測(cè)無(wú)明顯區(qū)別，但在中、高風(fēng)速下葉片發(fā)生了較大變形，尤其對(duì)于15 MW這種大柔性葉片，幾何精確梁理論與模態(tài)疊加法相比，在葉片結(jié)構(gòu)分析中考慮了葉片的扭轉(zhuǎn)變形和彎扭耦合效應(yīng)，因此更適合于求解帶有幾何非線性效應(yīng)的葉片。

圖5 5 MW和15 MW機(jī)組在不同風(fēng)速下的功率曲線Fig. 5 Power performance curves of the 5 MW and 15 MW wind turbines under different wind speeds

圖6給出了5 MW和15 MW機(jī)組在不同風(fēng)速下的風(fēng)輪推力。從圖中可以看出：兩個(gè)機(jī)組隨著風(fēng)速增加，風(fēng)輪推力逐漸增大，均在額定風(fēng)速下達(dá)到最大推力；對(duì)于5 MW機(jī)組，非線性方法下的風(fēng)輪推力要略小于線性結(jié)果，差距較小，總的來(lái)說(shuō)一致性較好；對(duì)于15 MW機(jī)組，最大推力出現(xiàn)在風(fēng)速11.5 m/s附近，與功率曲線完全對(duì)應(yīng)。在額定風(fēng)速附近，計(jì)算結(jié)果有明顯差異，非線性結(jié)果要小于線性計(jì)算結(jié)果。額定風(fēng)速之后，線性與非線性結(jié)果一致性較好。出現(xiàn)這種差異可能是由于15 MW的117 m葉片大幅度的變形使風(fēng)力機(jī)風(fēng)輪實(shí)際的掃掠面積變小，翼型的氣動(dòng)性能偏離了最優(yōu)狀態(tài)。

圖6 5 MW和15 MW機(jī)組在不同風(fēng)速下的風(fēng)輪推力Fig. 6 Rotor thrust of the 5 MW and 15 MW wind turbines under different wind speeds

葉尖揮舞變形量如圖7所示。對(duì)于5 MW機(jī)組，葉尖的最大揮舞變形發(fā)生在額定風(fēng)速下，此時(shí)風(fēng)力機(jī)剛達(dá)到滿發(fā)，在額定風(fēng)速附近線性與非線性結(jié)果差值為0.246 m。在低風(fēng)速范圍內(nèi)，線性與非線性結(jié)果吻合良好，而在中高風(fēng)速下，非線性結(jié)果要略小于線性結(jié)果，并且隨著風(fēng)速增大差值也越大，出現(xiàn)這種現(xiàn)象的原因是實(shí)際變形過(guò)程中，由于葉片產(chǎn)生的位移而導(dǎo)致長(zhǎng)度減小，從而發(fā)生了彎曲方向的剛度強(qiáng)化，可見(jiàn)非線性分析方法更加貼合實(shí)際情況。而15 MW機(jī)組的變化趨勢(shì)相同，線性結(jié)果與非線性結(jié)果在額定風(fēng)速附近差值最大為3.4 m，非線性結(jié)果在低風(fēng)速內(nèi)與線性結(jié)果吻合較好，而在中高風(fēng)速下，非線性結(jié)果與線性結(jié)果差異明顯，說(shuō)明在額定風(fēng)速附近和高風(fēng)速下，大功率機(jī)組的葉片顯示出了強(qiáng)非線性，普通的模態(tài)疊加法已經(jīng)失效，這時(shí)對(duì)幾何非線性的考慮尤為重要。

圖7 5 MW和15 MW機(jī)組在不同風(fēng)速下的葉尖揮舞位移Fig. 7 Flapwise tip deflection of the 5 MW and 15 MW wind turbines under different wind speeds

圖8 給出了5 MW和15 MW機(jī)組在不同風(fēng)速下的葉根揮舞彎矩。由圖可見(jiàn)，葉根揮舞彎矩在額定風(fēng)速附近達(dá)到最大值。對(duì)于5 MW機(jī)組，線性與非線性結(jié)果吻合良好，在額定風(fēng)速附近非線性結(jié)果略小于線性結(jié)果，差值為0.15 MN；對(duì)于15 MW機(jī)組，非線性結(jié)果在低于額定風(fēng)速范圍內(nèi)要小于線性結(jié)果，差值較大，并在額定風(fēng)速附近差值達(dá)到最大值13.1 MN，其余風(fēng)速下吻合較好。從5 MW到15 MW，數(shù)值偏差增加了21.23%。葉根的揮舞力矩主要由氣動(dòng)載荷決定，而葉片的扭轉(zhuǎn)變形與攻角緊密相關(guān)，進(jìn)而影響了氣動(dòng)載荷，基于歐拉伯努利梁理論的線性模態(tài)疊加法僅考慮了葉片的彎曲自由度，對(duì)于61.5 m葉片的葉根彎矩具有較高的計(jì)算精度，但是針對(duì)大功率級(jí)風(fēng)力機(jī)組的葉片，沒(méi)有考慮扭轉(zhuǎn)自由度就會(huì)導(dǎo)致在額定風(fēng)速下產(chǎn)生較大差異。

圖8 5 MW和15 MW機(jī)組在不同風(fēng)速下的葉根揮舞彎矩Fig. 8 Flapwise root moment of the 5 MW and 15 MW wind turbines under different wind speeds

圖9給出了5 MW和15 MW機(jī)組在不同風(fēng)速下的葉根擺振彎矩。從圖中可以看出，對(duì)于5 MW機(jī)組，葉根擺振力矩在額定風(fēng)速下達(dá)到了最大值，非線性結(jié)果在高風(fēng)速下要大于線性結(jié)果，差值隨著風(fēng)速增大而增大；對(duì)于15 MW機(jī)組，非線性結(jié)果在額定風(fēng)速到切出風(fēng)速范圍內(nèi)整體明顯大于線性結(jié)果。葉根擺振彎矩主要由從整體坐標(biāo)系到葉片局部坐標(biāo)系的重力分量決定。葉片槳距角和扭轉(zhuǎn)變形的大小直接影響了這兩個(gè)坐標(biāo)系中坐標(biāo)的轉(zhuǎn)換和重力分量的大小，達(dá)到額定風(fēng)速葉片變槳后兩種方法得到的槳距角差異較大，所以影響了重力分量的大小，進(jìn)而影響了中高風(fēng)速下葉根擺振彎矩的大小。

圖9 5 MW和15 MW機(jī)組在不同風(fēng)速下的葉根擺振彎矩Fig. 9 Edgewise root moment of the 5 MW and 15 MW wind turbines under different wind speeds

兩個(gè)機(jī)組的響應(yīng)在線性方法與非線性方法的最大差值如表5所示，ΔP、ΔT、ΔTx、ΔMx、ΔMy分別為功率、推力、葉尖揮舞位移、葉根擺振彎矩、葉根揮舞彎矩在不同風(fēng)速下的最大差值。對(duì)于功率、載荷和位移來(lái)說(shuō)，15 MW機(jī)組兩種方法下的最大差值要比5 MW機(jī)組明顯偏大，并且對(duì)于葉根擺振彎矩來(lái)說(shuō)，非線性計(jì)算結(jié)果更大。

表5 兩個(gè)機(jī)組的響應(yīng)在兩種方法下的最大差值Table 5 Maximum difference between the response of the two wind turbines under two methods

2.3 湍流風(fēng)工況

湍流風(fēng)場(chǎng)的建立基于IEC標(biāo)準(zhǔn)中給出的Kaimal模型，根據(jù)IEC 61400-1，對(duì)于中性和穩(wěn)定大氣，功率譜模型定義如下：

式中：f為頻率；k為速度分量方向的指數(shù)；Sk為單側(cè)速度分量譜； σk為速度分量標(biāo)準(zhǔn)差；Lk為速度分量積分標(biāo)度參數(shù)；Vhub為輪轂高度處的平均風(fēng)速。

湍流風(fēng)的設(shè)定采用正常湍流風(fēng)（NTM）?；诮o出的Kaimal速度譜模型，由頻域的速度分布來(lái)產(chǎn)生時(shí)域的風(fēng)速數(shù)據(jù)，在二維矩形區(qū)域建立一定湍流強(qiáng)度的脈動(dòng)風(fēng)速場(chǎng)，并作為氣動(dòng)模型的來(lái)流風(fēng)輸入，進(jìn)行湍流風(fēng)況的風(fēng)力機(jī)氣動(dòng)彈性仿真。表6對(duì)輪轂高度處風(fēng)速為10 m/s的湍流風(fēng)場(chǎng)進(jìn)行了描述，生成了湍流強(qiáng)度為18.34%的剪切湍流風(fēng)，風(fēng)速數(shù)據(jù)分布在二維平面961（31×31）個(gè)網(wǎng)格點(diǎn)上，網(wǎng)格平面高度和寬度均為260 m，可覆蓋整個(gè)風(fēng)輪及部分塔架。

表6 湍流風(fēng)場(chǎng)參數(shù)Table 6 Parameters of the turbulent wind field

分別采用線性分析與非線性分析方法，計(jì)算了湍流風(fēng)速10 m/s時(shí)兩個(gè)風(fēng)力機(jī)組的動(dòng)態(tài)響應(yīng)，并進(jìn)行了對(duì)比分析。有效模擬時(shí)長(zhǎng)為200 s，選取后100 s運(yùn)行時(shí)段的數(shù)據(jù)進(jìn)行分析，并將位移和載荷的時(shí)域變化進(jìn)行快速傅立葉變換，得到其頻域特征。

圖10（a）～圖10（d）分別描述了兩個(gè)機(jī)組的葉片在風(fēng)速10 m/s、湍流度18.34%條件下的葉尖揮舞位移在時(shí)域和頻域的情況。對(duì)于5 MW機(jī)組而言，線性與非線性結(jié)果差別并不明顯，圖10（c）中，在頻率為0.18 Hz、0.36 Hz處存在一些峰值，并且為該風(fēng)況下頻率0.183 Hz時(shí)的值的倍數(shù)。根據(jù)表3和表4，0.74 Hz處對(duì)應(yīng)的是一階揮舞模態(tài)，在高頻區(qū)非線性響應(yīng)運(yùn)動(dòng)能量要略大于線性響應(yīng)。而對(duì)于15 MW機(jī)組，線性與非線性結(jié)果之間的差別要明顯大于5 MW機(jī)組的，與穩(wěn)態(tài)風(fēng)況中額定風(fēng)速附近下的結(jié)果較吻合。其對(duì)應(yīng)轉(zhuǎn)頻為0.118 Hz，在0.119 Hz處峰值對(duì)應(yīng)此轉(zhuǎn)頻，0.579 Hz處對(duì)應(yīng)的是一階揮舞模態(tài)，并且在低頻區(qū)線性響應(yīng)能量更高，在高頻區(qū)非線性響應(yīng)能量更高。

圖10 10 m/s湍流風(fēng)條件下的葉尖揮舞位移Fig. 10 Flapwise tip deflection of the wind turbines under a turbulent inflow of 10 m/s

圖11（a）～圖11（d）描述了葉根揮舞彎矩在時(shí)域和頻域的變化情況，其隨時(shí)間的變化趨勢(shì)與葉尖揮舞位移的非常相似。5 MW機(jī)組的線性結(jié)果與非線性結(jié)果差別同樣不明顯，同樣在0.74 Hz處對(duì)應(yīng)的是一階揮舞模態(tài)，在高頻區(qū)非線性響應(yīng)略大于線性響應(yīng)。對(duì)于15 MW機(jī)組，無(wú)論是均值還是振幅，非線性結(jié)果都要小于線性結(jié)果，差別偏大。在0.579 Hz處對(duì)應(yīng)的是一階揮舞模態(tài)，在低頻區(qū)同樣線性響應(yīng)能量更高，高頻區(qū)葉根揮舞彎矩的線性與非線性值吻合較好。

圖11 10 m/s湍流風(fēng)下的葉根揮舞彎矩Fig. 11 Flapwise root moment of the wind turbines under a turbulent inflow of 10 m/s

從圖12（a）～圖12（d）看出，葉根擺振彎矩隨時(shí)間的變化呈明顯的周期性，基本圍繞均值正負(fù)波動(dòng)。對(duì)于5 MW和15 MW風(fēng)力機(jī)，線性與非線性時(shí)域結(jié)果均差別較小。在頻域分析中，5 MW機(jī)組在1.09 Hz處對(duì)應(yīng)的是一階擺振模態(tài)，并在兩種方法下捕捉到的峰值位置在高頻區(qū)有明顯差異，線性方法在對(duì)5 MW機(jī)組的峰值預(yù)測(cè)上在高頻區(qū)延遲了42.6%，根據(jù)文獻(xiàn)[22]，線性方法在高頻區(qū)捕捉到的峰值位置有明顯延遲，與本文結(jié)果相符。而高頻區(qū)非線性響應(yīng)的湍流能量大于線性響應(yīng)；線性方法預(yù)測(cè)的一階擺振模態(tài)頻率在15 MW機(jī)組上高估了3.2%，在低頻區(qū)非線性響應(yīng)的湍流能量要大于線性響應(yīng)的，高頻區(qū)線性響應(yīng)大于非線性響應(yīng)。

圖12 10 m/s湍流風(fēng)下的葉根擺振彎矩Fig. 12 Edgewise root moment of wind turbines under a turbulent inflow of 10 m/s

3 結(jié)論

本文計(jì)算了懸臂梁的純彎曲位移，并與文獻(xiàn)中的理論值進(jìn)行對(duì)比，驗(yàn)證了幾何精確梁理論求解帶有幾何非線性效應(yīng)的梁的精度。然后計(jì)算了NREL 5 MW與IEA 15 MW風(fēng)力機(jī)在線性與非線性分析下的動(dòng)態(tài)響應(yīng)，具體研究結(jié)論如下：

1）幾何精確梁理論考慮了彎扭耦合下葉片產(chǎn)生的額外扭轉(zhuǎn)，減小的槳距角相當(dāng)于彌補(bǔ)了額外的扭轉(zhuǎn)，從而影響了控制器的設(shè)定點(diǎn)。

2）對(duì)于揮舞方向的位移與載荷，15 MW風(fēng)力機(jī)線性與非線性計(jì)算結(jié)果在中高風(fēng)速下差別明顯，葉尖位移最大相差了22.5%，而對(duì)于擺振方向的載荷在過(guò)高風(fēng)速下才出現(xiàn)較為明顯的差別。

3）在考慮幾何非線性后，15 MW風(fēng)力機(jī)葉尖揮舞位移和葉根揮舞彎矩最大減小了22.5%和23.1%。因此對(duì)于百米級(jí)的大柔性葉片，在保證結(jié)構(gòu)安全的同時(shí)，可以進(jìn)一步降低葉片質(zhì)量。

4）對(duì)于揮舞方向的位移和載荷，與非線性方法相比，線性方法高估了低頻區(qū)的響應(yīng)，低估了高頻區(qū)的響應(yīng)。

5）對(duì)于葉根擺振彎矩，線性方法在對(duì)5 MW機(jī)組的峰值預(yù)測(cè)上在高頻區(qū)大幅度延遲，且15 MW機(jī)組中非線性方法下的頻率更加貼近一階擺振模態(tài)頻率。

綜上，對(duì)于5 MW機(jī)組等剛度較大的風(fēng)力機(jī)葉片而言，基于歐拉梁理論的模態(tài)疊加法仍然適用，其葉片質(zhì)心偏移量很小，不會(huì)引起響應(yīng)發(fā)生較大變化，且風(fēng)激發(fā)的最低模態(tài)主要是彎曲模態(tài)，位移較小，可以被經(jīng)典梁理論中的線性項(xiàng)準(zhǔn)確捕捉；但對(duì)于15 MW等大功率機(jī)組的百米級(jí)葉片，葉片柔性大、變形大，并具有彎扭耦合特征，模態(tài)疊加法已不再適用，且功率的變化可能需要優(yōu)化新的控制策略來(lái)改善風(fēng)力機(jī)氣動(dòng)性能。相對(duì)于經(jīng)典梁理論，本文基于幾何精確梁理論所建立的非線性葉片氣動(dòng)彈性響應(yīng)計(jì)算方法精度更高，更適用于各向異性復(fù)合材料的百米級(jí)葉片的數(shù)值求解。隨著風(fēng)力機(jī)葉片柔性化發(fā)展，對(duì)于更大功率的風(fēng)力機(jī)和更大的柔性葉片，應(yīng)考慮幾何非線性對(duì)風(fēng)力機(jī)葉片氣動(dòng)彈性響應(yīng)的影響，以準(zhǔn)確評(píng)估風(fēng)力機(jī)設(shè)計(jì)的安全性和穩(wěn)定性。