百米級大柔性風電葉片非線性氣彈響應分析

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

(1. 南京航空航天大學 江蘇省風力機設計高技術研究重點實驗室，南京 210016；2. 南京航空航天大學 民航學院，南京 210016)

0 引言

隨著風力機復合材料葉片尺寸的增加，幾何非線性效應、截面面內和面外翹曲等非經典的效應，對葉片結構動力學動態響應產生顯著影響。然而由于彈性耦合影響，復合材料結構的分析要比各向同性結構的更難[1]。傳統的模態疊加法仍然適用于計算小功率風力機的彈性變形。然而對于大型風力機而言，在復雜工況下可能會發生大的變形，在這時，傳統的線性方法無法再準確地預測葉片的動態響應。

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

本文主要將幾何非線性分析方法應用在兩種不同功率的大型風力機中的大柔性葉片，首先采用幾何精確梁理論計算懸臂梁變形，并與理論值進行對比，驗證結構計算方法的可靠性與準確性。然后通過葉素動量理論與幾何精確梁方法耦合計算，實現大型風力機5 MW和15 MW的動態響應計算。

1 理論描述

1.1 葉素動量理論

葉素動量方法實際上是葉素理論和動量理論的結合。氣流在流管內流動滿足動量定理，但是氣流在流經葉片時會受到擾動，從而導致了氣流的切向和軸向速度發生變化，通常引入切向和軸向誘導因子來反應氣流在通過風輪平面時速度的損失。葉素理論將葉片沿展向離散為有限數量的葉素，葉素隨著風輪旋轉形成了一個圓環，沿展向對升、阻力積分便可求得氣動推力和扭矩。葉片一個葉素微元的受力如圖1所示。

圖1 葉素上的來流速度與氣動力示意圖Fig. 1 Inflow velocity and aerodynamic forces on the blade

W為入流合速度，速度合成關系為：

式中 ?為入流角，V0為入流速度， ?為風輪轉速，a為軸向誘導因子，a′為切向誘導因子。

葉素升力、阻力表達式為：

式中：ρ為空氣密度，c為二維翼型弦長，CL、CD為二維翼型升力、阻力系數。

通過將葉段升、阻力分解到風輪旋轉平面及垂直于風輪旋轉平面的方向，得到了切向力Ftq和法向力Fth：

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

考慮葉片數量B，微元上受到的推力和轉矩通過葉素理論描述為：

根據動量理論和葉素理論下的推力和轉矩公式相等，可得軸向誘導因子a和 切向誘導因子a′的表達式：

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

1.2 幾何精確梁理論

幾何精確梁理論以受初始彎扭的梁能承受大的位移和旋轉的能力為特點。通過一個三維橫截面分析，六個自由度的所有耦合效應，包括拉伸、剪切、扭轉、彎曲、扭轉翹曲、面內翹曲等，都能被幾何精確梁理論涵蓋。“幾何精確”指的是在公式中沒有對初始幾何形狀和變形后的幾何形狀進行近似[16]，如圖2所示。

圖2 梁變形狀態示意圖Fig. 2 Schematic of the beam in deformed states

幾何精確梁理論的非線性運動控制方程如下：

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

基于小應變假設，動量與速度、應變與截面力之間的關系為：

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

一維應變與曲率定義為：

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

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

幾何精確梁理論的時間積分采用廣義α時間積分器來計算。由非線性控制運動方程定義的廣義α時間積分系統需要非線性系統在每個時間步長進行求解，相鄰兩個步數之間的差別在一定范圍內時判定為收斂。幾何精確梁理論應用的是能量停止準則，用每一次迭代時的內能增量與初始內能增量相比來判斷是否收斂，該準則提供了當位移和力接近其平衡值時的度量：

式中： |·|表 示絕對值； ?U為位移向量的增量；R為外部施加的節點載荷向量；F為對應于內部單元應力的節點力向量；εE為預設的能量容差。變量左側的上標表示時間值，說明處于動態分析中，右側上標表示Newton-Raphson迭代次數。

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

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

將風力機葉片簡化為梁單元，輸入相應結構參數后，對節點自由度通過Legendre譜有限元進行數值實現，采用梯形求積法 ，用單個單元對風力機葉片進行建模，用Wiener-Milenkovic參數表示三維旋轉，得到葉片各個節點的線位移、角位移，其次對非線性運動控制方程采用Newton-Raphson求解，線性化處理后，采用廣義α時間積分器判斷方程是否收斂。

1.3 氣動彈性響應計算

輸入風模型后，在一個時間步長內，BEM計算葉片的氣動載荷，然后計算葉片氣動力、重力、離心力合力。通過輸入的結構數據構建質量、剛度矩陣，建立葉片動力學方程，通過Newton-Raphson方法求解葉片響應，將葉片當前狀態反饋到氣動模型中，根據葉片變形后當地的來流條件和攻角確定升、阻力系數，葉片氣動外形更新后進行下一個時間步長的計算。氣動彈性響應計算流程圖見圖3。

圖3 氣彈響應計算流程圖Fig. 3 Flow chart of the aeroelastic response calculation

2 計算結果與分析

2.1 懸臂梁驗證

通過與文獻[18]中理論值進行對比，來驗證本文幾何精確梁計算模型的可靠性與準確性。文獻[18]中采用的是一個帶有預彎的懸臂梁。梁的橫截面是正方形，彈性模量為68.9 GPa。懸臂梁變形時的狀態如圖4所示，預彎梁與x軸夾角為45°，向葉尖分別施加1335 N和2670 N的力。

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

表1和表2分別是施加1335 N和2670 N的力時與文獻中的葉尖位移結果對比。從表中可以看出，本文采用的幾何精確梁算法與文獻[18]中的解析解吻合良好，說明幾何精確梁理論對于帶有預彎的梁位移的預測有較高的精度。對于風力機葉片來說，文獻[19]中已經證明了幾何精確梁模型的結果與2.3 MW風力機葉片的實驗值更加吻合。

表1 施加1 335 N時葉尖位移對比Table 1 Comparison of the blade tip displacements under a force of 1 335 N

表2 施加2 670 N時葉尖位移對比Table 2 Comparison of the blade tip displacements under a force of 2 670 N

2.2 穩態風工況

本文研究大型水平軸風力機柔性葉片，選用NREL 5 MW[20]和IEA 15 MW[21]風力機，風力機的葉片主要參數見表3和表4。

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

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

分別采用線性模態疊加法和幾何精確梁方法，對均勻來流條件下的5 MW和15 MW風力機組的功率和推力模擬進行對比分析，在單個風速下的功率值和推力值均采用風力機穩定運行周期的平均值。圖5～圖12中，Linear代表模態疊加法，Nonlinear代表幾何精確梁方法。

圖5給出了5 MW和15 MW機組在不同風速下的功率曲線。從圖5（a）中可以看出，對于5 MW風力機，額定風速在11.4 m/s附近。對于功率來說，在低于額定風速下，非線性結果均略小于線性結果，但總體差距不明顯。從圖5（b）中可以看出，對于15 MW風力機組，額定風速在10.56 m/s附近；非線性方法下的額定風速在11.5 m/s附近，并且對于功率來說，低于額定風速下的結果也均小于線性結果，但其差值比5 MW機組的要大，分析認為這是由于葉片的扭轉變形可能會造成輸出功率的降低，因為帶有扭轉變形的模型具有較低的載荷并且發電量也較低。越接近11.5 m/s風速，兩種方法的功率差值越大，這說明，模態疊加法和幾何精確梁理論兩種方法對低風速下小變形葉片的風力機發電功率的預測無明顯區別，但在中、高風速下葉片發生了較大變形，尤其對于15 MW這種大柔性葉片，幾何精確梁理論與模態疊加法相比，在葉片結構分析中考慮了葉片的扭轉變形和彎扭耦合效應，因此更適合于求解帶有幾何非線性效應的葉片。

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

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

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

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

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

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

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

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

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

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

表5 兩個機組的響應在兩種方法下的最大差值Table 5 Maximum difference between the response of the two wind turbines under two methods

2.3 湍流風工況

湍流風場的建立基于IEC標準中給出的Kaimal模型，根據IEC 61400-1，對于中性和穩定大氣，功率譜模型定義如下：

式中：f為頻率；k為速度分量方向的指數；Sk為單側速度分量譜； σk為速度分量標準差；Lk為速度分量積分標度參數；Vhub為輪轂高度處的平均風速。

湍流風的設定采用正常湍流風（NTM）。基于給出的Kaimal速度譜模型，由頻域的速度分布來產生時域的風速數據，在二維矩形區域建立一定湍流強度的脈動風速場，并作為氣動模型的來流風輸入，進行湍流風況的風力機氣動彈性仿真。表6對輪轂高度處風速為10 m/s的湍流風場進行了描述，生成了湍流強度為18.34%的剪切湍流風，風速數據分布在二維平面961（31×31）個網格點上，網格平面高度和寬度均為260 m，可覆蓋整個風輪及部分塔架。

表6 湍流風場參數Table 6 Parameters of the turbulent wind field

分別采用線性分析與非線性分析方法，計算了湍流風速10 m/s時兩個風力機組的動態響應，并進行了對比分析。有效模擬時長為200 s，選取后100 s運行時段的數據進行分析，并將位移和載荷的時域變化進行快速傅立葉變換，得到其頻域特征。

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

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

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

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

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

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

3 結論

本文計算了懸臂梁的純彎曲位移，并與文獻中的理論值進行對比，驗證了幾何精確梁理論求解帶有幾何非線性效應的梁的精度。然后計算了NREL 5 MW與IEA 15 MW風力機在線性與非線性分析下的動態響應，具體研究結論如下：

1）幾何精確梁理論考慮了彎扭耦合下葉片產生的額外扭轉，減小的槳距角相當于彌補了額外的扭轉，從而影響了控制器的設定點。

2）對于揮舞方向的位移與載荷，15 MW風力機線性與非線性計算結果在中高風速下差別明顯，葉尖位移最大相差了22.5%，而對于擺振方向的載荷在過高風速下才出現較為明顯的差別。

3）在考慮幾何非線性后，15 MW風力機葉尖揮舞位移和葉根揮舞彎矩最大減小了22.5%和23.1%。因此對于百米級的大柔性葉片，在保證結構安全的同時，可以進一步降低葉片質量。

4）對于揮舞方向的位移和載荷，與非線性方法相比，線性方法高估了低頻區的響應，低估了高頻區的響應。

5）對于葉根擺振彎矩，線性方法在對5 MW機組的峰值預測上在高頻區大幅度延遲，且15 MW機組中非線性方法下的頻率更加貼近一階擺振模態頻率。

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