## Solving Steady and Unsteady 2-D Heat Conduction Problem for a Square Plate

Objective:

1.  To use iterative solvers - Jacobi Method, Gauss-Seidel Method and Successive Over Relaxation(SOR) Method, to solve steady and unsteady 2D conduction equation
2. To use Explicit Method for unsteady state and see the effect of time step on the solution and analyse CFL condition
3. To use higher value of time step to check its stability on unsteady implicit solutions solved by the iterative solvers.

The staedy state equation for 2D heat conduction is by:

alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))+f=0

In the above equation alpha is thermal diffusivity, T is temperature and f is the source term.

For the convenience source term is neglected so the final equation is:

alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0

Finite difference equation using FTCS(Forward time Central Scheme) is given by:

(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(/_\x)^2+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(/_\y)^2=0

Unsteady state equation for 2D heat conduction is :

(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))+f

t here is time

Again here f is neglected so equation is :

(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))

Finite difference equation is:

(T_(i,j)^(n+1)-T_(i,j)^n)/(/_\t)=(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(/_\x)^2+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(/_\y)^2

let k1=(alpha/_\t)/(/_\x)^2 and k2=(alpha/_\t)/(/_\y)^2

T_(i,j)^(n+1)-T_(i,j)^n=k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))

For Explicit Approach the above equation is :

T_(i,j)^(n+1)-T_(i,j)^n=[k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))]^n

T_(i,j)^(n+1)=T_(i,j)^n+[k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))]^n

where n is the time index and i,j are the space indices in x-direction and y directions respectively

For Implicit Approach the equation is:

T_(i,j)^(n+1)-T_(i,j)^n=[k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))]^(n+1)

T_(i,j)^(n+1)=T_(i,j)^n+[k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))]^(n+1)

On rearranging it becomes:

T_(i,j)^(n+1)+2k1T_(i,j)^(n+1)+2k2T_(i,j)^(n+1)=T_(i,j)^n+[k1(T_(i+1,j)+T_(i-1,j))+k2(T_(i,j+1)+T_(i,j-1))]^(n+1)

T_(i,j)^(n+1)=(T_(i,j)^n+[k1(T_(i+1,j)+T_(i-1,j))+k2(T_(i,j+1)+T_(i,j-1))]^(n+1))/(1+2k1+2k2)

n is the time index here and i,j are space indices

Main Program for Steady State Solutions

%% To solve Steady 2D Conduction problem
% d^2T/dx^2+d^2T/dy^2=0

clc
close all
clear

nx=10;             % Number of grid points in X-direction
ny=10;             % Number of grid points in Y-direction
T=300*ones(nx,ny); % Room Temp = 300K Taken
% Creating grids along X and Y directions
x=linspace(0,1,nx);
y=linspace(0,1,ny);
% Reversing the direction of y vector
y_rev=flip(y);

% Creating Mesh
[xn,yn]=meshgrid(x,y_rev);

% Assigning BC
T_lft=400;
T_top=600;
T_rgt=800;
T_bot=900;
T(2:end-1,1)=T_lft;
T(2:end-1,end)=T_rgt;
T(1,2:end-1)=T_top;
T(end,2:end-1)=T_bot;
T(1,1)=(T_lft+T_top)/2;
T(end,1)=(T_lft+T_bot)/2;
T(1,end)=(T_rgt+T_top)/2;
T(end,end)=(T_rgt+T_bot)/2;

% Initializing values
tol=1e-4;   % Convergence Criterion
T_old=T;
T_initial=T;

for iterative_solver=1 :3

if iterative_solver==1
err_jac=1;
iter_jac=0;
while(err_jac>tol)

for j=2:ny-1
for i=2:nx-1
T(i,j)=(T_old(i-1,j)+T_old(i+1,j)+T_old(i,j-1)+T_old(i,j+1))/4;   % since dx=dy
end
end
err_jac=max(max(abs(T-T_old)));
iter_jac=iter_jac+1;
figure(1)
[c,d]=contourf(xn,yn,T);
colormap(jet);
colorbar;
txt=sprintf('Jacobi Method \nIteration Number= %d',iter_jac);
title(txt);
xlabel('X-direction along square plate\rightarrow');
ylabel('Y-direction along square plate\rightarrow');
clabel(c,d);
pause(0.003);

T_old=T; % Updating old Temperature values
end
end

if iterative_solver==2
err_gs=1;
iter_gs=0;
while(err_gs>tol)

for j=2:ny-1
for i=2:nx-1

T(i,j)=(T(i-1,j)+T_old(i+1,j)+T(i,j-1)+T_old(i,j+1))/4;
end
end

err_gs=max(max(abs(T-T_old)));
iter_gs=iter_gs+1;
figure(2)
[c,d]=contourf(xn,yn,T);
colormap(jet);
colorbar;
txt=sprintf('Gauss-Seidel Method \nIteration Number= %d',iter_gs);
title(txt);
xlabel('X-direction along square plate\rightarrow');
ylabel('Y-direction along square plate\rightarrow');
clabel(c,d);
pause(0.003);

T_old=T; % Updating old Temperature Values
end
end

if iterative_solver==3
% Scaling Factor Omega for SOR Method
omega=0.5:0.5:2;

for k=1:length(omega)
err_sor=1;
iter_sor=zeros(1,4);
T=T_initial;
T_old=T_initial;
while(err_sor>tol)

for j=2:ny-1
for i=2:nx-1

T_gs=(T(i-1,j)+T_old(i+1,j)+T(i,j-1)+T_old(i,j+1))/4;
T(i,j)=(1-omega(k))*T(i,j)+omega(k)*(T_gs);
end
end

err_sor=max(max(abs(T-T_old)));
iter_sor(k)=iter_sor(k)+1;
figure(k+2)
[c,d]=contourf(xn,yn,T);
colormap(jet);
colorbar;
txt=sprintf('Succesive Over Relaxation(SOR) Method \n Iteration Number= %d\n Omega value=%0.1f',iter_sor(k),omega(k));
title(txt);
xlabel('X-direction along square plate\rightarrow');
ylabel('Y-direction along square plate\rightarrow');
clabel(c,d);
pause(0.0003);

T_old=T; % Updating Old Temperature values

end
end
end

% Initializing T values;
T_old=T_initial;
T=T_initial;
end



In the program number of grid points taken is 10 in both X and Y directions along the square plate domain. For SOR method omega (scaling factor) taken as array of (0.5,1,1.5 and 2).

In steady state conditions we are using iterative solvers to converge the solution such that the error falls below the tolerance level.

1. From Jacobi Method

2. From Gauss-Seidel Method

3. From Successive Over Relaxation(SOR) Method

a) For Omega value equal to 0.5

b) For Omega value equal to 1

c) For Omega value equal to 1.5

d) For Omega value equal to 2

Plots Analysis of Steady State condition

• All the above plots show same results except for the result from SOR Method (Omega Value=2).
• The number of iterations in case of Jacobi method is 207,for Gauss-Seidel it is 111 which is less than the Jacobi method.So we can say that Gauss-Seidel method is more efficient than Jacobi method.
• In case of SOR method it is observed that for Omega value< 1(Under-Relaxation), the solution took more number of iterations to converge than at Omega value =1. Moreover the SOR method with omega value =1 is same as Gauss-Seidel method and so their number of interations are also same i.e 111.
• SOR method further proves to be more efficient than Gauss-Seidel method at Omega value>1(Over-Relaxation). For Omega value=1.5 it gives optimum solution with least number of iterations to reach the final solution. However at Omega value>1.5 the number of iterations again increases and at Omega value=2 the solutions blows up.

Main Program for Unsteady state Solutions

%% To solve Unsteady 2D conduction problem
% dT/dt=alpha(d^2T/dx^2+d^2T/dy^2)

clc
close all
clear
% Initializing Parameters
nx=10;
ny=10;
T=300*ones(nx,ny); % Room Temp = 300K Taken
time=0.3;
alpha=1;
x=linspace(0,1,nx);
y=linspace(0,1,ny);

dx=x(2)-x(1);
dy=y(2)-y(1);
% Time step size 1e-2 and 1e-3
dt=[1e-3 1e-2];
%% CFL calculation

CFL=alpha*(dt/dx^2+dt/dy^2);

% Number of time steps and constants k1 and k2 for dt=1e-3 and 1e-2
nt=1+(time./dt);
k1=dt/dx^2;
k2=dt/dx^2;
% Reversing the direction of y vector
y_rev=flip(y);

% Creating Mesh
[xn,yn]=meshgrid(x,y_rev);

% Assigning BC
T_lft=400;
T_top=600;
T_rgt=800;
T_bot=900;
T(2:end-1,1)=T_lft;
T(2:end-1,end)=T_rgt;
T(1,2:end-1)=T_top;
T(end,2:end-1)=T_bot;
T(1,1)=(T_lft+T_top)/2;
T(end,1)=(T_lft+T_bot)/2;
T(1,end)=(T_rgt+T_top)/2;
T(end,end)=(T_rgt+T_bot)/2;

% Initializing Variables

tol=1e-4; %  Convergence Criteria
T_old=T;
T_initial=T;
T_prev_dt=T;

%% ######## Explicit Direct Method #######

% Loop for each time step 1e-2 and 1e-3
for m=1:length(dt)

iter_explicit=0;
% Time Loop
for t=1:nt(m)
% Space Loop
for j=2:ny-1
for i=2:nx-1
term1=1-2*k1(m)-2*k2(m);
H=k1(m)*(T_old(i-1,j)+T_old(i+1,j));
V=k2(m)*(T_old(i,j-1)+T_old(i,j+1));
T(i,j)=T_old(i,j)*term1+H+V;
end
end
iter_explicit=iter_explicit+1;

figure(m)
[c,d]=contourf(xn,yn,T);
colormap(jet);
colorbar;
txt=sprintf('Explicit Method \nIteration Number= %d\n\\Deltat=%0.4f',iter_explicit,dt(m));
title(txt);
xlabel('X-direction along square plate\rightarrow');
ylabel('Y-direction along square plate\rightarrow');
clabel(c,d);
pause(0.003);
% Updating Temperature values at the end of Time loop
T_old=T;
end
end

%% ############  Implicit Method   ###########

% Specific dt values taken=1e-2 and corresponding nt to check whether solution blows or not
dt_imp= 1e-2;
nt_imp=1+time/dt_imp;

% Corresponding k1=alpha*dt/dx^2  and k2=alpha*dt/dy^2
k1=(alpha*dt_imp)/dx^2;
k2=(alpha*dt_imp)/dy^2;

for iterative_solver=1:3

% Initialising Temp at the beginning
T=T_initial;
T_old=T_initial;
T_prev_dt=T_initial;

if iterative_solver==1

%% Jacobi Method ##########

iter_jac=0;
% Time Loop
for t=1:nt_imp
err_jac=1;
% Convergence Loop
while(err_jac>tol)
% Space Loop
for j=2:ny-1
for i=2:nx-1
term1=1+2*k1+2*k2;
H=k1*(T_old(i-1,j)+T_old(i+1,j));
V=k2*(T_old(i,j-1)+T_old(i,j+1));
T(i,j)=(T_prev_dt(i,j)+H+V)/term1;
end
end
iter_jac=iter_jac+1;
err_jac=max(max(abs(T-T_old)));

figure(3)
[c,d]=contourf(xn,yn,T);
colormap(jet);
colorbar;
txt=sprintf('Jacobi Method \nIteration Number= %d\n\\Deltat=%0.3f',iter_jac,dt_imp);
title(txt);
xlabel('X-direction along square plate\rightarrow');
ylabel('Y-direction along square plate\rightarrow');
clabel(c,d);
pause(0.0003);
% Updating Temperature values at the end of Convergence loop
T_old=T;
end
% Updating previous time step Temperature values
T_prev_dt=T;

end
end

if iterative_solver==2
%% Gauss-Seidel Method ###

iter_gs=0;
% Time Loop
for t=1:nt_imp
err_gs=1;
% Convergence Loop
while(err_gs>tol)

for j=2:ny-1
for i=2:nx-1
term1=1+2*k1+2*k2;
H=k1*(T(i-1,j)+T_old(i+1,j));
V=k2*(T(i,j-1)+T_old(i,j+1));
T(i,j)=(T_prev_dt(i,j)+H+V)/term1;
end
end

iter_gs=iter_gs+1;
err_gs=max(max(abs(T-T_old)));

figure(4)
[c,d]=contourf(xn,yn,T);
colormap(jet);
colorbar;
txt=sprintf('Gauss-Seidel Method \nIteration Number= %d\n\\Deltat=%0.3f',iter_gs,dt_imp);
title(txt);
xlabel('X-direction along square plate\rightarrow');
ylabel('Y-direction along square plate\rightarrow');
clabel(c,d);
pause(0.0003);
% Updating Temperature values at the end of Convergence loop
T_old=T;
end
% Updating previous time step Temperature values
T_prev_dt=T;
end
end

if iterative_solver==3
%% SOR Method ####

% Scaling Factor Omega for SOR Method
omega=0.8:0.1:1.2;

for k=1:length(omega)

% Initializing for each omega value
iter_sor=zeros(1,5);
T=T_initial;
T_old=T_initial;
T_prev_dt=T_initial;
% Time Loop
for t=1:nt_imp

err_sor=1;
% Convergence Loop
while(err_sor>tol)
% Space Loop
for j=2:ny-1
for i=2:nx-1
term1=1+2*k1+2*k2;
H=k1*(T(i-1,j)+T_old(i+1,j));
V=k2*(T(i,j-1)+T_old(i,j+1));
T_gs=(T_prev_dt(i,j)+H+V)/term1;
T(i,j)=(1-omega(k))*T_prev_dt(i,j)+omega(k)*T_gs;
end
end

err_sor=max(max(abs(T-T_old)));

iter_sor(k)=iter_sor(k)+1;
figure(k+4)
[c,d]=contourf(xn,yn,T);
colormap(jet);
colorbar;
txt=sprintf('Successive Over Relaxation(SOR) Method \nIteration Number= %d\n\\Deltat=%0.3f\n Omega value=%0.1f',iter_sor(k),dt_imp,omega(k));

title(txt);
xlabel('X-direction along square plate\rightarrow');
ylabel('Y-direction along square plate\rightarrow');
clabel(c,d);

pause(0.0003);
T_old=T; % Updating Temperature values at the end of Convergence Loop
end
T_prev_dt=T; % Updating previous time step Temperature
end
end
end
end


In the program number of grid points taken is 10 in both X and Y directions along the square plate domain. For Explicit approach time step size used is 1e-2 and 1e-3 to check how the time step affects the solutions.

For Implicit approach time step size used is 1e-2 to check the stability of Implicit method at higher time step and also iterative solvers are used to converge the solutions such that the error falls below the set tolerance criteria.

For Implicit SOR method, Omega (scaling factor) is taken as array of (0.8 to 1.2).

1) From Explicit Method at time-step size=1e-3

2) From Explicit Method at time-step size=1e-2

Plots Analysis for Explicit Scheme

• The Explicit method does not use iterative solvers to converge and reach the solution infact it is just one equation which is solved directly. There is no tolerance criteria to converge to hence compared to Implicit Method it is faster and requires less computational time.
• At time step size =1e-3, the solution converges to the same steady-state solution.
• At time step size =1e-2, the solution blows up. This is due to the CFL condition which is given by- alpha((/_\t)/(/_\x)^2+(/_\t)/(/_\y)^2)<=1/2
• Here /_\x=/_\y. So CFL condition is actually alpha(/_\t)/(/_\x)^2<=1/4
• The CFL value is calculated for corresponding time step size using matlab and they are:

CFL value for 1e-2 is 1.62 which is greater than 0.25 and so the solution blows up but for 1e-3 =0.162 which is less than 0.25 so solution is stable.

3) From Implicit Scheme

a) Jacobi Method

b) Gauss-Seidel Method

c) SOR at omega value=0.8

d) SOR at omega value=0.9

e) SOR at omega value=1

f) SOR at omega value=1.1

g) SOR at omega value=1.2

Plots Analysis for Implicit Schemes

• Compared to steady state or explicit methods, implicit shemes shows more number of iterations because the coupled equations formed in implicit shemes are very large and so it requires more computational time to show converging results.
• All implicit scheme methods worked at /_\t=0.01 for which explicit scheme crashed. This proves that even at higher value of time step size implicit method is unconditionally stable.
• Gauss-Seidel Implicit scheme is more efficient than Jacobi or SOR implicit methods with number of iterations=587
• SOR method under relax condition (Omega value<1) although shows less iterations but the solutions are not same as steady state or explicit condition. The solutions from them are not correct.

SOR method with omega value>1 shows more number of iterations than SOR method with omega value=1. So SOR method does not seem to work well for implicit schemes. Although for some condition it works well when time step size is very small.

### Simulation of Flow through Pipe using Axisymmetry Element in OpenFOAM Md Imran Ansari · 2020-01-16 07:23:42

Aim  :    To simulate flow through a pipe upto the fully developed region Objective :  Take a wedge about the axis of the pipe and use axisymmetry condition to get the desired output Take length of the pipe to be entry length and simulate flow a Read more

### Review of Interpolation Schemes and need for a Flux Limiter Md Imran Ansari · 2020-01-06 09:42:38

In Finite Volume Method(FVM) the values of the variable are stored at the computational nodes (central position) of the control volume. The approximations of surface and volume integrals require values of the variables at the centre of the faces of the control volume. S Read more

### Simulating flow through Backward Facing Step using IcoFoam solver with different Grading Factor for mesh Md Imran Ansari · 2019-12-31 08:50:45

Aim:  To simulate the flow through a backward facing step using icoFoam solver   Objective: To create create the geometry as given and then develop mesh for that Use different scale factor to refine mesh near the walls Solve the flow through the domain Read more

### Quasi-One Dimensional Isentropic Subsonic-Supersonic flow through Convergent Divergent Nozzle using MacCormack s Technique Md Imran Ansari · 2019-12-24 17:10:08

Aim: To solve quasi one-dimensional isentropic flow through convergent-divergent nozzle using Mac-Cormack\'s technique Objective: To use conservation and non-conservation forms of partial differential equations to          construct fin Read more

### Fourth Order Approximation of Second Derivative of a function using Central Right sided and Left sided Difference Shemes Md Imran Ansari · 2019-11-28 11:56:29

Fourth Order Approximation of Second Derivative of a Function    1. Central Difference Scheme In Central difference scheme the number of nodes required to get fourth order approximation for a second derivative of a function is =5 ( A+O-1) where A is order o Read more