Solving the steady and unsteady 2D heat conduction problem using explicit and impicit methods in MATLAB.
The general 2D transient heat conduction equation is given as below -
`(delT)/(delt) = alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))`
where `alpha`is the thermal diffusivity of the material (Assumed to be copper)
and T is the temperature in Kelvin.
The above equation needs to be solved numerically using the Jacobi, Gauss Siedel and Successive Over Relaxation (SOR) iterative techniques.
Case 1 - Steady State 2D Equation
The steady state equation is given as -
`(del^2T)/(delx^2) + (del^2T)/(dely^2) = 0`
The steps followed in solving the problem numerically are described below -
Step 1 - Numerical discretization of the governing equation using Taylor's Series method. Since the steady state equation contains only second order PDEs, the central differencing schems has been utilized.
Upon discretization and rearranging the terms, we get the following algebraic form of the PDE.
`T_(i,j) = k*((T_(i+1,j) + T_(i-1,j))/(Deltax)^2 + (T_(i,j+1) + T_(i,j-1))/(Deltay)^2)`
where
`k = (Deltax*Deltay)^2/(2*((Deltax)^2 + (Deltay)^2)`.
This discretized form of the PDE now has been solved using MATLAB.
Step 2 - Defining the spatial discretization parameters, the domain and other variables to be used in the code.
Step 3 - Creating separate functions for each of the iterative techniques.
Jacobi -
`T_(i,j) = k*((Told_(i+1,j) + Told_(i-1,j))/(Deltax)^2 + (Told_(i,j+1) + Told_(i,j-1))/(Deltay)^2)`
Gauss Siedel
`T_(i,j) = k*((T_(i+1,j) + T_(i-1,j))/(Deltax)^2 + (T_(i,j+1) + T_(i,j-1))/(Deltay)^2)`
SOR(Gauss Siedel based)
`T_(i,j) = (1-alpha)T_(i,j) + alpha(k*((T_(i+1,j) + T_(i-1,j))/(Deltax)^2 + (T_(i,j+1) + T_(i,j-1))/(Deltay)^2))`
SOR(Jacobi based)
`T_(i,j) = (1-alpha)*Told_(i,j) + alpha(k*((Told_(i+1,j) + Told_(i-1,j))/(Deltax)^2 + (Told_(i,j+1) + Told_(i,j-1))/(Deltay)^2))`
Case 2 - Transient 2D Equation
The equation to be solved is given as
`(delT)/(delt) = alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))`
This equation can be solved both explicitly and implicitly.
The steps followed are the same as in the steady state case. The temporal discretization is one using the time marching method.
The explicit and implicit forms using the Jacobi iterative technique has been derived below. The algebraic form using the Gauss Siedel and SOR methods can also be similarly derived.
Explicit form using the Jacobi iteration Technique -
`(delT)/(delt) = alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))`
`(T_(i,j)^(n+1) -T_(i,j)^n) /(Deltat) =alpha*((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)`
For simplicity, we assing the following notation to the frid points.
`T_(i,j) = T_P`
`T_(i-1,j) = T_L`
`T_(i+1,j) = T_R`
`T_(i,j-1) = T_B`
`T_(i,j+1) = T_T`
We need to solve the above equation for the next time step i.e. n+1. Thus, after rearranging the terms, we get
`T_P^(n+1) =Told_P^n*(1-k_1-k_2) + k_1*(Told_L^n+Told_R^n) + k_2*(Told_T^n+Told_B^n)`
his is the final algebraic form of the PDE to be solved.
Implicit form using the Jacobi iteration technique -
The temperature values at the current time step are simulataneously solved at all the node points.
Thus the discretized is similar. However, instead of usnig the Temperature values at the point P from the last time step, we use the current values.
`(T_(i,j)^(n+1) -T_(i,j)^n) /(Deltat) =alpha*((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)`
Upon simplication and rearranging the terms, the equation becomes -
`T_P^(n+1) =(Told_P^n + k_1*(Told_L^n+Told_R^n) + k_2*(Told_T^n+Told_B^n))/(1+2*k_1+2*k_2)`
his is the final algebraic form of the PDE to be solved.
Approach used while developing the solver in MATLAB
To keep the code clean and easy to understand, the steady state and transient solvers have been developed separately. In addition to this, the jacobi, gauss siedel and sor iterative techniqes have bee developed as separate functions which are being called from the main algorithm where the initial, boundary and operating conditions have been defined.
Thus the entire solution consists of the following -
1. 2 Main Algorithms( 1 each for steady state and transient equation)
2. 12 additional functions (4 each for explicit and implicit approaches for solving transient equation and 4 for the steady state case)
All the codes and the results have been attached below.
MATLAB CODES -
Steady State Main Algorithm
%Week5 Mid Term Project1
%Steady 2D Heat Conduction equation
clear all
close all
clc
tic
%% Spatial Discretization paratmers
Lx = 1; %Domain Length along the x axis in meters
Ly = 1; %Domain Length along the y axis in meters
nx = 20; %Nodes along the x axis
ny = 20; %Nodes along the y axis
x = linspace(0,Lx,nx); %Mesh points along the x axis
y = linspace(0,Ly,ny); %Mesh points along the y axis
dx =x(2) - x(1); %Spatial step size along x axis
dy =y(2) - y(1); %Spatial step size along y axis
alpha = 1.11e-4; %Copper thermal diffusivity (m2/s). Wikipedia
%% Defining temp initial matrix
T = 350*ones(ny,nx); %Also serves as initial condition in transient problem
T(:,1) = 400; %Temp in Kelvin at left end
T(:,nx) = 800; %Temp in Kelvin at right end
T(1,:) = 600; %Temp in Kelvin at top end
T(ny,:) = 900; %Temp in Kelvin at bottom end
T_init = T;
%Constants to be used in the equations
k = (dx*dy)^2/(2*(dx^2 + dy^2));
k1 = k/dx^2;
k2 = k/dy^2;
%% Iterative solvers
% 1 - Explicit Jacobi; 2 - Explicit gauss siedel;
% 3 - SOR(Jac); 4 - SOR(GS);
iterative_solver = 4;
tol = 1e-4; % Error Tolerance
err = 1; % Initial Error
if iterative_solver == 1;
%Explicit Jacobi
j_iter = 1;
[T,j_iter] = jacobi_SS(T,nx,ny,k1,k2,err,tol,j_iter);
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Steady State Heat conduction Equation using Explicit Jacobi Method',...
sprintf('Initial internal temp = %0d K',T_init(3,3)),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',j_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
elseif iterative_solver == 2;
%Explicit Gauss Siedel
g_iter = 1;
[T,g_iter] = gauss_siedel_SS(T,nx,ny,k1,k2,err,tol,g_iter);
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Steady State Heat conduction Equation using Explicit Gauss Siedel Method',...
sprintf('Initial internal temp = %0d K',T_init(3,3)),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',g_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
elseif iterative_solver == 3
%Explicit JacobiSOR
sorj_iter = 1;
alpha_sor = 0.8;
[T,sorj_iter] = sor_jacobi_SS(T,nx,ny,k1,k2,err,tol,sorj_iter,alpha_sor);
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Steady State Heat conduction Equation using Explicit SOR(Jacobi) Method',...
sprintf('Initial internal temp = %0d K and Over/Under relaxation factor = %0.2f',T_init(3,3),alpha_sor),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',sorj_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
else
%Explicit GaussSiedelSOR
sorgs_iter = 1;
alpha_sor = 1.5;
[T,sorgs_iter] = sor_gs_SS(T,nx,ny,k1,k2,err,tol,sorgs_iter,alpha_sor);
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Steady State Heat conduction Equation using Explicit SOR(Jacobi) Method',...
sprintf('Initial internal temp = %0d K and Over/Under relaxation factor = %0.2f',T_init(3,3),alpha_sor),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',sorgs_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
endJacobi Steady state solver
function [T,j_iter]= jacobi_SS(T,nx,ny,k1,k2,err,tol,j_iter)
T_old = T; %Creating a copy of T for each iteration
Tjac = T; %Value at next time step
jac_iter = j_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = T_old(i+1,j)+T_old(i-1,j);
term2 = T_old(i,j+1)+T_old(i,j-1);
Tjac(i,j) = k1*term1 + k2*term2;
end
end
err = max(max(abs(Tjac-T_old))); %Single max statement gives max of each column.
T_old = Tjac;
jac_iter = jac_iter+1;
j_iter = jac_iter;
end
T = Tjac;Gauss Siedel Steady State Solver
function [T,g_iter]= gauss_siedel_SS(T,nx,ny,k1,k2,err,tol,g_iter)
T_old = T; %Creating a copy of T for each iteration
Tgs = T;
gs_iter = g_iter; %Iteration counter
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = Tgs(i+1,j)+Tgs(i-1,j); % Updated from current time step
term2 = Tgs(i,j+1)+Tgs(i,j-1); % Updated from current time step
Tgs(i,j) = k1*term1 + k2*term2; %Next iteration value
end
end
err = max(max(abs(Tgs-T_old))); %Single max statement gives max of each column.
T_old = Tgs;
gs_iter = gs_iter+1;
g_iter = gs_iter;
end
T = Tgs;
Steady State SOR (Jacobi) Solver
function [T,sorj_iter]= sor_jacobi_SS(T,nx,ny,k1,k2,err,tol,sorj_iter,alpha_sor)
T_old = T; %Creating a copy of T for each iteration
Tsorjac = T;
sorjac_iter = sorj_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = T_old(i+1,j)+T_old(i-1,j);
term2 = T_old(i,j+1)+T_old(i,j-1);
Tsorjac(i,j) = (1-alpha_sor)*T_old(i,j) + alpha_sor*(k1*term1 + k2*term2);
end
end
err = max(max(abs(Tsorjac-T_old))); %Single max statement gives max of each column.
T_old = Tsorjac;
sorjac_iter = sorjac_iter+1;
sorj_iter = sorjac_iter;
end
T = Tsorjac;Steady State SOR (Gauss Siedel) Solver
function [T,sorgs_iter]= sor_gs_SS(T,nx,ny,k1,k2,err,tol,sorgs_iter,alpha_sor)
T_old = T; %Creating a copy of T for each iteration
Tsorgs = T;
sorgauss_iter = sorgs_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = Tsorgs(i+1,j)+Tsorgs(i-1,j);
term2 = Tsorgs(i,j+1)+Tsorgs(i,j-1);
Tsorgs(i,j) = (1-alpha_sor)*Tsorgs(i,j) + alpha_sor*(k1*term1 + k2*term2);
end
end
err = max(max(abs(Tsorgs-T_old))); %Single max statement gives max of each column.
T_old = Tsorgs;
sorgauss_iter = sorgauss_iter+1;
sorgs_iter = sorgauss_iter;
end
T = Tsorgs;Transient Equation Main Algorithm
%Week5 Mid Term Project1
%Transient 2D Heat Conduction equation
clear all
close all
clc
tic
%% Spatial Discretization paratmers
Lx = 1; %Domain Length along the x axis in meters
Ly = 1; %Domain Length along the y axis in meters
nx = 20; %Nodes along the x axis
ny = 20; %Nodes along the y axis
x = linspace(0,Lx,nx); %Mesh points along the x axis
y = linspace(0,Ly,ny); %Mesh points along the y axis
dx =x(2) - x(1); %Spatial step size along x axis
dy =y(2) - y(1); %Spatial step size along y axis
alpha = 1.11e-4; %Copper thermal diffusivity (m2/s). Wikipedia
%% Temporal discretization parameters
nt = 5000; %Number of time steps
dt = 1e-2; %Time step size
%% Defining temp initial and BCs
T = 350*ones(ny,nx); %Also serves as initial condition
T(:,1) = 400; %Temp in Kelvin at left end
T(:,nx) = 800; %Temp in Kelvin at right end
T(1,:) = 600; %Temp in Kelvin at top end
T(ny,:) = 900; %Temp in Kelvin at bottom end
T_init = T; %Copy of initial temp matrix
%Extra terms defined to be used in discretized equation
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);
%% Iterative solvers
% 1 - Explicit Jacobi; 2 - Explicit gauss siedel;
% 3 - SOR(Jac); 4 - SOR(GS);
% 5 - Implicit Jacobi; 6 - Implicit gauss siedel;
% 7 - SOR(Imp Jac); 8 - SOR(Imp GS);
iterative_solver = 8;
tol = 1e-4; % Error Tolerance
err = 1; % Initial Error
if iterative_solver == 1;
%Explicit Jacobi
j_iter = 1;
for time = 1:nt
[T,j_iter] = jacobi(T,nx,ny,k1,k2,err,tol,j_iter);
end
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Transient Heat conduction Equation using Explicit Jacobi Method',...
sprintf('Time = %0d seconds with step size = %0.4f and initial internal temp = %0d K',nt*dt,dt,T_init(3,3)),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',j_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
elseif iterative_solver == 2;
%Explicit Gauss Siedel
g_iter = 1;
for time = 1:nt
[T,g_iter] = gauss_siedel(T,nx,ny,k1,k2,err,tol,g_iter);
end
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Transient Heat conduction Equation using Explicit Gauss Siedel Method',...
sprintf('Time = %0d seconds with step size = %0.4f and initial internal temp = %0d K',nt*dt,dt,T_init(3,3)),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',g_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
elseif iterative_solver == 3
%Explicit SORJacobi
sorj_iter = 1;
alpha_sor = 1.5;
for time = 1:nt
[T,sorj_iter] = sor_jacobi(T,nx,ny,k1,k2,err,tol,sorj_iter,alpha_sor);
end
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Transient Heat conduction Equation using Explicit SOR(Jacobi) Method',...
sprintf('Time = %0d seconds with step size = %0.4f, initial internal temp = %0d K and Over relaxation factor = %0.2f',nt*dt,dt,T_init(3,3),alpha_sor),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',sorj_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
elseif iterative_solver == 4
%Explicit SORGaussSiedel
sorgs_iter = 1;
alpha_sor = 1.5;
for time = 1:nt
[T,sorgs_iter] = sor_gs(T,nx,ny,k1,k2,err,tol,sorgs_iter,alpha_sor);
end
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Transient Heat conduction Equation using Explicit SOR(Gauss Siedel) Method',...
sprintf('Time = %0d seconds with step size = %0.4f, initial internal temp = %0d K and Over relaxation factor = %0.2f',nt*dt,dt,T_init(3,3),alpha_sor),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',sorgs_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
elseif iterative_solver == 5;
%Implicit Jacobi
jimp_iter = 1;
for time = 1:nt
[T,jimp_iter] = jacobi_imp(T,nx,ny,k1,k2,err,tol,jimp_iter);
end
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Transient Heat conduction Equation using Implicit Jacobi Method',...
sprintf('Time = %0d seconds with step size = %0.4f and initial internal temp = %0d K',nt*dt,dt,T_init(3,3)),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',jimp_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
elseif iterative_solver == 6;
%Implicit Gauss Siedel
gimp_iter = 1;
for time = 1:nt
[T,gimp_iter] = gauss_siedel_imp(T,nx,ny,k1,k2,err,tol,gimp_iter);
end
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Transient Heat conduction Equation using Implicit Gauss Siedel Method',...
sprintf('Time = %0d seconds with step size = %0.4f and initial internal temp = %0d K',nt*dt,dt,T_init(3,3)),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',gimp_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
elseif iterative_solver == 7
%Implicit SORJacobi
sorjimp_iter = 1;
alpha_sor = 1.5;
for time = 1:nt
[T,sorjimp_iter] = sor_jacobi_imp(T,nx,ny,k1,k2,err,tol,sorjimp_iter,alpha_sor);
end
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Transient Heat conduction Equation using Implicit SOR(Jacobi) Method',...
sprintf('Time = %0d seconds with step size = %0.4f, initial internal temp = %0d K and Over/Under relaxation factor = %0.2f',nt*dt,dt,T_init(3,3),alpha_sor),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',sorjimp_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
else
%Implicit SORGaussSiedel
sorgsimp_iter = 1;
alpha_sor = 1.5;
for time = 1:nt
[T,sorgsimp_iter] = sor_gs_imp(T,nx,ny,k1,k2,err,tol,sorgsimp_iter,alpha_sor);
end
contourf(x,y,T,'ShowText','on');
colorbar;
title({'2D Transient Heat conduction Equation using Implicit SOR(Gauss Siedel) Method',...
sprintf('Time = %0d seconds with step size = %0.4f, initial internal temp = %0d K and Over relaxation factor = %0.2f',nt*dt,dt,T_init(3,3),alpha_sor),...
sprintf('Computational Time = %0.3f seconds', toc),...
sprintf('Iterations = %0d',sorgsimp_iter)});
xlabel('Length of domain along x axis(m)');
ylabel('Length of domain along y axis(m)');
end
Functions being called in the transient algorithm
function [T,j_iter]= jacobi(T,nx,ny,k1,k2,err,tol,j_iter)
T_old = T; %Creating a copy of T for each iteration
T_prev = T; %value at previous time step
Tjac = T; %Value at next time step
jac_iter = j_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1-(2*k1)-(2*k2);
term2 = T_prev(i+1,j)+T_prev(i-1,j);
term3 = T_prev(i,j+1)+T_prev(i,j-1);
Tjac(i,j) = (T_prev(i,j)*term1) + (k1*term2) + (k2*term3);
end
end
err = max(max(abs(Tjac-T_old))); %Single max statement gives max of each column.
T_old = Tjac;
jac_iter = jac_iter+1;
j_iter = jac_iter;
end
T = Tjac;
function [T,g_iter]= gauss_siedel(T,nx,ny,k1,k2,err,tol,g_iter)
T_old = T; %Creating a copy of T for each iteration
T_prev = T; %Value at previous time step
Tgs = T;
gs_iter = g_iter; %Iteration counter
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1-(2*k1)-(2*k2);
term2 = T_old(i+1,j)+T_old(i-1,j); % Updated from current time step
term3 = T_old(i,j+1)+T_old(i,j-1); % Updated from current time step
Tgs(i,j) = (T_prev(i,j)*term1) + (k1*term2) + (k2*term3); %Next time step value
end
end
err = max(max(abs(Tgs-T_old))); %Single max statement gives max of each column.
T_old = Tgs;
gs_iter = gs_iter+1;
g_iter = gs_iter;
end
T = Tgs;
function [T,sorj_iter]= sor_jacobi(T,nx,ny,k1,k2,err,tol,sorj_iter,alpha_sor)
T_old = T; %Creating a copy of T for each iteration
T_prev = T; %Value at previous time step
Tsorjac = T; %Value at current time step
sorjac_iter = sorj_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1-(2*k1)-(2*k2);
term2 = T_prev(i+1,j)+T_prev(i-1,j);
term3 = T_prev(i,j+1)+T_prev(i,j-1);
Tsorjac(i,j) = (1-alpha_sor)*T_old(i,j) + alpha_sor*((T_prev(i,j)*term1) + (k1*term2) + (k2*term3));
end
end
err = max(max(abs(Tsorjac-T_old))); %Single max statement gives max of each column.
T_old = Tsorjac;
sorjac_iter = sorjac_iter+1;
sorj_iter = sorjac_iter;
end
T = Tsorjac;
function [T,sorgs_iter]= sor_gs(T,nx,ny,k1,k2,err,tol,sorgs_iter,alpha_sor)
T_old = T; %Creating a copy of T for each iteration
T_prev = T; %Creating a copy of T for each time step
Tsorgs = T;
sorgauss_iter = sorgs_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1-(2*k1)-(2*k2);
term2 = Tsorgs(i+1,j)+Tsorgs(i-1,j);
term3 = Tsorgs(i,j+1)+Tsorgs(i,j-1);
Tsorgs(i,j) = (1-alpha_sor)*T_prev(i,j) + alpha_sor*((T_prev(i,j)*term1) + (k1*term2) + (k2*term3));
end
end
err = max(max(abs(Tsorgs-T_old))); %Single max statement gives max of each column.
T_old = Tsorgs;
sorgauss_iter = sorgauss_iter+1;
sorgs_iter = sorgauss_iter;
end
T = Tsorgs;
function [T,jimp_iter]= jacobi_imp(T,nx,ny,k1,k2,err,tol,jimp_iter)
T_old = T;
T_prev = T; %Value at previous time step
Tjac_imp = T; %Value at current time step
jac_iter = jimp_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1+(2*k1)+(2*k2);
term2 = T_old(i+1,j)+T_old(i-1,j);
term3 = T_old(i,j+1)+T_old(i,j-1);
Tjac_imp(i,j) = (T_prev(i,j) + (k1*term2) + (k2*term3))/term1;
end
end
err = max(max(abs(Tjac_imp-T_old))); %Single max statement gives max of each column.
T_old = Tjac_imp;
jac_iter = jac_iter+1;
jimp_iter = jac_iter;
end
T = Tjac_imp;
function [T,gimp_iter]= gauss_siedel_imp(T,nx,ny,k1,k2,err,tol,gimp_iter)
T_old = T;
T_prev = T; %Value at previous time step
Tgs_imp = T; %Value at current time step
gs_iter = gimp_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1+(2*k1)+(2*k2);
term2 = Tgs_imp(i+1,j)+Tgs_imp(i-1,j);
term3 = Tgs_imp(i,j+1)+Tgs_imp(i,j-1);
Tgs_imp(i,j) = (T_prev(i,j) + (k1*term2) + (k2*term3))/term1;
end
end
err = max(max(abs(Tgs_imp-T_old))); %Single max statement gives max of each column.
T_old = Tgs_imp;
gs_iter = gs_iter+1;
gimp_iter = gs_iter;
end
T = Tgs_imp;
function [T,sorjimp_iter]= sor_jacobi_imp(T,nx,ny,k1,k2,err,tol,sorjimp_iter,alpha_sor)
T_old = T; %Creating a copy of T for each iteration
T_prev = T; %Creating a copy of T for each time step
Tsorjac_imp = T;
sorjacimp_iter = sorjimp_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1+(2*k1)+(2*k2);
term2 = T_old(i+1,j)+T_old(i-1,j);
term3 = T_old(i,j+1)+T_old(i,j-1);
Tsorjac_imp(i,j) = (1-alpha_sor)*T_old(i,j) + alpha_sor*((T_prev(i,j) + (k1*term2) + (k2*term3))/term1);
end
end
err = max(max(abs(Tsorjac_imp-T_old))); %Single max statement gives max of each column.
T_old = Tsorjac_imp;
sorjacimp_iter = sorjacimp_iter+1;
sorjimp_iter = sorjacimp_iter;
end
T = Tsorjac_imp;
function [T,sorgsimp_iter]= sor_gs_imp(T,nx,ny,k1,k2,err,tol,sorgsimp_iter,alpha_sor)
T_old = T; %Creating a copy of T for each iteration
T_prev = T; %Creating a copy of T for each time step
Tsorgs_imp = T;
sorgs_imp_iter = sorgsimp_iter;
while (err>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1+(2*k1)+(2*k2);
term2 = Tsorgs_imp(i+1,j)+Tsorgs_imp(i-1,j);
term3 = Tsorgs_imp(i,j+1)+Tsorgs_imp(i,j-1);
Tsorgs_imp(i,j) = (1-alpha_sor)*T_old(i,j) + alpha_sor*((T_prev(i,j) + (k1*term2) + (k2*term3))/term1);
end
end
err = max(max(abs(Tsorgs_imp-T_old))); %Single max statement gives max of each column.
T_old = Tsorgs_imp;
sorgs_imp_iter = sorgs_imp_iter+1;
sorgsimp_iter = sorgs_imp_iter;
end
T = Tsorgs_imp;Results and Discussion -
The results obtained after the simulation pf the problem have been shared below.
2D Transient
From the above results, it is evident that all the methods of solving the 2D transient equation yield the same results. The only differenced are in the computational time and number of iterations required to achieve convergence.
The same has been mentioned in the below table.
| Iterative Technique | Computational Time (Sec) | Iterations |
| Jacobi Explicit | 1.014 | 10001 |
| Gauss Siedel Explicit | 0.977 | 11932 |
| SOR Explicit Jacobi | 1.015 | 10001 |
| SOR Explicit Gauss Siedel | 0.997 | 11753 |
| Implicit Jacobi | 0.984 | 11925 |
| Implicit Gauss Siedel | 0.968 | 10925 |
| Implicit SOR Jacobi | 1.73 | 58088 |
| Implicit SOR Gauss Siedel | 1.78 | 58380 |
The above table highlights the differences in the different techniques. It can be observed that the explicit methods require fewer iterations to converge. However, the explicit methods might not converge all the time. For this case, all the methods have converged to yield a stable solution.
The implicit SOR methods take longer computaitonal times as multiple equation have to be solved simultaneously.
2D Steady State -
| Iterative Technique | Computational Time(Sec) | Iterations |
| Jacobi | 0.682 | 816 |
| Gauss Siedel | 0.704 | 437 |
| SOR Jacobi | 0.8 | 1000 |
| SOR Gauss Siedel | 0.78 | 155 |
All the iterative techniques also yield converged results. Howver, it needs to be highlighted that the SOR (jacobi) method doesn't converge above the relaxtion factoer of 1. Thus the relaxation factor utilized for this case is 0.8, which leads to conversion in 1000 iterations.