风力机matlab设计程序.docx

makedata%根据profili导出到翼型性能数据Excel表格生成翼型的结构体clear airfoil;for n=2:100 %sheet nairfoil(n-1).Re=22; %找到sheet n 截面翼型雷诺数所在行数nRe%try num,txt, = xlsread(yxdata.xlsx,n,A1:I5000);catchbreak;endnstr=strfind(txt,Re = );%找到每一个雷诺数翼型的起始记录位置 k=1; %nRe的变量 for i=1:length(nstr)%行数从第一行到最后一行开始判断 j=nstri; %将第i行的值赋给临时变量j if j %如果j存在则将行数给nReairfoil(n-1).nRe(k)=i; k=k+1;endairfoil(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).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行的值赋给临时变量jend %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);catchend %读入截面翼型拟合各雷诺数性能曲线和其它数据%for i=1:(length(airfoil(n-1).nRe)-1) %Retemp=(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);tempn=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); airfoil(n-1).xCd(:,i) airfoil(n-1).SxCd(:,i) = polyfit(airfoil(n-1).Alf(:,i),airfoil(n-1).Cd(:,i),6);endendsaveairfoildataqdclc;clear;filename=name;load(filename)%load xcloadairfoilData airfoilpi=3.141592653;qR=287.64;k=1.4;fq=0.12;u=1.698e-05;%pi 圆周率；qR气体常数；k 等商指数；fq风切指数；u 动力粘度；Pr=1200000;Ve=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)2Vr=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(风轮直径 %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，马赫数 %f,Vh,C,Ma);disp(str)if Ma0.3disp(放心！不可压);elsedisp(可压请修正！);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; %迭代次数%粗略计算雷诺数初始化迭代%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; %Ftfh=(2*acos(exp(-1*B*(ri(i)-R0)/(2*R0*sin(Phi(i,1)/pi; %FhF(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(i);if a(i,1) 0.3539 g1=cgm*Cn(i,1)/(4*F(i,1)*sin(Phi(i,1)*sin(Phi(i,1)*4. *a(i,1)*(1-a(i,1)/(0.6+0.61*a(i,1)+0.79*a(i,1)*a(i,1);else g1=cgm*Cn(i,1)/(4*F(i,1)*sin(Phi(i,1)*sin(Phi(i,1);end g2=cgm*Ct(i,1)/(4*F(i,1)*sin(Phi(i,1)*cos(Phi(i,1);tempa=a(i,1);tempb=b(i,1);a(i,1)=g1/(1+g1);b(i,1)=g2/(1-g2);kk=kk+1;ifkk10000breakendenddd(i,j)=kk; %迭代次end j=j+1;end%开始修型计算%tc=c; cs=4; %cse =1 多项式修订，2 smooth 修正，3 手动修正，4 功率修正，h=figure;switchcscase 1tempc=polyfit(ri,c,6); c=polyval(tempc,ri);bx;case 2 c=smooth(c); %xc 自动修改弦长bxcase 3 c=c-xc;bxcase 4Pe=0; c=smooth(c)while abs(Pe-Pr)0.001*Pr c=c-0.001 M=0.5*B*rou*Wi.*Wi.*c.*Ct.*ri.*dr; % 各段扭矩Fn=0.5*B*rou*Wi.*Wi.*c.*Cn.*dr; %各段的轴向力 SM=sum(M);SF=sum(Fn); P=omiR*SMPe=P*DJ_eta*CD_eta %SM 总扭矩; SF 总轴向推力；P 功率 PF=0.5*rou*Vr3*pi*R2; %风能Cp=P/PF %风能利用系数plot(ri,c,-r);hold onplot(ri,tc);endendzAlf=zAlf*pi/180;theta=Phi-zAlf;thetatheta=theta*180/pi; %Phi 倾角度，theta 安装角度，zAlf攻角Thick=Thick.*c;filename = name1,iname nameend,iname qd num2str(iname);save(filename)save tempfilename = name1,iname nameend,iname qd num2str(iname);print(h, -dbmp, filename);xnclc,clear;t1=clock;load temploadAirfoilData airfoil%功率控制计算%Cpmax=0.43;Vci=3.5;Vst=0.5;Vco=19; %Vci切入风速；Vst风速变化步长；Vco切出风速Vf=(Vci:Vst:Vco); jsc=0;%总循环次数统计 % 风速变化范围% Phi=Phi*pi/3180; % 弧度7% theta=theta*pi/180; % 弧度for j=1:length(Vf); % jj从切入风速开始 Vi(j,1)=Vf(j)*(Hhub/10)fq; %轮毂处的变风速omiRj(j,1)=Vi(j,1)*namR/R; %不同风速下的转速min_omiR=min_n/nz/CD*omiR; %最小叶尖速比max_omiR=max_n/nz/CD*omiR; %最大叶尖速比 %调整叶尖速比的范围%ifomiRj(j,1)max_omiRomiRj(j,1)=max_omiR;end %namRj(j,1)=omiRj(j,1)*R/Vi(j,1); %主轴角速度xCn(1:n,j)=0;xCt(1:n,j)=0; F(1:n,j)=0;xa(1:n,j)=a;xb(1:n,j)=b;xCp(j,1)=0;ij=0;temp_Cp=0;i_Cp=1;Pj(j,1)=1;% pitch(j,1)=0*pi/180;ifVf(j)Ve+0.5pitch(j,1)=pitch(j-1,1);elsepitch(j,1)=0*pi/180;end %当调整pitch到Cp=cPmax结束%ifVf(j)Ve mark1=1;mark2=0;end while (mark1|(abs(xCp(j,1)-Cpmax)=0.1)&(abs(Pj(j,1)-Pr)=0.01*Pr)|isnan(Pj(j,1)|mark2) %相差0.05是结束 % while xCp(j,1)-temp_Cp0 %相差0.005是结束 for i=1:n % 第i处截面xnamri(i,j)=namRj(j,1)*ri(i)/R; %不通风速下各截面的周速比 xWi(i,j)=(1-xa(i,j)*(1-xa(i,j)*Vi(j,1)2+(1+xb(i,j)*(1+xb(i,j)*(omiRj(j,1)*ri(i)2)0.5;%断面上的合成速度。xRe(i,j)=(rou*xWi(i,j)* c(i)/u;%断面的雷诺数 %匹配翼型性能数据% marke1=0; for jj=1:length(airfoil) %翼型 if(strcmp(namei,1,airfoil(jj).name) %找到名字相同翼型 temp=find(abs(xRe(i,j)-airfoil(1).Re)=min(abs(xRe(i,j)-airfoil(1).Re);%将雷诺数相近的翼型数据赋给断面ifisempty(temp) marke1=1;continue;endxSctRe(:,i,j)=airfoil(jj).Re(:,temp);sCl:,i=airfoil(jj).Cl(:,temp);sCd:,i=airfoil(jj).Cd(:,temp);sAlf:,i=airfoil(jj).Alf(:,temp);break;endend Alf=1; %清空AlfCl=1; Cd=1;len=length(sAlf:,i);Alf(1:len,1)=sAlf:,i;Cl(1:len,1)=sCl:,i;Cd(1:len,1)=sCd:,i;tempAlf=-20; tempi=1; while tempi= length(Alf) %修正翼型的数据中错误数据if Alf(tempi,1)0.001|abs(xb(i,j)-tempb)0.001jsc=jsc+1; Phi(i,j)=atan(1-xa(i,j)/(1+xb(i,j)*xnamri(i,j); f1=(2*acos(exp(-1*B*(R-ri(i)/(2*R*sin(Phi(i,j)/pi; %Ft f2=(2*acos(exp(-1*B*(ri(i)-R0)/(2*R0*sin(Phi(i,j)/pi; %FhF(i,j)=f1*f2;xALf(i,j)= (Phi(i,j)-theta(i,1)-pitch(j,1)*180/pi; %xAlf攻角; Phi 倾角; theta 安装角; pitch 变桨角nCl(i,j)=interp1(Alf(:,1),Cl(:,1),xALf(i,j),pchip); %这是引起误差的关键因素nCd(i,j)=interp1(Alf(:,1),Cd(:,1),xALf(i,j),pchip);xCn(i,j)=nCl(i,j)*cos(Phi(i,j)+nCd(i,j)*sin(Phi(i,j);xCt(i,j)=nCl(i,j)*sin(Phi(i,j)-nCd(i,j)*cos(Phi(i,j);cgm(i,j)=B*c(i,1)/(2*pi*ri(i);ifxa(i,j) 0.3539 g1=cgm(i,j)*xCn(i,j)/(4*F(i,j)*sin(Phi(i,j)*sin(Phi(i,j)*4. *xa(i,j)*(1-xa(i,j)/(0.6+0.61*xa(i,j)+0.79*xa(i,j)*xa(i,j);else g1=cgm(i,j)*xCn(i,j)/(4*F(i,j)*sin(Phi(i,j)*sin(Phi(i,j);end g2=cgm(i,j)*xCt(i,j)/(4*F(i,j)*sin(Phi(i,j)*cos(Phi(i,j);tempa=xa(i,j);tempb=xb(i,j);xa(i,j)=g1/(1+g1);xb(i,j)=g2/(1-g2);endendxM(:,j)=0.5*B*rou*xWi(:,j).*xWi(:,j).*c.*xCt(:,j).*ri*dr; % 各段扭矩Fn(:,j)=0.5*B*rou*xWi(:,j).*xWi(:,j).*c.*xCn(:,j).*dr; %各段的轴向力xSM(j,1)=sum(xM(:,j);xSF(j,1)=sum(Fn(:,j);xP(j,1)= omiRj(j,1)*xSM(j,1); %SM 总扭矩; SF 总轴向推力；P 功率Pj(j,1)=xP(j,1)*DJ_eta*CD_eta;disp