大型風力機葉片凈空計算方法研究

董笑宇 曲春明 葛文澎 夏 冰 苗得勝

明陽智慧能源集團股份公司

0 引言

隨著風能技術的不斷發展，風電機組越來越大型化，葉輪直徑也越來越大，在現有加工工藝和材料技術未完成革命性突破前，葉片越來越注重輕量化設計，從而使葉片更加柔性化［1］。輕型化、柔性化的葉片在實際運行中發生幾何非線性大變形的幾率增大，會出現葉片與塔筒相撞的掃塔現象，嚴重時葉片會斷裂，造成巨大的經濟損失。為保障風電機組大型化對葉片結構強度的要求及風電機組的安全運行，對風機葉片的凈空研究十分重要。

葉片凈空計算方法很多，風力機載荷計算方面主要分為葉素動量理論BEM、計算流體力學CFD、和升力線方法等；葉片變形方面，按幾何結構可分為梁單元法、殼單元法和實體單元法，按方法原理可分為多體動力學MED 和有限元結構力學FEM。Liu 等人［2］利用計算流體力學和多體運動系統耦合方法計算了海上風力機組運行狀態下的動態結構響應，Yang 等人［3］研究了在旋轉輪轂上歐拉-伯努利梁的柔性運動，推導了一組描述梁軸向、橫向和旋轉運動全耦合的非線性積分微分方程。

本文比較了幾種葉片凈空計算方法，首先從風力機空氣動力計算方法進行比較，再從葉片結構響應方法比較，最后從工程實用角度出發給出推薦的計算方法。

1 葉片模型

本文以某三葉片海上機組為研究對象，選取方位角為180°的葉片（葉輪呈“Y”形時葉尖向下的葉片）為計算模型，機組及葉片參數見表1，葉片模型見圖1。

圖1 機組葉片模型

表1 機組及葉片參數

有研究表明，當機組處于以下兩種狀態時，葉片凈空距離最?。孩亠L輪達到額定轉速、環境風速接近額定風速、葉片尚未變槳；②機組運行存在極端負切變，葉片載荷過大。鑒于海上風切變較小，本文選取第一種狀態作為分析工況。

本研究葉片采用纖維增強型復合材料，根據葉片弦長、扭角、預彎、厚度、主梁寬度及鋪層情況，FOCUS軟件可輸出葉片各截面的揮舞剛度、擺振剛度、軸向剛度和扭轉剛度，同時輸出用于結構計算的截面質量屬性。本文采用已知的氣動外形參數和FOCUS輸出結構參數作為葉片凈空計算的數據來源。

2 載荷計算方法

風力機氣動性能及載荷計算的方法有計算流體力學CFD、葉素動量理論BEM、升力線、升力面、面元法及區域求解法等［4-5］，其中BEM 理論和CFD 計算方法較普遍，以下為兩種方法的對比和結果分析。

2.1 葉素動量理論BEM

葉素動量理論BEM 由英國科學家Glauert 提出，將動量理論與葉素理論相結合求解軸向誘導因子a和周向誘導因子a′，進而計算風輪葉片上的力和力矩。葉素動量理論建立后，一直用于風力機氣動特性的研究，正確性已得到充分驗證［6］，該理論計算步驟如下［7-8］：

（1）a、a′賦初值，作為迭代第一步；

（2）根據公式（1）計算入流角：

其中：Ω 為葉輪轉速，r為流體微團所在的葉輪展向位置，u為流體沿流線方向的運動速度，下標∞為流體位于上游無窮遠處；

（3）利用公式（2）計算攻角：

θ為局部扭角和局部變槳角之和；

（4）根據翼型氣動特性曲線得到葉素的升力系數CL和阻力系數CD；

（5）通過軸向和周向誘導因子公式（3）、（4）計算a和a′：

其中，B為葉片數，F為考慮葉根和葉尖流動損失的總損失系數，H為推力系數的修正；

（6）根據更新后的a和a′判斷迭代過程的收斂性；

（7）由此得到葉輪的推力系數和功率系數，見公式（5）（6）：

ρ為微團密度，T為葉輪上的推力，M為葉輪上的轉矩。

2.2 計算流體力學CFD

計算流體力學CFD 就是通過計算機和數值方法求解流體力學控制方程，獲取各條件下流動的數據和作用在繞流物體上的力、力矩、流動圖像和熱量等。

流動的數值計算一般有三個方法：直接數值模擬法DNS，大渦模擬法LES 和雷諾平均方程RANS，但解決流體動力學及空氣動力學的流動和工程實際問題的主要是雷諾平均法。雷諾平均法的基本思想是將Navier-Stokes 方程時均化，用時均值來代表瞬時值。該方法可節省大量時間與資源，因此在工程中被廣泛應用［9］。

2.3 方法驗證

BEM 方法在無黏假設下的一維動量理論與實際三維黏性流動之間存在差異，為此，CFD 設置了整機、風輪、單葉片三種計算模型，圖2、圖3 為葉片方位角為180°時BEM和CFD載荷計算結果。

圖2 葉片推力

圖3 葉片力矩

由圖2 可知，兩種計算方法得到的沿葉展方向的推力整體變化趨勢一致，均先增大后急劇減小，但單葉片模型推力分布區別較大，BEM 在葉片中段位置較大，在葉尖位置較小，且推力沿葉展方向更光滑，即推力值沿葉展方向增長且增長率較平穩均勻，CFD 單葉片整體推力較BEM 存在-2.09%差值；由圖3可見，兩種計算方法得到的沿葉展方向的力矩整體變化趨勢一致，但CFD 葉根與葉尖位置比BEM 大，CFD 單葉片模型力矩較BEM 轉矩存在8.03%差值。

CFD 整機、三葉片風輪和單葉片模型沿葉片展向推力分布一致，單葉片推力最大，整機最小。與風輪模型相比，整機模型存在一定的塔影效應，使位于塔筒前的葉片所受推力減??；與單葉片模型相比，三葉片風輪模型葉片推力更小，這是因為當葉片以較高線速度旋轉時，每個葉片都會對軌跡后方葉片所處流場造成一定的影響，因此三個葉片氣動性能均存在一定程度的損失。

CFD 整機、三葉片風輪和單葉片模型力矩結果中，塔筒與葉片對氣動性能干擾作用更加明顯，尤其是葉片間流場的相互干擾，很大程度降低了葉片力矩。

圖4 為整機模型流場圖。由圖可見，由于塔筒的存在使葉片與塔筒之間流場方向變得紊亂且流速變慢。圖5、圖6 分別為葉片展向50%和80%位置翼型周向壓力分布圖，取單葉片模型與整機模型為對象，弦長方向位置0 處為前緣，1 處為后緣。由圖可見，單葉片與整機模型翼型周向壓力分布規律一致，壓力面前緣壓力最大，后緣稍低，中段為負壓；吸力面均為負壓，前緣壓力最小，由前向后逐漸升高。

圖5 葉片展向50%位置翼型周向壓力分布

圖6 葉片展向80%位置翼型周向壓力分布

此外，從圖5中還可看到，單葉片壓力面前緣壓力較整機大，吸力面則小，因此單葉片模型前緣壓差更大。該規律在圖6 展向80%位置處依然可得，即在更靠近葉尖處，單葉片模型前緣壓差更大，兩種模型在壓力面壓差最大達155 Pa，吸力面達-970 Pa，說明靠近葉尖位置受塔筒和其他葉片的干擾更明顯。單葉片模型前緣位置壓差更大說明單葉片比整機所受推力更大。

3 結構變形計算方法

現階段對風力機葉片振動及變形的研究大多采用梁模型或多體動力系統。本研究采用傳統有限元方法FEM 進行葉片靜力學分析：將葉片簡化為梁模型，通過調整葉片截面參數使梁模型獲得更準確的揮舞、擺振及拉伸方向剛度等參數，同時采用多體動力學方法開展葉片結構分析，即將葉片離散為多個彈簧與質量的組成。

3.1 FEM方法

FEM為將所求解連續區域的彈性體離散為有限個單元，各個單元通過相鄰節點連接，使復雜邊界連續體無限自由度問題簡化為單元節點有限自由度問題，從而得到整個區域的近似解，如果采用的單元滿足收斂要求，那么單元數越多越趨近于精確解。

風力機葉片具有細長結構的特點，在分析葉片結構時，葉片通常被假設為梁模型(1D-FEM)，在有限元中，梁單元模型的基本假設是原截面受彎曲力后依然是平面。

本文采用的線性梁單元B31 是考慮剪切變形的Timoshenko 梁單元［10］，受力前截面和中性面是垂直的，彎曲后不再垂直。該單元適用于模擬剪切變形起重要作用的深梁和模擬剪切變形不太重要的細長梁。

3.2 MBD方法

MBD 方法為將柔性葉片離散為若干個有一定質量與形狀的體，體與體之間用旋轉鉸、彈簧及阻尼器連接，組成彈簧質量多體系統，通過系統的位移表征葉片的運動和變形。

為準確模擬柔性葉片彈性變形及力學特性，根據最小勢能原理和伯努利梁假設確定彈簧剛度系數。多體系統動力學拉格朗日數學模型可通過虛功形式的動力學普遍方程建立，即多體系統在虛加慣性力、真實主動力與理想約束力的作用下處于靜平衡，聯立運動約束方程和受力平衡方程對靜力學進行計算［11-12］。

3.3 方法驗證

首先對梁模型與多體模型進行模態分析，驗證兩種模型質量剛度分布的準確性，模態分析結果見表2。

根據模態分析結果可看出，兩種模型的揮舞頻率較接近，梁模型模態二階揮舞頻率比多體動力系統模態結果誤差更大，相差0.9%。模態結果誤差的原因是模型沒有采用葉片的實際外形（以近似葉片模型的方法），故葉片整體質量與剛度等效存在一定的誤差，但與葉片實際測量值誤差均在2%以內，因此認為模型質量剛度分布符合精度要求。

為對比多體動力學方法與有限元方法計算葉片變形情況，采用BEM方法得到的載荷值對梁模型和多體系統進行靜力分析，葉尖變形結果見表3。

表3 葉尖變形結果

此機組葉片最大凈空距離為34 m，由表3 可見，MBD 和FEM 兩種靜力分析方法結果均未超過34 m，且十分接近，僅相差0.48，與最大凈空距離相比不超過20.59%，采用FEM 梁單元分析葉片變形稍大于MBD結果。

圖7 為1D-FEM 方法計算得到的葉片變形圖。由圖可見，取葉片展向50%與80%兩個位置a、b，原模型a 處葉片切向夾角θa為1.5°，加載變形后切向夾角θa′為-8.5°，角度變化10°，位移2.92 m，原模型b 處葉片切向夾角θb為7.0°，加載后b 處切向夾角θb′為-20.3°，角度變化27.3°，位移14.12 m。

圖7 葉片變形圖

4 葉片變形敏感因素分析

為分析葉片變形影響因素，本文從材料剛度與載荷角度出發，采用BEM 耦合FEM 的方法計算葉片空氣動力特性和結構響應，并將BEM載荷計算結果按比例作為葉片靜力分析的輸入。載荷對葉尖變形分析結果見表4，變形曲線見圖8。

圖8 葉尖變形隨載荷變化曲線

表4 葉尖變形隨載荷變化結果

由圖8 可見，葉尖變形隨載荷增大呈單調遞增正比關系，當載荷比值為1.24 時，葉片將發生掃塔現象。將材料剛度（包括彈性模量和剪切模量）按比例取值，分析剛度對葉尖變形影響，結果見表5，變形曲線見圖9。

圖9 葉尖變形隨剛度變化曲線

表5 葉尖變形隨剛度變化結果

由圖9可見，葉尖變形隨剛度增大而減小，當剛度比值為0.83時，機組將發生掃塔現象。

5 結論

本文以某108 m 長葉片為模型，對不同葉片凈空計算方法進行了研究和對比，計算方法包括氣動性能計算方法與結構響應計算方法，結論如下：

1）不同計算方法得到沿葉展方向推力與力矩均呈先增大后急劇減小的趨勢。相比CFD，BEM 載荷結果在葉片中段更大，沿葉展方向均勻增長，在葉尖位置更小。由于葉片氣動外形，導致CFD 推力在根部及中段較小，在葉尖有明顯增加。

2）塔筒效應和葉片間尾流相互影響均會導致葉片氣動性能的損失。

3）MBD 二階模態分析誤差較FEM 低0.9%，相同載荷下，MBD 葉尖變形27.00 m，相較FEM 變形少0.48 m。

4）載荷增大會引起葉尖變形正比增大，凈空減小，而剛度增大則會使葉尖變形減小，凈空增大，且變形隨剛度變化逐漸降低。

5）CFD 方法由于考慮了葉片實際外形而精確度更高，但計算耗時更長，BEM 方法與CFD 方法差值較小，但計算效率更高；結構分析方面，MBD與FEM 在計算效率方面差距不大，但模態分析MBD 比FEM 準確性更高一些，因此采用BEMMBD 計算葉片凈空可快速得到準確結果，較適合于實際工程。