风力机MATLAB设计程序.docx

1、风力机MATLAB设计程序makedata%根据profili导出到翼型性能数据Excel表格生成翼型的结构体clear airfoil;for n=2:100 %sheet n airfoil(n-1).Re=22; %找到sheet n 截面翼型雷诺数所在行数nRe% try num,txt, = xlsread(yxdata.xlsx,n,A1:I5000); catch break; end nstr=strfind(txt,Re = );%找到每一个雷诺数翼型的起始记录位置 k=1; %nRe的变量 for i=1:length(nstr)%行数 从第一行到最后一行开始判断 j=ns

2、tri; %将第i行的值赋给临时变量j if j %如果j存在则将行数给nRe airfoil(n-1).nRe(k)=i; k=k+1; end airfoil(n-1).nRe(k)=length(num)+5; end %找到翼型相近的雷诺数下的性能数据和name% k=length(airfoil(n-1).nRe); airfoil(n-1).Re(k)=0; airfoil(n-1).name=txtairfoil(n-1).nRe(1)(1:(nstrairfoil(n-1).nRe(1)-4); %将第i行的值赋给临时变量j for i=1:length(airfoil(n-1

3、).nRe) %行数 从第一行到最后一行开始判断 airfoil(n-1).Re(i)=str2double(txtairfoil(n-1).nRe(i)(nstrairfoil(n-1).nRe(i)+5):length(txtairfoil(n-1).nRe(i); %将第i行的值赋给临时变量j end % try wnum, = xlsread(yxdata.xlsx,n,H1:I500); lth=length(wnum); airfoil(n-1).x(1:lth,1)=wnum(1:lth,1); airfoil(n-1).y(1:lth,1)=wnum(1:lth,2); cat

4、ch end %读入截面翼型拟合各雷诺数性能曲线和其它数据% for i=1:(length(airfoil(n-1).nRe)-1) %Re temp=(airfoil(n-1).nRe(i):(airfoil(n-1).nRe(i+1)-5); lth=length(temp); airfoil(n-1).Alf(1:lth,i)=num(temp,1); airfoil(n-1).Cl(1:lth,i)=num(temp,2); airfoil(n-1).Cd(1:lth,i)=num(temp,3); airfoil(n-1).ClCd(1:lth,i)=num(temp,4); te

5、mpn=find(airfoil(n-1).ClCd(:,i)=max(airfoil(n-1).ClCd(:,i); airfoil(n-1).zAlf(i)=airfoil(n-1).Alf(tempn,i); airfoil(n-1).zCl(i)=airfoil(n-1).Cl(tempn,i); airfoil(n-1).zCd(i)=airfoil(n-1).Cd(tempn,i); airfoil(n-1).xCl(:,i) airfoil(n-1).SxCl(:,i) = polyfit(airfoil(n-1).Alf(:,i),airfoil(n-1).Cl(:,i),6)

6、; airfoil(n-1).xCd(:,i) airfoil(n-1).SxCd(:,i) = polyfit(airfoil(n-1).Alf(:,i),airfoil(n-1).Cd(:,i),6); endendsave airfoildataqdclc;clear;filename=name;load(filename)%load xcload airfoilData airfoilpi=3.141592653;qR=287.64;k=1.4;fq=0.12;u=1.698e-05;%pi 圆周率 ；qR气体常数；k 等商指数；fq 风切指数；u 动力粘度；Pr=1200000;Ve

7、=8.5;Pa=85.8;T=15;B=3;DJ_eta=0.95;CD_eta=0.95;%Pr 额定功率；Vr 额定风速；Pa 风场平均压强；T 平均气温；B 叶片数%DJ_eta 电机效率；CD_eta传动效率；%tempV=70;Cp=0.43;n=30;namR=9;BL1=0.15;BL2=0.05;%Cp 风能利用系数；n 等分段数；namR ;叶尖速比;BL1 叶根园比例；BL1 轮毂园比例；min_n=900;max_n=1950;e_n=1620;%发电机的转速范围iname=1;%开始迭代计算轮毂高度%Hhub=95;temp=0;while abs(Hhub-temp)

8、2 Vr=Ve*(Hhub/10)fq;rou=Pa*1000/(273+T)*qR); %Vr 设计风速；rou 空气密度 D=(8*Pr/(Cp*DJ_eta*CD_eta*rou*Vr3*pi)0.5; D1=floor(D);%取比圆整风轮直径 向上取 对Cp和功率的大小又决定性作用。对 R=D1/2; temp=Hhub; %D 风轮直径；R 风轮半径 Hhub=ceil(0.85*D1);R0=BL1*R;Lb=R-R0;dr=Lb/n;Rhub=BL2*R;% 轮毂半径 %Hhub 圆整轮毂高度； R0 叶根园半径；Lb 叶片有效长度；dr 每段长度； %str=sprintf(

9、风轮直径 %f，圆整风轮直径 %f,风轮半径R %f, 轮毂高度 %f,D,D1,R,Hhub); %disp(str)end%结束迭代计算轮毂高度%判断是否可压%C=(qR*k*(T+273)0.5;%当地声速Vz=Ve*(Hhub+R)/10)fq;Vh=Vz*namR*cos(5*pi/180);Ma=Vh/C;%Vz 最高处的风速；Vh 风轮上的最高风速；Ma 马赫数str=sprintf(namR= %f ，圆整轮毂高度 %f，设计风速 %f， 最高处的风速 %f,namR, Hhub,Vr,Vz);disp(str)str=sprintf(风轮上的最高风速 %f， 当地声速 %f，

10、马赫数 %f,Vh,C,Ma);disp(str)if Ma0.3 disp(放心！不可压);else disp(可压请修正！);end%粗略计算雷诺数初始化迭代%omiR=Vr*namR/R;%主轴角速度nz=60*omiR/(2*pi);%主轴转速CD=e_n/nz;c(1:n,1)=2;%初始化弦长i=(1:n);ri=R0+(i-0.5)*dr;%断面半径namri=namR*ri/R;%周速比Phi(1:n,1)=0;Cn(1:n,1)=0;Ct(1:n,1)=0; F(1:n,1)=0;a(1:n,1)=0.7;b(1:n,1)=0; dd(1:n,1)=0;j=1;dn=4; %

11、迭代次数%粗略计算雷诺数初始化迭代%while j0.0001|abs(b(i,1)-tempb)0.0001 Phi(i,1)=atan(1-a(i,1)/(1+b(i,1)*namri(i); ft=(2*acos(exp(-1*B*(R-ri(i)/(2*R*sin(Phi(i,1)/pi; %Ft fh=(2*acos(exp(-1*B*(ri(i)-R0)/(2*R0*sin(Phi(i,1)/pi; %Fh F(i,1)=ft*fh; c(i,1)=8*pi*ri(i)*a(i,1)*F(i,1)*(1-a(i,1)*F(i,1)*sin(Phi(i,1)*sin(Phi(i,1)/. (1-a(i,1)*(1-a(i,1)*B*zCl(i)*cos(Phi(i,1); Cn(i,1)=zCl(i)*cos(Phi(i,1)+zCd(i)*sin(Phi(i,1); Ct(i,1)=zCl(i)*sin(Phi(i,1)-zCd(i)*cos(Phi(i,1); cgm=B*c(i,1)/(2*pi*ri