Nested for loop for showing cooling of different thicknesses in a plot
2 views (last 30 days)
Show older comments
Kenneth Bisgaard Cristensen
on 17 Mar 2021
Commented: Mathieu NOE
on 17 Mar 2021
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
0 Comments
Accepted Answer
Mathieu NOE
on 17 Mar 2021
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
4 Comments
More Answers (0)
See Also
Categories
Find more on Function Creation 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!