Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

OBJECTIVE:

To tackle Steady and Unsteady state 2D Heat Conduction Equation by utilizing the point iterative methods like Jacobi, Gauss Seidel and Successive over-unwinding for both verifiable and express plans utilizing MATLAB.

 

Problem statement:

1. Computation area is to be a square.
2. Assume nx = ny [No. of points on x-direction is equal to no. of points in y-direction]
3. We assume nx = ny = 10
4. Boundary Conditions are
   Left      = 400k
   Right    = 600k
   Top       = 800k 
   Bottom = 900k
5. Apart from the boundary conditions, all other nodes are given an initial temperature of 300k.
6. Various methods to be employed in solving the equation are given below:

 

1. Various States: Steady-state & Unsteady-state Condition
2. Various Schemes: Jacobi, Gauss-seidel & Successive over-relaxation
3. Various Methods to be used in Unsteady-state: IMPLICIT & EXPLICIT

By utilizing Point iterative strategy, address the 2D Heat Conduction Equation and plot the Temperature Contour for every one of the three iterative techniques in particular Jacobi, Gauss-Siedel and SOR for Steady and Unsteady conditions under IMPLICIT plan likewise moreover tackling the EXPLICIT plan for Unsteady Condition.

 

THEORY:

The intensity condition is a partial differential equation that portrays how the dispersion of some quantity (like heat) develops over the long run in a solid medium, as it suddenly moves from where it is higher towards where it is lower.

For a function  u(x,y,z,t)  of three spatial variables (x,y,z)vthe heat equation is                                   

Where α is a real coefficient called the diffusivity of the medium.

 

Types of Heat Conduction schemes

This transfer of Heat energy may occur under two conditions

• Steady-state condition
• Unsteady/Transient state condition

Unsteady state conditions are a precursor to steady state conditions. No system exists initially under steady state conditions. Some time must pass, after heat transfer is initiated, before the system reaches steady state. During that period of transition the system is under unsteady state conditions.

-------------------------

 

Here,

TL =  Left boundary

TT = Top boundary

TP = Middle/Point boundary

TB = Bottom boundary

TR = Right boundary

 

Governing Equations

In order to solve the 2D Heat Conduction problem, the governing equations for respective conditions with respect to cartesian coordinates should be discretized to enter that in a code which is given under MAIN PROGRAM below.

 

MAIN PROGRAM:

Beneath there is a succession of seven individual projects, beginning with IMPLICIT Scheme of Steady and Unsteady-state conditions for three iterative strategies Jacobi, Gauss-siedel and SOR followed by EXPLICIT Scheme for Unsteady/Transient Conditions.

 

IMPLICIT Scheme - Steady-state Condition Program:

Steady-state Equation discretization:

In steady-state conduction, the amount of heat entering any region of an object is equal to the amount of heat coming out (if this were not so, the temperature would be rising or falling, as thermal energy was tapped or trapped in a region).

 

 

 

JACOBI ITERATION

%Jacobi Iteration

close all
clear all
clc

%Solving the Steady state 2D heat conduction (Jacobi Method)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 1;
tolerance = 1e-4;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1) = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_intial = T;
T_old = T;

%Calculation of 2D steady heat conduction EQUATION by JACOBI method
iterative_calc = 0;

while(error > tolerance)

for i = 2:nx - 1

for j = 2:ny - 1
T(i,j) = 0.25*(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;
iterative_calc=iterative_calc + 1;


end

%Plotting
figure
contourf(flipud(T),'ShowText','ON')
colorbar
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Jacobi Iterations = %d', iterative_calc));

 

 

OUTPUT

 

 

 

 

 GAUSS SEIDEL

%Gauss Seidel 

close all
clear all
clc

%Solving the Steady state 2D heat conduction (Gauss Method)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 1;
tolerance = 1e-4;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1) = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_intial = T;
T_old = T;

%Calculation of 2D steady heat conduction EQUATION by Gauss-Seidel method
iterative_calc = 0;

while(error > tolerance)

for i = 2:nx - 1

for j = 2:ny - 1
T(i,j) = 0.25*(T(i-1,j) + T_old(i+1,j) + T(i,j-1) + T_old(i,j+1));

end

end

error = max(max(abs(T_old - T)));
T_old = T;
iterative_calc=iterative_calc + 1;


end

%Plotting
figure
contourf(flipud(T),'ShowText','ON')
colorbar
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Gauss Iterations = %d', iterative_calc));

 

OUTPUT

 

 

 

SUCCESSIVE OVER-RELAXATION

%Successive over-relaxation

close all
clear all
clc

%Solving the Steady state 2D heat conduction (Successive over-relaxation Method)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 1;
tolerance = 1e-4;
r = 1.2;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1) = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_intial = T;
T_old = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method

iterative_calc = 0;

while(error > tolerance)

for i = 2:nx - 1

for j = 2:ny - 1
T(i,j) = r*(0.25*(T(i-1,j) + T_old(i+1,j) + T(i,j-1) + T_old(i,j+1))) + (1 - r)*T_old(i,j);
end

end

error = max(max(abs(T_old - T)));
T_old = T;
iterative_calc=iterative_calc + 1;

end


%Plotting
figure
contourf(flipud(T),'ShowText','ON')
colorbar
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Successive over-relaxation = %d', iterative_calc));

 

OUTPUT

 

 

 

 Unsteady/Transient-state Condition Program:

Unsteady-state Equation discretization:

Non-steady-state situations appear after an imposed change in temperature at a boundary as well as inside of an object, as a result of a new source or sink of heat suddenly introduced within an object, causing temperatures near the source or sink to change in time.

 

 

 

 

JACOBI ITERATION

%Jacobi Iteration

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Jacobi Method(IMPLICIT SCHEME)

%Input Parameters
nx = 10;
ny = nx;
nt = 1600;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
dt=1e-3

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1) = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_old = T;
T_intial = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

Jacobi_iteration = 1;
for k = 1:nt
error = 9e9;
while(error > tolerance)

for i = 2:nx - 1

for j = 2:ny - 1

Term_1 = 1/(1 + (2*k1) + (2*k2));
Term_2 = k1*Term_1;
Term_3 = k2*Term_1;
H = (T_old(i-1,j)) + (T_old(i+1,j));
V = (T_old(i,j+1)) + (T_old(i,j-1));

T(i,j) = (T_intial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);

end

end

error = max(max(abs(T_old - T)));
T_old = T;
Jacobi_iteration = Jacobi_iteration + 1;

end
T_intial = T;
end
%Plotting
figure
contourf(flipud(T),'ShowText','ON')
colorbar
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Jacobi Iterations = %d', Jacobi_iteration));

 

 

OUTPUT

 

 

 

GAUSS-SIEDEL ITERATION

%Gauss-Siedel Iteration

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Gauss Seidel Method(IMPLICIT SCHEME)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1600;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
dt = 1e-3;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1) = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_old = T;
T_intial = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

Gauss_Seidel_iteration = 1;
for k = 1:nt
error = 9e9;

while(error > tolerance)

for i = 2:nx - 1

for j = 2:ny - 1

Term_1 = 1/(1 + (2*k1) + (2*k2));
Term_2 = k1*Term_1;
Term_3 = k2*Term_1;
H = (T(i-1,j)) + (T_old(i+1,j));
V = (T(i,j+1)) + (T_old(i,j-1));

T(i,j) = (T_intial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);

end

end

error = max(max(abs(T_old - T)));
T_old = T;
Gauss_Seidel_iteration = Gauss_Seidel_iteration + 1;

end
T_intial = T;

end
%Plotting
figure
contourf(flipud(T),'ShowText','ON')
colorbar
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of GaussSeidel_iteration = %d', Gauss_Seidel_iteration));

 

OUTPUT

 

 

SUCCESSIVE-OVER RRELAXATION

%Successive over-relaxation

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Successive over-relaxation Method(IMPLICIT SCHEME)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1600;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
dt = 1e-3;
omega = 1.1;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1) = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_old = T;
T_intial = T;
T_end = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

Successive_over_relaxation_iteration = 1;
for k = 1:nt
error = 9e9;

while(error > tolerance)

for i = 2:nx - 1

for j = 2:ny - 1

Term_1 = 1/(1 + (2*k1) + (2*k2));
Term_2 = k1*Term_1;
Term_3 = k2*Term_1;
H = (T(i-1,j)) + (T_old(i+1,j));
V = (T(i,j+1)) + (T_old(i,j-1));

T(i,j) = (T_intial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);
T_end(i,j) = ((1-omega)*T_old(i,j)) + (T(i,j)*omega);

end

end

error = max(max(abs(T_old - T)));
T_old = T_end;
Successive_over_relaxation_iteration = Successive_over_relaxation_iteration + 1;

end

T_intial = T_end;

end
%Plotting
figure
contourf(flipud(T),'ShowText','ON')
colorbar
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Successive-over relaxation iteration = %d',Successive_over_relaxation_iteration));

 

OUTPUT

 

 

 

 

EXPLICIT Scheme - Unsteady/Transient state Condition Program

Unsteady-state Equation discretization(EXPLICIT):

 

 

 

 EXPLICIT TERM

%Explicit Scheme
close all
clear all
clc

%Solving the Unsteady state 2D heat conduction (EXPLICIT SCHEME)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1600;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
dt = 1e-3;
omega = 1.1;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1) = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_intial = T;
T_old = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
iterative_solver = 1;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

if iterative_solver == 1
unsteady_Explicit = 1;
for k = 1:nt
error = 1;

while(error > tolerance)

for i = 2:nx - 1

for j = 2:ny - 1
Term_1 = T_old(i,j);
Term_2 = k1*(T_old(i+1,j) - 2*T_old(i,j) + (T_old(i-1,j)));
Term_3 = k2*(T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1));
T(i,j) = Term_1 + Term_2 + Term_3;
end

end

error = max(max(abs(T_old - T)));
T_old = T;
unsteady_Explicit = unsteady_Explicit + 1;

end
end
end

%Plotting
figure
contourf(flipud(T),'ShowText','ON')
colorbar
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Unsteady Iterations(EXPLICIT) = %d',unsteady_Explicit ));

 

 OUTPUT

 

 

 CONCLUSION

From the above outputs of various methods & comparing their No. of iterations, we can conclude that

• Jacobi > Gauss-seidel > Successive over-relaxation according to the No. of iterations computed for all methods & Schemes.
  
  
• When comparing Steady state and Unsteady state, it is proven that Steady state is the most successful one as it takes least number of iterations and least time to compute.
  
  
• SOR iteration method is the most effective iteration method as it took least No. of iterations for coverging solutions in both the schemes (EXPLICIT & IMPLICIT) compared to Jacobi & Gauss-siedel.