《风力机matlab设计程序》由会员分享，可在线阅读，更多相关《风力机matlab设计程序（21页珍藏版）》请在金锄头文库上搜索。
1、makedata%根据 profili 导出到翼型性能数据 Excel 表格生成翼型的结构体clear airfoil;for n=2:100 %sheet nairfoil(n-1).Re=22;%找到 sheet n 截面翼型雷诺数所在行数 nRe%trynum,txt, = xlsread(yxdata.xlsx,n,A1:I5000);catchbreak;endnstr=strfind(txt,Re = );%找到每一个雷诺数翼型的起始记录位置k=1; %nRe 的变量for i=1:length(nstr)%行数从第一行到最后一行开始判断j=nstri; %将第 i 行的值赋给临时
2、变量 jif 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 行的值赋给临时变量 jfor i=1:length(airfoil(n-1).nRe) %行数从第一行到最后一行
3、开始判断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%trywnum, = 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%读入截面翼型拟合各雷诺数性能曲线和其它数
4、据%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(
5、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
6、(:,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.9
7、5;%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);
8、%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);%dis
9、p(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.0001|abs(b(i,1)-tempb)0.
10、0001Phi(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(
11、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.3539g1=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);elseg1=cgm*Cn(i,1)/(4*F(i,1)*sin(Phi(i,1)*sin(Phi(i,1);endg2=cgm*Ct(i
12、,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; %迭代次endj=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 2c=smooth(c
13、); %xc 自动修改弦长bxcase 3c=c-xc;bxcase 4Pe=0;c=smooth(c)while abs(Pe-Pr)0.001*Prc=c-0.001M=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 onplo
14、t(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/1
