基于升力線理論的風(fēng)力機(jī)氣動性能分析方法*

李德源，何正舉，郭坤翔

(廣東工業(yè)大學(xué) 機(jī)電工程學(xué)院，廣州 510006)

80年代Jeng提出升力線算法[9]，將葉片用一根環(huán)量連續(xù)變化的升力線來代替，根據(jù)霍姆亥茲渦量守恒定理，從升力線上拽出渦絲張成螺旋尾渦面，即由附著渦(升力線)和尾渦構(gòu)成流動模型.誘導(dǎo)速度由Biot-Savart定理求解，由渦絲相互之間的誘導(dǎo)考慮有限葉片的效應(yīng)，不需引入任何修正，是一種簡潔有效的算法.

近年來，風(fēng)力機(jī)組呈大型化發(fā)展趨勢，氣動載荷、重力、慣性載荷及操縱載荷隨著葉片長度的增加急劇上升[10]，需要提高氣動模型的三維計算能力.葉片近壁區(qū)域流場的展向流動和尾渦系統(tǒng)的三維運(yùn)動等因素將使基于二維流動假設(shè)的BEM模型面臨困難，需要對該模型引入人工修正模型，如葉根葉尖損失修正模型和三維無粘延遲失速模型等，才能得出符合實際的結(jié)果.而對升力線模型來說，雖然也需要根據(jù)二維翼型的升阻力曲線使其計算封閉，但通過Biot-Savart公式計算三維誘導(dǎo)速度在一定程度上提高了該模型的三維流場分析能力，尤其是尾渦模型的引入，大大地增加了模型對尾流分析的可靠性.

基于以上分析，本文采用螺旋尾渦升力線模型來計算葉片氣動性能參數(shù).通過對附著渦分布、控制點的誘導(dǎo)速度以及迭代法求解升力系數(shù)等問題進(jìn)行研究和分析，建立了基于螺旋尾渦升力線模型的風(fēng)力機(jī)氣動性能數(shù)值分析算法，并在FORTRAN平臺中創(chuàng)建了分析程序.算例分析了NREL實驗室公布的5 MW風(fēng)力機(jī)[11]的氣動性能，并與NREL發(fā)布的結(jié)果進(jìn)行了對比驗證.

1 風(fēng)力機(jī)葉片升力線模型

由伯努利的流體機(jī)械能守恒原理可知，由于流過翼型壓力面和吸力面的流速不同導(dǎo)致作用在壓力面和吸力面產(chǎn)生壓力差，從而使翼型產(chǎn)生升力.根據(jù)Kutta-Joukowski定理，翼型的升力L可以由自由來流風(fēng)速w與翼型環(huán)量Γ之間的相互作用來表示，即

dL=ρΓwdr

(1)

式中：ρ為當(dāng)?shù)乜諝饷芏?；dr為葉片沿展向方向的微段.根據(jù)這一思想，升力線模型利用一根渦量沿展向變化的線渦替代葉片，表征葉片與流場之間的相互作用，示意圖如圖1所示.

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

在升力線模型中，渦量沿展向分布的線渦代表了三維流場中真實三維葉片的附著環(huán)量分布，該線渦稱為集中附著渦.此外，由Helmholtz第二定理可知，在無粘流體環(huán)境中，渦量并不能在流體內(nèi)部終止，而只能延伸至流體邊界或構(gòu)成環(huán).因此，附著渦沿展向的變化量并不能消失而必然脫落流向下游，形成尾流渦.升力線模型以一根集中渦表征葉片的空氣動力作用，對于該集中渦來說，只能感受到力的作用而不能感受到力矩的作用，因此需要將該集中渦放置在翼型的某一點，使得該點的俯仰力矩系數(shù)不隨攻角的改變而變化，而這樣的點恰好為翼型的氣動中心[12].

升力線模型通過計算附著渦線上參考點位置處的誘導(dǎo)速度，獲取葉片在展向位置某一處的有效攻角，并通過插值獲得該處翼型的升阻力系數(shù)，從而實現(xiàn)當(dāng)前翼型的升阻力計算；再通過在整個葉片展向上進(jìn)行積分，實現(xiàn)葉片的扭矩和功率等氣動性能參數(shù)的計算.

1.1 誘導(dǎo)速度計算

直線線渦對空間內(nèi)目標(biāo)點位置(即控制點)處產(chǎn)生的誘導(dǎo)速度可用Biot-Savart公式描述，即

(2)

式中：wt為葉片控制點處的誘導(dǎo)速度；γ為流場中分布的渦量；r為葉片半徑向量；ds為線渦微段矢量.

由于本文采用的是剛性尾渦模型，該渦線為螺旋狀，螺旋尾渦向葉片下游延伸至無限遠(yuǎn)處，形成誘導(dǎo)速度場，坐標(biāo)系如圖2所示.

將螺旋尾渦離散為若干段直線線渦，由式(2)可推出由螺旋尾渦引起的葉片展向位置處的誘導(dǎo)速度，其表達(dá)式為

圖2 螺旋尾渦示意圖Fig.2 Schematic diagram of spiral wake vortex

(3)

式中：s為輪轂中心指向微端dη的向量；wi為葉片展向位置處的誘導(dǎo)速度；r′為輪轂中心指向目標(biāo)位置點的向量；dη為第k個葉片展向位置r處發(fā)散出來的螺旋線渦微段矢量；Γ為附著渦的環(huán)量.將dwi沿著坐標(biāo)軸分解成分量形式，可得

(4)

根據(jù)剛性尾渦理論可知，dwx=0.若將誘導(dǎo)速度wi沿著垂直和平行于相對風(fēng)速w′方向進(jìn)行分解(如圖3所示)，則有

dwn=dwzcosφ-dwysinφ

(5)

dwt=dwzsinφ+dwycosφ

(6)

(7)

式中：φ為入流角；V0為來流風(fēng)速；Ω為葉片轉(zhuǎn)速.

圖3 誘導(dǎo)速度分解Fig.3 Decomposition of induced velocity

誘導(dǎo)速度總量可以通過式(5)、(6)從葉片輪轂到葉尖處積分得到.Moriya證明了在葉片上任意位置均有wt=0，設(shè)葉片上展向無量綱位置參數(shù)為

(8)

式中，R為葉片半徑.葉片展向位置任意位置ξ′處，由所有尾渦引起的總誘導(dǎo)速度為

(9)

當(dāng)θ=θk=0時，積分在點ξ=ξ′處的值變?yōu)闊o窮大，產(chǎn)生奇點，而且奇點附近的積分值也無法確定，為此，本文引進(jìn)一個額外的參數(shù)I，即誘導(dǎo)因子，其定義為

(10)

式中：wn2和wn1分別為渦量強(qiáng)度相同的螺旋線渦和直線渦產(chǎn)生的誘導(dǎo)速度；I為連續(xù)函數(shù).

由文獻(xiàn)[9]可知，誘導(dǎo)因子I是關(guān)于ξ、ξ′和λ0的函數(shù)，即

(11)

當(dāng)θ=θk=0時，積分在點ξ=ξ′處的值同樣會變?yōu)闊o窮大，實際上可以近似地認(rèn)為在這個點處dwn2=dwn1，即I=1，且由于誘導(dǎo)因子為連續(xù)函數(shù)，可以通過在奇點附近進(jìn)行樣條插值來確定奇點以及其附近的誘導(dǎo)因子的值.

當(dāng)葉片上各點處的誘導(dǎo)因子確定之后，其誘導(dǎo)速度也可確定，即

(12)

誘導(dǎo)速度與各攻角的關(guān)系如圖4所示，其中，圖4中各變量的表達(dá)式為

(13)

(14)

αe=αi+αg

(15)

(16)

根據(jù)圖4可以得到

(17)

圖4 誘導(dǎo)速度與各攻角的關(guān)系Fig.4 Relationship among induced velocity and attack angles

設(shè)環(huán)量為傅里葉正弦級數(shù)，并且滿足邊界條件Γ(ξhub)=Γ(1)=0，則環(huán)量函數(shù)為

(18)

式中，無窮級數(shù)求和可近似用有限項求和來計算，項數(shù)m越多，環(huán)量Γ(ξ)的計算越精確，一般m取為計算控制點的數(shù)目.

根據(jù)Kutta-Joukovski定理可知，環(huán)量與升力系數(shù)之間的關(guān)系為

(19)

式中：c為翼型弦長；CL為升力系數(shù).對于每一個有效攻角，均可以通過二維翼型數(shù)據(jù)表并通過線性插值得到一個與之相對應(yīng)的升力系數(shù)，即CL=CL(αe).

對于初始計算，合成風(fēng)速w即為來流風(fēng)速V0，根據(jù)式(19)可得

(20)

在小攻角范圍內(nèi)，升力系數(shù)與有效攻角之間近似為線性關(guān)系，即

CL=a0αe+b0

(21)

式中，a0、b0為常數(shù).綜合式(12)、(17)、(18)、(21)可得

(22)

將方程(22)應(yīng)用于升力線各計算點ξ′處，即可得到M個方程，求解該方程組可以得到Am的值，將Am代入式(18)即可得到環(huán)量Γ的值；將式(21)代入式(19)，可得有效攻角與環(huán)量之間的關(guān)系，即

(23)

由環(huán)量值即可求出有效攻角αe.

1.2 迭代法求解有效攻角

迭代結(jié)束之后，可以計算出葉片的軸向力F、力矩Q以及功率P等葉片氣動性能參數(shù)，其表達(dá)式為

(24)

(25)

(26)

式中，CD為阻力系數(shù)，對于每一個有效攻角，都可以通過二維翼型數(shù)據(jù)表并通過線性插值得到一個與之相對應(yīng)的阻力系數(shù)CD.本文中ψ取為0°.

2 算例分析

根據(jù)以上理論分析與數(shù)值求解分析方法，在FORTRAN平臺上建立了完整的氣動性能數(shù)值計算程序，實現(xiàn)輸入已知參數(shù)即可求得多葉片風(fēng)輪的氣動性能參數(shù).

為了驗證建立的升力線模型的準(zhǔn)確性，以NREL實驗室公布的5 MW風(fēng)力機(jī)葉片數(shù)據(jù)為研究案例，對其氣動性能進(jìn)行計算，并將計算結(jié)果與NREL實驗室提供的結(jié)果(基于葉素動量理論進(jìn)行的計算)進(jìn)行對比.

為了保證計算效率和計算結(jié)果的可靠性，離散控制點數(shù)和位置的選取，一般可根據(jù)葉素動量理論計算控制點的選取方法，理論上控制點越多，相應(yīng)的計算精度越高，但一般不宜取得太多，否則增大了計算量且對精度沒有太大提高[9].本文計算中，根據(jù)周傳捷等[13]的研究可知，將葉片離散為13段，既能保證計算效率，又能得到較為滿意的結(jié)果.由于葉片靠近輪轂中心的前兩段為圓筒狀，其升阻力系數(shù)在任意攻角下均為恒定值，故其升阻力系數(shù)與攻角的函數(shù)曲線斜率為0，在進(jìn)行式(12)的計算時分母為0，會導(dǎo)致結(jié)果變?yōu)闊o窮，使得計算無法進(jìn)行，因此，可以將其剝離出來單獨(dú)討論，由于升阻力系數(shù)均為已知，故可通過式(24)、(25)直接計算出前兩段圓筒的軸向力以及力矩.

基于以上分析，假定前兩段圓筒形的環(huán)量均為0，則可將葉片第三段起始截面假想為輪轂半徑ξhub，因此，葉片的扭角和弦長分布如圖5所示.

圖5 葉片扭角與弦長分布Fig.5 Blade twist angle and chord length distribution

根據(jù)NREL公布的5 MW風(fēng)力機(jī)葉片二維翼型數(shù)據(jù)表[11]，各翼型的升力系數(shù)曲線的斜率以及截距如表1所示.

表1 葉片各截面位置參數(shù)Tab.1 Parameters of blade sections

將上述已知參數(shù)作為輸入文件導(dǎo)入用FORTRAN編寫的升力線模型程序，計算出該5 MW風(fēng)力機(jī)在不同風(fēng)速下的氣動力特性.圖6為不同風(fēng)速下葉片展向位置的環(huán)量分布圖.

圖6 不同風(fēng)速下沿葉片展向方向的環(huán)量分布Fig.6 Distribution of circular rector along blade span direction at different wind velocities

根據(jù)文獻(xiàn)[11]可知，當(dāng)風(fēng)速低于額定風(fēng)速時，無需對槳距角進(jìn)行調(diào)整，而在更高的風(fēng)速下，為了使風(fēng)力機(jī)產(chǎn)生穩(wěn)定的機(jī)械功率，在恒定的轉(zhuǎn)速12.1 r/min下需對葉片的槳距角進(jìn)行相應(yīng)調(diào)整，其相應(yīng)風(fēng)速下的變槳距角為文獻(xiàn)[11]通過風(fēng)洞試驗所得.表2為不同風(fēng)速下變槳距后利用升力線模型所求得的風(fēng)力機(jī)軸向力、功率以及扭矩.

將不同風(fēng)速下采用升力線理論計算所得結(jié)果中的推力、功率和扭矩，與NREL實驗室公布的5 MW風(fēng)力機(jī)采用引進(jìn)修正模型的葉素動量理論計算所得的相應(yīng)參數(shù)[11]進(jìn)行對比，結(jié)果如圖7所示.

由圖7可以看出，用升力線理論計算得到的風(fēng)力機(jī)的功率、力矩、軸向力與文獻(xiàn)[11]基于修正的葉素動量理論計算結(jié)果比較接近.通過NREL實驗室公布的5 MW風(fēng)力機(jī)葉片計算案例可以看出，風(fēng)力機(jī)葉片氣動性能參數(shù)的升力線模型具有較高的準(zhǔn)確性.

3 結(jié) 論

本文采用基于螺旋尾渦的升力線模型來研究風(fēng)力機(jī)三維流場的空氣動力學(xué)性能，研究了升力線模型中附著渦分布和誘導(dǎo)速度的計算，并通過對升力系數(shù)進(jìn)行迭代收斂來修正環(huán)量，從而確定準(zhǔn)確的攻角和升阻力系數(shù)，進(jìn)而計算出風(fēng)力機(jī)的各項氣動性能參數(shù).值得注意的是，誘導(dǎo)速度的求解考慮了風(fēng)力機(jī)各個葉片之間的相互誘導(dǎo)效應(yīng)，因此，相比于葉素動量理論而言，升力線理論可以極大提高風(fēng)力機(jī)的三維計算能力.

基于升力線模型的理論基礎(chǔ)，本文建立了風(fēng)力機(jī)葉片的升力線數(shù)值計算模型，編譯了適用于多個葉片風(fēng)輪的氣動性能參數(shù)分析程序，并與NREL實驗室公布的基于葉素動量理論的5 MW風(fēng)力機(jī)氣動性能參數(shù)進(jìn)行對比，結(jié)果表明，基于螺旋剛性尾渦模型的升力線算法通過結(jié)合翼型的升阻力特性可以有效計算風(fēng)力機(jī)的氣動性能參數(shù)，而且無需引入應(yīng)用葉素動量理論時需要的大量人為修正模型即可以滿足工程要求的精度.對于葉素動量理論而言，人為引入模型的不同可能會導(dǎo)致計算結(jié)果存在一定偏差，而升力線模型由于無需引入各種修正模型，減少了繁瑣的模型引入，相對于葉素動量理論和CFD理論而言，其計算效率也會提高，同時也可以在一定程度上保證計算的穩(wěn)定性.由于升力線理論考慮到了葉片展向之間的流動，其不僅適用于直葉片，對于積疊線彎曲的、三維流動更強(qiáng)的葉片(如后掠葉片)設(shè)計與氣動性能計算同樣有效，并能與葉片結(jié)構(gòu)模型結(jié)合進(jìn)行高效而準(zhǔn)確的氣彈耦合分析.

表2 不同風(fēng)速下風(fēng)輪氣動性能計算結(jié)果Tab.2 Calculated results of aerodynamic loads of wind turbine at different wind velocities

圖7 不同風(fēng)速下兩種理論計算結(jié)果的比較Fig.7 Comparison between two theoretical calculation results at different wind velocities