趙 騫， 厲 偉， 姚興佳， 邵一川， 郭慶鼎

(1. 沈陽工業(yè)大學(xué) 電氣工程學(xué)院， 沈陽 110870； 2. 沈陽大學(xué) 信息工程學(xué)院， 沈陽 110004)

本文對風(fēng)力機(jī)葉片攻角α及風(fēng)輪錐角β進(jìn)行優(yōu)化，并對α與β之間的相關(guān)性做進(jìn)一步研究，為進(jìn)一步設(shè)計新型更高效風(fēng)力機(jī)提供最佳的優(yōu)化方案.采用兩套優(yōu)化方案對NACA0012翼型的1.5 MW風(fēng)力機(jī)α、β角進(jìn)行實(shí)例優(yōu)化，方案1對α、β進(jìn)行獨(dú)立優(yōu)化：首先應(yīng)用傳統(tǒng)葉素動量(BEM)理論對葉片攻角α進(jìn)行優(yōu)化，確定α的最佳工作點(diǎn)，再在該最佳工作點(diǎn)處，應(yīng)用對β錐角修正所得的β-BEM理論進(jìn)一步對錐角β進(jìn)行優(yōu)化，計算優(yōu)化后風(fēng)能利用系數(shù)Cp1.方案2對α、β角進(jìn)行統(tǒng)一優(yōu)化：應(yīng)用β-BEM理論同時對α、β角進(jìn)行優(yōu)化，計算優(yōu)化后風(fēng)能利用系數(shù)Cp2.通過對Cp1與Cp2進(jìn)行對比分析，研究α、β之間相關(guān)性，并確定最佳優(yōu)化方案.

1 傳統(tǒng)葉素動量理論

1.1 貝茨理論模型

圖1為Betz理論的氣流圖.在風(fēng)輪上游遠(yuǎn)端處，風(fēng)速為υ0，氣流截面積為S0；風(fēng)輪處，風(fēng)速為υ，氣流截面積為S；風(fēng)輪下游遠(yuǎn)端處，風(fēng)速為υ1，氣流截面積為S1.由Betz理論[6]可知，當(dāng)滿足式(1)，即

(1)

(2)

式中，ρ為空氣密度.風(fēng)力機(jī)工作于低風(fēng)速環(huán)境下，此時氣流視為不可壓縮流體，式(1)結(jié)合連續(xù)性方程S0υ0=Sυ=S1υ1得

S0∶S∶S1=1∶2∶3

(3)

式(3)說明風(fēng)輪上、中、下游的氣流截面積逐漸增大，表明在風(fēng)輪處氣流除具有平行于風(fēng)輪旋轉(zhuǎn)軸的軸向速度外，同時具有垂直于旋轉(zhuǎn)軸的徑向速度分量，這為進(jìn)一步引入風(fēng)輪錐角β提供了依據(jù).

圖1 貝茨理論風(fēng)輪氣流圖Fig.1 Airflow of wind wheel in Betz theory

1.2 傳統(tǒng)葉素動量理論

葉素動量理論一直是研究風(fēng)力機(jī)氣動特性的有效方法，相關(guān)研究在文獻(xiàn)中多有報道[7-9].傳統(tǒng)的葉素動量理論采用零錐角模型(β=0°)，即風(fēng)輪旋轉(zhuǎn)面與旋轉(zhuǎn)軸垂直，此時風(fēng)輪處氣流的徑向速度分量與旋轉(zhuǎn)面平行，對風(fēng)輪不產(chǎn)生轉(zhuǎn)矩作用，因此只考慮氣流軸向速度對風(fēng)輪產(chǎn)生的轉(zhuǎn)矩作用.

由動量理論可知，氣流作用在風(fēng)輪根部距離r處時葉素的軸向推力dFn和轉(zhuǎn)矩dM分別為

(4)

dM=4πρωυ0a′(1-a)r3dr

(5)

式中：a為軸向誘導(dǎo)因子；a′為周向誘導(dǎo)因子；ω為轉(zhuǎn)速.同時，由葉素理論可得到

(6)

(7)

式中：Nb為葉片數(shù)；c為葉素弦長；CL、CD為風(fēng)力機(jī)翼型升力、阻力系數(shù)；φ為入流角，其表達(dá)式為

(8)

式中，λ為當(dāng)?shù)厮俣缺?，其表達(dá)式為

(9)

對于某確定的翼型，CL、CD主要由攻角α決定，其表達(dá)式為

α=φ-θ

(10)

式中，θ為節(jié)距角.

分別將式(4)、(6)及式(5)、(7)建立等量關(guān)系，并引入普朗特葉尖損失因子f，可計算出軸向誘導(dǎo)因子a及周向誘導(dǎo)因子a′為

(11)

(12)

式中，σ為實(shí)度，其表達(dá)式為

(13)

進(jìn)一步結(jié)合式(2)、(7)計算出風(fēng)能利用系數(shù)為

Cp=σλ[(1-a)2+λ2(1+a′)2]·

(CLsinφ-CDcosφ)

(14)

給定當(dāng)?shù)厮俣缺圈思肮?jié)距角θ，通過迭代式(8)～(13)計算出a、a′的收斂值，并進(jìn)一步由式(14)計算出風(fēng)能利用系數(shù)Cp.

2 基于錐角β的BEM修正理論

傳統(tǒng)BEM理論模型中沒有涉及到錐角β，因此，若要對β實(shí)現(xiàn)優(yōu)化，須對傳統(tǒng)BEM理論進(jìn)行β修正.圖2為將β引入到傳統(tǒng)葉素動量理論后的氣流分布示意圖，此時風(fēng)輪旋轉(zhuǎn)面為圓錐面，氣流在與旋轉(zhuǎn)軸垂直平面內(nèi)的徑向速度分量與風(fēng)輪旋轉(zhuǎn)圓錐面不平行，因此，徑向速度分量亦將對風(fēng)輪產(chǎn)生轉(zhuǎn)矩作用，進(jìn)而轉(zhuǎn)換為風(fēng)輪輸出功率.

圖2 β錐角時風(fēng)輪氣流分布Fig.2 Airflow distribution of wind wheel with cone angle β

如圖2所示，取與葉根距離為r的風(fēng)輪葉素，其距離旋轉(zhuǎn)軸的垂直距離為

r⊥=rcosβ

(15)

氣流速度與風(fēng)輪旋轉(zhuǎn)錐面垂直的速度分量υ⊥為

(16)

式中：γ為氣流發(fā)散角，表示氣流流線與旋轉(zhuǎn)軸夾角；ac為對軸向誘導(dǎo)因子a的修正值.

將修正式(15)、(16)及β代入傳統(tǒng)BEM理論表達(dá)式(4)～(14)中，得到修正后β-BEM理論表達(dá)式為

(17)

dMc=4πρωυ0a′c(1-ac)r3cos4βdr

(18)

CDsinφc)cosβdr

(19)

CDcosφc)rcos2βdr

(20)

(21)

(22)

αc=φc-θ

(23)

(24)

(25)

(26)

(CLsinφc-CDcosφc)

(27)

β-BEM理論表達(dá)式(17)～(27)在風(fēng)輪錐角β=0°時與修正前表達(dá)式(4)～(14)完全一致，說明β-BEM理論表達(dá)式具有較強(qiáng)的通用性.

3 α、β獨(dú)立優(yōu)化

3.1 NACA0012翼型風(fēng)力機(jī)的α優(yōu)化

為實(shí)現(xiàn)方案1中對α的優(yōu)化，采用BEM理論對NACA0012翼型的1.5 MW風(fēng)力機(jī)CP值進(jìn)行數(shù)值計算.該翼型為NACA經(jīng)典翼型，低風(fēng)速時，該翼型的特性曲線如圖3所示，在小攻角段，CL與α幾乎成正比關(guān)系，在12°攻角時翼型發(fā)生失速，通常情況下，翼型均工作在失速點(diǎn)之前.

圖3 NACA0012翼型特性曲線Fig.3 Characteristic curves of NACA0012 airfoil

選取翼型中部葉素作為計算區(qū)域時，不涉及葉尖損失，因此葉尖損失因子f取1.將λ、θ做為循環(huán)變量，對于每一組λ、θ，通過式(8)～(13)迭代求解各參數(shù)，進(jìn)一步代入式(14)計算CP，最終得到CP隨λ、θ變化的曲面圖.使用C++語言編寫該迭代程序，對6萬組λ、θ點(diǎn)的CP值進(jìn)行計算，近似得到了CP與λ、θ之間的連續(xù)變化關(guān)系如圖4所示.曲面最高點(diǎn)處λopt=4.02，θopt=4.0°，CP-opt=0.482，攻角αopt=5.9°，該點(diǎn)即為最佳工作點(diǎn)，αopt即為最佳攻角(opt下標(biāo)表示最佳工作點(diǎn)).

圖4 CP、λ、θ曲面圖Fig.4 Surface diagram of CP，λ and θ

進(jìn)一步將采用最佳工作點(diǎn)與采用最大升阻比工作點(diǎn)的CP值做比較.在圖4中分別提取攻角α等于最佳攻角(α=αopt=5.9°)和最大升阻比攻角(α=αmldr=4.2°)的各數(shù)據(jù)點(diǎn)，構(gòu)成兩條等攻角曲線，分別稱之為最佳攻角曲線和最大升阻比攻角曲線，兩者對比圖如圖5所示.

圖5 最佳攻角與最大升阻比攻角時CP對比Fig.5 Comparison of CP at optimum attack angle andattack angle with maximum lift-drag ratio

由圖5可見，在其它因素相同的條件下，采用最佳攻角比采用最大升阻比攻角的CP值顯著提高.兩曲線CP值均在θ=4°時達(dá)到極大值，兩極值點(diǎn)分別稱為最大升阻比工作點(diǎn)和最佳工作點(diǎn)，對應(yīng)CP值分別為0.462 0和0.482 1.可見最佳工作點(diǎn)比采用最大升阻比工作點(diǎn)時，CP提高了4.3%.

計算結(jié)果同時表明，采用最佳工作點(diǎn)時，當(dāng)?shù)厮俣缺圈擞稍捎米畲笊璞裙ぷ鼽c(diǎn)時的4.7降為4.0，減小15%.這說明在不改變風(fēng)輪其它運(yùn)行參數(shù)的條件下，轉(zhuǎn)速降低15%，可使風(fēng)力機(jī)由最大升阻比工作點(diǎn)調(diào)整到最佳工作點(diǎn).

3.2 NACA0012翼型風(fēng)力機(jī)的β優(yōu)化

圖6 α、β獨(dú)立優(yōu)化曲線簇Fig.6 Curve cluster of after independent optimization of α and β

4 α、β的統(tǒng)一優(yōu)化

為實(shí)現(xiàn)優(yōu)化方案2中對α、β的統(tǒng)一優(yōu)化，將β、γ、λc、θ四參數(shù)在迭代程序中進(jìn)行統(tǒng)一考慮.程序流程圖如圖7所示.

圖7 α、β統(tǒng)一優(yōu)化程序流程圖Fig.7 Flow chart of procedure for unified optimization of α and β

圖8 α、β統(tǒng)一優(yōu)化曲線簇Fig.8 Curve cluster of after unified optimization of α and β

5 結(jié) 論

分別采用對α、β的獨(dú)立優(yōu)化方案和統(tǒng)一優(yōu)化方案對NACA0012翼型的1.5 MW風(fēng)力機(jī)α、β角進(jìn)行實(shí)例優(yōu)化，結(jié)果表明：

1) 采用α最佳工作點(diǎn)攻角與采用最大升阻比工作點(diǎn)攻角相比，在其它因素相同的情況下，風(fēng)力機(jī)Cp提高4.3%；在α最佳工作點(diǎn)處繼續(xù)對β優(yōu)化，Cp再提高4.8%.

α、β具有弱相關(guān)性，理論上宜采用統(tǒng)一優(yōu)化方案，但由于Cp提升率差異不大，故從工程方便角度考慮，仍可以采取獨(dú)立優(yōu)化方案.