风力机matlab设计程序Word文档下载推荐.docx
《风力机matlab设计程序Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《风力机matlab设计程序Word文档下载推荐.docx（13页珍藏版）》请在冰豆网上搜索。
（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（:
saveairfoildata
qd
clc;
clear;
filename='
name'
;
load（filename）
%loadxc
loadairfoilDataairfoil
pi=3.3;
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）
%%%%%%%%%%%%%%%%%%%%%%%%%结束迭代计算轮毂高度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%判断是否可压%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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）
风轮上的最高风速%f，当地声速%f，马赫数%f'
Vh,C,Ma）;
ifMa<
0.3
disp（'
放心！
！
不可压'
else
可压请修正！
'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%粗略计算雷诺数初始化迭代%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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:
Ct（1:
F（1:
a（1:
n,1）=0.7;
b（1:
dd（1:
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;
%断面的雷诺数
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（:
zCl（i,:
）=airfoil（jj）.zCl（:
zCd（i,:
）=airfoil（jj）.zCd（:
Thick（i,:
）=max（airfoil
（1）.y）-min（airfoil
（1）.y）;
%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（i））;
ifa（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））;
g1=cgm*Cn（i,1）/（4*F（i,1）*sin（Phi（i,1））*sin（Phi（i,1）））;
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;
ifkk>
10000
break
dd（i,j）=kk;
%迭代次
j=j+1;
%%%%%%%%%%%%%%%%%%%%开始修型计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tc=c;
cs=4;
%cse=1多项式修订，2smooth修正，3手动修正，4功率修正，
h=figure;
switchcs
case1
tempc=polyfit（ri,c,6）;
c=polyval（tempc,ri）;
bx;
case2
c=smooth（c）;
%xc自动修改弦长
bx
case3
c=c-xc;
case4
Pe=0;
c=smooth（c）
whileabs（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*SM
Pe=P*DJ_eta*CD_eta
%SM总扭矩;
SF总轴向推力；
P功率
PF=0.5*rou*Vr^3*pi*R^2;
%风能
Cp=P/PF%风能利用系数
plot（ri,c,'
-r'
holdon
plot（ri,tc）;
zAlf=zAlf*pi/180;
theta=Phi-zAlf;
thetatheta=theta*180/pi;
%Phi倾角度，theta安装角度，zAlf攻角
Thick=Thick.*c;
filename=[name{1,iname}name{end,iname}'
qd'
num2str（iname）];
save（filename）
savetemp
print（h,'
-dbmp'
filename）;
xn
clc,clear;
t1=clock;
loadtemp
loadAirfoilDataairfoil
%%%%%%%%%%%%%%%%%%%%%功率控制计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
%弧度
forj=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）<
min_omiR
omiRj（j,1）=min_omiR;
elseifomiRj（j,1）>
max_omiR
omiRj（j,1）=max_omiR;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
namRj（j,1）=omiRj（j,1）*R/Vi（j,1）;
%主轴角速度
xCn（1:
n,j）=0;
xCt（1:
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.5
pitch（j,1）=pitch（j-1,1）;
pitch（j,1）=0*pi/180;
%%%%%%%%%%%%%当调整pitch到Cp=cPmax结束%%%%%%%%%%%%%%%%%%%%%%%%%%
ifVf（j）<
=Ve
mark1=0;
mark2=1;
elseifVf（j）>
Ve
mark1=1;
mark2=0;
while（mark1||（abs（xCp（j,1）-Cpmax）>
=0.1））&
&
（（（abs（Pj（j,1）-Pr）>
=0.01*Pr）||isnan（Pj（j,1）））||mark2）%相差0.05是结束
%whilexCp（j,1）-temp_Cp>
0%相差0.005是结束
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;
if（strcmp（name{i,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;
xSctRe（:
i,j）=airfoil（jj）.Re（:
sCl{:
i}=airfoil（jj）.Cl（:
sCd{:
i}=airfoil（jj）.Cd（:
sAlf{:
i}=airfoil（jj）.Alf（:
Alf=1;
%清空Alf
Cl=1;
Cd=1;
len=length（sAlf{:
i}）;
Alf（1:
len,1）=sAlf{:
i};
Cl（1:
len,1）=sCl{:
Cd（1:
len,1）=sCd{:
tempAlf=-20;
tempi=1;
whiletempi<
=length（Alf）%修正翼型的数据中错误数据
ifAlf（tempi,1）<
tempAlf
Alf（tempi）='
Cl（tempi）='
Cd（tempi）='
tempi=tempi-1;
ifisnan（Alf（tempi,1））
tempAlf=Alf（tempi,1）;
tempi=tempi+1;
tempa=0;
tempb=0;
%%%%%%%%%%%%%%%%%%%%结束匹配翼型性能数据%%%%%%%%%%%%%%
whileabs（xa（i,j）-tempa）>
0.001||abs（xb（i,j）-tempb）>
0.001
jsc=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;
f2=（2*acos（exp（-1*B*（ri（i）-R0）/（2*R0*sin（Phi（i,j））））））/pi;
F（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（:
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）>
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））;
g1=cgm（i,j）*xCn（i,j）/（4*F（i,j）*sin（Phi（i,j））*sin（Phi（i,j）））;
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）;
xM（:
j）=0.5*B*rou*xWi（:
j）.*xWi（:
j）.*c.*xCt（:
j）.*ri*dr;
Fn（:
j）.*c.*xCn（:
j）.*dr;
xSM（j,1）=sum（xM（:
j））;
xSF（j,1）=sum（Fn（:
xP（j,1）=omiRj（j,1）*xSM（j,1）;
Pj（j,1）=xP（j,1）*DJ_eta*CD_eta;
disp（['
Vf（j）='
num2str（Vf（j））]）
Pj（j,1）='
num2str（Pj（j,1））]）
xPF（j,1）=0.5*rou*Vi（j,1）^3*pi*R^2;
Cp_temp=xCp（j,1）;
xCp（j,1）=xP（j,1）/xPF（j,1）;
%%%%%%%%%%%%%%
pitch（j,1）=pitch（j,1）+0.005*pi/180;
pitch=
升级会员 