风力机课程设计一种简单叶片的设计过程Hank一、设计任务1）基于叶素和动量理论设计设计水平轴风力机叶片；2）绘制风力机叶片弦长随叶片展向长度的变化曲线；3）绘制风力机叶片扭角随叶片展向长度的变化曲线；4）绘制设计风力机的性能曲线；5）绘制设计叶片的图纸，以及各位置的叶片翼型结构图纸；6）编写设计说明书，并附上必要的计算公式二、常用符号Pd 设计功率U∞设计风速W 合成流度ω风轮旋转角速度ρ空气密度c 弦长λ叶尖速比D 风轮直径R 风轮半径σ叶片弦长实度rφ叶素倾角β叶素安装角α叶素攻角Cl 升力系数Cd 阻力系数M 总力矩P 总功率Cp 风能利用系数a 轴向诱导因子b 切向诱导因子三、 基础理论切向速度： (1-1) 轴向速度：(1-2)叶素合成流速:(1-3) 入流角：(1)sin U a Wφ∞-=(1-4a)(1)cos 2r b aωφ+=(1-4b) 攻角：φα=-β(1-5) 垂直分力系数：sin cos y l d C C C φφ=- (1-6) 水平分力系数：cos sin x l d C C C φφ=+(1-7) 迭代方程：22214sin 4sin r r x y a C C a φφσσ=(-)-(1-8)14sin cos r yb C b φφσ=+(1-9)叶片弦长实度： 22r N c N cR R σ==ππμ(1-10)叶素单位圆环扇面：24()(1)M U r b a r r ωδ∞δ=πρ-(1-11)风能利用系数表达式： 3212p PC U R ρ∞=π(1-12)翼型与尖速比的关系：2r l C σλ=(1-13)0(1b)y V r ω=+0(1)x V U a ∞=-W ==风轮直径设计：D=(1-14)四、 设计流程图设计参数：10/U m s ∞= 6λ= 100d P kW =, 给定 p C = 0.44图 1-1a 叶片设计流程图图 1-1b 叶片设计流程图五、 计算说明书1. 估算D 、R 、ω根据公式（1-14），带入数据可得：21.157D m =/210.578R D m ==5.672/U rad s Rλω∞==2/6054.164/min n r πω==2. 翼型及其气动特性的确定选择翼型NACA4415，其最佳攻角及此时的l C 、d C 如下：06α=1.1119l C = 0.0109d C =表1-1 Polars of NACA4415拟合出l C 、d C 关于α的曲线如下32= -0.00016097 - 0.00060958 + 0.11374+0.48717l C ααα(1-15)6352= -5.85710 + 8.147910 - 0.000333840.0081917 d C ααα--⨯⨯+(1-16)3. 弦长分布的确定根据经验公式有：169l c C N πλμ=对其进行修正如下：0.416/9l c R C N πμλμ=(0.8)c ≤ (1-17) 3/0.0540.6050.0303c R μμ=-++(0.8)c ≤(1-18)M/μ分布曲线4.迭代求解a、b、β、d d过程如下：图1-2 方程迭代流程图相关曲线及表格如下：0.10.20.30.40.50.60.70.80.91μ速度诱导因子图 1-3 速度诱导因子 -μ0.10.20.30.40.50.60.70.80.91μc /R图 1-4 c/R -μ0.10.20.30.40.50.60.70.80.91μβ (︒)图 1-5 β -μ0.10.20.30.40.50.60.70.80.914μd M /d μ图 1-6dMd μ- μ表1-2 叶片参数表结果： 在给定的设计参数下 0.4426p C =3267.055168.99150.7149.315βμμμ=-+-+(1-19)25860.4458571285.2dMd μμμ=-+-(1-20)5. 求解不同λ值下的p C给定上面设计的叶片和来流风速U ∞，带入不同的λ可求出响应的p C ，进而得到曲线λ-p C ：5.56 6.577.580.350.40.450.50.55λC p图1-7λ - p C六、 附录1. 叶片立体图2.翼型截面(c = 200.00mm)3.程序代码M/μ分布曲线a.迭代求解a、b、β、d dclcclear% 风轮直径和转速设计Pd = 100e3; %设计功率ro = 1.293; %空气密度Cp = 0.44; %风能利用系数U = 10; %设计风速lemta = 6; %叶尖速比D = ((8*Pd)/(pi*ro*Cp*U^3))^.5 %风轮直径R = D/2Omega = lemta*U/R %设计角速度n = 60*U*lemta/(pi*D) %设计转速%最佳升阻比下的参数Cl = 1.1119;Cd = 0.0109;alfa = 6;K3 = [-0.054 0 0.0605 0.0313];%[-0.051 0 0.0605 0.0313]i = 1;for miu = 0.1:0.02:1.0a = 0.3;b = 0.1;C = lemta*miu;if(miu>=0.06 & miu<=0.8)c = 16*pi*R/(9*Cl*3*C*lemta)*miu^0.4;endif(miu>0 &miu<0.06)c = (-2*miu^2 - 0.5*miu + 0.25)*R;endif(miu>0.8 & miu<1)c = polyval(K3,miu)*R; %(-0.051*miu^3 + 0.0605*miu + 0.0313)endsd = 3*c/(2*pi*miu*R);count = 1;while(count<1000)A = 1-a;B = 1+b;phi = atan(A/C/B);sf = sin(phi);cf = cos(phi);Cx = Cl*cf+Cd*sf;Cy = Cl*sf-Cd*cf;E = sd/4/sf^2;r1 = E*(Cx-E*Cy*Cy);a_new = r1/(1+r1);r2 = E*tan(phi)*Cy;b_new = r2/(1-r2);exp = (abs(b/b_new-1)<1e-2) & (abs(a/a_new-1)<1e-2);if(exp)break;elseb = b_new;a = a_new;endcount = count + 1;endx(i) = miu;ax(i) = a;bx(i) = b;cx(i) = c;A = 1-a;C = lemta*miu;B = 1+b;phi = atan(A/C/B);beta(i) = phi*180/pi - alfa;sf = sin(phi);F = 0.5*ro*U*U*pi*R*R*R*C*miu;M = 8*b*A*miu;N = A*3*c*Cd*B/(sf*pi*R);M_u(i) = F*(M-N);i = i + 1;endfigure(1)hold ongrid onxlabel('\mu')ylabel('速度诱导因子');plot(x,ax,'r',x,bx,'b');figure(2)hold ongrid onxlabel('\mu');ylabel('c/R');plot(x,cx/R);figure(3)hold ongrid onxlabel('\mu');ylabel('\beta (\circ)');AB = polyfit(x,beta,3);z = polyval(AB,x);plot(x,beta,'*',x,z)figure(4)A = polyfit(x,M_u,2)y = polyval(A,x);hold ongrid onxlabel('\mu');ylabel('dM/d\mu');plot(x,M_u,'+',x,y)B = [A(1)*2/3, A(2)*1/2, A(3),0] z = polyval(B,[0.0,1.0])P = (z(2)-z(1))*OmegaCp = P/(0.5*ro*U*U*U*pi*R*R)。