根据叶型坐标计算叶型参数的代码
导读我们工作时经常有这样的需求，已经有了叶型数据，想知道它的中弧线形状、厚度分布、进出口金属角等参数。我曾经编写过这样一个程序，现在分享出来，希望对大家的工作有用。
程序简介
二维叶型坐标就是叶片的叶盆和叶背的XY坐标。在不同的软件中表达不同，比如MISES中用的是从尾缘开始绕过前缘再到尾缘的顺序来列出坐标，也有些软件中是把叶盆和叶背单列的方式。这里给出的代码前一部分是先把MISES的坐标变成叶盆和叶背单列的方式，然后再提取中弧线和计算厚度分布。
更多的说明请参见程序中的说明语句。
MISES是叶轮机械的二维分析和设计软件包，不属于本文关心内容，这里就不多介绍了，实际上读者也不需要关心。本文给出的代码稍加改动就可以适用于任何形式的叶型坐标。
不废话，下面给出程序源代码。
程序源代码
% 从叶型数据坐标点计算叶型的几何参数
% 输入叶型数据文件格式为MISES生成的叶型blade.in
% 输入坐标点顺序为逆时针，尾缘→叶背→前缘→叶盆→尾缘，不含尾缘小圆
% 输出的叶型参数包括：中弧线坐标，厚度分布，金属角分布，曲率分布
clear;
npoint=200; % 前缘到尾缘共200点
nfit1=4; % 光顺金属角分布时的拟合多项式次数
fin=fopen('blade.in','r'); % 输入叶型数据
bladexy=textscan(fin,'%f %f','HeaderLines',2); % 前两行不读
bladex=bladexy{1}; % 叶型x坐标
bladey=bladexy{2}; % 叶型y坐标
fid=fopen('blade1.out','w'); % 输出叶型坐标
fid1=fopen('bpara1.dat','w'); % 输出叶型参数
% 将坐标表达成从前缘到后缘的排列方式
nall=length(bladex);
[minx,nmin]=min(bladex); % 暂以x最小点为叶盆和叶背分界
np=nmin; % 叶盆点数
ns=nall-nmin; % 叶背点数
for i=1:np
bpx(i)=bladex(np-i 1); % 叶盆坐标
bpy(i)=bladey(np-i 1);
end
for i=1:ns
bsx(i)=bladex(i np); % 叶背坐标
bsy(i)=bladey(i np);
end
bpx(1)=bsx(1);
bpy(1)=bsy(1);
% 将叶盆和叶背插值成相同多的点数并输出
bpxi=linspace(min(bpx),max(bpx),npoint);
bsxi=linspace(min(bsx),max(bsx),npoint);
bpyi=interp1(bpx,bpy,bpxi,'pchip');
bsyi=interp1(bsx,bsy,bsxi,'pchip');
for i=1:npoint
fprintf(fid,'%7.4f %7.4f %7.4f %7.4f\n',bpxi(i),bpyi(i),bsxi(i),bsyi(i));
end
% 对于叶片各段进行局部旋转，得到当地中弧线和厚度
xp=bpxi;
yp=bpyi;
xs=bsxi;
ys=bsyi;
xc0=(xp xs)/2; % 这只是近似的中弧线
yc0=(yp ys)/2;
for i=1:npoint
if i==1
theta0(i)=atand((yc0(i 1)-yc0(i))/(xc0(i 1)-xc0(i)));
elseif i==npoint
theta0(i)=atand((yc0(i)-yc0(i-1))/(xc0(i)-xc0(i-1)));
else
theta0(i)=atand((yc0(i 1)-yc0(i-1))/(xc0(i 1)-xc0(i-1)));
end
for j=1:npoint
xp1(j)=xp(j)*cosd(theta0(i)) yp(j)*sind(theta0(i)); % 中弧线旋转到水平
yp1(j)=-xp(j)*sind(theta0(i)) yp(j)*cosd(theta0(i));
xs1(j)=xs(j)*cosd(theta0(i)) ys(j)*sind(theta0(i));
ys1(j)=-xs(j)*sind(theta0(i)) ys(j)*cosd(theta0(i));
end
xc1=(xp1 xs1)/2;
yp1i=interp1(xp1,yp1,xc1,'pchip'); % 插值到同样x座标
ys1i=interp1(xs1,ys1,xc1,'pchip');
thick0=ys1i-yp1i; % 厚度
thick(i)=thick0(i);
yc1=(ys1i yp1i)/2;
xc(i)=xc1(i)*cosd(-theta0(i)) yc1(i)*sind(-theta0(i)); % 旋转回去得到中弧线
yc(i)=-xc1(i)*sind(-theta0(i)) yc1(i)*cosd(-theta0(i));
end
yc(1)=(yc(3)-yc(2))/(xc(3)-xc(2))*(xc(1)-xc(2)) yc(2); % 中弧线起点
yc(npoint)=(yc(npoint-2)-yc(npoint-1))/(xc(npoint-2)-xc(npoint-1))*(xc(npoint)-xc(npoint-1)) yc(npoint-1); % 中弧线终点
thick=abs(thick); % 把厚度变为正值
% 计算叶片金属角theta沿轴向（x向）的分布
for i=1:npoint
if i==1
theta(i)=atand((yc(i 1)-yc(i))/(xc(i 1)-xc(i)));
elseif i==npoint
theta(i)=atand((yc(i)-yc(i-1))/(xc(i)-xc(i-1)));
else
theta(i)=atand((yc(i 1)-yc(i-1))/(xc(i 1)-xc(i-1)));
end
end
p1=polyfit(xc,theta,nfit1); % 前缘附近金属角计算误差大，拟合光顺
thetai=polyval(p1,xc);
% 计算中弧线曲率沿轴向（x向）的分布
for i=1:npoint
if i==1
delta(i)=thetai(i 1)-thetai(i);
chord=sqrt((yc(i 1)-yc(i))^2 (xc(i 1)-xc(i))^2);
curve(i)=-sind(0.5*delta(i))/(0.5*chord);
elseif i==npoint
delta(i)=thetai(i)-thetai(i-1);
chord=sqrt((yc(i)-yc(i-1))^2 (xc(i)-xc(i-1))^2);
curve(i)=-sind(0.5*delta(i))/(0.5*chord);
else
delta(i)=thetai(i 1)-thetai(i-1); % 弯角
chord=sqrt((yc(i 1)-yc(i-1))^2 (xc(i 1)-xc(i-1))^2)/2; % 当地弦长
curve(i)=-sind(0.5*delta(i))/(0.5*chord); % 曲率
end
end
% 输出叶型参数
for i=1:npoint
fprintf(fid1,'%7.4f %7.4f %7.4f %7.4f %7.4f\n',xc(i),yc(i),thick(i),thetai(i),curve(i));
end
% 画图
figure(1)
subplot(2,2,1);
plot(xc0,yc0,'b.',xc,yc,'k-',xp,yp,'k-',xs,ys,'k-');
title('meanline');
axis equal;
grid on;
subplot(2,2,2);
plot(xc,thick,'b.-');
title('thick vs x');
grid on;
subplot(2,2,3);
plot(xc,theta,'b.-',xc,thetai,'r.-');
title('theta vs x');
grid on;
subplot(2,2,4);
plot(xc,curve,'b.-');
title('curve vs x');
grid on;
fclose(fin);
fclose(fid);
fclose(fid1);
输入数据文件balde.in
t01
1.1106 0.2126 0.6000 0.6000 0.8501
0.997921 0.522033
0.997266 0.521930
0.995794 0.521699
0.993683 0.521367
0.991010 0.520946
0.987821 0.520443
0.984151 0.519863
0.980025 0.519209
0.975464 0.518485
0.970485 0.517692
0.965101 0.516832
0.959326 0.515906
0.953171 0.514916
0.946644 0.513861
0.939755 0.512744
0.932512 0.511564
0.924920 0.510320
0.916987 0.509014
0.908719 0.507644
0.900120 0.506211
0.891195 0.504713
0.881949 0.503149
0.872386 0.501519
0.862510 0.499820
0.852325 0.498052
0.841832 0.496211
0.831036 0.494295
0.819939 0.492302
0.808543 0.490228
0.796851 0.488069
0.784864 0.485821
0.772586 0.483479
0.760016 0.481038
0.747157 0.478492
0.734009 0.475834
0.720575 0.473057
0.706855 0.470151
0.692849 0.467109
0.678559 0.463919
0.663985 0.460571
0.649128 0.457052
0.633987 0.453349
0.618565 0.449448
0.602861 0.445331
0.586875 0.440982
0.570608 0.436382
0.554062 0.431510
0.537237 0.426343
0.520135 0.420859
0.502757 0.415032
0.485107 0.408832
0.467435 0.402326
0.449996 0.395594
0.432797 0.388634
0.415841 0.381444
0.399137 0.374023
0.382689 0.366372
0.366505 0.358492
0.350591 0.350388
0.334954 0.342064
0.319600 0.333524
0.304535 0.324777
0.289767 0.315830
0.275302 0.306693
0.261144 0.297377
0.247299 0.287895
0.233770 0.278260
0.220565 0.268486
0.207687 0.258587
0.195144 0.248575
0.182940 0.238465
0.171083 0.228271
0.159580 0.218009
0.148438 0.207700
0.137666 0.197364
0.127269 0.187025
0.117254 0.176707
0.107623 0.166436
0.098380 0.156236
0.089527 0.146135
0.081065 0.136160
0.072995 0.126336
0.065315 0.116691
0.058024 0.107252
0.051120 0.098045
0.044602 0.089098
0.038467 0.080438
0.032710 0.072091
0.027331 0.064086
0.022325 0.056448
0.017691 0.049207
0.013426 0.042391
0.009530 0.036030
0.006003 0.030157
0.002848 0.024807
0.000070 0.020017
-0.002321 0.015833
-0.004309 0.012309
-0.005869 0.009516
-0.006952 0.007562
-0.007431 0.006691
-0.007361 0.006647
-0.008024 0.005406
-0.008507 0.004086
-0.008801 0.002710
-0.008901 0.001307
-0.008805 -0.000096
-0.008513 -0.001472
-0.008033 -0.002794
-0.007373 -0.004036
-0.006547 -0.005174
-0.005570 -0.006186
-0.004462 -0.007052
-0.003243 -0.007755
-0.001939 -0.008281
-0.000574 -0.008621
0.000825 -0.008766
0.002231 -0.008716
0.003615 -0.008470
0.004952 -0.008033
0.006216 -0.007415
0.007380 -0.006626
0.007431 -0.006691
0.008245 -0.006125
0.010072 -0.004850
0.012682 -0.003018
0.015975 -0.000692
0.019882 0.002092
0.024355 0.005307
0.029349 0.008933
0.034832 0.012955
0.040770 0.017357
0.047136 0.022127
0.053904 0.027252
0.061050 0.032719
0.068553 0.038514
0.076393 0.044623
0.084551 0.051032
0.093010 0.057725
0.101753 0.064684
0.110767 0.071892
0.120039 0.079330
0.129556 0.086978
0.139308 0.094816
0.149288 0.102825
0.159486 0.110982
0.169897 0.119267
0.180517 0.127659
0.191342 0.136140
0.202370 0.144691
0.213601 0.153294
0.225037 0.161936
0.236681 0.170603
0.248532 0.179284
0.260593 0.187966
0.272864 0.196638
0.285346 0.205291
0.298041 0.213915
0.310947 0.222503
0.324067 0.231050
0.337399 0.239551
0.350947 0.248001
0.364710 0.256396
0.378694 0.264731
0.392898 0.273001
0.407328 0.281203
0.421986 0.289334
0.436875 0.297389
0.451998 0.305366
0.467359 0.313262
0.482959 0.321075
0.498803 0.328802
0.514893 0.336441
0.531004 0.343888
0.546909 0.351048
0.562608 0.357934
0.578099 0.364559
0.593380 0.370936
0.608450 0.377076
0.623306 0.382990
0.637946 0.388687
0.652365 0.394178
0.666562 0.399471
0.680533 0.404575
0.694274 0.409499
0.707782 0.414249
0.721054 0.418832
0.734086 0.423256
0.746874 0.427527
0.759414 0.431650
0.771704 0.435630
0.783739 0.439474
0.795515 0.443184
0.807028 0.446766
0.818276 0.450224
0.829253 0.453561
0.839957 0.456781
0.850382 0.459885
0.860525 0.462878
0.870381 0.465761
0.879946 0.468536
0.889215 0.471205
0.898184 0.473770
0.906847 0.476231
0.915199 0.478590
0.923236 0.480848
0.930949 0.483003
0.938335 0.485058
0.945385 0.487011
0.952092 0.488862
0.958448 0.490610
0.964445 0.492253
0.970073 0.493792
0.975320 0.495223
0.980175 0.496543
0.984622 0.497751
0.988646 0.498842
0.992226 0.499810
0.995336 0.500651
0.997944 0.501355
1.000004 0.501911
1.001440 0.502298
1.002079 0.502470
运行结果
图：
就这样。
声明：原创文章，欢迎留言与我讨论，如需转载留言
请问有没有什么书或者理论涉及这部分的知识？
谢谢老师
