拉压弯曲线图及柱状图2022/12
更新于2023.03.13
{% folding child:codeblock open:flase color:yellow 点击展开代码详情 %}
%% 基础信息
% 反复修改的代码简直屎山！
clc;clear
close all
% color1=[77,133,189]/255;
% color2=[247,144,61]/255;
% color3=[89,169,90]/255;
% 马卡龙配色
color1=[147 213 220]/255;
color2=[242 202 201]/255;
color3=[198 223 200]/255;
% color1=[255 230 215]/255;
% color2=[225 240 219]/255;
% color3=[222 234 248]/255;
% color1=[240 181 125]/255;
% color2=[211 225 174]/255;
% color3=[113 171 182]/255;
%
% color1=[ 76 139 192]/255;
% color2=[ 232 49 51]/255;
% color3=[ 94 183 91]/255;
% 大论文使用的配色
% color1=[ 110 158 206]/255;
% color2=[ 118 186 128]/255;
% color3=[ 230 146 143]/255;
FontSize=13;
color3s={color1;color2;color3};
colors={color1;color1;color1;color2;color2;color2;color3;color3;color3};
colors9=[
199 92 100
240 181 125
211 225 174
113 171 182
75 90 161
76 139 192
232 49 51
163 96 173
94 183 91]/255;
ki=[ 1 0 0
0.75 0.25 0
0.5 0.5 0
0.25 0.75 0
0 1 0
0 0.75 0.25
0 0.5 0.5
0 0.25 0.75
0 0 1];
bar_colors=ki*[color1;color3;color2];
legends=["v=240mm/s T=210°C h=0.16mm"
"v=250mm/s T=215°C h=0.20mm"
"v=260mm/s T=220°C h=0.24mm"
"v=240mm/s T=215°C h=0.24mm"
"v=250mm/s T=220°C h=0.16mm"
"v=260mm/s T=210°C h=0.20mm"
"v=240mm/s T=220°C h=0.20mm"
"v=250mm/s T=210°C h=0.24mm"
"v=260mm/s T=215°C h=0.16mm"]';
legends9=[legends;repmat("",5,9)];
%%
hexc=[
'32aeec'
'E3863c'
'E9262e'];
r=hex2dec(hexc(:,[1,2]));
g=hex2dec(hexc(:,[3 4]));
b=hex2dec(hexc(:,[5 6]));
rgb=[r,g,b];
%% 拉伸
Tensile_x=cell(9,6); % read strain to Tensile_x
Tensile_y=cell(9,6); % read stress to Tensile_y
% 读取原始数据并储存
% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);
% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";
% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];
% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
T_L0=50; % 标距L0
T_S0=4*13; % 横截面积S0
for Gi=1:9
for ti=1:6
filesname=sprintf("/Users/mawenqian/Documents/硕士/VSL/科研/小论文1/基本试件打印参数正交试验/G%d/N%d.txt",Gi,ti);
readdata = readtable(filesname, opts);
% 转换为输出类型
readdata = table2array(readdata);
Tensile_x{Gi,ti}=readdata(:,5)/T_L0*100; % 应变 6sh/L^2
Tensile_y{Gi,ti}=readdata(:,1)/T_S0; % 应力 3Fl/2bh^2
clear readdata filesname
end
end
clear opts
save Tensile_raw_data Tensile_x Tensile_y
%%
% 通过下方原始图像人工判断获得，亦可绘制后展示。
Tensile_start_line=[
133 137 137 144 145 139
151 151 147 157 160 155
176 178 177 179 191 174
168 205 143 194 162 192
157 162 163 154 222 170
196 242 175 206 246 174
143 141 145 130 138 146
133 129 135 126 130 127
145 150 163 145 138 142];
Tensile_end_line=[
2022 1764 3568 2073 2596 3338
2222 1912 2363 3740 1927 1832
2613 1821 2390 2353 2041 2040
4000 3000 2449 2667 2557 1937
2000 1943 2000 1776 2000 2000
1719 1719 1719 1719 1719 1719
2312 1975 2537 2161 1730 1979
1407 1238 1345 1260 1593 1684
1656 2419 2000 2000 2000 2000];
Tensile_maxy_line=[
740 709 744 730 755 716
1066 1103 1064 1167 1058 1157
834 709 735 779 779 723
1377 1441 1622 1512 1698 1399
1043 1074 1069 1140 1072 1036
1214 1142 1164 1214 1142 1164
1178 1099 1213 1252 1191 1182
1097 1091 1081 1035 1100 1063
853 903 920 896 891 922];
% ttt=1;
% for Gi=[ 7 8]
% for ti=1:6
% subplot(4,6,ttt)
%
%
% x=Tensile_x{Gi,ti};
% y=Tensile_y{Gi,ti};
% [~,Tensile_maxy_line(Gi,ti)]=max(y);
% plot(y,MarkerIndices=Tensile_end_line(Gi,ti),MarkerSize=6,Marker="*",MarkerFaceColor='r')
% title(sprintf("%d-%d",Gi,ti))
% ttt=ttt+1;
% end
%
% end
%%
% % 绘制原始数据图像总图3x3
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=1:9
% subplot(3,3,Gi)
% hold on
% for ti=1:6
% x=Tensile_x{Gi,ti}(Tensile_start_line(Gi,ti):Tensile_end_line(Gi,ti));x=x-x(1);
% y=Tensile_y{Gi,ti}(Tensile_start_line(Gi,ti):Tensile_end_line(Gi,ti));y=y-y(1);
% Tensile_x{Gi,ti}=x;
% Tensile_y{Gi,ti}=y;
% plot(Tensile_y{Gi,ti},LineWidth=1.5)
% % plot(Tensile_x{Gi,ti},Tensile_y{Gi,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[Tensile_start_line(Gi,ti)])
% % plot(Tensile_y{Gi,ti},LineWidth=1.5)
% end
% xlabel("Strain (%)")
% ylabel("Stress (MPa)")
% % xlim([0 0.1])
% % ylim([0,70])
% % legend(Location="northwest",Box="off");
% box('on');
% title(['Tensile ',num2str(Gi)])
% set(gca,'FontName','Arial','FontSize',FontSize)
% end
%%
%
% % 分组绘制原始图像9张单独图（应力）
% for Gi=[1 2 3 7 8]
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 40 28]);
% for ti=1:6
% subplot(2,3,ti)
% % [~,maxy]=max(Tensile_y{Gi,ti}(1:Tensile_end_line(Gi,ti)));
% % Tensil_maxy_line(Gi,ti)=maxy;
% plot(Tensile_y{Gi,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[Tensile_start_line(Gi,ti),Tensil_maxy_line(Gi,ti),Tensile_end_line(Gi,ti)])
% xlabel("N")
% ylabel("Stress (MPa)")
% % xlim([0 10])
% % ylim([0,70])
% % legend(Location="northwest",Box="off");
% box('on');
%
% title(['Tensile ',num2str(Gi), '-',num2str(ti)])
% set(gca,'FontName','Arial','FontSize',FontSize)
% end
% end
% 绘制受拉总图
%%
Tensile_maxy=zeros(9,6);
clear fh1
fh1=figure();
set(gcf,'unit','centimeters','Position',[15 15 15 15]);
axes('parent',fh1,'unit','centimeters','position',[1.5 1.5 12.5 12.5]);
hold on
LineStyles=repmat(["-",":","--"],1,3);
for Gi=1:9
for ti=1:6
x=Tensile_x{Gi,ti}(Tensile_start_line(Gi,ti):Tensile_end_line(Gi,ti));x=x-x(1);
y=Tensile_y{Gi,ti}(Tensile_start_line(Gi,ti):Tensile_end_line(Gi,ti));y=y-y(1);
Tensile_maxy(Gi,ti)=max(y);
plot(x,y,LineWidth=1.5,Color=colors{Gi},LineStyle=LineStyles(Gi))
clear x y
end
xlabel("Strain (%)")
ylabel("Stress (MPa)")
xlim([0 8])
ylim([0,60])
legend(legends9,Location="northeast",Box="off");
box('on');
set(gca,'FontName','Arial','FontSize',FontSize)
end
%% 压缩1屈服强度
compressive1_x=cell(9,6); % read strain to compressive1_x
compressive1_y=cell(9,6); % read stress to compressive1_y
Com1_L0=25.4; % 标距L0
Com1_S0=12.7*12.7*pi*0.25; % 横截面积S0
% 读取原始数据并储存
% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 2);
% 指定工作表和范围
opts.Sheet = "Sheet1";
opts.DataRange = "A2:B2000";
% 指定列名称和类型
opts.VariableNames = ["LoadValue", "PositionValue", "PlayTime", "ExtendValue"];
opts.VariableTypes = ["double", "double", "double", "double"];
for Gi=19:27
for ti=1:6
filesname=sprintf("/Users/mawenqian/Documents/硕士/VSL/科研/小论文1/基本试件打印参数正交试验/G%d/N%d.xls",Gi,ti);
readdata = readtable(filesname, opts);
% 转换为输出类型
readdata = table2array(readdata);
compressive1_x{Gi-18,ti}=readdata(:,2)/Com1_L0*100; % 应变
compressive1_y{Gi-18,ti}=readdata(:,1)/Com1_S0; % 应力
clear readdata filesname
end
end
clear opts
% 通过下方原始图像人工判断获得，亦可绘制后展示。
compressive1_start_line=[
73 78 74 77 12 1
77 87 79 78 75 77
82 82 80 82 86 83
369 70 64 72 72 74
71 67 66 63 64 66
70 70 68 67 69 69
62 6 6 5 4 4
3 3 3 181 1 7
65 12 5 4 5 64];
compressive1_end_line=[
471 413 432 390 303 322
422 432 413 380 380 378
471 442 480 390 303 322
827 529 493 497 503 482
511 526 453 475 496 496
513 509 528 486 473 494
528 516 475 414 463 470
463 454 519 605 463 457
619 497 496 365 420 478];
compressive1_maxy_line=[310 307 318 314 234 219
354 377 344 374 352 345
321 325 333 324 303 322
728 412 381 403 404 404
357 350 340 344 345 347
349 330 332 341 329 339
331 295 293 279 286 314
306 285 300 499 313 293
436 274 300 289 298 323];
% % 绘制原始数据图像总图3x3
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=19:27
% subplot(3,3,Gi-18)
% hold on
% for ti=1:6
% [~,maxy]=max(compressive1_y{Gi-18,ti}(1:compressive1_end_line(Gi-18,ti)));
% plot(compressive1_x{Gi-18,ti},compressive1_y{Gi-18,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[maxy])
% % plot(compressive1_y{Gi-18,ti},LineWidth=1.5)
% end
% xlabel("Strain (%)")
% ylabel("Stress (MPa)")
% % xlim([0 10])
% % ylim([0,70])
% % legend(Location="northwest",Box="off");
% box('on');
% title(['compressive1 ',num2str(Gi-18)])
% set(gca,'FontName','Arial','FontSize',FontSize)
% end
% % 分组绘制原始图像9张单独图（应力）
% for Gi=27:-1:19
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 40 28]);
% for ti=1:6
% subplot(2,3,ti)
% [~,maxy]=max(compressive1_y{Gi-18,ti}(1:compressive1_end_line(Gi-18,ti)));
% compressive1_maxy_line(Gi,ti)=maxy;
% plot(compressive1_y{Gi-18,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[compressive1_start_line(Gi-18,ti),compressive1_maxy_line(Gi,ti),compressive1_end_line(Gi-18,ti)])
% xlabel("N")
% ylabel("Stress (MPa)")
% % xlim([0 10])
% % ylim([0,70])
% % legend(Location="northwest",Box="off");
% box('on');
%
% title(['compressive1 ',num2str(Gi-18), '-',num2str(ti)])
% set(gca,'FontName','Arial','FontSize',FontSize)
% end
% end
% 绘制受压总图
compressive1_maxy=zeros(9,6);
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=19:27
% Gi=Gi-18;
% subplot(3,3,Gi)
%
% hold on
% for ti=1:6
% compressive1_maxy(Gi,ti)=compressive1_y{Gi,ti}(compressive1_maxy_line(Gi,ti))-compressive1_x{Gi,ti}(compressive1_start_line(Gi,ti));
% x=compressive1_x{Gi,ti}(compressive1_start_line(Gi,ti):compressive1_end_line(Gi,ti));x=x-x(1);
% y=compressive1_y{Gi,ti}(compressive1_start_line(Gi,ti):compressive1_end_line(Gi,ti));y=y-y(1);
% compressive1_maxy(Gi,ti)=max(y);
% plot(x,y,LineWidth=1.5,Color=colors{Gi})
% end
% xlabel("Strain (%)")
% ylabel("Stress (MPa)")
% xlim([0 8])
% ylim([0,60])
%
% legend(legends(Gi),Location="north",Box="off");
% box('on');
% title(['Compressive ',num2str(Gi)])
% set(gca,'FontName','Arial','FontSize',FontSize)
% end
clear fh1
fh1=figure();
set(gcf,'unit','centimeters','Position',[15 15 15 15]);
axes('parent',fh1,'unit','centimeters','position',[1.5 1.5 12.5 12.5]);
hold on
LineStyles=repmat(["-",":","--"],1,3);
for Gi=19:27
Gi=Gi-18;
hold on
for ti=1:6
x=compressive1_x{Gi,ti}(compressive1_start_line(Gi,ti):compressive1_end_line(Gi,ti));x=x-x(1);
y=compressive1_y{Gi,ti}(compressive1_start_line(Gi,ti):compressive1_end_line(Gi,ti));y=y-y(1);
compressive1_maxy(Gi,ti)=max(y);
plot(x,y,LineWidth=1.5,Color=colors{Gi},LineStyle=LineStyles(Gi))
end
xlabel("Strain (%)")
ylabel("Stress (MPa)")
xlim([0 8])
ylim([0,60])
legend(legends9,Location="southeast",Box="off");
box('on');
% title('Compressive')
set(gca,'FontName','Arial','FontSize',FontSize)
end
%% 弯曲
Flexural_x=cell(9,6); % read strain to Flexural_x
Flexural_y=cell(9,6); % read stress to Flexural_y
% 读取原始数据并储存
% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);
% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";
% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];
% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
L0=6*3.2/52^2; % 标距L0
S0=3*52/(2*12.7*3.2^2); % 横截面积S0
for Gi=10:18
for ti=1:6
filesname=sprintf("/Users/mawenqian/Documents/硕士/VSL/科研/小论文1/基本试件打印参数正交试验/G%d/N%d.txt",Gi,ti);
readdata = readtable(filesname, opts);
% 转换为输出类型
readdata = table2array(readdata);
Flexural_x{Gi-9,ti}=readdata(:,4).*L0*100; % 应变 6sh/L^2
Flexural_y{Gi-9,ti}=readdata(:,1).*S0; % 应力 3Fl/2bh^2
clear readdata filesname
end
end
clear opts
% 通过下方原始图像人工判断获得，亦可绘制后展示。
Flexural_start_line=[
467 518 522 779 368 484
519 337 270 390 415 252
517 357 717 515 553 654
770 553 453 344 498 318
498 521 569 667 510 417
360 431 407 364 472 358
308 493 677 617 648 471
379 529 641 781 915 439
693 701 745 696 734 744];
Flexural_maxy_line=[
4554 4543 4534 4729 4459 4457
4391 4252 4080 4298 4264 4114
4496 4306 4691 4533 4453 4458
4431 4205 4134 3991 4242 4116
4529 4633 4710 4794 4624 4509
4078 4153 4095 4148 4172 4140
4059 4267 4356 4334 4289 4146
4165 4235 4366 4562 4611 4207
4725 4725 4725 4725 4725 4725];
% % 绘制原始数据图像总图3x3
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=10:18
% subplot(3,3,Gi-9)
% hold on
% for ti=1:6
% plot(Flexural_x{Gi-9,ti},Flexural_y{Gi-9,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[Flexural_start_line(Gi-9,ti)])
% % plot(Flexural_y{Gi-9,ti},LineWidth=1.5)
% end
% xlabel("Strain (%)")
% ylabel("Stress (MPa)")
% % xlim([0 10])
% % ylim([0,70])
% % legend(Location="northwest",Box="off");
% box('on');
% title(['Flexural ',num2str(Gi)])
% set(gca,'FontName','Arial','FontSize',FontSize)
% end
%
% % 分组绘制原始图像9张单独图（应力）
% for Gi=18:-1:10
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 40 28]);
% for ti=1:6
% subplot(2,3,ti)
% [~,maxy]=max(Flexural_y{Gi-9,ti});
% Flexural_maxy_line(Gi-9,ti)=maxy;
% plot(Flexural_y{Gi-9,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[Flexural_start_line(Gi-9,ti),Flexural_maxy_line(Gi-9,ti)])
% xlabel("N")
% ylabel("Stress (MPa)")
% % xlim([0 10])
% % ylim([0,70])
% % legend(Location="northwest",Box="off");
% box('on');
%
% title(['Flexural ',num2str(Gi-9), '-',num2str(ti)])
% set(gca,'FontName','Arial','FontSize',FontSize)
% end
% end
% 绘制受弯总图
Flexural_maxy=zeros(9,6);
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=10:18
% Gi=Gi-9;
% subplot(3,3,Gi)
%
% hold on
% for ti=1:6
% Flexural_maxy(Gi,ti)=Flexural_y{Gi,ti}(Flexural_maxy_line(Gi,ti))-Flexural_x{Gi,ti}(Flexural_start_line(Gi,ti));
% x=Flexural_x{Gi,ti}(Flexural_start_line(Gi,ti):end);x=x-x(1);
% y=Flexural_y{Gi,ti}(Flexural_start_line(Gi,ti):end);y=y-y(1);
% Flexural_maxy(Gi,ti)=max(y);
% plot(x,y,LineWidth=1.5,Color=colors{Gi})
% end
% xlabel("Strain (%)")
% ylabel("Stress (MPa)")
% xlim([0 8])
% ylim([0,70])
%
% legend(legends(Gi),Location="north",Box="off");
% box('on');
% title(['Flexural ',num2str(Gi)])
% set(gca,'FontName','Arial','FontSize',FontSize)
% end
%
%
clear fh1
fh1=figure();
set(gcf,'unit','centimeters','Position',[15 15 15 15]);
axes('parent',fh1,'unit','centimeters','position',[1.5 1.5 12.5 12.5]);
hold on
LineStyles=repmat(["-",":","--"],1,3);
for Gi=10:18
Gi=Gi-9;
hold on
for ti=1:6
x=Flexural_x{Gi,ti}(Flexural_start_line(Gi,ti):end);x=x-x(1);
y=Flexural_y{Gi,ti}(Flexural_start_line(Gi,ti):end);y=y-y(1);
Flexural_maxy(Gi,ti)=max(y);
plot(x,y,LineWidth=1.5,Color=colors{Gi},LineStyle=LineStyles(Gi))
end
end
xlabel("Strain (%)")
ylabel("Stress (MPa)")
xlim([0 8])
ylim([0,70])
legend(legends9,Location="southeast",Box="off");
box('on');
% title(['Flexural'])
set(gca,'FontName','Arial','FontSize',FontSize)
%% 汇总柱状图
% close all
com_bar=figure();
set(gcf,'unit','centimeters','Position',[5 5 28 28]);
axes('parent',com_bar,'unit','centimeters','position',[1.5 1.5 25 25],'fontsize',20);
x=1:9;
all_mean=reshape(mean([Tensile_maxy;compressive1_maxy;Flexural_maxy],2),9,3)';
all_mean_std=reshape(std([Tensile_maxy;compressive1_maxy;Flexural_maxy],0,2),9,3)';
bar_fig= bar(x,all_mean,FaceColor="flat",BarWidth=0.8);
hold on
for typei=1:3
set(bar_fig(typei),'FaceColor',color3s{typei});
errorbar(x+0.225*(typei-2),all_mean(typei,:),all_mean_std(typei,:),'LineStyle','none',LineWidth=1.5,Color=[0.5 0.5 0.5]);
xtips = bar_fig(typei).XEndPoints;
ytips = bar_fig(typei).YEndPoints+all_mean_std(typei,:);
label = string(vpa(bar_fig(typei).YData,3));
text(xtips,ytips,label,'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',12)
end
legend(["Tensile Yield Strength" "Compressive Yield Strength" "Flexural Yield Strength" ],Location="northwest",box="on")
% xlabel("Type")
ylabel("Stress (MPa)")
set(gca,'FontName','Arial','FontSize',14)
%% 单独柱状图
figure();
ti=1;
set(gcf,'unit','centimeters','Position',[3 3 15 15]);
axes('unit','centimeters','position',[1.5 1.5 12.5 12.5]);
bar_fig_i= bar(x,all_mean(ti,:),BarWidth=0.8,FaceColor=color2);
hold on
errorbar(x,all_mean(ti,:),all_mean_std(ti,:),'LineStyle','none',LineWidth=1.5,Color=[0.5 0.5 0.5]);
xtips = bar_fig_i.XEndPoints;
ytips = bar_fig_i.YEndPoints+all_mean_std(ti,:);
label = string(vpa(bar_fig_i.YData,3));
text(xtips,ytips,label,'HorizontalAlignment','center','VerticalAlignment','bottom')
ylabel("Stress (MPa)")
set(gca,'FontName','Arial','FontSize',FontSize)
%% 保存计算获得的应力矩阵
all_sigma=[Tensile_maxy(:),compressive1_maxy(:),Flexural_maxy(:)];
save sigma_mat all_sigma all_mean all_mean_std Tensile_maxy compressive1_maxy Flexural_maxy
%% 加载应力矩阵
load('sigma_mat.mat')
A=repmat([1 1 1 2 2 2 3 3 3]',6,1);
B=repmat([1 2 3 1 2 3 1 2 3]',6,1);
C=repmat([1 2 3 2 3 1 3 1 2]',6,1);
D=repmat([1 2 3 3 1 2 2 3 1]',6,1);
save taguchi9_mat A B C D % 保存相关变量
%% SN方差分析
sn=-10*log10(1./all_sigma.^2);
anova_soure_tbl= [{'A','B','C','D','Tensile','Compressive','Flexural'};[num2cell([A B C D]), num2cell(sn)]];
all_anova=cell(3,1);
for ni=1:3
[~,tbli]=anovan(sn(:,ni),{A,B,C,D},"varnames",{'A','B','C','D'});
tbli{1,1}='Factor';
percent_source=cell2mat(tbli(2:6,5));
percent=[{'Contribution'};num2cell(round(100*[percent_source;sum(percent_source)]./sum(percent_source),2))];
all_anova{ni}=[tbli(:,[1 2 3 5 6 7]),percent];
end
save anovan_mat all_anova anova_soure_tbl
%% 信噪比分析（采用越大越好计算公式）
all_sigma2square=1./all_sigma.^2;
all_sn=cell(3,1);
% all_ki=cell(3,1);
for ni=1:3
[tbl_A,grp_A]=grpstats(all_sigma2square(:,ni),A,["mean","gname"]);
[tbl_B,grp_B]=grpstats(all_sigma2square(:,ni),B,["mean","gname"]);
[tbl_C,grp_C]=grpstats(all_sigma2square(:,ni),C,["mean","gname"]);
[tbl_D,grp_D]=grpstats(all_sigma2square(:,ni),D,["mean","gname"]);
sni=-10*log10([tbl_A,tbl_B,tbl_C,tbl_D]);
sni_delta=max(sni)-min(sni);
[~,sni_sort]=sort(sni_delta,'descend'); % 降序排列
[~,sni_sort]=sort(sni_sort); % 获得排位
all_sn{ni}=cell(size(sni)+[3 1]);
all_sn{ni}(1,:)={'Level' 'A','B','C','D'};
all_sn{ni}(2:end,1)={'1' '2' '3' 'Delta' 'Rank'};
all_sn{ni}(2:end,2:end)=[num2cell([sni;sni_delta]);num2cell(sni_sort)];
% 直观分析计算过程（不采用）
% tbl_A=grpstats(all_sigma(:,ni),A,@(x)sum(x,1));
% tbl_B=grpstats(all_sigma(:,ni),B,@(x)sum(x,1));
% tbl_C=grpstats(all_sigma(:,ni),C,@(x)sum(x,1));
% tbl_D=grpstats(all_sigma(:,ni),D,@(x)sum(x,1));
% all_ki{ni}=[tbl_A,tbl_B,tbl_C,tbl_D];
end
save all_sn_mat all_sn
%% 绘制信噪比SN曲线
close all
which_sigma=1; % 1:Tensile 2:Compressive1 3:Flexural
y1=cell2mat(all_sn{which_sigma}(2:4,2));
y2=cell2mat(all_sn{which_sigma}(2:4,3));
y3=cell2mat(all_sn{which_sigma}(2:4,4));
y4=cell2mat(all_sn{which_sigma}(2:4,5));
% y1=[0 0 0 ];
% y2=[0 0 0 ];
% y3=[0 0 0 ];
% y4=[0 0 0 ];
plot([y1;y2;y3;y4])
ylim_range=[18 32]; % set ylim range from the upper figure
FontSize=13;
fh1=figure();
set(gcf,'unit','centimeters','Position',[2 2 23 11]);
axes('parent',fh1,'unit','centimeters','position',[2 2 5 8],'fontsize',10);
plot([1 2 3],y1,'linewidth',3,'Marker','.','MarkerSize',30);
xlim([0.5 3.5])
ylim(ylim_range)
set(gca,'XTick',[1 2 3],'yTickLabel',num2str(get(gca,'yTick')','%.1f'));
xlabel("Filament","FontSize",FontSize,'FontName','Arial')
ylabel("Mean of SN Radio","FontSize",FontSize,'FontName','Arial')
box off
axes('parent',fh1,'unit','centimeters','position',[7 2 5 8],'fontsize',10);
plot([240 250 260],y2,'linewidth',3,'Marker','.','MarkerSize',30);
xlim([235 265])
ylim(ylim_range)
set(gca,'XTick',[240 250 260],'YTickLabel','')
xlabel("Speed(mm/s)","FontSize",FontSize,'FontName','Arial')
box off
axes('parent',fh1,'unit','centimeters','position',[12 2 5 8],'fontsize',10);
plot([210 215 220],y3,'linewidth',3,'Marker','.','MarkerSize',30)
xlim([205 225])
ylim(ylim_range)
set(gca,'XTick',[210 215 220],'YTickLabel','')
xlabel("Temperature(\circC)","FontSize",FontSize,'FontName','Arial')
box off
axes('parent',fh1,'unit','centimeters','position',[17 2 5 8],'fontsize',10);
plot([0.16 0.20 0.24],y4,'linewidth',3,'Marker','.','MarkerSize',30)
xlim([0.14 0.26])
ylim(ylim_range)
set(gca,'XTick',[0.16 0.20 0.24],'YTickLabel','')
xlabel("Layer thickness(mm)","FontSize",FontSize,'FontName','Arial')
box off
ax2 = axes('parent',fh1,'unit','centimeters','position',[2 2 20 8],...
'XAxisLocation','top',...
'YAxisLocation','right',...
'Color','none',...
'XColor','k','YColor','k');
set(ax2,'YTick', []);
set(ax2,'XTick', []);
%% 线性回归预测
X=[ones(size(A)) A B C D];
all_regress=cell(3,1);
for ni=1:3
[b,~,~,~,stats] = regress(all_sigma(:,ni),X);
all_regress{ni}=[{'' 'A' 'B' 'C' 'D' 'R^2' 'p-value'}',num2cell([b;stats([1 3])'])];
end
save all_regress_mat all_regress
{% endfolding %}
拉压弯曲线绘图2022/10/27
{% folding child:codeblock open:flase color:yellow 点击展开代码详情 %}
clc;clear
close all
% color1=[77,133,189]/255;
% color2=[247,144,61]/255;
% color3=[89,169,90]/255;
color1=[147 213 220]/255;
color2=[242 202 201]/255;
color3=[198 223 200]/255;
FontSize=12;
%% 拉压弯标准对比图
% 试验时间：2022年10月19-24日
% 3个子图
fh1=figure(1);
% Tensile
set(gcf,'unit','centimeters','Position',[5 5 33 13]);
axes('parent',fh1,'unit','centimeters','position',[3,3,8,8],'fontsize',10);
hold on
% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);
% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";
% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];
% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% 标距L0
L0=[ones(1,10).*50,ones(1,5)*7.62];
% 横截面积S
S=[ones(1,5).*40,ones(1,5).*4*13,ones(1,5).*4*3.18];
% 曲线颜色
% colorrgb=[ones(5,1)*color1;ones(5,1)*color2;ones(5,1)*color3];
colorrgb=[repmat(color1,5,1);repmat(color2,5,1);repmat(color3,5,1)];
maxf_T=zeros(5,3);
for filesnum=1:15
filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/Tensile/N%d.txt",filesnum);
readdata = readtable(filesname, opts);
% 转换为输出类型
readdata = table2array(readdata);
x=readdata(:,4)./L0(filesnum)*100; % 应变
y=readdata(:,1)./S(filesnum); % 应力
plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
maxf_T(filesnum)=max(y);
end
ave_stress_T=sum(maxf_T,1)/5;
legendtext=["ISO 527 1B",repmat([""],1,4) ,"ASTM D638 Type I",repmat([""],1,4),"ASTM D638 Type V"];
xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
title('Tensile')
set(gca,'FontName','Arial','FontSize',FontSize)
clear opts
clear colorrgb
% Compressive
axes('parent',fh1,'unit','centimeters','position',[13,3,8,8],'fontsize',10);
hold on
% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 4);
% 指定工作表和范围
opts.Sheet = "Sheet1";
opts.DataRange = "A2:D5860";
% 指定列名称和类型
opts.VariableNames = ["LoadValue", "PositionValue", "PlayTime", "ExtendValue"];
opts.VariableTypes = ["double", "double", "double", "double"];
% 标距L0
L0=ones(1,10).*50.8;
% 横截面积S
S=[ones(1,5).*12.7*12.7,ones(1,5).*12.7*12.7*pi*0.25];
colorrgb=[repmat(color2,5,1);repmat(color3,5,1)];
maxf_C=zeros(5,2);
for filesnum=1:10
filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/Compressive/N%d.xls",filesnum);
readdata = readtable(filesname, opts, "UseExcel", false);
% 转换为输出类型
readdata = table2array(readdata);
x=readdata(:,2)./L0(filesnum)*100; % 应变
y=readdata(:,1)./S(filesnum); % 应力
plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
maxf_C(filesnum)=max(y);
end
ave_stress_C=sum(maxf_C,1)/5;
legendtext=["ASTM D695 Prise",repmat([""],1,4) ,"ASTM D695 Cylinder"];
xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
title('Compressive')
set(gca,'FontName','Arial','FontSize',FontSize)
clear opts
clear colorrgb
% Flexural
axes('parent',fh1,'unit','centimeters','position',[23,3,8,8],'fontsize',10);
hold on
% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);
% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";
% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];
% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% 标距L0
L0=[ones(1,5).*6*4/56^2,ones(1,5)*6*3.2/52^2];
% 横截面积S
S=[ones(1,5).*3*56/(2*10*4^2),ones(1,5).*3*52/(2*12.7*3.2^2)];
% 曲线颜色
% colorrgb=[ones(5,1)*color1;ones(5,1)*color2;ones(5,1)*color3];
colorrgb=[repmat(color1,5,1);repmat(color2,5,1)];
maxf_F=zeros(5,2);
for filesnum=1:10
filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/Flexural/N-%d.txt",filesnum);
readdata = readtable(filesname, opts);
% 转换为输出类型
readdata = table2array(readdata);
x=readdata(:,4).*L0(filesnum)*100; % 应变 6sh/L^2
y=readdata(:,1).*S(filesnum); % 应力 3Fl/2bh^2
plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
maxf_F(filesnum)=max(y);
end
ave_stress_F=sum(maxf_F,1)/5;
legendtext=["ISO 178",repmat([""],1,4) ,"ASTM D790"];
xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
title('Flexural')
set(gca,'FontName','Arial','FontSize',FontSize)
clear opts
clear colorrgb
%% ASTM标准下压、弯不同打印方向对比图
% 试验时间：2022年10月19-24日
% 3个子图
fh2=figure(2);
% Compressive-1
set(gcf,'unit','centimeters','Position',[5 5 33 13]);
axes('parent',fh2,'unit','centimeters','position',[3,3,8,8],'fontsize',10);
hold on
% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 4);
% 指定工作表和范围
opts.Sheet = "Sheet1";
opts.DataRange = "A2:D5860";
% 指定列名称和类型
opts.VariableNames = ["LoadValue", "PositionValue", "PlayTime", "ExtendValue"];
opts.VariableTypes = ["double", "double", "double", "double"];
% 标距L0
L0=[ones(1,10).*50.8];
% 横截面积S
S=[ones(1,10).*12.7*12.7];
colorrgb=[repmat(color2,5,1);repmat(color3,5,1)];
maxf_CN5=zeros(5,2);
for filesnum=1:10
filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/compressive-diff-directions/N5-%d.xls",filesnum);
readdata = readtable(filesname, opts, "UseExcel", false);
% 转换为输出类型
readdata = table2array(readdata);
x=readdata(:,2)./L0(filesnum)*100; % 应变
y=readdata(:,1)./S(filesnum); % 应力
plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
maxf_CN5(filesnum)=max(y);
end
ave_stress_CN5=sum(maxf_CN5,1)/5;
legendtext=["水平打印",repmat([""],1,4) ,"竖直打印"];
xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
title('ASTM D695 Prise')
set(gca,'FontName','Arial','FontSize',FontSize)
clear opts
clear colorrgb
% Compressive-2
axes('parent',fh2,'unit','centimeters','position',[13,3,8,8],'fontsize',10);
hold on
% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 4);
% 指定工作表和范围
opts.Sheet = "Sheet1";
opts.DataRange = "A2:D5860";
% 指定列名称和类型
opts.VariableNames = ["LoadValue", "PositionValue", "PlayTime", "ExtendValue"];
opts.VariableTypes = ["double", "double", "double", "double"];
% 标距L0
L0=[ones(1,10).*50.8];
% 横截面积S
S=[ones(1,10).*12.7*12.7*pi*0.25];
colorrgb=[repmat(color2,5,1);repmat(color3,5,1)];
maxf_CN6=zeros(5,2);
for filesnum=1:10
filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/compressive-diff-directions/N6-%d.xls",filesnum);
readdata = readtable(filesname, opts, "UseExcel", false);
% 转换为输出类型
readdata = table2array(readdata);
x=readdata(:,2)./L0(filesnum)*100; % 应变
y=readdata(:,1)./S(filesnum); % 应力
plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
maxf_CN6(filesnum)=max(y);
end
ave_stress_CN6=sum(maxf_CN6,1)/5;
legendtext=["水平打印",repmat([""],1,4) ,"竖直打印"];
xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
title('ASTM D695 Cylinder')
set(gca,'FontName','Arial','FontSize',FontSize)
clear opts
clear colorrgb
% Flexural
axes('parent',fh2,'unit','centimeters','position',[23,3,8,8],'fontsize',10);
hold on
% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);
% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";
% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];
% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% 标距L0
L0=[ones(1,5).*6*4/56^2,ones(1,5)*6*3.2/52^2];
% 横截面积S
S=[ones(1,5).*3*56/(2*10*4^2),ones(1,5).*3*52/(2*12.7*3.2^2)];
% 曲线颜色
% colorrgb=[ones(5,1)*color1;ones(5,1)*color2;ones(5,1)*color3];
colorrgb=[repmat(color2,5,1);repmat(color3,5,1)];
maxf_FN8=zeros(5,2);
for filesnum=1:10
filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/Flexural-diff-direction/N-%d.txt",filesnum);
readdata = readtable(filesname, opts);
% 转换为输出类型
readdata = table2array(readdata);
x=readdata(:,4).*L0(filesnum)*100; % 应变 6sh/L^2
y=readdata(:,1).*S(filesnum); % 应力 3Fl/2bh^2
plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
maxf_FN8(filesnum)=max(y);
end
ave_stress_FN8=sum(maxf_FN8,1)/5;
legendtext=["水平打印",repmat([""],1,4) ,"竖直打印"];
xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
title('ASTM D790')
set(gca,'FontName','Arial','FontSize',FontSize)
clear opts
clear colorrgb
{% endfolding %}
拉伸曲线2022/06/05
{% folding child:codeblock open:flase color:yellow 点击展开代码详情 %}
clc;clear
close all
out_data_file=[...
4 5 6
7 8 9
10 11 12
13 14 15
16 19 20
21 22 23
27 28 29
41 30 31
32 33 34
24 25 26
35 36 37
38 39 40];
group_name={'ITIL-折3';'ITIL-折6';'ITIL-折9';'ITIL-直-3';'ITIL-直-6';'ITIL-直-9';'Zebra-3';'Zebra-6';'Zebra-9';'Theta-15';'复合ITIL-直-9';'复合ITIL-折-9'};
data_size=size(out_data_file);
stress=zeros(data_size);
maxf=[];
ave_xys={};
h1=figure(1);
colorrgb=[77,133,189;247,144,61;89,169,90]/255;
for rows=1:data_size(1)
subplot(4,3,rows)
hold on
legendtext=[];
xys={};
for cols=1:data_size(2)
readdata = readmatrix(sprintf("/Users/wenqian/Documents/硕士/VSL/2022.5/PLA- PLA拉伸试验/20220510.is_tens_Exports/20220510_%d.csv",out_data_file(rows,cols)));
[x,y,stress(rows,cols)]=filter_date(readdata(:,2),readdata(:,3));
% N->kN
y=1000*y;
stress(rows,cols)=stress(rows,cols)*1000;
plot(x,y,LineWidth=1.5,Color=colorrgb(cols,:))
legendtext=[legendtext,sprintf("%s-%d",cell2mat(group_name(rows)),cols)];
maxf(rows,cols)=max(y);
xys=[xys;{x y}];
end
[ax,ay]=average_date(xys);
ave_xys=[ave_xys;{ax ay}];
plot(ax,smooth(ay),LineWidth=2.5);
legendtext=[legendtext,"平均曲线"];
xlabel("Displacement (mm)")
ylabel("Force (N)")
xlim([0 15])
ylim([0,1250])
legend(legendtext,Location="southeast",Box="off");
box('on');
title(sprintf('(%s)',char(96+rows)),'position',[7.5,-600])
set(gca,'FontName','Arial','FontSize',14)
end
ave_force=mean(maxf,2);
eb=std(stress,0,2);
fig2_date_index={[1 2 3] [4 5 6] [7 8 9 ] [11 12]};
draw_from_index(ave_xys,group_name,fig2_date_index,2,4)
draw_from_index_bar(ave_force,eb,group_name,fig2_date_index,2,4)
fig3_date_index={[1 4 7] [2 5 8] [3 6 9 11 12 ]};
draw_from_index(ave_xys,group_name,fig3_date_index,2,4)
draw_from_index_bar(ave_force,eb,group_name,fig3_date_index,2,4)
function draw_from_index(fig_date_xys,fig_titles,fig_date_index,m,n)
fig_size=length(fig_date_index);
h=figure();
for i=1:fig_size
subplot(m,n,i)
hold on
for j=1:length(fig_date_index{i})
subfig_index=fig_date_index{i};
index_num=subfig_index(j);
x=fig_date_xys{index_num,1};
y=fig_date_xys{index_num,2};
plot(x,y,LineWidth=1.5)
legend_title=fig_titles(subfig_index);
end
xlabel("Displacement (mm)")
ylabel("Force (N)")
xlim([0 15])
ylim([0,1200])
legend(legend_title,Location="southeast",Box="off");
box('on');
title(sprintf('(%s)',char(96+i)),'position',[7.5,-300])
set(gca,'FontName','Arial','FontSize',14)
end
end
function draw_from_index_bar(bar_date,ebbar_date,fig_titles,fig_date_index,m,n)
fig_size=length(fig_date_index);
h=figure();
for i=1:fig_size
subplot(m,n,i)
hold on
subfig_index=fig_date_index{i};
x=categorical(fig_titles(subfig_index));
x=reordercats(x,fig_titles(subfig_index));
y=bar_date(subfig_index);
eb=ebbar_date(subfig_index);
bar_fig= bar(x,y,1.5*length(y)/10,FaceColor="flat");
errorbar(x,y,eb,'LineStyle','none',LineWidth=1.5);
% legend_title=fig_titles(subfig_index);
% xlabel("Type")
ylabel("Force (N)")
% xlim([0 60])
ylim([0,1300])
% legend(legend_title,Location="southeast",Box="off");
date=y;
colorrgb=[77,133,189;247,144,61;89,169,90;77,133,189;247,144,61;89,169,90]/255;
for num=1:size(date)
bar_fig.CData(num,:)=colorrgb(num,:);
text_date=string(round(date(num),2));
text(num,date(num)+eb(num)+0.1,text_date,'VerticalAlignment','bottom','HorizontalAlignment','center')
end
box('on');
title(sprintf('(%s)',char(96+i)),Position=[length(y)/2+0.5 -200])
set(gca,'FontName','Arial','FontSize',14)
end
end
{% endfolding %}