clear
clc
%%样本数据获取
filename=ls('D:\RAE2822\2822base\RAE 2822.dat');%翼型数据文件需和代码放同一个路径
delimiterIn=' ';
headerlinesIn=1;%不要第一行的翼型名称 从第二行开始读取翼型坐标
data_cell=importdata(filename,delimiterIn,headerlinesIn);
bound = find(data_cell.data(:,1)==0);
numend = size(data_cell.data,1);
for i=1:bound%从第一个坐标点到最前缘点的行数 即翼型上表面曲线坐标 取决于翼型dat文件
a(i)=num2cell(data_cell.data(i,1));%将翼型横坐标x赋给矩阵a, 不知道为啥要转换，不然运行不了，下面得再转回来
b(i)=data_cell.data(i,2);%将翼型纵坐标y赋给矩阵b
end
for j=bound:numend%从最前缘点到最后一行的行数 即翼型下表面曲线坐标 取决于翼型dat文件
aa(j)=data_cell.data(j,1);%将翼型横坐标x赋给矩阵aa
bb(j)=data_cell.data(j,2);%将翼型纵坐标y赋给矩阵bb
end
for ii=1:numend
aaa(ii)=data_cell.data(ii,1);%将翼型所有的横坐标x赋给矩阵aaa
bbb(ii)=data_cell.data(ii,2);
end
a_=cell2mat(a);
x_up = a_';
x_low = aa(bound:numend)';
y_up = b';
y_low = bb(bound:numend)';
xall = [x_up;x_low];
y_old = [y_up;y_low];
%%
%上下翼面BPO阶数设定
m = 6; %上翼面
n = 6; %下翼面
%翼型后缘厚度
dzu = 0.00002;
dzl = 0;
%前缘半径
r = 0.008496;
%上下翼面参数拟合结果
wu = CSTfit_weight(m, x_up, y_up,dzu);
wl = CSTfit_weight(n, x_low, y_low,-dzl);
%Class类函数确定
N1 = 0.5;
N2 = 1;
%生成新的上下翼面数据
ynewup = ClassShape(wu,x_up,N1,N2,r,dz);
ynewlow = ClassShape(wl,x_low,N1,N2,r,-dz);
yall = [ynewup;ynewlow];
%%
%拟合结果对比
%show
figure(1)
one = scatter(xall,y_old,'r+');
hold on;
two = plot(x_up,ynewup,'g-');
hold on;
three = plot(x_low,ynewlow,'b-');
xlabel('x')
ylabel('y')
legend('原始翼型','CST-BPO6上表面','CST-BPO6下表面')
%总体标准差
sum=0;
for i=1:size(xall,1)
sum=sum+(yall(i)-y_old(i))^2;
end
s = (sum/(size(xall,1)-1))^0.5;
%error 残差
for i=1:size(xall,1)
e(i)=yall(i)-y_old(i);
end
x_show=[x_up;-x_low];
figure(2)
plot(x_show,e,'r-')
xlabel('x')
ylabel('残差')
string = {'伯恩斯坦阶数BPO = 6';'总体标准差s =';s};
title(string)
评论0