Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
steady state analysis;
code for solving 2D steady state laplace equation is given by:
clear all
close all
clc
nx=10;
ny=10;
tol=1e-5;
error=9e9;
l=1;
x=linspace(0,l,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
t=300*ones(nx,ny);
t(1,1)=650;
t(1,nx)=850;
t(ny,1)=500;
t(ny,nx)=700;
t(1,2:nx-1)=900;
t(ny,2:nx-1)=600;
t(2:ny-1,1)=400;
t(2:ny-1,nx)=800;
k=2*((dx^2)+(dy^2))/((dx^2)*(dy^2));
t_initial=t;
told=t;
w=1.2;
[xx,yy]=meshgrid(x,y);
iter_sol=input('enter the iterative solver number=')
% jacobi_iter
jacobi_iter=1;
gs_iter=1;
sor_iter=1;
if iter_sol==1
tol=1e-4;
error=9e9;
% jacobi_iter
% convergence loop
while (error>tol)
% space loop
for i=2:nx-1
for j=2:ny-1
% t(i,j)=(1/k)((told(i-1,j)+told(i+1,j))/(dx^2)+(told(i,j+1)+told(i,j-1))/(dy^2))
ti=(told(i-1,j)+told(i+1,j))/(dx^2);
tj=(told(i,j+1)+told(i,j-1))/(dy^2);
t(i,j)=(1/k)*(ti+tj);
end
end
error=max(max(abs(told-t)));
% update
told=t;
jacobi_iter=jacobi_iter+1;
[c,d]=contourf(xx,yy,t);
colormap(jet);
clabel(c,d);
colorbar;
title('jacobi iterations=',jacobi_iter)
end
end
% gauss_seidal
if iter_sol==2
tol=1e-4;
error=9e9;
while error>tol
% time loop
% space loop
for i=2:nx-1
for j=2:ny-1
% t(i,j)=(1/k)((told(i-1,j)+told(i+1,j))/(dx^2)+(told(i,j+1)+told(i,j-1))/(dy^2))
ti=(t(i-1,j)+told(i+1,j))/(dx^2);
tj=(t(i,j+1)+t(i,j-1))/(dy^2);
t(i,j)=(1/k)*(ti+tj);
end
end
error=max(max(abs(told-t)));
% update
told=t;
gs_iter=gs_iter+1;
[c,d]=contourf(xx,yy,t);
colormap(jet);
clabel(c,d);
colorbar;
title('gauss iterations=',gs_iter)
end
end
% sor_method
if iter_sol==3
tol=1e-4;
error=9e9;
while error>tol
% space loop
for i=2:nx-1
for j=2:ny-1
% t(i,j)=(1/k)((told(i-1,j)+told(i+1,j))/(dx^2)+(told(i,j+1)+told(i,j-1))/(dy^2))
ti=(t(i-1,j)+told(i+1,j))/(dx^2);
tj=(told(i,j+1)+t(i,j-1))/(dy^2);
t_c(i,j)=(1/k)*(ti+tj);
t(i,j)=(1-w)*told(i,j)+w*(t_c(i,j));
end
end
error=max(max(abs(told-t)));
% update
told=t;
sor_iter=sor_iter+1;
[c,d]=contourf(xx,yy,t);
colormap(jet);
clabel(c,d);
colorbar;
title('sor iterations=',sor_iter)
end
end
we assume the inner node starting temperature as 300k so that solution will converge in a faster manner. the equation employed here is a 2D laplace equation which governs the behavior of object at steady state.
steady state is the condition of the body after attaining thermal equilibrium.
we see that in sor mthod the solution converges at lesser number of iterations compared to jacobi or gauss seidal, as we considered overrelaxation factor greater than 1.
un-steady state analysis;
code for solving 2D unsteady state equation is given by:
clear all
close all
clc
nx=10;
ny=10;
alpha=1.5
tol=1e-4;
error=9e9;
l=1;
x=linspace(0,l,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=0.0017;
nt=140;
t=300*ones(nx,ny);
t(1,1)=650;
t(1,nx)=850;
t(ny,1)=500;
t(ny,nx)=700;
t(1,2:nx-1)=900;
t(ny,2:nx-1)=600;
t(2:ny-1,1)=400;
t(2:ny-1,nx)=800;
t_initial=t;
told=t;
w=1.2;
[xx,yy]=meshgrid(x,y);
iter_sol=input('enter the condition or iterative solver number=')
if iter_sol==1
tic;
t_explicit=explicit_transient(t,alpha,dt,dx,dy,nx,ny,nt,told,xx,yy);
time_exp=toc;
end
if iter_sol==2
tic;
t_imp_jac=imp_jac(t,alpha,dt,dx,dy,nx,ny,nt,told,tol,xx,yy);
time_imp_j=toc;
end
if iter_sol==3
tic;
t_imp_gs=imp_gs(t,alpha,dt,dx,dy,nx,ny,nt,told,tol,xx,yy);
time_imp_gs=toc;
end
if iter_sol==4
tic;
t_imp_sor=imp_sor(t,alpha,dt,dx,dy,nx,ny,nt,told,tol,xx,yy);
time_imp_sor=toc;
end
subfunctions:
function t=explicit_transient(t,alpha,dt,dx,dy,nx,ny,nt,told,xx,yy)
k1=alpha*dt/(dx^2);
k2=alpha*dt/(dy^2);
for k=1:nt
for i=2:nx-1
for j=2:ny-1
t(i,j)=told(i,j)+k1*(told(i-1,j)-2*t(i,j)+told(i+1,j))+k2*(told(i,j-1)-2*t(i,j)+told(i,j+1));
end
end
told=t;
[c,d]=contourf(xx,yy,t);
colormap(jet);
clabel(c,d);
colorbar;
title('transient explicit method');
end
function t=imp_jac(t,alpha,dt,dx,dy,nx,ny,nt,told,tol,xx,yy)
k1=alpha*dt/(dx^2);
k2=alpha*dt/(dy^2);
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
tp=t;
jac_iter=1;
% time loop
for k=1:nt
error=9e9;
while (error>tol)
% space loop
for i=2:nx-1
for j=2:ny-1
t(i,j)=term1*tp(i,j)+term2*(told(i-1,j)+told(i+1,j))+term3*(told(i,j-1)+told(i,j+1));
end
end
error=max(max(abs(told-t)));
told=t;
jac_iter=jac_iter+1;
end
tp=t;
[c,d]=contourf(xx,yy,t);
title(sprintf('jacobi iterations=%d',jac_iter));
colormap(jet);
clabel(c,d);
colorbar;
pause(0.03)
end
function t=imp_gs(t,alpha,dt,dx,dy,nx,ny,nt,told,tol,xx,yy)
k1=alpha*dt/(dx^2);
k2=alpha*dt/(dy^2);
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
tp=t;
gs_iter=1;
% time loop
for k=1:nt
error=9e9;
while (error>tol)
% space loop
for i=2:nx-1
for j=2:ny-1
t(i,j)=term1*tp(i,j)+term2*(t(i-1,j)+told(i+1,j))+term3*(t(i,j-1)+told(i,j+1));
end
end
error=max(max(abs(told-t)));
told=t;
gs_iter=gs_iter+1;
end
tp=t;
[c,d]=contourf(xx,yy,t);
title(sprintf('gauss seidal iterations=%d',gs_iter));
colormap(jet);
clabel(c,d);
colorbar;
end
function t=imp_sor(t,alpha,dt,dx,dy,nx,ny,nt,told,tol,xx,yy)
k1=alpha*dt/(dx^2);
k2=alpha*dt/(dy^2);
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
w=1.2;
tp=t;
sor_iter=1;
% time loop
for k=1:nt
error=9e9;
while (error>tol)
% space loop
for i=2:nx-1
for j=2:ny-1
t(i,j)=(1-w)*told(i,j)+w*(term1*tp(i,j)+term2*(t(i-1,j)+told(i+1,j))+term3*(t(i,j-1)+told(i,j+1)));
end
end
error=max(max(abs(told-t)));
told=t;
sor_iter=sor_iter+1;
end
tp=t;
[c,d]=contourf(xx,yy,t);
title(sprintf('sor iterations=%d',sor_iter));
colormap(jet);
clabel(c,d);
colorbar;
end
with given boundary conditions steady state is obtained with 140 time steps.
we employed iteration methods only in implicit methods
we see that in sor mthod the solution converges at lesser number of iterations compared to jacobi or gauss seidal, as we considered overrelaxation factor greater than 1.