基于渦旋升力線模型的風機后掠葉片性能優化*

李德源， 劉佳宇， 夏鴻建， 張海波

(廣東工業大學 機電工程學院， 廣州 510006)

目前，大功率兆瓦級風力機的空氣動力學和結構特性優化設計已經成為了研究重點.而風輪作為風力機的主要組成部分，其成本占到整個風力機的20%左右，同時風輪也是風力機捕獲風能的關鍵部分，因此，國內外學者從不同的方面，對于風力機葉片進行優化設計研究，以獲取良好的氣動特性和結構特性.研究者們提出在保持風輪直徑不變或增加較少的前提下，將葉片適當后掠的概念，也稱為幾何自適應[1].研究表明：作為一種自適應的氣動降載布局設計，后掠型葉片在相同結構載荷下，可增加風輪直徑，從而增大捕風面積，增加發電量，降低度電成本；在保證相同發電量下，可有效減少結構載荷，降低風輪與相關零部件的制造成本，提高結構可靠性，降低維護成本，因此，這種葉片具備較好的使用價值和應用前景[2].

后掠葉片在氣動載荷作用下發生沿順槳方向的扭轉，雖然降低了葉片上的氣動力，但是也可能降低葉片的功率輸出，因此，在設計后掠自適應葉片時，進行一體化設計優化是必需的[3]，這其中就包括各運行風況下，后掠葉片槳距角的優化設計.丁勇鋼等[4]通過對風機控制系統的優化，實現發電量的最大化和葉片軸向推力的最小化.

后掠葉片在結構和空氣動力學方面表現出較強非線性，需要采用恰當的非線性氣彈耦合模型表達空氣動力與葉片彈性變形以及振動時的相互作用.為此，本文通過將柔性葉片離散為由運動副和力元聯接的多體系統[5-6]，應用計算多體系統動力學理論建立葉片非線性氣彈耦合方程[7].在葉片氣動載荷分析方面，后掠葉片表現出顯著的三維流動，螺旋尾渦升力線模型將葉片用環量連續變化的升力線代替，具有較好的三維流場分析能力[8].

采用“超級單元”方法[9]將后掠葉片離散成若干個由運動副和力元聯接的剛體，應用Roberson-Wittenburg遞推建模方法，結合螺旋尾渦升力線模型建立葉片非線性氣彈耦合方程.以變槳角為優化設計變量，最大年發電量作為目標函數，將風力機額定功率、最大葉尖揮舞位移和最大軸向推力設為約束條件，對變槳角進行尋優并對結果進行對比分析.

1 后掠葉片非線性氣彈耦合模型

后掠葉片具有與直葉片不同的氣彈特性，有必要建立恰當的后掠葉片氣彈耦合模型.本文在后掠葉片結構非線性模型的基礎上，結合升力線理論建立葉片氣彈耦合模型.

1.1 后掠葉片結構非線性模型

為了準確描述葉片的結構變形，采用“超級單元”建立葉片模型.其中，每個超級單元由4個剛體(B1～B4)組成，各剛體通過運動副、彈簧和阻尼器連接[5-6]，用來描述葉片的揮舞、擺振和扭轉變形，其模型如圖1所示.

圖1 超級單元模型Fig.1 Super-element model

應用多體系統動力學中的R-W方法，取各鉸的轉動方向為廣義坐標，則從柔性葉片的拓撲構形可得其廣義坐標陣為

q(t)=(q1(t)，q2(t)，…，qs(t))

(1)

(2)

1.2 基于升力線理論的后掠葉片氣動載荷計算

由于后掠葉片具有顯著的展向流動和積疊線彎曲特點，采用升力線模型模擬后掠葉片流場的三維流動.

升力線模型的核心思想是將風力機葉片三維流場中的渦量分布簡化為一根渦量沿葉片展向變化的線渦，表示葉片與流場之間的相互作用.該線渦代表了葉片三維流場中的環量分布，因而該線渦也稱作附著渦.由Helmholtz第二定理可知，在無粘流體環境中，渦量并不能在流體內部終止，而只能延伸至流體邊界或構成封閉區間，附著渦沿葉展方向變化的渦量會從附著渦上脫落形成尾渦，渦方法模型如圖2所示.

圖2 渦方法模型示意圖Fig.2 Schematic model diagram of vortex method

從附著渦上脫落的螺旋尾渦流向葉片下游無限遠處，形成誘導速度場，本文采用剛性尾渦模型，尾渦坐標系如圖3所示.

圖3 剛性尾渦及其坐標系示意圖Fig.3 Schematic diagram of rigid wake vortex and its coordinate system

文獻[10]指出，誘導速度Wn垂直于未受擾動的相對風速W′，并給出了附著渦任意位置e′(已知慣性坐標(x′，y′，z′))處由所有尾渦引起的誘導速度，即

(3)

式中：Γ(ξ)為螺旋尾渦的環量；R為風輪半徑；ξ、ξ′為葉展方向無量綱位置參數；ξhub為葉展方向輪轂位置處的無量綱參數；λ0=V0/(RΩcosΛ)；V0為來流風速；Λ為當地后掠角；Ω為風輪轉速；ψ為風輪錐角；θ為螺旋尾渦的方位角；θk為控制點所在的葉片與發出尾渦的葉片之間的夾角；N為葉片數；K為控制點所在葉片序號.

附著渦線為空間曲線，渦量從輪轂沿著葉片徑向流向葉尖，形成誘導速度場，附著渦坐標系如圖4所示.

葉片展向位置處受到附著渦上某一直線渦量微段dn′的誘導速度為

(4)

圖4 附著渦線及其坐標系示意圖Fig.4 Schematic diagram of attached vortex line and its coordinate system

式中：dWi2為附著渦誘導速度；dwx2、dwy2、dwz2分別為dWi2在x、y、z軸的分量；dn′為附著渦上的微段向量；s′為輪轂中心指向微段dn′的向量；dn′和s′與積疊線形狀有關.附著渦誘導速度垂直于相對風速，經過葉片積分得到所有葉片附著渦對控制點的誘導速度為

(5)

式中：φ為入流角；Vh和Vb分別為葉片的揮舞、擺振速度.當N=K時，即|s′-e′|=0，式(4)中分子和分母均為0，這不符合物理實際.實際上，由于粘性的作用會使得靠近渦線的中心速度趨于有限值，形成所謂的粘性渦核[8].基于粘性渦核模型建立修正的Biot-Savart公式，在計算誘導速度時引入粘性修正因子，即

(6)

式中，Re為渦核的有效半徑.修正的Biot-Savart公式為

dwy2sinφ)dξ

(7)

控制點上總的誘導速度為

Wn=Wn1+Wn2

(8)

誘導速度表達式(3)和(7)中環量項Γ(ξ)未知.假設附著渦的環量分布為Fourier級數，并且滿足邊界條件Γ(ξhub)=Γ(1)=0，設環量分布函數為

(9)

一般m取9即可[10].根據Kutta-Joukovski定理可知，環量與升力系數之間的關系為

Γ=0.5cWCL

(10)

式中：c為對應翼型的弦長；CL為升力系數；W為截面位置受擾動的相對風速.

聯立式(3)和(7)～(10)，代入m個控制點的已知數據，即可得到m個方程，方程組中只有Am未知，求解該方程組可得到Am的初值，利用迭代法[10]可以得到數值精度更高的有效攻角αe.迭代結束得到有效攻角αe和誘導攻角αi，通過二維翼型數據表線性插值得到升力系數CL和阻力系數CD后，沿著葉片徑向積分求得軸向推力F和轉動功率P，即

(11)

式中：ρ為空氣密度；φ′為實際入流角.

1.3 后掠葉片氣彈耦合模型驗證算例

為驗證后掠葉片非線性氣彈耦合模型的準確性和可行性，選用NREL實驗室的5 MW風力機直葉片數據[11]和后掠4 m葉片數據[12]進行模型驗證.

采用本文建立的氣彈耦合模型分別分析了直葉片和后掠葉片風輪在5～25 m/s風速下的氣動性能與結構響應，并與文獻[12]采用葉素動量理論(BEM)氣動模型的分析結果進行對比，結果如圖5所示.

從圖5可以看出，無論直、后掠葉片，應用本文所建立的氣彈耦合模型能夠比較準確地進行氣動性能計算，結果與文獻[12]的結果基本一致.兩種分析模型結果隨風速增大相差更大一些，這主要是因為對于大型柔性葉片，存在著沿葉片展向的流動，后掠葉片這種流動更為顯著一些，使用葉素動量理論計算氣動性能時，忽略了氣流在葉片展向的流動影響，之后可通過葉根和葉尖損失等經驗修正模型進行彌補.而升力線模型考慮了葉片的展向流動與尾流對整個葉片的影響，無需引入經驗修正模型，考慮的影響因素和分析效率更高.在結構上由于本文所建立的“超級單元”模型和文獻[12]鐵木辛柯梁的模型差異，帶來計算模擬結果的誤差.

圖5 升力線理論和葉素動量理論的穩態計算結果Fig.5 Steady-state calculation results for lift line and blade element momentum theories

算例驗證結果表明，相較于葉素動量理論，基于升力線所建立的氣彈耦合模型具有更好的三維流場分析能力，驗證了該方法的有效性與準確性.

2 優化設計模型

2.1 后掠葉片外形設計與運行條件

文獻[12]后掠葉片積疊線形狀計算表達式為

(12)

此處分別取a=8，b=12，即葉尖后掠4 m，R0為葉片長度，R0=61.5 m，則根據式(12)積疊線形狀為

(13)

根據文獻[12]取后掠葉片的切入風速為5 m/s、切出風速為25 m/s.在較低風速時葉片保持尖速比8.7運行，隨著風速的增加，葉片轉速從6 r/min開始增大，當轉速提高到12.1 r/min時不再保持額定尖速比，在額定風速達到11.4 m/s后，風輪轉速保持為12.1 r/min.

2.2 目標函數

年發電量(AEP)與風場風速的威布爾分布有關[13]，理論年發電量的計算公式為

(14)

式中：H為年小時，此處取8 766 h；vi和vi+1為風速；P(vi)和P(vi+1)為風速vi和vi+1對應的功率；B和u分別為風速威布爾分布的尺度因子和形狀因子.本文中計算AEP時u取2，B取7.334 5.

2.3 約束條件

為了確保傳遞到發電機的軸功率不超過其最大容量，需要約束風力機風輪在各風速下的轉動功率不能超過額定功率5.3 MW[11]，即

Pmax≤5.3 MW

(15)

為保證葉片不會因為揮舞方向變形過大引起結構破壞，變槳角優化后的風機葉片最大葉尖揮舞位移小于NREL 5 MW直葉片的最大葉尖揮舞位移[12]，即

δmax≤5.6 m

(16)

風力機風輪的工作壽命、結構強度均與最大軸向推力有關.變槳角優化后的風機葉片最大軸向推力小于NREL 5 MW直葉片的最大軸向推力[12]，即

Fmax≤0.779 MN

(17)

2.4 優化變量與優化方法

采用一維搜索優化方法，以變槳角為優化變量，尋求各個風速下，變槳角與功率的擬合函數關系，并求出擬合函數最大功率所對應的變槳角，計算氣動性能和結構響應.

3 優化結果分析

按照前面所建立的優化設計方法，對后掠量為4 m的風機葉片在各個風速下的變槳角進行優化，分析中采用升力線模型.圖6為優化的后掠葉片與直葉片、后掠葉片[12]在5～25 m/s風速范圍內氣動性能和結構響應的結果對比(仿真環境為勻流風，不考慮塔影效應和重力影響，設計尖速比為8.7，轉速范圍為6～12.1 r/min，優化變量為變槳角).

表1為直葉片、后掠葉片與變槳角優化的后掠葉片在風速11 m/s時最大葉尖揮舞位移、最大軸向推力.變槳角優化的后掠葉片與原后掠葉片最大葉尖揮舞位移小于5.6 m，最大軸向推力小于0.779 MN，保證了葉片不會因為揮舞方向變形過大引起結構破壞且不會對風輪的工作壽命和結構強度產生較大影響.

圖6a結果表明，在6～11 m/s風速下，變槳角優化的后掠葉片輸出功率顯著提升；在12～25 m/s風速下，各葉片的輸出功率已經達到額定功率5.3 MW，功率不再增加.通過式(14)計算變槳角優化的后掠葉片、原后掠葉片和直葉片的年發電量，其值分別為14 393 942、13 257 036、13 378 251 kW·h，變槳角優化的后掠葉片年發電量相較于其余葉片的年發電量有明顯提高，與原后掠葉片相比提高約8.58%、與直葉片相比提高約7.59%.

圖6b、c結果表明，在6～11 m/s風速下，后掠葉片的軸向推力、葉尖揮舞位移有所增加；在12～25 m/s風速下，軸向推力、葉尖揮舞位移顯著減少.通過優化后掠葉片的變槳角，使葉片在高風速下維持5.3 MW的同時，減少了對風輪和塔架的載荷，提高結構可靠性.

圖6d結果表明，原后掠葉片在5～11 m/s風速下變槳角一直減小，這是由于葉片積疊線后掠，在相同條件下后掠葉片切向力要乘上后掠角余弦值，使功率(轉動力矩)比直葉片要??；要讓后掠葉片獲得與直葉片相近的功率(轉動力矩)則需要增大氣動載荷，所以后掠葉片變槳距角比直葉片小，在11 m/s風速下取得最小值.在12～25 m/s風速下變槳角一直增大，使得功率維持在5.3 MW附近，變槳角優化的后掠葉片在7～11 m/s風速下，變槳角進一步減小，增大了有效攻角，提高了升力系數CL，從而提高了功率，增大軸向推力和葉尖揮舞位移.

圖6 氣動性能和結構響應計算結果Fig.6 Calculation results for aerodynamic performance and structural response

表1 11 m/s風速下最大葉尖揮舞位移、最大軸向推力Tab.1 Maximum flapwise tip deflection and axial thrust at wind speed of 11 m/s

4 結 論

本文建立葉片非線性氣彈耦合方程，并驗證模型的有效性和準確性.以5 MW后掠4 m葉片為原型，變槳角為優化變量，年發電量作為目標函數，將風力機額定功率、最大軸向推力、最大葉尖揮舞位移設為約束條件，通過后掠葉片非線性氣彈耦合程序，計算葉片的氣動性能和結構響應，得到結論如下：

1) 利用基于升力線的葉片非線性氣彈耦合程序，計算后掠葉片的氣動性能和結構響應具有較高的效率和精度，為風力機的優化設計提供理論依據.

2) 對變槳角進行優化，所得到的槳距角隨風速變化曲線可以作為風機運行時槳距角變化的實際控制運行規律曲線.

3) 通過該優化模型，變槳角優化的后掠葉片相較于原后掠葉片和直葉片，在年發電量上分別提高了8.58%、7.59%；該優化模型也為風力機葉片的優化設計提供了參考.