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

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

Answer 1

[Meg Noah]

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

[mehmet salihi]

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]

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.