Nested for loop for showing cooling of different thicknesses in a plot

Hi Community,

I can't get my plot to show all the diffrent thicnesses of steel from 5 mm to 25 mm with a 5 mm with intevals:

(for thick=.005:.005:.025; % pipe thickness)

Would like some advice on how to integrate this into my plot.

Thanks in advance.

%--------------------------------------------------------------------------
k=55; %steel thermal conductivity (W/m/K)
dens=8000; %steel density (kg/m^3)
cp=500; % steel specific heat capacity (J/kg/K)
h_n=200;%Liquid Nitrogen convective heat transfer coefficient (W/m^2/K)
t_n=-150; %Nitrogen temperature (C)
t_a=20; % air temperature (C)
h_a=20; % not accurate (W/m^2/K)
t_i=t_a; % initial temperature (C)
t_target=-25; % (C)
%--------------------------------------------------------------------------
len=.5; % full length
for thick=.005:.005:.025; % pipe thickness 
%--------------------------------------------------------------------------
del=.005;
dt=.4;
t_t=7200;
%--------------------------------------------------------------------------
results=[];
for b=.05:.01:.3
    l_LN=b/2;
    l_a=(len-b)/2;
    n1=round(l_LN/del+1);
    n=round((l_LN+l_a)/del)+1;
    m=round(thick/del)+1;
    nt=round(t_t/dt)+1;
    %--------------------------------------------------------------------------
    %boundary nodes
    jb=sort([1,2:n1,n1+1:n-1,n,n+1:n:(m-2)*n+1,(m-1)*n+1,m*n,2*n:n:(m-1)*n,(m-1)*n+2:m*n-1]);
    %--------------------------------------------------------------------------
    alfa=k/dens/cp;%thermal diffusitivity (m^2/s)
    t=alfa*dt/del^2;%dimensionless time
    %--------------------------------------------------------------------------
    %check for stability
    max_time_step=del^2/4/alfa/(1+h_n*del/k);
    if dt>=max_time_step
        disp(['Unstable! time step = ',num2str(dt), ' > max permissible time step =', num2str(max_time_step)])
        return
    end
    t2t_target=t_t;
    u=ones(nt,m*n)*t_i;
    for i=1:nt-1
        %--------------top surface
        %top left
        tinf=t_n;hlk=h_n*del/k;
        u(i+1,1)=(1-4*t-2*t*hlk)*u(i,1)+2*t*(u(i,2)+u(i,n+1)+hlk*tinf);
        
        %LN part except top left
        for j=2:n1
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j+n)+2*hlk*tinf);
        end
        %air part except top right
        tinf=t_a;hlk=h_a*del/k;
        for j=n1+1:n-1
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j+n)+2*hlk*tinf);
        end
        %top right
        u(i+1,n)=(1-4*t-4*t*hlk)*u(i,n)+2*t*(u(i,n-1)+u(i,2*n)+2*hlk*tinf);
        %---------------left surface
        for j=n+1:n:(m-2)*n+1%left surface except top left and bottom left
            u(i+1,j)=(1-4*t)*u(i,j)+t*(u(i,j-n)+2*u(i,j+1)+u(i,j+n));
        end
        j=(m-1)*n+1; % bottom left
        u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+2*t*(u(i,j+1)+u(i,j-n)+hlk*tinf);
        
        %------------------------right surface
        %bottom right
        u(i+1,m*n)=(1-4*t-4*t*hlk)*u(i,m*n)+2*t*(u(i,m*n-1)+u(i,(m-1)*n)+2*hlk*tinf);
        
        for j=2*n:n:(m-1)*n %right surface except top right and bottom right
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-n)+u(i,j+n)+2*u(i,j-1)+2*hlk*tinf);
        end
        %---------------------------bottom surface
        for j=(m-1)*n+2:m*n-1 %bottom surface except bottom left and bottom right
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j-n)+2*hlk*tinf);
        end
        %---------------------------internal
        for j=1:m*n
            if ~ismember(j,jb)
                u(i+1,j)=(1-4*t)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+u(i,j-n)+u(i,j+n));
            end
        end
        %--------------------------------------------------------------------------
        if u(i+1,(m-1)*n+1)<=t_target
            disp(['Time to reach @ T = ', num2str(t_target), ' C = ' num2str((i+1)*dt), 's'])
            u(i+2:nt,:)=[];
            break
        end
    end
results=[results;b,(i+1)*dt]
end
end
figure
plot(results(:,1),results(:,2))
title(['L = ',num2str(len),' m; t = ',num2str(thick),' m' ])
xlabel('b (m)')
ylabel('Time to T_{target}  (s)')
grid minor

Answer 1

[Mathieu NOE]

hello Kenneth

simply moved the figure display before the last end

clc
clearvars
close all

%--------------------------------------------------------------------------
k=55; %steel thermal conductivity (W/m/K)
dens=8000; %steel density (kg/m^3)
cp=500; % steel specific heat capacity (J/kg/K)
h_n=200;%Liquid Nitrogen convective heat transfer coefficient (W/m^2/K)
t_n=-150; %Nitrogen temperature (C)
t_a=20; % air temperature (C)
h_a=20; % not accurate (W/m^2/K)
t_i=t_a; % initial temperature (C)
t_target=-25; % (C)

ck = 0;
%--------------------------------------------------------------------------
len=.5; % full length
for thick=.005:.005:.025; % pipe thickness 
%--------------------------------------------------------------------------
del=.005;
dt=.4;
t_t=7200;
%--------------------------------------------------------------------------
results=[];
for b=.05:.01:.3
    l_LN=b/2;
    l_a=(len-b)/2;
    n1=round(l_LN/del+1);
    n=round((l_LN+l_a)/del)+1;
    m=round(thick/del)+1;
    nt=round(t_t/dt)+1;
    %--------------------------------------------------------------------------
    %boundary nodes
    jb=sort([1,2:n1,n1+1:n-1,n,n+1:n:(m-2)*n+1,(m-1)*n+1,m*n,2*n:n:(m-1)*n,(m-1)*n+2:m*n-1]);
    %--------------------------------------------------------------------------
    alfa=k/dens/cp;%thermal diffusitivity (m^2/s)
    t=alfa*dt/del^2;%dimensionless time
    %--------------------------------------------------------------------------
    %check for stability
    max_time_step=del^2/4/alfa/(1+h_n*del/k);
    if dt>=max_time_step
        disp(['Unstable! time step = ',num2str(dt), ' > max permissible time step =', num2str(max_time_step)])
        return
    end
    t2t_target=t_t;
    u=ones(nt,m*n)*t_i;
    for i=1:nt-1
        %--------------top surface
        %top left
        tinf=t_n;hlk=h_n*del/k;
        u(i+1,1)=(1-4*t-2*t*hlk)*u(i,1)+2*t*(u(i,2)+u(i,n+1)+hlk*tinf);
        
        %LN part except top left
        for j=2:n1
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j+n)+2*hlk*tinf);
        end
        %air part except top right
        tinf=t_a;hlk=h_a*del/k;
        for j=n1+1:n-1
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j+n)+2*hlk*tinf);
        end
        %top right
        u(i+1,n)=(1-4*t-4*t*hlk)*u(i,n)+2*t*(u(i,n-1)+u(i,2*n)+2*hlk*tinf);
        %---------------left surface
        for j=n+1:n:(m-2)*n+1%left surface except top left and bottom left
            u(i+1,j)=(1-4*t)*u(i,j)+t*(u(i,j-n)+2*u(i,j+1)+u(i,j+n));
        end
        j=(m-1)*n+1; % bottom left
        u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+2*t*(u(i,j+1)+u(i,j-n)+hlk*tinf);
        
        %------------------------right surface
        %bottom right
        u(i+1,m*n)=(1-4*t-4*t*hlk)*u(i,m*n)+2*t*(u(i,m*n-1)+u(i,(m-1)*n)+2*hlk*tinf);
        
        for j=2*n:n:(m-1)*n %right surface except top right and bottom right
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-n)+u(i,j+n)+2*u(i,j-1)+2*hlk*tinf);
        end
        %---------------------------bottom surface
        for j=(m-1)*n+2:m*n-1 %bottom surface except bottom left and bottom right
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j-n)+2*hlk*tinf);
        end
        %---------------------------internal
        for j=1:m*n
            if ~ismember(j,jb)
                u(i+1,j)=(1-4*t)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+u(i,j-n)+u(i,j+n));
            end
        end
        %--------------------------------------------------------------------------
        if u(i+1,(m-1)*n+1)<=t_target
            disp(['Time to reach @ T = ', num2str(t_target), ' C = ' num2str((i+1)*dt), 's'])
            u(i+2:nt,:)=[];
            break
        end
    end
results=[results;b,(i+1)*dt];



end

hold on
plot(results(:,1),results(:,2))
title('Time to Reach vs. thickness')
ck = ck+1;
legend_string{ck} = ['L = ',num2str(len),' m; t = ',num2str(thick),' m' ];
legend(legend_string')
xlabel('b (m)')
ylabel('Time to T_{target}  (s)')
grid minor

end
hold off

[Kenneth Bisgaard Cristensen]

Hi Matheiu,

Thanks, totally missed that end. Now it works perfectly.

Just a follow-up quastion, what do I do if I need the loop to be 5, 10, 25, 45 and 65 mm?

[Mathieu NOE]

Kenneth

this is the modified code

thickness is now defined as a vector containing your new data (before the loop)

thick_list = [5, 10, 25, 45, 65]*1e-3; % thickness (in m) 

also I modified the legend : it displays thickness in mm instead of m.

all the best,

%--------------------------------------------------------------------------
k=55; %steel thermal conductivity (W/m/K)
dens=8000; %steel density (kg/m^3)
cp=500; % steel specific heat capacity (J/kg/K)
h_n=200;%Liquid Nitrogen convective heat transfer coefficient (W/m^2/K)
t_n=-150; %Nitrogen temperature (C)
t_a=20; % air temperature (C)
h_a=20; % not accurate (W/m^2/K)
t_i=t_a; % initial temperature (C)
t_target=-25; % (C)

ck = 0;
%--------------------------------------------------------------------------
len=.5; % full length
thick_list = [5, 10, 25, 45, 65]*1e-3; % thickness (in m) 

for ithick= 1:numel(thick_list) % pipe thickness 
    thick = thick_list(ithick);
%--------------------------------------------------------------------------
del=.005;
dt=.4;
t_t=7200;
%--------------------------------------------------------------------------
results=[];
for b=.05:.01:.3
    l_LN=b/2;
    l_a=(len-b)/2;
    n1=round(l_LN/del+1);
    n=round((l_LN+l_a)/del)+1;
    m=round(thick/del)+1;
    nt=round(t_t/dt)+1;
    %--------------------------------------------------------------------------
    %boundary nodes
    jb=sort([1,2:n1,n1+1:n-1,n,n+1:n:(m-2)*n+1,(m-1)*n+1,m*n,2*n:n:(m-1)*n,(m-1)*n+2:m*n-1]);
    %--------------------------------------------------------------------------
    alfa=k/dens/cp;%thermal diffusitivity (m^2/s)
    t=alfa*dt/del^2;%dimensionless time
    %--------------------------------------------------------------------------
    %check for stability
    max_time_step=del^2/4/alfa/(1+h_n*del/k);
    if dt>=max_time_step
        disp(['Unstable! time step = ',num2str(dt), ' > max permissible time step =', num2str(max_time_step)])
        return
    end
    t2t_target=t_t;
    u=ones(nt,m*n)*t_i;
    for i=1:nt-1
        %--------------top surface
        %top left
        tinf=t_n;hlk=h_n*del/k;
        u(i+1,1)=(1-4*t-2*t*hlk)*u(i,1)+2*t*(u(i,2)+u(i,n+1)+hlk*tinf);
        
        %LN part except top left
        for j=2:n1
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j+n)+2*hlk*tinf);
        end
        %air part except top right
        tinf=t_a;hlk=h_a*del/k;
        for j=n1+1:n-1
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j+n)+2*hlk*tinf);
        end
        %top right
        u(i+1,n)=(1-4*t-4*t*hlk)*u(i,n)+2*t*(u(i,n-1)+u(i,2*n)+2*hlk*tinf);
        %---------------left surface
        for j=n+1:n:(m-2)*n+1%left surface except top left and bottom left
            u(i+1,j)=(1-4*t)*u(i,j)+t*(u(i,j-n)+2*u(i,j+1)+u(i,j+n));
        end
        j=(m-1)*n+1; % bottom left
        u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+2*t*(u(i,j+1)+u(i,j-n)+hlk*tinf);
        
        %------------------------right surface
        %bottom right
        u(i+1,m*n)=(1-4*t-4*t*hlk)*u(i,m*n)+2*t*(u(i,m*n-1)+u(i,(m-1)*n)+2*hlk*tinf);
        
        for j=2*n:n:(m-1)*n %right surface except top right and bottom right
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-n)+u(i,j+n)+2*u(i,j-1)+2*hlk*tinf);
        end
        %---------------------------bottom surface
        for j=(m-1)*n+2:m*n-1 %bottom surface except bottom left and bottom right
            u(i+1,j)=(1-4*t-2*t*hlk)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+2*u(i,j-n)+2*hlk*tinf);
        end
        %---------------------------internal
        for j=1:m*n
            if ~ismember(j,jb)
                u(i+1,j)=(1-4*t)*u(i,j)+t*(u(i,j-1)+u(i,j+1)+u(i,j-n)+u(i,j+n));
            end
        end
        %--------------------------------------------------------------------------
        if u(i+1,(m-1)*n+1)<=t_target
            disp(['Time to reach @ T = ', num2str(t_target), ' C = ' num2str((i+1)*dt), 's'])
            u(i+2:nt,:)=[];
            break
        end
    end
results=[results;b,(i+1)*dt];

end

hold on
plot(results(:,1),results(:,2))
title('Time to Reach vs. thickness')
ck = ck+1;
legend_string{ck} = ['L = ',num2str(len),' m; t = ',num2str(thick*1000),' mm' ];
legend(legend_string')
xlabel('b (m)')
ylabel('Time to T_{target}  (s)')
grid minor

end
hold off

[Kenneth Bisgaard Cristensen]

Hi Mathieu,

Thanks, that is great, really appreciate the help.

[Mathieu NOE]

Glad I could help

All the best,