适合高亚音速和跨音速的压气机叶型生成源代码
导读这个程序是基于Wennerstrom的书“Design of highly loaded axial-flow fans and compressors”中建议的叶型编写的，比较适合于高亚音速和跨音速，比如航空发动机中的高压压气机设计。基于这个叶型设计的压气机是可以有较高性能的，不过目前最先进的压气机中采用的预压缩叶型和可控扩散叶型要更好一点，当然设计难度也大。
叶型简介
这里偷一下懒，把原书中的几页直接放出来，给大家参考，我就不用整理了，哈哈。
需要说明的是，书中的中弧线形状是要根据气动设计给出的，而我编的程序用的是自己定的一种中弧线。有需要的人可以把我程序里中弧线部分替换成适合你们的情况即可。
程序源代码
% Wennerstrom叶型生成程序
% 基于Wennerstrom的书“Design of highly loaded axial-flow fans and compressors”
% 这种叶型适用于高亚音速或有一点跨音速的轴流压气机
% 用户可以自定进出口金属角，最大厚度及位置，前尾缘厚度等参数。
% 造型时吸力面在上压力面在下，最后根据需要转换
% 同时输出成适合MISES的叶型
clear;
Gridswitch=-1.0; % Gridswitch>0对应H型网格，Gridswitch<0对应O型网格
Rsign=1.0; % 1.0对应吸力面朝上，-1.0对应吸力面朝下
Npoint=101; % 从前缘到尾缘的叶型数据点数
Ncirc=10; % 前尾缘小圆点数的一半
Metal1=-45; % 进口金属角(应给正值，保证吸力面向上)
Metal2=0; % 出口金属角
Xthick=0.40; % 最大厚度位置（最大厚度位置/弦长）
Tmax=0.060; % 最大厚度（最大厚度/弦长）
T1=0.010; % 前缘厚度（厚度/弦长）
T2=0.010; % 尾缘厚度（厚度/弦长）
Nscale1=1.7; % 用于调整中弧线前部点疏密
Nscale2=1.7; % 用于调整中弧线后部点疏密
Pitch=0.75; % 栅距，只对之后MISES的计算有用
fid1=fopen('meanline&thick.dat','w'); % 中弧线和厚度
fid2=fopen('profile1.dat','w'); % 输出叶型
fid3=fopen('blade.t01','w'); % 输出叶型给MISES
% 中弧线的x坐标（前尾缘加密）
xc00=0:1/(Npoint-1):1;
for i=1:Npoint
if xc00(i)<=0.5
xc(i)=(2*xc00(i))^Nscale1/2;
else
xc(i)=1-(2*(1-xc00(i)))^Nscale2/2;
end
end
% 叶型中弧线的y坐标（弯角沿x线性变化）
Metal(1)=Metal1;
yc(1)=0;
for i=2:Npoint
Metal(i)=Metal(i-1) (Metal2-Metal1)*(xc(i)-xc(i-1))/(xc(end)-xc(1));
yc(i)=yc(i-1) (xc(i)-xc(i-1))*tand(Metal(i));
end
% 叶型厚度分布
if Xthick>=0.5
a=(T1/2-Tmax/2)/(2*Xthick^3);
b=0;
c=3*(Tmax/2-T1/2)/(2*Xthick);
d=T1/2;
e=3*(Tmax/2-T1/2)/(2*Xthick^2*(1-Xthick))-(Tmax/2-T2/2)/(1-Xthick)^3;
f=3*(T1/2-Tmax/2)/(2*Xthick^2);
g=0;
h=Tmax/2;
else
a=(Tmax/2-T1/2)/Xthick^3-3*(Tmax/2-T2/2)/(2*Xthick*(1-Xthick)^2);
b=3*(Tmax/2-T2/2)/(1-Xthick)^2-3*(Tmax/2-T1/2)/Xthick^2;
c=3*(Tmax/2-T1/2)/Xthick-3*Xthick*(Tmax/2-T2/2)/(2*(1-Xthick)^2);
d=T1/2;
e=(Tmax/2-T2/2)/(2*(1-Xthick)^3);
f=3*(T2/2-Tmax/2)/(2*(1-Xthick)^2);
g=0;
h=Tmax/2;
end
for i=1:Npoint
if xc(i)<=Xthick
thick(i)=a*xc(i)^3 b*xc(i)^2 c*xc(i) d;
% 前半段的三次曲线
else
thick(i)=e*(xc(i)-Xthick)^3 f*(xc(i)-Xthick)^2 g*(xc(i)-Xthick) h;
% 后半段的三次曲线
end
fprintf(fid1,'%7.4f %7.4f %7.4f %7.2f\n',xc(i),yc(i),thick(i),Metal(i));
% 输出中弧线和厚度
end
% 生成叶型
xs=xc-thick.*sind(Metal);
ys=yc thick.*cosd(Metal);
xp=xc thick.*sind(Metal);
yp=yc-thick.*cosd(Metal);
% 前缘小圆
St1=(Tmax/2-T1/2)/Xthick;
St2=(Tmax/2-T2/2)/(1-Xthick);
rr1=(sqrt((xc(1)-xp(1))^2 (yc(1)-yp(1))^2) sqrt((xc(1)-xs(1))^2 (yc(1)-ys(1))^2))/2;
delta1=atand(St1); % 此处的delta永远为正
r1=rr1/cosd(delta1);
dxc=rr1*tand(delta1);
xc0=xc(1) dxc;
yc0=(yc(2)-yc(1))/(xc(2)-xc(1))*(xc0-xc(1)) yc(1);
if xc0==xs(1) % Theta都转换成0-360°
Theta_start=90;
elseif xs(1)>xc0
Theta_start=atand((ys(1)-yc0)/(xs(1)-xc0));
else
Theta_start=atand((ys(1)-yc0)/(xs(1)-xc0)) 180;
end
if xc0==xp(1)
Theta_end=270;
elseif xp(1)>xc0
Theta_end=atand((yp(1)-yc0)/(xp(1)-xc0)) 360;
else
Theta_end=atand((yp(1)-yc0)/(xp(1)-xc0)) 180;
end
if Gridswitch<0 % 以前缘点为起点，适应O形网格
Theta_p=linspace((Theta_start Theta_end)/2,Theta_end,Ncirc 1);
Theta_s=linspace((Theta_start Theta_end)/2,Theta_start,Ncirc 1);
else % 以轴向最前点为起点，适应H形网格
Theta_p=linspace(180,Theta_end,Ncirc 1);
Theta_s=linspace(180,Theta_start,Ncirc 1);
end
xple=r1*cosd(Theta_p) xc0;
yple=r1*sind(Theta_p) yc0;
xsle=r1*cosd(Theta_s) xc0;
ysle=r1*sind(Theta_s) yc0;
% 尾缘小圆
rr2=(sqrt((xc(end)-xp(end))^2 (yc(end)-yp(end))^2) sqrt((xc(end)-xs(end))^2 (yc(end)-ys(end))^2))/2;
delta2=atand(St2); % 此处的delta永远为正
r2=rr2/cosd(delta2);
dxc=rr2*tand(delta2);
xcn=xc(end)-dxc;
ycn=(yc(end-1)-yc(end))/(xc(end-1)-xc(end))*(xcn-xc(end)) yc(end);
if xcn==xp(end) % Theta都转换成-180°-180°
Theta_start=-90;
elseif xp(end)>xcn
Theta_start=atand((yp(end)-ycn)/(xp(end)-xcn));
else
Theta_start=atand((yp(end)-ycn)/(xp(end)-xcn))-180;
end
if xcn==xs(end)
Theta_end=90;
elseif xs(end)>xcn
Theta_end=atand((ys(end)-ycn)/(xs(end)-xcn));
else
Theta_end=atand((ys(end)-ycn)/(xs(end)-xcn)) 180;
end
if Gridswitch<0
Theta_p=linspace(Theta_start,(Theta_start Theta_end)/2,Ncirc 1);
Theta_s=linspace(Theta_end,(Theta_start Theta_end)/2,Ncirc 1);
else
Theta_p=linspace(Theta_start,0,Ncirc 1); % 以轴向最后点为起点，适应Denton的H形网格
Theta_s=linspace(Theta_end,0,Ncirc 1);
end
xpte=r2*cosd(Theta_p) xc(end);
ypte=r2*sind(Theta_p) yc(end);
xste=r2*cosd(Theta_s) xc(end);
yste=r2*sind(Theta_s) yc(end);
% 生成带前尾缘小圆的叶型
xp1=[xple(1:Ncirc) xp xpte(2:Ncirc 1)];
yp1=[yple(1:Ncirc) yp ypte(2:Ncirc 1)]*Rsign; % 改变叶片旋转方向
xs1=[xsle(1:Ncirc) xs xste(2:Ncirc 1)];
ys1=[ysle(1:Ncirc) ys yste(2:Ncirc 1)]*Rsign;
for i=1:length(xp1)
fprintf(fid2,'%7.4f %7.4f %7.4f %7.4f\n',xp1(i),yp1(i),xs1(i),ys1(i));
% 输出叶型坐标
end
% 输出只有前缘带小圆，且顺序为TE->SS->LE->PS->TE的叶型，适合MISES
fprintf(fid3,'t01\n');
fprintf(fid3,'%7.3f %7.3f %7.3f %7.3f %7.3f\n',tand(Metal1),tand(Metal2),0.6,0.6,Pitch);
for i=1:length(xs)
fprintf(fid3,'%10.6f %10.6f\n',xs(end-i 1),ys(end-i 1));
end
for i=1:length(xsle)
fprintf(fid3,'%10.6f %10.6f\n',xsle(end-i 1),ysle(end-i 1));
end
for i=2:length(xple)
fprintf(fid3,'%10.6f %10.6f\n',xple(i),yple(i));
end
for i=1:length(xp)
fprintf(fid3,'%10.6f %10.6f\n',xp(i),yp(i));
end
% 画图
subplot(2,2,1);
plot(xc,yc,'b.-');
title('meanline');
axis equal;
grid on;
subplot(2,2,2);
plot(xc,Metal,'b.-');
title('metal angle');
grid on;
subplot(2,2,3);
plot(xc,thick,'b.-');
title('thickness');
grid on;
subplot(2,2,4);
plot(xp1,yp1,'b.-',xs1,ys1,'r.-');
title('profile');
axis equal;
grid on;
fclose(fid1);
fclose(fid2);
fclose(fid3);
运行结果
进口角：-45°
出口角：0°
最大厚度位置：40%
最大厚度：6%
前缘厚度：1%
尾缘厚度：1%
进口角：50°
出口角：30°
最大厚度位置：55%
最大厚度：6%
前缘厚度：1%
尾缘厚度：1%
进口角：50°
出口角：-50°
最大厚度位置：50%
最大厚度：15%
前缘厚度：3%
尾缘厚度：3%
进口角：30°
出口角：-70°
最大厚度位置：30%
最大厚度：20%
前缘厚度：8%
尾缘厚度：3%
从后面这两个叶型看，这个程序用来做涡轮设计似乎也是可以的，不过效果肯定不会太好。
