风力机MATLAB设计程序.docx
《风力机MATLAB设计程序.docx》由会员分享，可在线阅读，更多相关《风力机MATLAB设计程序.docx（14页珍藏版）》请在冰豆网上搜索。
风力机MATLAB设计程序
makedata
%根据profili导出到翼型性能数据Excel表格生成翼型的结构体
clearairfoil;
forn=2:
100%sheetn
airfoil（n-1）.Re=22;
%%%%找到sheetn截面翼型雷诺数所在行数nRe%%%%
try
[num,txt,~]=xlsread（'yxdata.xlsx',n,'A1:
I5000'）;
catch
break;
end
nstr=strfind（txt,'Re='）;%找到每一个雷诺数翼型的起始记录位置
k=1;%nRe的变量
fori=1:
length（nstr）%行数从第一行到最后一行开始判断
j=nstr{i};%将第i行的值赋给临时变量j
ifj%如果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=txt{airfoil（n-1）.nRe
（1）}（1:
（nstr{airfoil（n-1）.nRe
（1）}-4））;%将第i行的值赋给临时变量j
fori=1:
length（airfoil（n-1）.nRe）%行数从第一行到最后一行开始判断
airfoil（n-1）.Re（i）=str2double（txt{airfoil（n-1）.nRe（i）}（（nstr{airfoil（n-1）.nRe（i）}+5）:
length（txt{airfoil（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）;
catch
end
%%%%读入截面翼型拟合各雷诺数性能曲线和其它数据%%%%
fori=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）;
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）;
end
end
saveairfoildata
qd
clc;clear;
filename='name';
load（filename）
%loadxc
loadairfoilDataairfoil
pi=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;
whileabs（Hhub-temp）>2
Vr=Ve*（Hhub/10）^fq;rou=Pa*1000/（（273+T）*qR）;
%Vr设计风速；rou空气密度
D=（8*Pr/（Cp*DJ_eta*CD_eta*rou*Vr^3*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）
ifMa<0.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;%迭代次数
%%%%%%%%%%%%%%%%%%%%%%%粗略计算雷诺数初始化迭代%%%%%%%%%%%%%%%%%%%%%%%%%
whilej<=dn
%%%%%%%%%%%%%%%%%%%%%%%多次迭代直到雷诺数基本不变%%%%%%%%%%%%%%%%%%%%%%%%%
Wi=（（1-a）.*（1-a）*Vr^2+（1+b）.*（1+b）.*（omiR*ri）.^2）.^0.5;%断面上的合成速度。
Re（:
j）=（rou*Wi.*c）/u;%断面的雷诺数
fori=1:
n%截面
%%%%%%%%%%%%%%%读入数A据%%%%%%%%%%%%%%%
%根部s818中部s830尖部s817
forjj=1:
length（airfoil）%翼型
if（strcmp（name{i,iname},airfoil（jj）.name））
temp=find（abs（Re（i,j）-airfoil
（1）.Re）==min（abs（Re（i,j）-airfoil
（1）.Re）））;%将雷诺数相近的翼型数据赋给断面
SctRe（:
i）=airfoil（jj）.Re（:
temp）;
zAlf（i,:
）=airfoil（jj）.zAlf（:
temp）;
zCl（i,:
）=airfoil（jj）.zCl（:
temp）;
zCd（i,:
）=airfoil（jj）.zCd（:
temp）;
Thick（i,:
）=max（airfoil
（1）.y）-min（airfoil
（1）.y）;
end
end
%zAlf攻角；zCl升力系数；zCd阻力系数；Thick翼型最大厚度
tempa=0.6;tempb=0.1;
%a轴向诱导；b周向诱导
kk=0;
whileabs（a（i,1）-tempa）>0.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
升级会员 