Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Solving Steady state and Transient state 2D Heat conduction
Please find the attached report and code.
Code:
close all
clc
% Declaring variables
nx=51;
ny=51;
lx=1;
ly=1;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=lx/(nx-1);
dy=ly/(ny-1);
tol=1e-4;
error=9e9;
[xx,yy]=meshgrid(x,y);
%------------ Type of analysis --------------
state=0; % 1--> Steady state 0--> Transient state
method=3;% 1--> Jacobi 2--> Gauss Seidel 3--> SOR
exim=1; % 0--> Explicit 1--> Implicit
SORfactor=1.5; % Relaxation factor
%-------------------------------------------
% Iteration counter
iter=0;
% Temperature matrix and Boundary conditions
T=ones(nx,ny)*298;
T(:,1)=400; % left wall
T(1,:)=600; % top wall
T(:,end)=800; % right wall
T(end,:)=900; % bottom wall
% Edge averaging
T(1,1) = (600+400)/2;
T(1, end) = (600+800)/2;
T(end, 1) = (400+900)/2;
T(end, end) = (800+900)/2;
% Temperature inial array initialization
T_prev=T;
T_initial=T;
T_view=T;
%--------------------------------------------
% Transient state Heat Conduction
%--------------------------------------------
if state==0
%--------------------------------------------
% Transient State Implicit
%--------------------------------------------
if exim==1
% Variables for Transient calculations
alpha=1.11e-4; % Copper
dt=0.9; % Such that CFLx + CFLy= 0.5 )
% Initializing k1, k2 & term
k1=(alpha*dt)/((dx)^2);
k2=(alpha*dt)/((dy)^2);
term=(1+(2*k1)+(2*k2));
if method==1 % Jacobi Method
for nt=1:1500 % Time loop
error=1e9; % error changes after every time loop
while error>tol % Convergence loop
for i=2:nx-1 % Space loop
for j=2:ny-1
Tp=T_initial(i,j);
Tr=T_prev(i+1,j);
Tl=T_prev(i-1,j);
Tt=T_prev(i,j-1);
Tb=T_prev(i,j+1);
A=(k1/term)*(Tr+Tl);
B=(k2/term)*(Tt+Tb);
T(i,j)=(Tp/term)+(A+B);
end
end
error=max(max(abs(T-T_prev)));
T_prev=T;
end
T_initial=T;
end
[C,h] = contourf(xx,yy,flipud(T));
colorbar;
colormap (jet);
clabel(C,h);
title_text=sprintf("Transient state Implicit Jacobi Method, Time steps:%d seconds",nt);
title(title_text);
end
if method==2 % Gauss Seidal Method
for nt=1:1500 % Time loop
error=1e9;
while error>tol % Convergence loop
for i=2:nx-1 % Space loop
for j=2:ny-1
Tp=T_initial(i,j);
Tr=T_prev(i+1,j);
Tl=T(i-1,j);
Tt=T(i,j-1);
Tb=T_prev(i,j+1);
A=k1/term*(Tr+Tl);
B=k2/term*(Tt+Tb);
T(i,j)=(Tp/term) + A + B;
end
end
error=max(max(abs(T-T_prev)));
T_prev=T;
end
T_initial=T;
end
[C,h] = contourf(xx,yy,flipud(T));
colorbar;
colormap (jet);
clabel(C,h);
title_text=sprintf("Transient state Implicit Gauss Seidal Method, Time steps:%d seconds",nt);
title(title_text);
end
if method==3 % SOR Method
for nt=1:1500 % Time loop
error=1e9;
while error>tol % Convergence loop
for i=2:nx-1 % Space loop
for j=2:ny-1
Tp1=T_initial(i,j);
Tp2=T_prev(i,j);
Tr=T_prev(i+1,j);
Tl=T(i-1,j);
Tt=T(i,j-1);
Tb=T_prev(i,j+1);
A=k1/term*(Tr+Tl);
B=k2/term*(Tt+Tb);
T(i,j)=(1-SORfactor)*Tp2 +...
(SORfactor)*((Tp1/term) + A + B);
end
end
error=max(max(abs(T-T_prev)));
T_prev=T;
end
T_initial=T;
end
[C,h] = contourf(xx,yy,flipud(T));
colorbar;
colormap (jet);
clabel(C,h);
title_text=sprintf("Transient state Implicit SOR Method, Time steps:%d seconds",nt);
title(title_text);
end
end
%--------------------------------------------
% Transient State Explicit
%--------------------------------------------
if exim==0
% Variables for Transient calculations
alpha=1.11e-4; % Copper
dt=0.9; % Such that CFLx + CFLy= 0.5 (Ideally CFL <=0.5 for explicit calculations)
% Initializing k1, k2 & term
k1=(alpha*dt)/((dx)^2);
k2=(alpha*dt)/((dy)^2);
term=(1-2*k1-2*k2);
for nt=1:2000 % Time loop
for i= 2:nx-1 % Space loop
for j= 2:ny-1
Tp=T_prev(i,j);
Tl=T_prev(i-1,j);
Tr=T_prev(i+1,j);
Tt=T_prev(i,j-1);
Tb=T_prev(i,j+1);
T(i,j)=(Tp*term) + (k1*(Tl+Tr)) + (k2*(Tt+Tb));
end
end
T_prev=T;
end
[C,h] = contourf(xx,yy,flipud(T));
colorbar;
colormap (jet);
clabel(C,h);
title_text=sprintf("Transient State Explicit Method, Time steps:%d seconds",nt);
title(title_text);
end
end
%--------------------------------------------
% Steady state
%--------------------------------------------
if state==1
% initializing k1 and k2
k1=dx^2/(2*(dx^2+dy^2));
k2=dy^2/(2*(dx^2+dy^2));
if method==1 % Jacobi
while error>tol
for i=2:nx-1
for j=2:ny-1
Tr=T_prev(i+1,j);
Tl=T_prev(i-1,j);
Tt=T_prev(i,j-1);
Tb=T_prev(i,j+1);
T(i,j)=k1*(Tl+Tr)+k2*(Tt+Tb);
end
end
error=max(max(abs(T-T_prev)));
T_prev=T;
iter=iter+1;
end
[C,h] = contourf(xx,yy,flipud(T));
colorbar;
colormap (jet);
clabel(C,h);
title_text=sprintf("Steady state Jacobi Method, Iterations:%d",iter);
title(title_text);
end
if method==2 % Gauss Seidel
while error>tol
for i=2:nx-1
for j=2:ny-1
Tr=T_prev(i+1,j);
Tl=T(i-1,j);
Tt=T(i,j-1);
Tb=T_prev(i,j+1);
T(i,j)=k1*(Tl+Tr)+k2*(Tt+Tb);
end
end
error=max(max(abs(T-T_prev)));
T_prev=T;
iter=iter+1;
end
[C,h] = contourf(xx,yy,flipud(T));
colorbar;
colormap (jet);
clabel(C,h);
title_text=sprintf("Steady state Gauss Seidel Method, Iterations:%d",iter);
title(title_text);
end
if method==3 % SOR method
while error>tol
for i=2:nx-1
for j=2:ny-1
Tp=T_prev(i,j);
Tr=T_prev(i+1,j);
Tl=T(i-1,j);
Tt=T(i,j-1);
Tb=T_prev(i,j+1);
T(i,j)=(1-SORfactor)*Tp+SORfactor*...
(k1*(Tl+Tr)+k2*(Tt+Tb));
end
end
error=max(max(abs(T-T_prev)));
T_prev=T;
iter=iter+1;
end
[C,h] = contourf(xx,yy,flipud(T));
colorbar;
colormap (jet);
clabel(C,h);
title_text=sprintf("Steady state SOR Method, Iterations:%d",iter);
title(title_text);
end
end
Plots:
Steady State, Jacobi Method
Steady State, Gauss Seidel Method
Steady State, SOR Method
Transient State Explicit
Transient State Implicit, Jacobi Method
Transient State Implicit, Gauss Seidel Method
Transient State Implicit, SOR Method
