Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Code for steady state
%%%%%%% 2D Heat Equation%%%%%%%%%%%%%
%%%%%%% Steady state %%%%%%%%%%%%%%%%
%%%%%%%% Arghajit Chakraborty %%%%%%%
close all
clear all
clc
%%%%%%
nx=21;
ny=21;
%%%%%% DEFINE BOUNDARY CONDITIONS %%%%%%%%
T=300*ones(nx,ny);
T(1:nx,1)=400;
T(1:nx,end)=800;
T(1,1:end)=600;
T(end,1:end)=900;
%%%%%%%DEFINE 'X' & 'Y'DOMAIN %%%%%%%%%%%
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
tol=1e-4;
Told=T;
error=9e9;
k=2*(dx^2+dy^2)/(dx^2*dy^2);
%%%%%%%%%%%%% STEADY STATE %%%%%%%%%%%%%%%%%%%%%%%%%%
iterative_solver=input('Input "1" for Jacobi method, "2" for Gauss Siedel &"3" for SOR')
%%%%%%%%%%%% JACOBI %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if iterative_solver==1
jacobi_iter=1;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term1=1/k*((Told(i-1,j)+Told(i+1,j))/dx^2);
term2=1/k*((Told(i,j-1)+Told(i,j+1))/dy^2);
T(i,j)=term1+term2;
end
end
error=max(max(abs(T-Told)));
Told=T;
jacobi_iter=jacobi_iter+1;
end
%Plot
plot_iterative(x,y,T,iterative_solver)
title_text=sprintf('Steady state 2D heat conduction using Jacobi iterative solver n No.of iterations=%d',jacobi_iter);
title(title_text)
end
%%%%%%%%%%%GAUSS SIEDEL%%%%%%%%%%%%%%%%%%%%%%%%
if iterative_solver==2
gs_iter=1;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term1=1/k*((T(i-1,j)+Told(i+1,j))/dx^2);
term2=1/k*((T(i,j-1)+Told(i,j+1))/dy^2);
T(i,j)=term1+term2;
end
end
error=max(max(abs(T-Told)));
Told=T;
gs_iter=gs_iter+1;
end
plot_iterative(x,y,T,iterative_solver)
title_text=sprintf('Steady state 2D heat conduction using Gauss Siedel iterative solver n No.of iterations=%d',gs_iter);
title(title_text)
end
%%%%%%%%%%%% SOR %%%%%%%%%%%%%%%%%%%%%%%%%%%
if iterative_solver==3
omega=1.2;
sor_iter=1;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term1=1/k*((T(i-1,j)+Told(i+1,j))/dx^2);
term2=1/k*((T(i,j-1)+Told(i,j+1))/dy^2);
T(i,j)=(1-omega)*Told(i,j)+omega*(term1+term2);
end
end
error=max(max(abs(T-Told)));
Told=T;
sor_iter=sor_iter+1;
end
%Plot
plot_iterative(x,y,T,iterative_solver)
title_text=sprintf('Steady state 2D heat conduction using SOR iterative solver n No.of iterations=%d',sor_iter);
title(title_text)
end
%%%%%Function for plotting%%%%%%%%%%%%%%%%%%
function plot_iterative(x,y,T,iterative_solver)
figure(iterative_solver);clf;
[C,h]=contourf(x,y,T);
colorbar;
colormap(jet);
clabel(C,h);
set(gca,'YDIR','reverse');
xlabel('x-axis')
ylabel('y-axis')
% title_text=sprintf('Steady state 2D heat conduction using Gauss Siedel iterative solver n No.of iterations=%d',gs_iter);
% title(title_text)Plots for steady state:
Code for transient state:
%%% 2D HEAT CONDUCTION EQUATION%%%%%%%%%%%%%%%%
%%%%%Transient state%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Arghajit Chakraborty%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
%input
nx=10; %no. of spatial nodes in x-direction
ny=nx; %no. of spatial nodes in y-direction
T=300*ones(nx,ny); %Initialize temperature array
T(:,1)=400; %Temperature boundary conditions
T(:,end)=800; %Temperature boundary conditions
T(1,:)=600; %Temperature boundary conditions
T(end,:)=900; %Temperature boundary conditions
x=linspace(0,1,nx); %define mesh in x-direction
y=linspace(0,1,ny); %define mesh in y-direction
alpha=1.4; %thermal diffusivity
tol=1e-4; %tolerance value for convergence
error=1; %initialize error vale
Told=T;
dx=x(2)-x(1);
dy=y(2)-y(1);
Tprev=T;
method=input('Input "1" for Explicit method or "2" for Implicit method');
if method==1
dt=0.001
k1=alpha*dt/dx^2;
k2=alpha*dt/dy^2;
for k=1:1400
for i=2:nx-1
for j=2:ny-1
term1=Told(i,j-1)-2*Told(i,j)+Told(i,j+1);
term2=Told(i+1,j)-2*Told(i,j)+Told(i-1,j);
T(i,j)=Told(i,j)+k1*term1+k2*term2;
end
end
if max(max(abs(Told-T)))tol
for i=2:nx-1
for j=2:ny-1
term1=(1+(2*k1)+(2*k2))^-1;
term2=k2*term1;
term3=k1*term1;
h=Told(i-1,j)+Told(i+1,j);
v=Told(i,j-1)+Told(i,j+1);
T(i,j)=(Tprev(i,j)*term1)+(term2*h)+(term3*v);
end
end
error=max(max(abs(Told-T)));
Told=T;
jacobi_iter=jacobi_iter+1;
end
Tprev=T;
end
plot_iterative(x,y,T,iterative_solver)
title_text=sprintf('Transient 2D heat conduction (Implicit Method)n using Jacobi iterative solver n No.of iterations=%d',jacobi_iter);
title(title_text)
end
%%%%%%GS%%%%%%%%%%%%%%
if iterative_solver==2
gs_iter=1
for nt=1:600
error=1;
while error>tol
for i=2:nx-1
for j=2:ny-1
term1=(1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
h=T(i-1,j)+Told(i+1,j);
v=T(i,j-1)+Told(i,j+1);
T(i,j)=(Tprev(i,j)*term1)+(term2*h)+(term3*v);
end
end
error=max(max(abs(Told-T)));
Told=T;
gs_iter=gs_iter+1;
end
Tprev=T;
end
plot_iterative(x,y,T,iterative_solver)
title_text=sprintf('Transient 2D heat conduction (Implicit Method)n using Gauss Siedel iterative solver n No.of iterations=%d',gs_iter);
title(title_text)
end
%%%%%%SOR%%%%%%%%%%%%%%
if iterative_solver==3
sor_iter=1
omega=1.3;
for nt=1:600
error=1;
while error>tol
for i=2:nx-1
for j=2:ny-1
term1=(1+(2*k1)+(2*k2))^-1;
term2=k2*term1;
term3=k1*term1;
h=T(i-1,j)+Told(i+1,j);
v=T(i,j-1)+Told(i,j+1);
T(i,j)=((1-omega)*Told(i,j))+(omega*((Tprev(i,j)*term1)+(term2*h)+(term3*v)));
end
end
error=max(max(abs(Told-T)));
Told=T;
sor_iter=sor_iter+1;
end
Tprev=T;
end
plot_iterative(x,y,T,iterative_solver)
title_text=sprintf('Transient 2D heat conduction (Implicit Method)n using SOR iterative solver n No.of iterations=%d',sor_iter);
title(title_text)
end
end
Plots for transient state:
Explicit:
Implicit:
