How can I plot on the same four figures on my code which changes while list numbers changing?

1 次查看（过去 30 天）
mehmet salihi
mehmet salihi 2022-1-1
编辑： mehmet salihi 2022-1-1
Hello Dear frinends.
I have 15 different configuration + base one. each configuration gives me four plots.
I am trying to plot the figures of each configuration to see all 16 cases in the same figure to see the difference. maybe the below photo will give you idea what I mean.
my For loop starting with k, also excel list is uploaded for the base and distribution in the same file.
You can see my base configuration after deleting for loop with k index.
close all
clear all
clc
%%% Fixed Parameters
B = 2; % Number of Blades
R = 10.06; % Blade Radius in meter
Rho = 1.225; % Density of air kg/m3
w = 7.501; % Rotational Speed rad/s
Theta_P = -3; % Tip Pitch Angle in degree
ac = 0.2; % Glauert correction
eps = 1e-5; % Error tolerance
%%% Airfoil Date
afoil = readmatrix('airfoil data5e5.xlsx');
Aoa = afoil(:,1);
Cl = afoil(:,2);
Cd = afoil(:,3);
Velocity =[5,6,7,8,9,11]; % Wind Speed m/s
MTotal=[];
%%% Torque vs Velocity Data
torquee = readmatrix('TorqueWindSpeedExperi.xlsx');
vupper=torquee(1:5,1); % Radial Distance in meter
vlower=torquee(6:10,1); % Chord length in meter
torquelower=torquee(1:5,2); % Twist in degree
torqueupper=torquee(6:10,2);
figure(1)
hold on
plot(vupper,torqueupper,'--r',vlower,torquelower,'--r')
grid on
%%% Given Blade Data
ns=26;
sections = readmatrix('chord and twist 15 distributions.xlsx');
for k=3:2:33
r=sections(1:26,2); % Radial Distance in meter
SS=sections(1:26,3); % Span Station %
Chord=sections(1:26,k+1); % Chord length in meter
Beta=sections(1:26,k+2); % Twist in degree
% disp(' ');disp(' BEM Theory ')
% disp('Velocity (m/s) Torque (J) Power (KW)')
% disp(' ')
for j=1:length(Velocity)
for i = 8:length(r)
sigma(i) = (B*Chord(i))/(2*pi*r(i));
a=0;
a_prime=0;
alist(i)=[0];
alistprime(i)=[0];
% disp(['Section ',num2str(i)])
for n=1:1000
Phi = atand(((1-a)*Velocity(j))/((1+a_prime)*w*r(i))); % Flow Angle
Theta = Beta(i) + Theta_P ;
aoa = Phi- Theta; % local angle of attack in degree
f = (B/2)* ((R-r(i))/(r(i)*sind(Phi)));
F = (2/pi)*acos(exp(-f)); % Prandtl correction factor
cl=interp1(Aoa,Cl,aoa,'linear','extrap');
cd=interp1(Aoa,Cd,aoa,'linear','extrap');
Cn = cl*cosd(Phi)+cd*sind(Phi);
Ct =cl*sind(Phi)-cd*cosd(Phi);
K=(4*F*(sind(Phi))^2)/(sigma(i)*Cn);
%%%%%%%%%% a and ac condition
a = 1/(((4*F*(sind(Phi)^2))/(sigma(i)*Cn))+1);
if a < ac || a==ac
a=a;
else
a = 0.5*(2+K*(1-2*ac)-sqrt(((K*(1-2*ac)+2)^2)+4*((K*(ac^2))-1)));
end
a_prime =1/(((4*F*(sind(Phi)*cosd(Phi)))/(sigma(i)*Ct))-1);
alist(i,n+1)=a;
alistprime(i,n+1)=a_prime;
% disp([" a ",a,' a_old ',alist(i,n)," a_prime ",a_prime,' a_Prime_old ',alistprime(i,n)])
%%% condition for epsilon or tolerance
if a-alist(i,n) < eps && a_prime-alistprime(i,n) < eps
break;
end
end
%%% V relative
v_rel(i) = sqrt( (Velocity(j)*(1-alist(end)))^2 + (w*r(i)*(1+alistprime(end)))^2 );
%%% The tangential force per length Pt
Pt(i) = 0.5*Rho*(v_rel(i)^2)*Chord(i)*Ct;
end
for i=8:length(Pt)-1
Ai(i) = (Pt(i+1)-Pt(i) ) / (r(i+1)-r(i));
Bi(i) = ((Pt(i)*r(i+1))-(Pt(i+1)*r(i)) ) / (r(i+1)-r(i));
Pt_total(i) = Ai(i)*r(i)+Bi(i);
end
% disp(['Pt_total at ',num2str(V0(j)),' m/s is ',num2str(sum(Pt_total))])
for i=8:length(Pt)-1
M(i) = ((1/3)*Ai(i)*((r(i+1)^3)-(r(i)^3))) + ((0.5)*Bi(i)*((r(i+1)^2)-(r(i)^2)));
end
% % total shaft torque and power
MTotal(j) = B* sum(M);
Power(j)=MTotal(j)*w*1e-3;
disp([' ',num2str(Velocity(j)),' ',...
num2str(MTotal(j)),' ',num2str(Power(j))])
end
plot(Velocity,MTotal,'--b','LineWidth',1.5)
legend('Experimental upper range ','Experimental lower range ','Baseline computation','Position',[.25, 0.8, .11, .11]); xlabel(' Wind Speed (m/s) '); ylabel(' Torque (Nm) ');
set(gca,'XTick',(4:1:12)); set(gca,'YTick',(0:300:1800))
ylim([0 1800]); xlim([4 12]);
for i =1
if Aoa(i)>-19
title(' Reynolds number = 500.000' )
else
title(' Reynolds number = 1.000.000' )
end
end
VelocityEx =[5,6,7,8,9,11,13,15,17,19,23,25]; % Wind Speed m/s
PowerEx=[2.34,4.03,5.87,7.68,9.62,11.12,9.19,8.93,6.75,6.62,7.62,9.05];
figure(2) % Power
hold on
plot(Velocity,Power,'-or',VelocityEx,PowerEx,'-oblue','LineWidth',1.5);grid on
set(gca,'XTick',(0:5:30)); set(gca,'YTick',(0:2:15))
ylim([0 15]); xlim([0 30]);
xlabel(' Wind Speed (m/s)'); ylabel(' Power(KW) ');
legend('BEM Theory ','Experimental Data','Position',[.25, 0.8, .11, .11]); xlabel(' Wind Speed (m/s) '); ylabel(' Torque (Nm) ');
for i =1
if Aoa(i)>-19
title(' Reynolds number = 500.000' )
else
title(' Reynolds number = 1.000.000' )
end
end
figure(3) % Twist vs Span Station in percent
plot(SS,Beta,'-*r')
set(gca,'XTick',(0:0.2:1)); set(gca,'YTick',(-5:5:25))
ylim([-5 25]); xlim([0 1.2]);
xlabel(' Span Station % '); ylabel(' Twist Angle (degrees) ');
title(' Baseline Twist Angle Distribution ')
grid on
figure(4) % Chord vs Span Station in percent
plot(SS,Chord,'-*b')
set(gca,'XTick',(0:0.2:1)); set(gca,'YTick',(0:0.2:1))
ylim([0 1]); xlim([0 1.2]);
xlabel(' Span Station % '); ylabel(' Chord length ');
title(' Baseline Chord Distribution ')
grid on
end

请先登录，再回答此问题。

采纳的回答

Meg Noah
Meg Noah 2022-1-1
There's something amiss with the xlsx reading, but that wasn't the question. This will update the plots, with a slight delay between time steps to allow for visualization.
close all
clear all
clc
%%% Fixed Parameters
B = 2; % Number of Blades
R = 10.06; % Blade Radius in meter
Rho = 1.225; % Density of air kg/m3
w = 7.501; % Rotational Speed rad/s
Theta_P = -3; % Tip Pitch Angle in degree
ac = 0.2; % Glauert correction
eps = 1e-5; % Error tolerance
%%% Airfoil Date
afoil = readmatrix('airfoil data5e5.xlsx');
Aoa = afoil(:,1);
Cl = afoil(:,2);
Cd = afoil(:,3);
Velocity =[5,6,7,8,9,11]; % Wind Speed m/s
MTotal=[];
%%% Torque vs Velocity Data
torquee = readmatrix('TorqueWindSpeedExperi.xlsx');
vupper=torquee(1:5,1); % Radial Distance in meter
vlower=torquee(6:10,1); % Chord length in meter
torquelower=torquee(1:5,2); % Twist in degree
torqueupper=torquee(6:10,2);
%%% Given Blade Data
ns=26;
sections = readmatrix('chord and twist 15 distributions.xlsx');
flagFirst = 1;
for k=3:2:33
r=sections(1:26,2); % Radial Distance in meter
SS=sections(1:26,3); % Span Station %
Chord=sections(1:26,k+1); % Chord length in meter
Beta=sections(1:26,k+2); % Twist in degree
% disp(' ');disp(' BEM Theory ')
% disp('Velocity (m/s) Torque (J) Power (KW)')
% disp(' ')
for j=1:length(Velocity)
for i = 8:length(r)
sigma(i) = (B*Chord(i))/(2*pi*r(i));
a=0;
a_prime=0;
alist(i)=[0];
alistprime(i)=[0];
% disp(['Section ',num2str(i)])
for n=1:1000
Phi = atand(((1-a)*Velocity(j))/((1+a_prime)*w*r(i))); % Flow Angle
Theta = Beta(i) + Theta_P ;
aoa = Phi- Theta; % local angle of attack in degree
f = (B/2)* ((R-r(i))/(r(i)*sind(Phi)));
F = (2/pi)*acos(exp(-f)); % Prandtl correction factor
cl=interp1(Aoa,Cl,aoa,'linear','extrap');
cd=interp1(Aoa,Cd,aoa,'linear','extrap');
Cn = cl*cosd(Phi)+cd*sind(Phi);
Ct =cl*sind(Phi)-cd*cosd(Phi);
K=(4*F*(sind(Phi))^2)/(sigma(i)*Cn);
%%%%%%%%%% a and ac condition
a = 1/(((4*F*(sind(Phi)^2))/(sigma(i)*Cn))+1);
if a < ac || a==ac
a=a;
else
a = 0.5*(2+K*(1-2*ac)-sqrt(((K*(1-2*ac)+2)^2)+4*((K*(ac^2))-1)));
end
a_prime =1/(((4*F*(sind(Phi)*cosd(Phi)))/(sigma(i)*Ct))-1);
alist(i,n+1)=a;
alistprime(i,n+1)=a_prime;
% disp([" a ",a,' a_old ',alist(i,n)," a_prime ",a_prime,' a_Prime_old ',alistprime(i,n)])
%%% condition for epsilon or tolerance
if a-alist(i,n) < eps && a_prime-alistprime(i,n) < eps
break;
end
end
%%% V relative
v_rel(i) = sqrt( (Velocity(j)*(1-alist(end)))^2 + (w*r(i)*(1+alistprime(end)))^2 );
%%% The tangential force per length Pt
Pt(i) = 0.5*Rho*(v_rel(i)^2)*Chord(i)*Ct;
end
for i=8:length(Pt)-1
Ai(i) = (Pt(i+1)-Pt(i) ) / (r(i+1)-r(i));
Bi(i) = ((Pt(i)*r(i+1))-(Pt(i+1)*r(i)) ) / (r(i+1)-r(i));
Pt_total(i) = Ai(i)*r(i)+Bi(i);
end
% disp(['Pt_total at ',num2str(V0(j)),' m/s is ',num2str(sum(Pt_total))])
for i=8:length(Pt)-1
M(i) = ((1/3)*Ai(i)*((r(i+1)^3)-(r(i)^3))) + ((0.5)*Bi(i)*((r(i+1)^2)-(r(i)^2)));
end
% % total shaft torque and power
MTotal(j) = B * sum(M);
Power(j)=MTotal(j)*w*1e-3;
disp([' ',num2str(Velocity(j)),' ',...
num2str(MTotal(j)),' ',num2str(Power(j))])
end
VelocityEx =[5,6,7,8,9,11,13,15,17,19,23,25]; % Wind Speed m/s
PowerEx=[2.34,4.03,5.87,7.68,9.62,11.12,9.19,8.93,6.75,6.62,7.62,9.05];
if (flagFirst == 1)
f1 = figure(1);
subplot(2,2,1);
hold on
h1A = plot(vupper,torqueupper,'--r','DisplayName','Experimental Upper Range');
h1B = plot(vlower,torquelower,'--r','DisplayName','Experimental Lower Range');
grid on
h1C = plot(Velocity,MTotal,'--b','LineWidth',1.5,'DisplayName','Baseline Computatation');
legend('Location','best');
xlabel(' Wind Speed (m/s) ');
ylabel(' Torque (Nm) ');
set(gca,'XTick',(4:1:12)); set(gca,'YTick',(0:300:1800))
ylim([0 1800]); xlim([4 12]);
for i =1
if Aoa(i)>-19
title(' Reynolds number = 500.000' )
else
title(' Reynolds number = 1.000.000' )
end
end
drawnow()
%f2 = figure(2); % Power
subplot(2,2,2);
hold on
h2A = plot(Velocity,Power,'-or','LineWidth',1.5,'DisplayName','BEM Theory');
h2B = plot(VelocityEx,PowerEx,'-oblue','LineWidth',1.5,'DisplayName','Experimental Data');
grid on
set(gca,'XTick',(0:5:30)); set(gca,'YTick',(0:2:15))
ylim([0 15]); xlim([0 30]);
xlabel(' Wind Speed (m/s)'); ylabel(' Power(KW) ');
legend('Location','best');
xlabel(' Wind Speed (m/s) '); ylabel(' Torque (Nm) ');
for i =1
if Aoa(i)>-19
title(' Reynolds number = 500.000' )
else
title(' Reynolds number = 1.000.000' )
end
end
drawnow()
%f3 = figure(3); % Twist vs Span Station in percent
subplot(2,2,3)
h3 = plot(SS,Beta,'-*r','DisplayName','\beta');
set(gca,'XTick',(0:0.2:1)); set(gca,'YTick',(-5:5:25))
ylim([-5 25]); xlim([0 1.2]);
xlabel(' Span Station % '); ylabel(' Twist Angle (degrees) ');
title(' Baseline Twist Angle Distribution ')
grid on
legend('location','best')
drawnow()
%f4 = figure(4); % Chord vs Span Station in percent
subplot(2,2,4)
h4 = plot(SS,Chord,'-*b','DisplayName','Chord');
set(gca,'XTick',(0:0.2:1)); set(gca,'YTick',(0:0.2:1))
ylim([0 1]); xlim([0 1.2]);
xlabel(' Span Station % '); ylabel(' Chord length ');
title(' Baseline Chord Distribution ')
grid on
legend('location','best')
drawnow()
else
h1C.XData = Velocity;
h1C.YData = MTotal;
drawnow()
h2A.XData = Velocity;
h2A.YData = Power;
drawnow()
h2B.XData = VelocityEx;
h2B.YData = PowerEx;
drawnow()
h3.XData = SS;
h3.YData = Beta;
drawnow()
h4.XData = SS;
h4.YData = Chord;
drawnow()
pause(0.1); % delay for viewing
end
flagFirst = 0;
end
  2 个评论
mehmet salihi
mehmet salihi 2022-1-1
I am sorry for that excel file, I fixed it.
it gives all four plots in the same figure and changing like simulation but the are not saved while coming plots coming over and over with different legends.
the code and fixed excel is updated. could you please check it one more time beacuse i am receiving the below figure at the end.
close all
clear all
clc
%%% Fixed Parameters
B = 2; % Number of Blades
R = 10.06; % Blade Radius in meter
Rho = 1.225; % Density of air kg/m3
w = 7.501; % Rotational Speed rad/s
Theta_P = -3; % Tip Pitch Angle in degree
ac = 0.2; % Glauert correction
eps = 1e-5; % Error tolerance
%%% Airfoil Date
afoil = readmatrix('airfoil data5e5.xlsx');
Aoa = afoil(:,1);
Cl = afoil(:,2);
Cd = afoil(:,3);
Velocity =[5,6,7,8,9,11]; % Wind Speed m/s
MTotal=[];
%%% Torque vs Velocity Data
torquee = readmatrix('TorqueWindSpeedExperi.xlsx');
vupper=torquee(1:5,1); % Radial Distance in meter
vlower=torquee(6:10,1); % Chord length in meter
torquelower=torquee(1:5,2); % Twist in degree
torqueupper=torquee(6:10,2);
%%% Given Blade Data
ns=26;
sections = readmatrix('chord and twist 15 distributions.xlsx');
flagFirst = 1;
for k=3:2:33
r=sections(1:26,2); % Radial Distance in meter
SS=sections(1:26,3); % Span Station %
Chord=sections(1:26,k+1); % Chord length in meter
Beta=sections(1:26,k+2); % Twist in degree
% disp(' ');disp(' BEM Theory ')
% disp('Velocity (m/s) Torque (J) Power (KW)')
% disp(' ')
for j=1:length(Velocity)
for i = 8:length(r)
sigma(i) = (B*Chord(i))/(2*pi*r(i));
a=0;
a_prime=0;
alist(i)=[0];
alistprime(i)=[0];
% disp(['Section ',num2str(i)])
for n=1:1000
Phi = atand(((1-a)*Velocity(j))/((1+a_prime)*w*r(i))); % Flow Angle
Theta = Beta(i) + Theta_P ;
aoa = Phi- Theta; % local angle of attack in degree
f = (B/2)* ((R-r(i))/(r(i)*sind(Phi)));
F = (2/pi)*acos(exp(-f)); % Prandtl correction factor
cl=interp1(Aoa,Cl,aoa,'linear','extrap');
cd=interp1(Aoa,Cd,aoa,'linear','extrap');
Cn = cl*cosd(Phi)+cd*sind(Phi);
Ct =cl*sind(Phi)-cd*cosd(Phi);
K=(4*F*(sind(Phi))^2)/(sigma(i)*Cn);
%%%%%%%%%% a and ac condition
a = 1/(((4*F*(sind(Phi)^2))/(sigma(i)*Cn))+1);
if a < ac || a==ac
a=a;
else
a = 0.5*(2+K*(1-2*ac)-sqrt(((K*(1-2*ac)+2)^2)+4*((K*(ac^2))-1)));
end
a_prime =1/(((4*F*(sind(Phi)*cosd(Phi)))/(sigma(i)*Ct))-1);
alist(i,n+1)=a;
alistprime(i,n+1)=a_prime;
% disp([" a ",a,' a_old ',alist(i,n)," a_prime ",a_prime,' a_Prime_old ',alistprime(i,n)])
%%% condition for epsilon or tolerance
if a-alist(i,n) < eps && a_prime-alistprime(i,n) < eps
break;
end
end
%%% V relative
v_rel(i) = sqrt( (Velocity(j)*(1-alist(end)))^2 + (w*r(i)*(1+alistprime(end)))^2 );
%%% The tangential force per length Pt
Pt(i) = 0.5*Rho*(v_rel(i)^2)*Chord(i)*Ct;
end
for i=8:length(Pt)-1
Ai(i) = (Pt(i+1)-Pt(i) ) / (r(i+1)-r(i));
Bi(i) = ((Pt(i)*r(i+1))-(Pt(i+1)*r(i)) ) / (r(i+1)-r(i));
Pt_total(i) = Ai(i)*r(i)+Bi(i);
end
% disp(['Pt_total at ',num2str(V0(j)),' m/s is ',num2str(sum(Pt_total))])
for i=8:length(Pt)-1
M(i) = ((1/3)*Ai(i)*((r(i+1)^3)-(r(i)^3))) + ((0.5)*Bi(i)*((r(i+1)^2)-(r(i)^2)));
end
% % total shaft torque and power
MTotal(j) = B * sum(M);
Power(j)=MTotal(j)*w*1e-3;
disp([' ',num2str(Velocity(j)),' ',...
num2str(MTotal(j)),' ',num2str(Power(j))])
end
VelocityEx =[5,6,7,8,9,11,13,15,17,19,23,25]; % Wind Speed m/s
PowerEx=[2.34,4.03,5.87,7.68,9.62,11.12,9.19,8.93,6.75,6.62,7.62,9.05];
if (flagFirst == 1)
f1 = figure(1);
subplot(2,2,1);
hold on
h1A = plot(vupper,torqueupper,'--r','DisplayName','Experimental Upper Range');
h1B = plot(vlower,torquelower,'--r','DisplayName','Experimental Lower Range');
grid on
h1C = plot(Velocity,MTotal,'--b','LineWidth',1.5,'DisplayName','Baseline Computatation');
legend('Location','best');
xlabel(' Wind Speed (m/s) ');
ylabel(' Torque (Nm) ');
set(gca,'XTick',(4:1:12)); set(gca,'YTick',(0:300:1800))
ylim([0 1800]); xlim([4 12]);
for i =1
if Aoa(i)>-19
title(' Reynolds number = 500.000' )
else
title(' Reynolds number = 1.000.000' )
end
end
drawnow()
%f2 = figure(2); % Power
subplot(2,2,2);
hold on
h2A = plot(Velocity,Power,'-or','LineWidth',1.5,'DisplayName','BEM Theory');
h2B = plot(VelocityEx,PowerEx,'-oblue','LineWidth',1.5,'DisplayName','Experimental Data');
grid on
set(gca,'XTick',(0:5:30)); set(gca,'YTick',(0:2:15))
ylim([0 15]); xlim([0 30]);
xlabel(' Wind Speed (m/s)'); ylabel(' Power(KW) ');
legend('Location','best');
xlabel(' Wind Speed (m/s) '); ylabel(' Torque (Nm) ');
for i =1
if Aoa(i)>-19
title(' Reynolds number = 500.000' )
else
title(' Reynolds number = 1.000.000' )
end
end
drawnow()
%f3 = figure(3); % Twist vs Span Station in percent
subplot(2,2,3)
h3 = plot(SS,Beta,'-*r','DisplayName','\beta');
set(gca,'XTick',(0:0.2:1)); set(gca,'YTick',(-5:5:25))
ylim([-5 25]); xlim([0 1.2]);
xlabel(' Span Station % '); ylabel(' Twist Angle (degrees) ');
title(' Baseline Twist Angle Distribution ')
grid on
legend('location','best')
drawnow()
%f4 = figure(4); % Chord vs Span Station in percent
subplot(2,2,4)
h4 = plot(SS,Chord,'-*b','DisplayName','Chord');
set(gca,'XTick',(0:0.2:1)); set(gca,'YTick',(0:0.2:1))
ylim([0 1]); xlim([0 1.2]);
xlabel(' Span Station % '); ylabel(' Chord length ');
title(' Baseline Chord Distribution ')
grid on
legend('location','best')
drawnow()
else
h1C.XData = Velocity;
h1C.YData = MTotal;
drawnow()
h2A.XData = Velocity;
h2A.YData = Power;
drawnow()
h2B.XData = VelocityEx;
h2B.YData = PowerEx;
drawnow()
h3.XData = SS;
h3.YData = Beta;
drawnow()
h4.XData = SS;
h4.YData = Chord;
drawnow()
pause(0.1); % delay for viewing
end
flagFirst = 0;
end
mehmet salihi
mehmet salihi 2022-1-1
编辑：mehmet salihi 2022-1-1
As an addition for my comment.
in figure 3, I want 16 lines shown for each configuration. all these 16 lines shown in the same figure.
the legends also should be writtenas ( baseline config1 confg2 config3 ------congfig16)
the same for other figures.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by