SOLVING THE STEADY AND UNSTEADY STATE 2-D CONDUCTION EQUATION USING DIFFERENT SOLVERS.
Simulate 2-D Steady-state and unsteady state conduction equation
2-D Conduction equation
`(partialT)/(partialt) = alpha.[(partial^2T)/(partialx^2) + (partial^2T)/(partialy^2)]`
STEADY STATE EQUATION
`[(partial^2T)/(partialx^2) + (partial^2T)/(partialy^2)] = 0`
Discretization
`T(i,j) = (T(i+1,j) + T(i-1,j))/(k.dx^2) + (T(i,j+1) + T(i,j-1))/(k.dy^2)`
We have employed Jacobi, Gauss and Sor techniques.
Steady_state solver
clear all
close all
clc
%Defining the domain
Lx = 1; % X-Dimension
Ly = 1; % Y-dimension
Nx = 80; % X nodes
Ny = Nx; % Y - nodes
%Length between nodes
dx = Lx/Nx-1;
dy = Ly/Ny-1;
%Creating space domain
x = linspace(0,Lx,Nx);
y = linspace(0,Ly,Ny);
[X,Y] = meshgrid(x,y);
% Defining the feild
t =298*ones(Nx,Ny); % defining the temprature feild
% Applying Boundary conditions
T_left = 400; % temprature(k)
T_right = 800;
T_top = 600;
T_bottom = 900;
% assigining these values to temprature feild
t(1,2:end-1) = T_top;
t(end,2:end-1) = T_bottom;
t(2:end-1,1) = T_left;
t(2:end-1,end) = T_right;
t(1,1) = (T_top + T_left)/2;
t(end,1) = (T_bottom + T_left)/2;
t(1,end) = (T_top + T_right)/2;
t(end,end) = (T_bottom + T_right)/2;
%Error tolerence
tol = 1e-4 ;
% assigning the each solver a value
jacobi = 1;
Gauss = 2;
Sor = 3;
%%solver selection counter
solver = 1;
%Jacobi
if solver == 1
t = steady_jacobi(x,y,t,tol,Nx,Ny,dx,dy);
else
%Gauss
if solver == 2
t = steady_gauss(x,y,t,tol,Nx,Ny,dx,dy);
else
%Sor_gauss
if solver == 3
t = steady_Sor_gauss(x,y,t,tol,Nx,Ny,dx,dy);
else
%Sor_jacobi
if solver == 4
t = steady_Sor_jacobi(x,y,t,tol,Nx,Ny,dx,dy);
end
end
end
end
1) Jacobi method
%Steady_state jacobi
function t = steady_jacobi(x,y,t,tol,Nx,Ny,dx,dy)
[ X,Y] = meshgrid(x,y);
jacobi_iteration = 0;
error =9e9;
t_old = t;
tic
k = 2*(dx^2 + dy^2)/(dx^2*dy^2);
while(error>tol)
for i = 2:Nx-1
for j = 2:Ny-1
value_1 = (t_old(i-1,j) + t_old(i+1,j))/(k*dx^2);
value_2 = (t_old(i,j-1) + t_old(i,j+1))/(k*dx^2);
t(i,j) = value_1 + value_2;
end
end
error = max(max(abs(t_old-t)));
t_old = t;
jacobi_iteration = jacobi_iteration+1;
end
time = toc;
%% Visualization
% Contour Map
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temprature distribution \\n jacobi method : %d iterations\',jacobi_iteration))
xlabel(\'x\')
ylabel(\'y\')
clabel(C,\'FontSize\',10)
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Time to solve
fprintf(\'steady state 2D heat conduction (jacobi method)\');
fprintf(\'Computational Time : %f s/n\',time);
fprintf(\'number of iterations : %f\', jacobi_iteration);
% Mesh-plot
figure (2)
mesh(X,Y,t)
grid on
title(sprintf(\'Computational Time : %f seconds\',time));
xlabel(\'x\')
ylabel(\'y\')
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
% Gradient plot
figure(3)
%Gradient
[Tx, Ty] = gradient(t);
qx = -k.*Tx;
qy = -k.*Ty;
imagesc(x,y,t);
set(gca,\'ydir\',\'norm\')% gca is get current axis
hold on
quiver(X,Y,qx,qy,\'w\')
% streamline visualization
[sx,sy] = meshgrid(0.2:0.2:1,0.2:0.2:1);
h = streamline(X,Y,qx,qy,sx,sy);
title(\'streamline of the temprature gradient\')
colormap(jet)
set(h,\'color\',\'red\')
hold on
end# Results
2-D Contour Plot
Number of Iterations = 10726
Computational Time = 0.978 sec
3-D Mesh plot
Temprature gradient across the domain
Quiver and streamline plot used for better visualization.
2) Gauss-sidel Method
% steady_gauss-sidel
function t = steady_gauss(x,y,t,tol,Nx,Ny,dx,dy)
[ X,Y] = meshgrid(x,y);
gauss_iteration = 0;
error =9e9;
t_old = t;
tic
k = 2*(dx^2 + dy^2)/(dx^2*dy^2);
while(error>tol)
for i = 2:Nx-1
for j = 2:Ny-1
value_1 = (t(i-1,j) + t(i+1,j))/(k*dx^2);
value_2 = (t(i,j-1) + t(i,j+1))/(k*dx^2);
t(i,j) = value_1 + value_2;
end
end
error = max(max(abs(t_old-t)));
t_old = t;
gauss_iteration = gauss_iteration+1;
end
time = toc;
%% Visualization block
% Contour-Map
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temprature distribution \\n Gauss-sidel method : %d iterations\',gauss_iteration))
xlabel(\'x\')
ylabel(\'y\')
clabel(C,\'FontSize\',10)
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Time to solve
fprintf(\'steady state 2D heat conduction (Gauss method)\');
fprintf(\'Computational Time : %f s/n\',time);
fprintf(\'number of iterations : %f\', gauss_iteration);
% Mesh-Plot
figure (2)
mesh(X,Y,t)
grid on
title(sprintf(\'Computational Time : %f seconds\',time));
xlabel(\'x\')
ylabel(\'y\')
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Gradient Plot
figure(3)
%Gradient
[Tx, Ty] = gradient(t);
qx = -k.*Tx;
qy = -k.*Ty;
imagesc(x,y,t);
set(gca,\'ydir\',\'norm\')% gca is get current axis
hold on
quiver(X,Y,qx,qy,\'w\')
% streamline visualization
[sx,sy] = meshgrid(0.2:0.2:1,0.2:0.2:1);
h = streamline(X,Y,qx,qy,sx,sy);
title(\'streamline of the temprature gradient\')
colormap(jet)
set(h,\'color\',\'red\')
hold on
end
#Results
2-D contour plot
Number of Iterations = 5812
Computational Time = 0.464 sec
3-D Mesh plot
Temprature gradient
3) SOR Method
a) Sor_gauss Method
% Steady SOR
function t = steady_Sor_gauss(x,y,t,tol,Nx,Ny,dx,dy)
[ X,Y] = meshgrid(x,y);
Sor_iteration = 0;
error =9e9;
t_old = t;
tic
k = 2*(dx^2 + dy^2)/(dx^2*dy^2);
alpha = 1.1;
A = (1/(k*dx^2));
while(error>tol)
for i = 2:Nx-1
for j = 2:Ny-1
t(i,j) = (1-alpha)*t_old(i,j) + alpha*A*(t(i-1,j) + t(i+1,j) + t(i,j-1) + t(i,j+1));
end
end
error = max(max(abs(t_old-t)));
t_old = t;
Sor_iteration = Sor_iteration+1;
end
time = toc;
%% Visualization
% Contour Map
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temprature distribution \\n SOR method : %d iterations\',Sor_iteration))
xlabel(\'x\')
ylabel(\'y\')
clabel(C,\'FontSize\',10)
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Time to solve
fprintf(\'steady state 2D heat conduction (SOR method)\');
fprintf(\'Computational Time : %f s/n\',time);
fprintf(\'number of iterations : %f\', Sor_iteration);
% Mesh-Plot
figure (2)
mesh(X,Y,t)
grid on
title(sprintf(\'Computational Time : %f seconds\',time));
xlabel(\'x\')
ylabel(\'y\')
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Gradient Plot
figure(3)
%Gradient
[Tx, Ty] = gradient(t);
qx = -k.*Tx;
qy = -k.*Ty;
imagesc(x,y,t);
set(gca,\'ydir\',\'norm\')% gca is get current axis
hold on
quiver(X,Y,qx,qy,\'w\')
% streamline visualization
[sx,sy] = meshgrid(0.2:0.2:1,0.2:0.2:1);
h = streamline(X,Y,qx,qy,sx,sy);
title(\'streamline of the temprature gradient\')
colormap(jet)
set(h,\'color\',\'red\')
hold on
end
#Results
2-D Contour Map
Number of iterations = 4860
Computational time = 0.414 sec
3-D Mesh Plot
Temprature gradient
b) Sor_Jacobi Method
Solution did not converge for over relaxed condition(alpha>1).
Convergence Achieved for under relaxed case (alpha<1). However the number of iterations(and therefore computational time) is more than the conventional Jacobi Method
% Steady SOR
function t = steady_Sor_jacobi(x,y,t,tol,Nx,Ny,dx,dy)
[ X,Y] = meshgrid(x,y);
Sor_iteration = 0;
error =9e9;
t_old = t;
tic
k = 2*(dx^2 + dy^2)/(dx^2*dy^2);
alpha = 0.9;
A = (1/(k*dx^2));
while(error>tol)
for i = 2:Nx-1
for j = 2:Ny-1
t(i,j) = (1-alpha)*t_old(i,j)+...
alpha*A*(t_old(i-1,j)+ ...
t_old(i+1,j) + t_old(i,j-1) + t_old(i,j+1));
end
end
error = max(max(abs(t_old-t)));
t_old = t;
Sor_iteration = Sor_iteration+1;
end
time = toc;
%% Visualization
%Contour Map
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temprature distribution \\n SOR method : %d iterations\',Sor_iteration))
xlabel(\'x\')
ylabel(\'y\')
clabel(C,\'FontSize\',10)
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Time to solve
fprintf(\'steady state 2D heat conduction (SOR method)\');
fprintf(\'Computational Time : %f s/n\',time);
fprintf(\'number of iterations : %f\', Sor_iteration);
%Mesh-Plot
figure (2)
mesh(X,Y,t)
grid on
title(sprintf(\'Computational Time : %f seconds\',time));
xlabel(\'x\')
ylabel(\'y\')
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Gradient Plot
figure(3)
%Gradient
[Tx, Ty] = gradient(t);
qx = -k.*Tx;
qy = -k.*Ty;
imagesc(x,y,t);
set(gca,\'ydir\',\'norm\')% gca is get current axis
hold on
quiver(X,Y,qx,qy,\'w\')
% streamline visualization
[sx,sy] = meshgrid(0.2:0.2:1,0.2:0.2:1);
h = streamline(X,Y,qx,qy,sx,sy);
title(\'streamline of the temprature gradient\')
colormap(jet)
set(h,\'color\',\'red\')
hold on
end
#Results
Number of iterations = 11770
Computational Time = 1.003 sec
2-D Contour
3-D Mesh
Temprature Gradient
CONCLUSION - We can observe that SOR_Gauss Method takes the least number of iterations and is computationally efficient
UNSTEADY STATE EQUATION
`(partialT)/(partialt) = alpha.[(partial^2T)/(partialx^2) + (partial^2T)/(partialy^2)]`
Discretization Of Implicit case
`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)])/(1+2*k_1+2*k_2)`
Discretization for Explicit case
`T_(i,j)^(n+1) = T_(i,j)^n+K_1*[T_(i+1,j)^(n)-2*T_(i,j)^n + T_(i-1,j)^(n+1)] + k_2*[T_(i,j+1)^(n)-2*T_(i,j)^n + T_(i,j-1)^(n)]`
Time step length(dt) is derived from the CFL number.
I have fixed the Simulation timr for each simulation for 10sec
Therefore the number of time steps = Simulation time/Length of each time step
Unsteady State Solver
% main run file for unsteady solvers
clear all
close all
clc
% Define the domain
Lx = 1; % x-dimension
Ly = 1; % y-Dimension
% Number of nodes nodes in x and y
Nx = 61;
Ny = Nx;
x = linspace(0,Lx,Nx);
y = linspace(0,Ly,Ny);
dx = Lx/(Nx-1); %element size
dy = dx;
[X,Y] = meshgrid(x,y);
%Temprature feild
t = 298*ones(Nx,Ny); % setting the temprature feild as 298
T_top = 600;
T_bottom = 900;
T_left = 400;
T_right = 800;
% Temprature at the four corners
t(1,2:end-1) = T_top;
t(end,2:end-1) = T_bottom;
t(2:end-1,1) = T_left;
t(2:end-1,end) = T_right;
t(1,1) = (T_top + T_left)/2;
t(end,1) = (T_bottom + T_left)/2;
t(1,end) = (T_top + T_right)/2;
t(end,end) = (T_bottom + T_right)/2;
%Parameters
alpha = 1; % thermal diffusivity
tol =1e-4; % tolerence
error = 9e9;% error
%Stability criteria or CFL number
% CFL = alpha*dt*z; %where z is defined below
CFL = 0.5 ; % k1 + k2
z = (1/(dx)^2 + 1/(dy)^2);
dt = CFL/(alpha*z) ;% time-step length
%dt=1e-4;
k1 = (alpha)*dt/(dx^2);
k2 = (alpha)*dt/(dy^2);
k = 1/(1+ 2*k1 + 2*k2);
%Assigning the each solver a value
Explicit_transient =1;
implicit_jacobi = 2;
implicit_Gauss = 3;
implicit_Sor = 4;
%Solver state
solver =4;
if solver == 1
t = unsteady_explicit(x,y,t,Nx,dt,k1,k2,Ny);
%gauss-sidel
else
if solver == 2
t = Implicit_jacobi(x,y,t,Nx,dt,k1,k2,Ny,k);
%sor
else
if solver == 3
t = Implicit_gauss(x,y,t,Nx,dt,k1,k2,Ny,k);
%SOR
else
if solver == 4
t = Implicit_sor(x,y,t,Nx,dt,k1,k2,Ny,k);
end
end
end
end
Explicit Solver
% unsteady_explicit_solver
function t = unsteady_explicit(x,y,t,Nx,dt,k1,k2,Ny)
[X,Y] = meshgrid(x,y);
explicit_iteration = 0;
t_old = t;
tic
Time = 10;% simulation time in seconds
Nt = (Time/dt)+dt; % No. of iterations
T_start = 0;
for n = 1:Nt
T_start = T_start + dt;% updation simulation time
for i = 2:Nx-1
for j = 2:Ny-1
A = t_old(i,j);
B = k1*(t_old(i-1,j)-2*t_old(i,j)+t_old(i+1,j));
C = k2*(t_old(i,j-1)-2*t_old(i,j)+t_old(i,j+1));
t(i,j) = A + B +C;
end
end
%error = max(max(abs(t_old-t)));
t_old = t;
explicit_iteration = explicit_iteration+1;
end
time = toc;
%% Visualization
% Temprature contour plot
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temprature distribution \\n Explicit : %d iterations\',explicit_iteration))
xlabel(\'x\')
ylabel(\'y\')
clabel(C,\'FontSize\',10)
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Time to solve
fprintf(\'Steady state 2D heat conduction (Explicit method)\');
fprintf(\'Computational Time : %f s/n\',time);
fprintf(\'number of iterations : %f\', explicit_iteration);
%Mesh plot
figure (2)
mesh(X,Y,t)
grid on
title(sprintf(\'Simulation Time = %f seconds\',T_start))
colormap(jet)
xlabel(\'x-dimension\')
ylabel(\'y-dimension\')
zlabel(\'temprature\')
set(gca,\'Ydir\',\'reverse\')
hold on;
%Gradient Plot
figure(3)
%Gradient
[Tx, Ty] = gradient(t);
qx = -k1.*Tx;
qy = -k2.*Ty;
imagesc(x,y,t);
colormap(jet)
set(gca,\'ydir\',\'reverse\')% gca is get current axis
hold on
quiver(X,Y,qx,qy,\'r\')
% streamline visualization
[sx,sy] = meshgrid(0.2:0.2:1,0.2:0.2:1);
h = streamline(X,Y,qx,qy,sx,sy);
title(sprintf(\'Streamline of the temprature gradient \\n Computational Time: %d time\',time))
set(h,\'color\',\'r\')
hold on
end#RESULTS
Number of iterations = 249640
Simulation Time = 10 sec
Computational time = 17.912 sec
2-D Contour plot
Mesh plot
Temprature gradient
2) Implicit_Jacobi
% unsteady_jacobi_solver
function t = Implicit_jacobi(x,y,t,Nx,dt,k1,k2,Ny,k)
[ X,Y] = meshgrid(x,y);
error = 9e9;
implicit_jacobi_iteration = 0;
t_old = t;
tol = 1e-4; %tolerence
term_1 = k1*k;
term_2 = k2*k;
Time = 10;% simulation time in seconds
Nt = (Time/dt) + dt; % No. of iterations
T_start = 0;
tic
for n = 1:Nt
T_start = T_start + dt; % update time
while(error>tol)
for i = 2:Nx-1
for j = 2:Ny-1
%implicit_jacobi
A = t_old(i,j);
B = (t_old(i-1,j)+t_old(i+1,j));
C = (t_old(i,j-1)+t_old(i,j+1));
t(i,j) = (A*k + B*term_1 + C*term_2);
end
end
error = max(max(abs(t_old-t)));
t_old = t;
implicit_jacobi_iteration = implicit_jacobi_iteration+1;
end
end
time = toc;
%% visualization
% Temprature contour plot
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temprature distribution \\n Implicit Jacobi : %d iterations\',implicit_jacobi_iteration))
xlabel(\'x\')
ylabel(\'y\')
clabel(C,\'FontSize\',10)
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Time to solve
fprintf(\'Steady state 2D heat conduction (Implicit Jacobi method)\');
fprintf(\'Computational Time : %f s/n\',time);
fprintf(\'number of iterations : %f\', implicit_jacobi_iteration);
% Mesh plot
figure (2)
mesh(X,Y,t)
grid on
title(sprintf(\'Simulation Time = %f seconds\',T_start))
colormap(jet)
xlabel(\'x-dimension\')
ylabel(\'y-dimension\')
zlabel(\'temprature\')
set(gca,\'Ydir\',\'reverse\')
hold on;
%Gradient Plot
figure(3)
%Gradient
[Tx, Ty] = gradient(t);
qx = -k1.*Tx;
qy = -k2.*Ty;
imagesc(x,y,t);
colormap(jet)
set(gca,\'ydir\',\'reverse\')% gca is get current axis
hold on
quiver(X,Y,qx,qy,\'r\')
% streamline visualization
[sx,sy] = meshgrid(0.2:0.2:1,0.2:0.2:1);
h = streamline(X,Y,qx,qy,sx,sy);
title(sprintf(\'Streamline of the temprature gradient \\n Computational Time: %d time\',time))
set(h,\'color\',\'r\')
hold on
end#Results
Number of iterations = 19702
simulation Time = 10 sec
Computational Time = 1.633 sec
2-D Contour Plot
3-D Mesh
Temprature Gradient
3) Implicit_Gauss
% unsteady_gauss_solver
function t = Implicit_gauss(x,y,t,Nx,dt,k1,k2,Ny,k)
[ X,Y] = meshgrid(x,y);
error = 9e9;
implicit_gauss_iteration = 0;
tol = 1e-4; %tolerence
term_1 = k1*k;
term_2 = k2*k;
Time = 10;% simulation time in seconds
Nt = (Time/dt) + dt; % No. of iterations
T_start = 0;
t_old = t; % Temprature update
tic
for n = 1:Nt
T_start = T_start + dt; % update time
while(error>tol)
for i = 2:Nx-1
for j = 2:Ny-1
%implicit_jacobi
A = t(i,j);
B = (t(i-1,j)+t_old(i+1,j));
C = (t(i,j-1)+t_old(i,j+1));
t(i,j) = (A*k + B*term_1 + C*term_2);
end
end
error = max(max(abs(t_old-t)));
t_old = t;
implicit_gauss_iteration = implicit_gauss_iteration+1;
end
end
time = toc;
%% visualization
% Temprature contour plot
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temprature distribution \\n Implicit Gauss : %d iterations\',implicit_gauss_iteration))
xlabel(\'x\')
ylabel(\'y\')
clabel(C,\'FontSize\',10)
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Time to solve
fprintf(\'Steady state 2D heat conduction (Implicit Gauss method)\');
fprintf(\'Computational Time : %f s/n\',time);
fprintf(\'number of iterations : %f\', implicit_gauss_iteration);
% Mesh plot
figure (2)
mesh(X,Y,t)
grid on
title(sprintf(\'Simulation Time = %f seconds\',T_start))
colormap(jet)
xlabel(\'x-dimension\')
ylabel(\'y-dimension\')
zlabel(\'temprature\')
set(gca,\'Ydir\',\'reverse\')
hold on;
%Gradient Plot
figure(3)
%Gradient
[Tx, Ty] = gradient(t);
qx = -k1.*Tx;
qy = -k2.*Ty;
imagesc(x,y,t);
colormap(jet)
set(gca,\'ydir\',\'reverse\')% gca is get current axis
hold on
quiver(X,Y,qx,qy,\'r\')
% streamline visualization
[sx,sy] = meshgrid(0.2:0.2:1,0.2:0.2:1);
h = streamline(X,Y,qx,qy,sx,sy);
title(sprintf(\'Streamline of the temprature gradient \\n Computational Time: %d time\',time))
set(h,\'color\',\'r\')
hold on
end#RESULTS
Number of iterations = 15334
Simulation Time = 10 sec
Computational Time = 1.22 sec
2-D Contour Plot
3-D Mesh
Temprature gradient
4) Unsteady_SOR
% unsteady_sor_solver
function t = Implicit_sor(x,y,t,Nx,dt,k1,k2,Ny,k)
[ X,Y] = meshgrid(x,y);
error = 9e9;
implicit_sor_iteration = 0;
tol = 1e-4; %tolerence
W = 1.8; % relaxation factor
term_1 = k1*k;
term_2 = k2*k;
Time = 10;% simulation time in seconds
Nt = (Time/dt) + dt; % No. of iterations
T_start = 0;
t_old = t; % Temprature update
tic
for n = 1:Nt
T_start = T_start + dt; % update time
while(error>tol)
for i = 2:Nx-1
for j = 2:Ny-1
%implicit_sor
A = t(i,j);
B = (t(i-1,j)+t_old(i+1,j));
C = (t(i,j-1)+t_old(i,j+1));
t(i,j) = (A*k + B*term_1 + C*term_2);
t(i,j) = (1-W)*(t_old(i,j)) + W*(t(i,j));
end
end
error = max(max(abs(t_old-t)));
t_old = t;
implicit_sor_iteration = implicit_sor_iteration+1;
end
end
time = toc;
%% visualization
% Temprature contour plot
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temprature distribution \\n Implicit SOR : %d iterations\',implicit_sor_iteration))
xlabel(\'x\')
ylabel(\'y\')
clabel(C,\'FontSize\',10)
set(gca,\'Ydir\',\'reverse\')
colormap(jet)
hold on;
%Time to solve
fprintf(\'Steady state 2D heat conduction (Implicit SOR method)\');
fprintf(\'Computational Time : %f s/n\',time);
fprintf(\'number of iterations : %f\', implicit_sor_iteration);
% Mesh plot
figure (2)
mesh(X,Y,t)
grid on
title(sprintf(\'Simulation Time = %f seconds\',T_start))
colormap(jet)
xlabel(\'x-dimension\')
ylabel(\'y-dimension\')
zlabel(\'temprature\')
set(gca,\'Ydir\',\'reverse\')
hold on;
%Gradient Plot
figure(3)
%Gradient
[Tx, Ty] = gradient(t);
qx = -k1.*Tx;
qy = -k2.*Ty;
imagesc(x,y,t);
colormap(jet)
set(gca,\'ydir\',\'reverse\')% gca is get current axis
hold on
quiver(X,Y,qx,qy,\'r\')
% streamline visualization
[sx,sy] = meshgrid(0.2:0.2:1,0.2:0.2:1);
h = streamline(X,Y,qx,qy,sx,sy);
title(sprintf(\'Streamline of the temprature gradient \\n Computational Time: %d time\',time))
set(h,\'color\',\'r\')
hold on
end#Result
Number of iterations = 6946
Simulation time = 10 sec
Computational Time = 0.6277 sec
2-D Contour Plot
3-D Mesh
Temprature Gradient
CONCLUSION
We can Conclude that SOR Method takes the Least number of iteration(therefore least time ) when solving the unsteady state equation implicitily.