Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Generalised 2-D heat conduction equation
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=1/alpha(delT)/(delt)`
Discretised form
`(T_(i+1,j)-2 T_(i,j)+ T_(i-1,j))/(Deltax^2) +(T_(i,j+1) -2 T_(i,j) +T_(i,j-1))/(Deltay^2)=1/alpha(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat)`
Given boundary conditions are:
- Top boundary, `T_T =600K`
- Bottom boundary, `T_B=900K`
- Left boundary, `T_L=400K`
- Right boundary, `T_R=800K`
Case 1: Steady state condition
`(T_(i+1,j)^n-2 T_(i,j)^(n+1)+ T_(i-1,j)^n)/(Deltax^2) +(T_(i,j+1)^n -2 T_(i,j)^(n+1) +T_(i,j-1)^n)/(Deltay^2)=0`
`T_(i,j)^(n+1)(2/(Deltax^2)+2/(Deltay^2))=(T_(i+1,j)^n+T_(i-1,j)^n)/(Deltax^2)+(T_(i,j+1)^n+T_(i,j-1)^n)/(Deltay^2)`
`T_(i,j)^(n+1)=(Deltay^2)/(2*(Deltax^2+Deltay^2))(T_(i+1,j)^n+T_(i-1,j)^n)+(Deltax^2)/(2*(Deltax^2+Deltay^2))(T_(i,j+1)^n+T_(i,j-1)^n)`
Let `K_1=(Deltay^2)/(2*(Deltax^2+Deltay^2)) and K_2=(Deltax^2)/(2*(Deltax^2+Deltay^2))`
`T_(i,j)^(n+1)=K_1(T_(i+1,j)^n+T_(i-1,j)^n)+K_2(T_(i,j+1)^n+T_(i,j-1)^n)` [Explicit Jacobi Method]
`T_(i,j)^(n+1)=K_1(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))+K_2(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))` [Implicit Gauss-Seidel Method]
%% case 1. steady state condition
clear all
close all
clc
nx=100; ny=100; % number of nodes in x and y direction
t=ones(nx,ny);
%% Boundary conditions
t(1,2:nx-1)=600;
t(nx,2:nx-1)=900;
t(2:nx-1,nx)=800;
t(2:nx-1,1)=400;
t(1,1)=500;
t(1,nx)=700;
t(nx,ny)=850;
t(nx,1)=650;
%% Initial Condition
t(2:nx-1,2:ny-1)=300;
%% defining grid
domain_size=1;
x=linspace(0,domain_size,nx); dx=domain_size/(nx-1);
y=linspace(0,domain_size,ny); dy=domain_size/(ny-1);
[xx,yy]=meshgrid(x,y);
%% Calculations
k1=dy^2/(2*(dx^2+dy^2)); k2=dx^2/(2*(dx^2+dy^2));
error=9e9;
tolerance=1e-4;
told=t; told2=t; told3=t;
tinitial=t;
jacobi_iteration=1; gs_iteration=1; sor_iteration=1;
t_jac=t;
%% Jacobi method
iterative_solver1=1;
if iterative_solver1==1
while (error>tolerance)
for i=2:nx-1
for j=2:ny-1
t_jac(i,j)=(k1*(told(i+1,j)+told(i-1,j))+k2*(told(i,j+1)+told(i,j-1)));
end
end
error=max(max(t_jac-told));
told=t_jac;
jacobi_iteration=jacobi_iteration+1;
end
end
%% Gauss seidel Method
iterative_solver2=1;
t_gs=t;
if iterative_solver2==1
error=9e9;
while error>tolerance
for i=2:nx-1
for j=2:ny-1
t_gs(i,j)=(k1*(told2(i+1,j)+t_gs(i-1,j))+k2*(told2(i,j+1)+t_gs(i,j-1)));
end
end
error=max(max(t_gs-told2));
told2=t_gs;
gs_iteration=gs_iteration+1;
end
end
%% SOR Method
iterative_solver3=1;
w=1.2;
t_sor=told3;
if iterative_solver3==1
error=9e9;
while error>tolerance
for i=2:nx-1
for j=2:ny-1
t_sor(i,j)=(1-w)*told3(i,j)+w*(k1*(told3(i+1,j)+t_sor(i-1,j))+k2*(told3(i,j+1)+t_sor(i,j-1)));
end
end
error=max(max(t_sor-told3));
told3=t_sor;
sor_iteration=sor_iteration+1;
end
end
%% Jacobi plot
subplot(3,1,1)
contourf(x,y,t_jac,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse'); % reverse Y direction
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('Jacobi Iterations=%d',jacobi_iteration))
%% Gauss seidel plot
subplot(3,1,2)
contourf(x,y,t_gs,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse'); % reverse Y direction
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('Gauss Seidel Iterations=%d',gs_iteration))
%% SOR plot
subplot(3,1,3)
contourf(x,y,t_sor,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse');
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('Successive over relaxation Iterations=%d',sor_iteration))
Case 2: Unsteady state condition
`(T_(i+1,j)^n-2 T_(i,j)^n+ T_(i-1,j)^n)/(Deltax^2) +(T_(i,j+1)^n -2 T_(i,j)^n +T_(i,j-1)^n)/(Deltay^2)=1/alpha(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat)`
`T_(i,j)^(n+1)=(1-(2alphaDeltat)/(Deltax^2)-(2alphaDeltat)/(Deltay^2))T_(i,j)^n+(alphaDeltat)/(Deltax^2)(T_(i+1,j)^n+T_(i-1,j)^n)+(alphaDeltat)/(Deltay^2)(T_(i,j+1)^n+T_(i,j-1)^n)`
`Let K_1 = (alphaDeltat)/(Deltax^2) and K_2= (alphaDeltat)/(Deltay^2)`
In the case of diffusion, Courant Number. `CFL_x=(alphaDeltat)/(Deltax^2) and CFL_y=(alphaDeltat)/(Deltay^2)`
Diffusion stability condition: `CFLle0.5`
`(alphaDeltat)/(Deltax^2)+(alphaDeltat)/(Deltay^2) le0.5`
`(alphaDeltat)/(Deltax^2)=(alphaDeltat)/(Deltay^2)=0.25` (for equisized grid)
`T_(i,j)^(n+1)=(1-2K_1-2K_2)T_(i,j)^n+K_1(T_(i+1,j)^n+T_(i-1,j)^n)+K_2(T_(i,j+1)^n+T_(i,j-1)^n)` [Explicit Jacobi Method]
%% case 3. Unsteady state condition
clear all
close all
clc
nx=100; ny=100; % number of nodes in x and y direction
t=ones(nx,ny);
%% Boundary conditions
t(1,2:nx-1)=600;
t(nx,2:nx-1)=900;
t(2:nx-1,nx)=800;
t(2:nx-1,1)=400;
t(1,1)=500;
t(1,nx)=700;
t(nx,ny)=850;
t(nx,1)=650;
%% Initial Condition
t(2:nx-1,2:ny-1)=300;
%% defining grid
domain_size=1;
x=linspace(0,domain_size,nx); dx=domain_size/(nx-1);
y=linspace(0,domain_size,ny); dy=domain_size/(ny-1);
[xx,yy]=meshgrid(x,y);
%% Calculations
error=9e9;
tolerance=1e-4;
told=t;
tinitial=t;
iteration=1;
alpha=0.88; dt=0.00001;
nt=8000;
k1=(alpha*dt)/(dx)^2; k2=(alpha*dt)/(dy)^2;
CFL=(alpha*dt)/(dx)^2;
jacobi_iteration=1;
%% Jacobi method
iterative_solver2=1;
if iterative_solver2==1
for k=1:nt
for i=2:nx-1
for j=2:ny-1
t(i,j)=k1*(told(i+1,j)+told(i-1,j))+k2*(told(i,j+1)+told(i,j-1))+(1-(2*k1)-(2*k2))*told(i,j);
end
end
error=max(max(t-told));
told=t;
jacobi_iteration=jacobi_iteration+1;
end
%% plotting
contourf(x,y,t,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse'); % reverse Y direction
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('Unsteady state heat conduction ( Jacobi Iterations=%d)',jacobi_iteration))
end
Case 4: Unsteady state condition (Implicit scheme)
`(T_(i+1,j)^(n+1)-2 T_(i,j)^(n+1)+ T_(i-1,j)^(n+1))/(Deltax^2) +(T_(i,j+1)^(n +1)-2 T_(i,j)^(n+1) +T_(i,j-1)^(n+1))/(Deltay^2)=1/alpha(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat)`
`(1+(2alphaDeltat)/(Deltax^2)+(2alphaDeltat)/(Deltay^2))T_(i,j)^(n+1)=T_(i,j)^n+(alphaDeltat)/(Deltax^2)(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))+(alphaDeltat)/(Deltay^2)(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))`
`(1+2K_1+2K_2)T_(i,j)^(n+1)=T_(i,j)^n+K_1(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))+K_2(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))`
`T_(i,j)^(n+1)=1/((1+2K_1+2K_2))[T_(i,j)^n+K_1(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))+K_2(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))]`
Here`K=1/(1+2K_1+2K_2`
`T_(i,j)^(n+1)=1/K[T_(i,j)^n+K_1(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))+K_2(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))]` [Implicit Gauss-Seidel Method]