一种高性能压气机叶型的生成源代码
导读
这个Matlab程序是我以前编的，不进行流场优化，而是直接用数学曲线生成叶型，参考的是前一篇文章中Wennerstrom建议的叶型生成方法。对生成的叶型进行过大量数值模拟，发现只要参数选择适当，性能还是非常不错的。最适合的是高压压气机静子的设计，尤其是高负荷叶型以及带吸气的叶型，也可以生成看起来不错的转子叶型，不过没有具体计算验证过。大家如果需要可以自由使用和改进。我本人没用这个发表过文章，如果有人用这个程序生成的叶型计算或做实验，用于发表文章，我觉着也不需要引用。
叶型简介
叶型用前后两段曲线构建，在中部满足曲线的连续性。中弧线是两段四次曲线，厚度分布也是两段四次曲线。尾缘是圆弧，前缘用一段四次曲线和一段三次曲线替换了圆弧。
之所以可以不考虑流动而直接造型，是因为对已有的，经过气动优化的叶型的几何进行分析，总结了一些特征如下：
下面是叶型生成方法的简介：
至于生成的叶型性能如何，这里只给出一个数值模拟的例子。下面的图中，ori代表某原始叶型，MISES代表经过气动优化的叶型，Poly代表本程序生成的叶型。可以看到，用本程序生成的多项式叶型的损失系数与经过气动优化后的叶型相当。而落后角要明显更小，也就是有更大的做功潜力。落后角小的主要原因是本程序生成的叶型更倾向于前加载，后部的弯角小。
Matlab程序源代码
分为两个程序，第一个程序生成的叶型前尾缘都是圆弧。事实证明对于可控扩散叶型，前缘形状影响很大，所以采用圆弧前缘的性能不会太好。第二个程序对第一个程序生成的叶型进行修改，把前缘变成多项式的曲线。
该叶型控制参数较多，用户可以自行尝试修改代码前几行中的各个控制参数，用来生成适合相应任务的高负荷叶型。
1. 基本叶型生成程序：
% 可控扩散叶型或吸附式叶型
% 两段式四次曲线中弧线,两段式四次曲线厚度分布
% 前后缘为圆弧
% 造型时吸力面在上压力面在下，最后根据需要转换
% 同时输出成适合MISES的叶型
% 2011-03-02
clear;
fid1=fopen('meanline&thick_cda.dat','w'); % 输出中弧线和厚度
fid2=fopen('profile_cda.dat','w'); % 输出叶型
fid3=fopen('blade.p5','w'); % 输出给MISES的叶型
Gridswitch=-1.0; % Gridswitch>0对应H型网格，Gridswitch<0对应O型网格
Rsign=1.0; % 吸力面向上为1.0，吸力面向下为-1.0
Npoint=101; % 叶型数据点数
Ncirc=4; % 前尾缘小圆点数/2
DeltaBeta=5; % 进口附加的金属角，让前缘叶型“低头”
Metal1=48 DeltaBeta; % 进口金属角(应给正值，保证吸力面向上)
Metal2=12; % 出口金属角
Xmean=0.08; % 中弧线拐点位置
K_m=0.18; % Xmean处已经发生的转角占总转角的比例
Ks1=5.00; % 前缘处s的斜率与前缘到尾缘s总斜率的倍数
Ks2=0.10; % 尾缘处s的斜率与前缘到尾缘s总斜率的倍数
Xthick=0.25; % 最大厚度位置
Tmax=0.08; % 最大厚度
T1=0.015; % 前缘厚度
T2=0.015; % 尾缘厚度
SS1=3.2; % 进口段厚度变化系数(0-4.0)，越大表示厚度增加越快
SS2=0.8; % 出口段厚度变化系数(0-2.0)，越小表示厚度减小越快
Nscale1=1.7; % 调整中弧线前部疏密
Nscale2=1.7; % 调整中弧线后部疏密
Pitch=0.8501; % For MISES only（栅距）
% 沿中弧线座标
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
% 中弧线参数
Sm1=tand(Metal1);
Sm2=tand(Metal2);
Sm0=tand(Metal1 (Metal2-Metal1)*K_m);
SSm1=Ks1*(Sm2-Sm1); % 进口角度正切值变化率
SSm2=Ks2*(Sm2-Sm1); % 出口角度正切值变化率
a1=1/2*(-SSm1*Xmean-Xmean^3*SSm2 Xmean^3*SSm1 Xmean^2*SSm2 3*Xmean^2*Sm1-3*Xmean^2*Sm2 Sm0-Sm1 ...
2*Xmean*Sm0-2*Xmean*Sm1)/(-1 Xmean)/Xmean^3;
a2=-1/2/Xmean^2*(-Xmean^3*SSm2 Xmean^3*SSm1 Xmean^2*SSm2 3*Xmean^2*Sm1-3*Xmean^2*Sm2 2*SSm1*Xmean^2-...
3*SSm1*Xmean 3*Sm0-3*Sm1)/(-1 Xmean);
a3=SSm1;
a4=Sm1;
a5=-1/2*(SSm1*Xmean-6*Xmean*Sm1-Xmean^3*SSm2 Xmean^3*SSm1 3*Xmean^2*Sm1 3*Xmean^2*SSm2 2*Xmean*Sm0-...
2*SSm1*Xmean^2-3*Xmean^2*Sm2 4*Xmean*Sm2-2*SSm2*Xmean-3*Sm0 3*Sm1)/Xmean/(-1 3*Xmean-3*Xmean^2 Xmean^3);
a6=-1/2*(-Xmean^3*SSm2 Xmean^3*SSm1 5*Xmean^2*SSm2-2*SSm1*Xmean^2 3*Xmean^2*Sm1-3*Xmean^2*Sm2-...
4*SSm2*Xmean SSm1*Xmean-6*Xmean*Sm1 6*Xmean*Sm2 3*Sm1-3*Sm0)/Xmean/(-1 Xmean)^2;
a7=-SSm2;
a8=Sm2;
% 厚度参数
St1=SS1*(Tmax/2-T1/2)/Xthick;
St2=SS2*(Tmax/2-T2/2)/(1-Xthick);
term1=(Tmax-T1)*(1-Xthick)^2;
term2=(Tmax-T2)*Xthick^2;
term3=St1*Xthick*(1-Xthick)^2;
term4=St2*(1-Xthick)*Xthick^2;
b1=(3/2*term1-3*term2-term3 3*term4)/((1-Xthick)^2*Xthick^4);
b2=(-4*term1 6*term2 3*term3-6*term4)/((1-Xthick)^2*Xthick^3);
b3=(3*term1-3*term2-3*term3 3*term4)/((1-Xthick)^2*Xthick^2);
b4=St1;
b5=T1/2;
b6=(-3/2*Tmax 3/2*T2 2*St2*(1-Xthick))/(1-Xthick)^4;
b7=(2*Tmax-2*T2-3*St2*(1-Xthick))/(1-Xthick)^3;
b8=0;
b9=St2;
b10=T2/2;
for i=1:Npoint
if xc(i)<=Xmean
sc(i)=a1*xc(i)^3 a2*xc(i)^2 a3*xc(i) a4;
yc(i)=1/4*a1*xc(i)^4 1/3*a2*xc(i)^3 1/2*a3*xc(i)^2 a4*xc(i);
else
sc(i)=a5*(1-xc(i))^3 a6*(1-xc(i))^2 a7*(1-xc(i)) a8;
yc(i)=-1/4*a5*(1-xc(i))^4-1/3*a6*(1-xc(i))^3-1/2*a7*(1-xc(i))^2-a8*(1-xc(i)) (1/4*a1*Xmean^4 1/3*a2*Xmean^3 ...
1/2*a3*Xmean^2 a4*Xmean)-(-1/4*a5*(1-Xmean)^4-1/3*a6*(1-Xmean)^3-1/2*a7*(1-Xmean)^2-a8*(1-Xmean));
end
if xc(i)<=Xthick
thick(i)=b1*xc(i)^4 b2*xc(i)^3 b3*xc(i)^2 b4*xc(i) b5;
else
thick(i)=b6*(1-xc(i))^4 b7*(1-xc(i))^3 b8*(1-xc(i))^2 b9*(1-xc(i)) b10;
end
fprintf(fid1,'%7.4f %7.4f %7.4f %7.4f\n',xc(i),yc(i),sc(i),thick(i));
end
% 生成叶型
xs=xc-thick.*sin(atan(sc));
ys=yc thick.*cos(atan(sc));
xp=xc thick.*sin(atan(sc));
yp=yc-thick.*cos(atan(sc));
% 前缘小圆
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
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
Theta_p=linspace(180,Theta_end,Ncirc 1); % 以轴向最前点为起点，适应H形网格
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); % 以轴向最后点为起点，适应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,'%8.4f %8.4f %8.4f %8.4f %8.4f\n',tand(Metal1-8),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) % for blunt LE, use 2:end
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
% 画图
plot(xp1,yp1,'b.-',xs1,ys1,'b.-'); hold on;
title('profile');
axis equal;
grid on;
fclose(fid1);
fclose(fid2);
fclose(fid3);
2. 前缘修改程序：
% 前缘生成（鹰嘴形）
% 压力面是三次曲线，吸力面是四次曲线
% 对吸力面的四次曲线，去掉后部的二阶导连续条件，在前缘附近再给一个一阶导条件
% 2011-04-08
clear;
fout=fopen('profile_cda_le.dat','w'); % 生成的叶型
fout1=fopen('blade.p56','w'); % 生成的叶型（MISES）
bladexy=load('profile_cda.dat'); % 原始叶型座标
xp=bladexy(:,1);
yp=bladexy(:,2);
xs=bladexy(:,3);
ys=bladexy(:,4);
Metal1=48 5;
Metal2=12;
Pitch=0.8501;
rc=0.0075; % 前缘半径
alpha=0.25; % 前缘点位移系数，一般不可大于0.4，否则吸力面曲线会有拐点
beta=2.5; % 压力面交点位置系数
gama=5.0; % 吸力面交点位置系数
delta=0.1; % 交点一二阶导计算用系数
beta2=0.5; % 在前缘点处控制吸力面曲线用的系数
ncirc=20; % 新叶型前缘点数/2
NN=4; % 原叶型尾缘点数/2
Xc=rc*alpha; % 前缘点向压力面移动的距离
Yc=0; % 前缘点座标，基本保持弦长不变。
Yp=rc*beta; % 压力面曲线与压力面前缘曲线的交点
Ys=rc*gama; % 吸力面曲线与吸力面前缘曲线的交点
Xs1=0*rc; % 在此点处给一个吸力面前缘曲线的斜率来控制其走向
Midp=floor(length(yp)/2); % 压力面前半段终点（要保证转成竖起来后没有重复的X座标）
Mids=floor(length(ys)/5); % 吸力面前半段终点（要保证转成竖起来后没有重复的X座标）
% 旋转叶型座标，使进口部分叶型“竖起来”
xp0=xp-xp(1);
yp0=yp-yp(1);
xs0=xs-xp(1);
ys0=ys-yp(1);
xxp=xp0*cosd(-90 Metal1) yp0*sind(-90 Metal1);
yyp=-xp0*sind(-90 Metal1) yp0*cosd(-90 Metal1);
xxs=xs0*cosd(-90 Metal1) ys0*sind(-90 Metal1);
yys=-xs0*sind(-90 Metal1) ys0*cosd(-90 Metal1);
% 求; Xp, Yp', Yp"; Xs, Ys', Ys"
% Xp,Yp,Xs,Ys分别为前缘曲线与叶型曲线的交点
interptype='pchip';
Xp=interp1(yyp(1:Midp),xxp(1:Midp),Yp,interptype); % 通过事先确定的Yp和Ys来插值得到Xp和Xs
Xs=interp1(yys(1:Mids),xxs(1:Mids),Ys,interptype);
Yp11=Yp-rc*delta; % 通过事先确定的DeltaY来确定Yp11
Xp11=interp1(yyp(1:Midp),xxp(1:Midp),Yp11,interptype); % 通过Yp11确定Xp11
DeltaX=Xp-Xp11; % 通过Xp11和Yp11确定DeltaX
Xs11=Xs-DeltaX;
Ys11=interp1(xxs(1:Mids),yys(1:Mids),Xs11,interptype); % 通过Xs11确定Ys11
Xp12=Xp DeltaX; % 通过DeltaX确定Xp12
Xs12=Xs DeltaX;
Yp12=interp1(xxp(1:Midp),yyp(1:Midp),Xp12,interptype);
Ys12=interp1(xxs(1:Mids),yys(1:Mids),Xs12,interptype);
Yp_prim=(Yp12-Yp11)/(2*DeltaX);
Ys_prim=(Ys12-Ys11)/(2*DeltaX);
Yp_prim2=(Yp12-2*Yp Yp11)/DeltaX^2;
Ys_prim2=(Ys12-2*Ys Ys11)/DeltaX^2;
AA_p=[Xc^3,Xc^2,Xc,1; 3*Xc^2,2*Xc,1,0; Xp^3,Xp^2,Xp,1; 3*Xp^2,2*Xp,1,0];
bb_p=[0; 0; Yp; Yp_prim]; % BC为：前缘点Y座标和一阶导等于0，交点处座标和一阶导连续
cc_p=AA_p\bb_p; % 解出系数
xple=linspace(Xs,Xp,ncirc*2);
yple=cc_p(1)*xple.^3 cc_p(2)*xple.^2 cc_p(3)*xple cc_p(4); % 压力面前缘曲线
% 利用压力面前缘原曲线左边延伸部分来控制吸力面前缘曲线
Ys1_left=interp1(xple,yple,Xs1-DeltaX,interptype); % 压力面前缘原曲线的左点
Ys1_right=interp1(xple,yple,Xs1 DeltaX,interptype); % 压力面前缘原曲线的右点
Ys1_prim=(Ys1_right-Ys1_left)/(2*DeltaX)*beta2; % Xs1处吸力面前缘曲线的斜率
AA_s=[Xc^4,Xc^3,Xc^2,Xc,1; 4*Xc^3,3*Xc^2,2*Xc,1,0; Xs^4,Xs^3,Xs^2,Xs,1; 4*Xs^3,3*Xs^2,2*Xs,1,0; 4*Xs1^3,3*Xs1^2,2*Xs1,1,0];
bb_s=[0; 0; Ys; Ys_prim; Ys1_prim];
cc_s=AA_s\bb_s;
xsle=linspace(Xs,Xc,ncirc);
ysle=cc_s(1)*xsle.^4 cc_s(2)*xsle.^3 cc_s(3)*xsle.^2 cc_s(4)*xsle cc_s(5);
xple=linspace(Xc,Xp,ncirc);
yple=cc_p(1)*xple.^3 cc_p(2)*xple.^2 cc_p(3)*xple cc_p(4); % 重新计算压力面前缘曲线，从前缘点开始算起
% 生成全部叶型座标
for i=1:ncirc
xxsle(i)=xsle(ncirc-i 1);
yysle(i)=ysle(ncirc-i 1);
end
ii=0;
for i=1:length(xxs) % 去掉前部小圆的数据
if yys(i)>=Ys
ii=ii 1;
xxs1(ii)=xxs(i);
yys1(ii)=yys(i);
end
end
ii=0;
for i=1:length(xxp) % 去掉前部小圆的数据
if yyp(i)>=Yp
ii=ii 1;
xxp1(ii)=xxp(i);
yyp1(ii)=yyp(i);
end
end
xxxs=[xxsle xxs1]; % 生成全叶型座标
yyys=[yysle yys1];
xxxp=[xple xxp1];
yyyp=[yple yyp1];
% 旋转回去
bxp=xxxp*cosd(90-Metal1) yyyp*sind(90-Metal1);
byp=-xxxp*sind(90-Metal1) yyyp*cosd(90-Metal1);
bxs=xxxs*cosd(90-Metal1) yyys*sind(90-Metal1);
bys=-xxxs*sind(90-Metal1) yyys*cosd(90-Metal1);
DX=bxp(end)-xp(end); % 与原叶型对齐
DY=byp(end)-yp(end);
bxp=bxp-DX;
byp=byp-DY;
bxs=bxs-DX;
bys=bys-DY;
plot(xs,ys,'b -',xp,yp,'b -',bxs,bys,'cx-',bxp,byp,'cx-'); % 画图与原叶型对比
axis equal;
xlim([-0.015 0.05]);
ylim([-0.015 0.035]);
grid on;
for i=1:length(bxp)
fprintf(fout,'%7.4f %7.4f\n',bxp(i),byp(i));
end
for i=1:length(bxs)
fprintf(fout,'%7.4f %7.4f\n',bxs(i),bys(i));
end
% 输出到MISES
fprintf(fout1,'t01_le_1\n');
fprintf(fout1,'%8.4f %8.4f %8.4f %8.4f %8.4f\n',tand(Metal1-5),tand(Metal2),0.6,0.6,Pitch);
for i=1 NN:length(bxs)
fprintf(fout1,'%10.6f %10.6f\n',bxs(end-i 1),bys(end-i 1));
end
for i=2:length(bxp)-NN
fprintf(fout1,'%10.6f %10.6f\n',bxp(i),byp(i));
end
fclose(fout);
fclose(fout1);
运行结果
进口角：48°
出口角：12°
最大厚度位置：25%
最大厚度：8%
前缘厚度：1.5%
尾缘厚度：1.5%
进口角：40°
出口角：0°
最大厚度位置：30%
最大厚度：8%
前缘厚度：1.5%
尾缘厚度：1.5%
进口角：55°
出口角：45°
最大厚度位置：50%
最大厚度：6%
前缘厚度：0.6%
尾缘厚度：1%
下面是第二个程序进行前缘修改后输出的叶型前缘图。有些文献上把这种前缘称为“鹰嘴形”前缘。这种前缘有两个好处，一是有利于消除吸力面的压力峰，这对吸力面的边界层尽量保持较长的层流很有好处。即使不考虑转捩，没有经过吸力峰的边界层也更饱满，而不容易产生分离泡。另一方面，这样的前缘也可以使叶型的可用攻角范围更大，对非设计工况适应更好。当然，采用多项式前缘对加工和检测的要求都高一些。
