赵 骞， 厉 伟， 姚兴佳， 邵一川， 郭庆鼎

(1. 沈阳工业大学 电气工程学院， 沈阳 110870； 2. 沈阳大学 信息工程学院， 沈阳 110004)

本文对风力机叶片攻角α及风轮锥角β进行优化，并对α与β之间的相关性做进一步研究，为进一步设计新型更高效风力机提供最佳的优化方案.采用两套优化方案对NACA0012翼型的1.5 MW风力机α、β角进行实例优化，方案1对α、β进行独立优化：首先应用传统叶素动量(BEM)理论对叶片攻角α进行优化，确定α的最佳工作点，再在该最佳工作点处，应用对β锥角修正所得的β-BEM理论进一步对锥角β进行优化，计算优化后风能利用系数Cp1.方案2对α、β角进行统一优化：应用β-BEM理论同时对α、β角进行优化，计算优化后风能利用系数Cp2.通过对Cp1与Cp2进行对比分析，研究α、β之间相关性，并确定最佳优化方案.

1 传统叶素动量理论

1.1 贝茨理论模型

图1为Betz理论的气流图.在风轮上游远端处，风速为υ0，气流截面积为S0；风轮处，风速为υ，气流截面积为S；风轮下游远端处，风速为υ1，气流截面积为S1.由Betz理论[6]可知，当满足式(1)，即

(1)

(2)

式中，ρ为空气密度.风力机工作于低风速环境下，此时气流视为不可压缩流体，式(1)结合连续性方程S0υ0=Sυ=S1υ1得

S0∶S∶S1=1∶2∶3

(3)

式(3)说明风轮上、中、下游的气流截面积逐渐增大，表明在风轮处气流除具有平行于风轮旋转轴的轴向速度外，同时具有垂直于旋转轴的径向速度分量，这为进一步引入风轮锥角β提供了依据.

图1 贝茨理论风轮气流图Fig.1 Airflow of wind wheel in Betz theory

1.2 传统叶素动量理论

叶素动量理论一直是研究风力机气动特性的有效方法，相关研究在文献中多有报道[7-9].传统的叶素动量理论采用零锥角模型(β=0°)，即风轮旋转面与旋转轴垂直，此时风轮处气流的径向速度分量与旋转面平行，对风轮不产生转矩作用，因此只考虑气流轴向速度对风轮产生的转矩作用.

由动量理论可知，气流作用在风轮根部距离r处时叶素的轴向推力dFn和转矩dM分别为

(4)

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

(5)

式中：a为轴向诱导因子；a′为周向诱导因子；ω为转速.同时，由叶素理论可得到

(6)

(7)

式中：Nb为叶片数；c为叶素弦长；CL、CD为风力机翼型升力、阻力系数；φ为入流角，其表达式为

(8)

式中，λ为当地速度比，其表达式为

(9)

对于某确定的翼型，CL、CD主要由攻角α决定，其表达式为

α=φ-θ

(10)

式中，θ为节距角.

分别将式(4)、(6)及式(5)、(7)建立等量关系，并引入普朗特叶尖损失因子f，可计算出轴向诱导因子a及周向诱导因子a′为

(11)

(12)

式中，σ为实度，其表达式为

(13)

进一步结合式(2)、(7)计算出风能利用系数为

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

(CLsinφ-CDcosφ)

(14)

给定当地速度比λ及节距角θ，通过迭代式(8)～(13)计算出a、a′的收敛值，并进一步由式(14)计算出风能利用系数Cp.

2 基于锥角β的BEM修正理论

传统BEM理论模型中没有涉及到锥角β，因此，若要对β实现优化，须对传统BEM理论进行β修正.图2为将β引入到传统叶素动量理论后的气流分布示意图，此时风轮旋转面为圆锥面，气流在与旋转轴垂直平面内的径向速度分量与风轮旋转圆锥面不平行，因此，径向速度分量亦将对风轮产生转矩作用，进而转换为风轮输出功率.

图2 β锥角时风轮气流分布Fig.2 Airflow distribution of wind wheel with cone angle β

如图2所示，取与叶根距离为r的风轮叶素，其距离旋转轴的垂直距离为

r⊥=rcosβ

(15)

气流速度与风轮旋转锥面垂直的速度分量υ⊥为

(16)

式中：γ为气流发散角，表示气流流线与旋转轴夹角；ac为对轴向诱导因子a的修正值.

将修正式(15)、(16)及β代入传统BEM理论表达式(4)～(14)中，得到修正后β-BEM理论表达式为

(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理论表达式(17)～(27)在风轮锥角β=0°时与修正前表达式(4)～(14)完全一致，说明β-BEM理论表达式具有较强的通用性.

3 α、β独立优化

3.1 NACA0012翼型风力机的α优化

为实现方案1中对α的优化，采用BEM理论对NACA0012翼型的1.5 MW风力机CP值进行数值计算.该翼型为NACA经典翼型，低风速时，该翼型的特性曲线如图3所示，在小攻角段，CL与α几乎成正比关系，在12°攻角时翼型发生失速，通常情况下，翼型均工作在失速点之前.

图3 NACA0012翼型特性曲线Fig.3 Characteristic curves of NACA0012 airfoil

选取翼型中部叶素作为计算区域时，不涉及叶尖损失，因此叶尖损失因子f取1.将λ、θ做为循环变量，对于每一组λ、θ，通过式(8)～(13)迭代求解各参数，进一步代入式(14)计算CP，最终得到CP随λ、θ变化的曲面图.使用C++语言编写该迭代程序，对6万组λ、θ点的CP值进行计算，近似得到了CP与λ、θ之间的连续变化关系如图4所示.曲面最高点处λopt=4.02，θopt=4.0°，CP-opt=0.482，攻角αopt=5.9°，该点即为最佳工作点，αopt即为最佳攻角(opt下标表示最佳工作点).

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

进一步将采用最佳工作点与采用最大升阻比工作点的CP值做比较.在图4中分别提取攻角α等于最佳攻角(α=αopt=5.9°)和最大升阻比攻角(α=αmldr=4.2°)的各数据点，构成两条等攻角曲线，分别称之为最佳攻角曲线和最大升阻比攻角曲线，两者对比图如图5所示.

图5 最佳攻角与最大升阻比攻角时CP对比Fig.5 Comparison of CP at optimum attack angle andattack angle with maximum lift-drag ratio

由图5可见，在其它因素相同的条件下，采用最佳攻角比采用最大升阻比攻角的CP值显著提高.两曲线CP值均在θ=4°时达到极大值，两极值点分别称为最大升阻比工作点和最佳工作点，对应CP值分别为0.462 0和0.482 1.可见最佳工作点比采用最大升阻比工作点时，CP提高了4.3%.

计算结果同时表明，采用最佳工作点时，当地速度比λ由原采用最大升阻比工作点时的4.7降为4.0，减小15%.这说明在不改变风轮其它运行参数的条件下，转速降低15%，可使风力机由最大升阻比工作点调整到最佳工作点.

3.2 NACA0012翼型风力机的β优化

图6 α、β独立优化曲线簇Fig.6 Curve cluster of after independent optimization of α and β

4 α、β的统一优化

为实现优化方案2中对α、β的统一优化，将β、γ、λc、θ四参数在迭代程序中进行统一考虑.程序流程图如图7所示.

图7 α、β统一优化程序流程图Fig.7 Flow chart of procedure for unified optimization of α and β

图8 α、β统一优化曲线簇Fig.8 Curve cluster of after unified optimization of α and β

5 结 论

分别采用对α、β的独立优化方案和统一优化方案对NACA0012翼型的1.5 MW风力机α、β角进行实例优化，结果表明：

1) 采用α最佳工作点攻角与采用最大升阻比工作点攻角相比，在其它因素相同的情况下，风力机Cp提高4.3%；在α最佳工作点处继续对β优化，Cp再提高4.8%.

α、β具有弱相关性，理论上宜采用统一优化方案，但由于Cp提升率差异不大，故从工程方便角度考虑，仍可以采取独立优化方案.