唐新姿 王效禹 袁可人 彭銳濤

(湘潭大學機械工程學院,湖南湘潭 411105)

引言

風力機常年運行于地表大氣層中,工作環(huán)境非常復雜.由于地表地貌,高度,氣候等差異,實際自然來流風速隨機波動.風速不確定性造成風力機氣動力和載荷分布不確定性,對風電穩(wěn)定性和可靠性產(chǎn)生嚴重影響,開展風力機復雜工況下氣動不確定性研究對于保障風電系統(tǒng)高效可靠運行具有重要意義.

目前風力機氣動載荷研究主要方法有數(shù)值計算方法、葉素動量(blade-element momentum,BEM)理論和風洞試驗等.其中,葉素動量理論是風力機氣動設計普遍采用的方法,該方法應用多基于確定性工況條件.為提高該方法在復雜工況的適應性,國內(nèi)外學者針對該方法進行了相關改進,包括氣動彈性修正[1]、動態(tài)失速修正[2]、三維效應修正[3-7]、阻力修正[8]和風切變修正[9]等,這些修正使得BEM 應用更為廣泛,但不確定工況問題仍然有待解決.針對風電不確定性,國內(nèi)外學者在預測模型[10],控制策略[11-15]和翼型[16-22]等方面開展了相關研究.考慮風速變化,Keshavarzzadeh 等[23]結(jié)合Von Karman 模型和非嵌入式混沌多項式將風速表示為時間、半徑和隨機變量的函數(shù),基于葉素動量理論和牛頓迭代法提出水平軸風力機葉片不確定性優(yōu)化方法.Murcia 等[24]構(gòu)建四維湍流空間,提出在多種不確定湍流參數(shù)下預測風力機氣動性能方法.Antonini 等[25]基于四種湍流模型并對數(shù)值模擬結(jié)果進行加權(quán)處理,揭示考慮風向不確定性的重要性.Toft 等[26]基于加速因子和風力機與測量位置之間距離建立風速不確定性概率模型以評估風力機結(jié)構(gòu)疲勞載荷.Horna 等[27]研究了風速和風向不確定性以及環(huán)境載荷不確定性對海上風力機地基疲勞可靠性的影響.Yan 等[28]提出一種包含槳矩角和風向參數(shù)的風力機功率曲線模型描述能量轉(zhuǎn)化不確定性.

復雜工況下風力機氣動力研究本身是具有挑戰(zhàn)性的研究熱點問題,不確定問題正受到國內(nèi)外學者重點關注.采用數(shù)值模擬或風洞試驗方法進行工況不確定條件下風力機氣動力研究將非常復雜且工作量巨大.本文提出一種基于修正葉素動量理論和非嵌入式概率配置點法的風力機不確定氣動計算方法,建立水平軸風力機不確定性空氣動力學響應模型;基于該模型,以NREL Phase VI S809 風力機為研究對象,量化不確定工況下隨機氣動邊界條件對風力機氣動力的影響程度;從氣動誘導上揭示不確定性在葉片展長方向上的傳播機制,為風力機設計與應用提供理論依據(jù)和重要參考.

1 不確定風力機氣動力學計算模型

1.1 修正葉素動量理論

目前描述風輪作用力與來流速度之間的關系主要有考慮尾渦旋轉(zhuǎn)的修正葉素動量理論.該理論認為,當氣流作用在風輪上產(chǎn)生轉(zhuǎn)矩時,氣流同樣會受到風輪反作用.風輪對氣流的影響用誘導因子表示,葉片截面速度關系如圖1 所示.

圖1 葉片截面速度示意圖Fig.1 Velocity of blade element

根據(jù)動量理論,得葉片扭矩Q和推力T為

式中,a為軸向誘導因子;b為周向誘導因子;r為葉素半徑,m;ρ 為空氣密度,kg/m3;V1為來流風速,m/s;? 為葉輪轉(zhuǎn)速,r/min;α 為攻角,(?);β 為槳矩角,(?).

根據(jù)葉素理論,同樣得到葉片扭矩Q和推力T為

式中,B為葉片數(shù);c為弦長,m;Cl為翼型升力系數(shù);Cd為翼型阻力系數(shù);W為相對風速,m/s;? 為入流角,(?).其中,相對風速、來流風速與葉輪轉(zhuǎn)速的關系為

動量理論和葉素理論結(jié)合得

式中,Ct為葉素切向力系數(shù);Cn為葉素法向力系數(shù);σ 為風輪實度,σ=Bc/(2πr).

當軸向誘導因子a>0.38 時,根據(jù)Wilson 和Prandtl 方法進行修正

其中F為葉尖和輪轂損失修正.

葉尖損失修正系數(shù)Ftip為

輪轂損失修正系數(shù)Fhub為

式中,Ftip為葉尖損失;Fhub為輪轂損失;R為葉片半徑,m;r為葉素半徑,m;rhub為輪轂半徑,m.

修正后軸向誘導因子為

修正后周向誘導因子為

通過式(12)和式(13)求解每個半徑位置上周向和軸向誘導因子描述風力機準三維流場變化,進而得到風力機氣動力.

1.2 非嵌入式概率配置點法

本文將風力機氣動力計算方法和非嵌入式概率配置點法 (no-instrusive probabilistic collocation,NIPRC)結(jié)合,進行不確定氣動力計算.NIPRC 是混沌多項式法對隨機變量進行多項式展開形式變換和隨機配置點法的結(jié)合.隨機變量?(x,t,θ)是關于時間t,空間x和隨機變量θ 的函數(shù),其表達式為

式中,?k(x,t)是隨機變量?(x,t,θ)在配置點θk的解,是由確定性計算獲得,hk表示與該點對應的拉格朗日差值多項式,是通過n個配置點的n-1 階多項式,該多項式表達式為

式中,hk(ξ(θi))=δki是和配置點對應的拉格朗日差值多項式,n為配置點個數(shù).

考慮到風速為自然隨機等特點,隨機性符合正態(tài)分布,因此本文配置點采用隨機空間? 上對?(x,t,θ)積分的高斯積分點,其表達式為

式中,xk是和權(quán)函數(shù)ρ(x)對應的正交多項式的零點,即高斯點;w(x)是各高斯點對應的求積系數(shù).

對于高斯分布隨機變量,其概率密度函數(shù)為

與其對應的Hermite 多項式的權(quán)函數(shù)為

對應高斯積分為Gaussian-Hermite 積分,高斯積分點即為n次Hermite 多項式的n個零點,對應的求積系數(shù)即為權(quán)重.其具體表達式為

在各個配置點進行確定性分析之后,根據(jù)統(tǒng)計學中均值和方差的定義,推導出均值和方差的計算公式,得到不確定風速對于氣動力的影響,具體表達式為

2 研究對象及方法驗證

NREL Phase VI 風力機[29]為三葉片失速調(diào)節(jié)風力機,葉片為扭曲錐形,基于S809 翼型設計,葉片半徑為5.03 m,額定功率為20 kW,錐角為3.4?,轉(zhuǎn)速為72 r/min.葉片輪轂長度為0.724 m,從0.724 m 處的輪轂到S809 翼型之間為線性過度.該風力機模型由美國國家航空航天局基于24.4 m×36.6 m 風洞進行空氣動力學實驗.

首先基于修正葉素動量理論建立確定性風力機氣動模型,編制MATLAB 程序進行風力機氣動分析.圖2 為功率和推力計算曲線和實驗曲線對比.由圖可知,功率和推力計算值與實際測量值基本吻合,驗證了本文所建立氣動模型的正確性.

根據(jù)自然風速統(tǒng)計規(guī)律,假設隨機風速服從高斯分布,采用非嵌入式概率配置點法結(jié)合修正葉素動量理論進行不確定性分析,計算流程如圖3 所示.首先輸入風速工況條件,然后利用NIPRC 方法得到風速配置點及其求積系數(shù),再應用BEM 方法求解風力機氣動參數(shù),最后計算輸出功率、推力、揮舞彎矩和擺振彎矩,及其均值和標準差.

圖2 輸出功率和推力計算值和實驗值對比Fig.2 Comparison of calculated and experimental output power and thrust

圖3 不確定性分析計算流程Fig.3 Flow chart of uncertainty calculation

平均風速為7 m/s,考慮到平坦地形大氣層常規(guī)湍流強度為5%～15%[30],取標準差為平均值10%.對應一階、二階和三階配置點和權(quán)重如表1 所示.

表1 平均風速為7 m/s 的不同階配置點及其響應的權(quán)重Table 1 Different order configuration points with average wind speed of 7 m/s and the weight of their responses

表2 給出了風速為7 m/s 時一階、二階和三階配置點功率和推力計算值對比.二階和三階配置點計算結(jié)果基本吻合,因此后續(xù)基于二階NIPRC 計算.

表2 不同階配置點功率和推力的計算值Table 2 Calculated power and thrust with different order configuration points

為驗證不確定分析計算模型,采用蒙特卡洛(Monte Carlo,MC)方法進行驗證.表3 給出了風速為7 m/s 時采用NIPRC 和MC 計算功率和推力的平均值和標準差對比.采用兩種方法計算所得結(jié)果基本一致,功率和推力最大誤差分別為3.10%和2.81%,驗證了不確定分析模型的正確性.

表3 NIPRC 和MC 功率和推力計算值對比Table 3 Calculated power and thrust based on NIPRC and MC

3 結(jié)果分析

取5～9 m/s 五個平均風速,標準差S tD分別為5%,10%和15%倍均值,進行多風速工況多種不確定度對比分析.

3.1 功率和推力

圖4 為不同平均風速及標準差下風力機輸出功率不確定性對比.相同平均風速不同標準差下輸出功率均值基本一致,標準差越大輸出功率波動幅度越大.平均風速為7 m/s 輸出功率波動達到最大,其中標準差為5%時最大波動幅度為13.44%,標準差為10%時最大波動幅度為24.94%,標準差為15%時最大波動幅度為35.11%.

圖4 不同風速及標準差下風力機輸出功率不確定性Fig.4 Uncertainty of wind turbine output power under different wind speeds and standard deviations

圖5 為不同平均風速及標準差下風力機推力不確定性對比.相同平均風速不同標準差下推力均值基本一致,標準差越大推力波動幅度越大.標準差為5%時,推力的波動幅度在風速為7 m/s 最大,最大波動幅度為8.00%;標準差為10%時,推力的波動幅度在風速為7 m/s 時最大,最大波動幅度為14.96%;標準差為15%時,推力的波動幅度在風速為7 m/s 時最大,最大波動幅度為22.02%.

圖5 不同風速及標準差下風力機推力不確定性Fig.5 Uncertainty of wind turbine thrust under different wind speeds and standard deviations

3.2 揮舞彎矩和擺振彎矩

圖6 為不同風速及標準差下風力機揮舞彎矩不確定性對比.不同標準差下?lián)]舞彎矩均值基本一致,標準差越大,揮舞彎矩波動越大.三種標準差下?lián)]舞彎矩波動均在風速為7 m/s 時最大,標準差為5%時最大波動幅度為7.20%,標準差為10%時最大波動幅度為13.64%,標準差為15%最大波動幅度為19.90%.

圖6 不同風速及標準差下葉根揮舞彎矩不確定性Fig.6 Uncertainty of the flapwise moment of the blade root under different wind speeds and standard deviations

為進一步分析不確定性影響,計算揮舞彎矩在葉片全長上的分布,由圖6 分析得到在風速為7 m/s時揮舞彎矩的波動最大,因此選擇7 m/s 計算.

圖7 為風速為7 m/s 時不同標準差下?lián)]舞彎矩在葉片全長上的分布.三種標準差下均值基本相同,標準差越大揮舞彎矩波動越大.標準差為5%時,揮舞彎矩在r/R為0.83 時波動最大,最大波動幅度為7.01%;標準差為10%時,揮舞彎矩在r/R為0.88 時波動最大,最大波動幅度為13.39%;標準差為15%時,揮舞彎矩在r/R為0.88 時波動最大,最大波動幅度為19.33%.風力機揮舞彎矩波動在葉根部分比較小,在葉尖部分比較大,主要因為彎矩計算受半徑影響,而載荷不確定性在葉片全長上變化不大.

圖7 不同標準差下?lián)]舞彎矩在葉片全長上的分布Fig.7 Distribution of flapwise moment on the full length of the blade under different standard deviations

圖8 不同風速及標準差下葉根擺振彎矩不確定性Fig.8 Uncertainty of the edgewise moment of the blade root under different wind speeds and standard deviations

圖8 為不同風速及標準差下風力機擺振彎矩不確定性對比.三種標準差下擺振彎矩均值基本一致,標準差越大擺振彎矩波動越大.三種標準差下擺振彎矩波動均在風速為7 m/s 時最大,標準差為5%時最大波動幅度為12.84%,標準差為10%時最大波動幅度為23.72%,標準差為15%時最大波動幅度為33.49%.

為進一步分析不確定性影響,計算擺振彎矩在葉片全長上的分布,由圖8 分析得到在風速為7 m/s時擺振彎矩的波動最大,因此選擇7 m/s 計算.

圖9 為風速為7 m/s 時不同標準差下擺振彎矩在葉片全長上的分布.三種標準差下均值基本相同,標準差越大擺振彎矩波動越大.三種標準差下擺振彎矩均在r/R為0.6～0.8 時波動最大,標準差為5%時最大波動幅度為13.32%,標準差為10%時最大波動幅度為25.16%,標準差為15%時最大波動幅度為35.18%.

圖9 不同標準差下擺振彎矩在葉片全長上的分布Fig.9 Distribution of edgewise moment on the full length of the blade under different standard deviations

3.3 周向和軸向誘導因子

為進一步分析不確定性影響,計算軸向誘導因子和周向誘導因子在葉片全長上分布.參考文獻[29]中給出了平均風速為6.7 m/s 時軸向誘導因子實際測量值,因此選擇6.7 m/s 進行計算.

圖10 為三種標準差下周向誘導因子不確定性在葉片全長分布對比.三種標準差下均值基本相同,在相同半徑位置標準差越大,周向誘導因子波動越大,且越靠近葉根位置周向誘導因子波動越大.三種標準差下,周向誘導因子均在r/R為0.18 時最大,標準差為5%時最大波動幅度為13.23%,標準差為10%時最大波動幅度為24.38%,標準差為15%時最大波動幅度為35.78%.表明風輪葉片中部到葉尖位置流動較為穩(wěn)定,受不確定性影響較小,而葉片根部靠近風輪中心受風速不確定性影響較大,可以在根部位置采取措施以改善不確定流動對風力機氣動性能的影響.

圖10 不同標準差下周向誘導因子不確定性分布Fig.10 Uncertainty distribution of tangential induction factor under different standard deviations

圖11 為平均風速為6.7 m/s 時,三種標準差下軸向誘導因子不確定性在葉片全長分布對比.在三種標準差下均值基本相同,相同半徑位置,風速標準差越大軸向誘導因子波動越大.標準差為5%時,軸向誘導因子的波動在r/R為0.23 最大,最大波動幅度為4.71%;標準差為10%時,軸向誘導因子的波動在r/R為0.38 時最大,最大波動幅度為7.42%;標準差為15%時,軸向誘導因子波動在r/R為0.43 時最大,最大波動幅度為11.64%.

圖11 不同標準差下軸向誘導因子不確定性分布Fig.11 Uncertainty distribution of axial induction factor under different standard deviations

4 結(jié)論

基于修正葉素動量理論和非嵌入式概率配置點法,建立風力機不確定氣動分析模型,獲得風速不確定性對風力機功率和氣動力影響程度,通過流動誘導葉片全長分布特性揭示不確定性在風力機葉片流向方向傳播機制.主要結(jié)論如下:

(1)風速波動對風力機功率和氣動力影響顯著.在討論風速范圍內(nèi),風速標準差由0.05 倍增大至0.15倍均值,功率和推力最大波動幅度分別由13.44%和8.00%增大至35.11%和22.02%;葉片揮舞彎矩和擺振彎矩最大波動幅度分別由7.20%和12.84%增大至19.90%和33.49%.

(2)風速標準差越大葉片彎矩波動越大.風速標準差為15%時,揮舞彎矩最大波動幅度為19.33%,擺振彎矩最大波動幅度為35.18%.風力機揮舞彎矩和擺振彎矩波動在葉根部分比較小在葉尖部分比較大,主要因為彎矩計算受半徑影響,而載荷不確定性在葉片全長上變化不大.

(3)風力機葉片周向氣動誘導因子不確定性從葉根到葉尖逐漸減小,在根部位置流動誘導受不確定因素影響最大.可以考慮在該部分采取流動控制措施降低葉片氣動性能對不確定因素敏感程度.