How can I plot on the same four figures on my code which changes while list numbers changing?
2 views (last 30 days)
Show older comments
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/849965/image.jpeg)
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
0 Comments
Accepted Answer
Meg Noah
on 1 Jan 2022
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 Comments
More Answers (0)
See Also
Categories
Find more on Graphics Performance in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!