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

Aim: Solving the iterative for 2D heat conduction equation of Steady-state & Transient State Analysis, forgiven boundary condition.

Objective:

Steady-state analysis & Transient State Analysis

Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

You will implement the following methods for solving implicit equations.

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Your absolute error criteria is 1e-4

Theory:

Steady-state analysis

A system or a process is in a steady state if the variables which define the behavior of the system or the process are unchanging in time.

heat conduction under steady conditions, the temperature of a body at any point does not change with time. Thus, it can be concluded that this type of heat conduction is a function of position only.
The formula for steady-state heat conduction without heat generation is given by the formula:

We make use of numerical methods making use of boundary elements. The solutions are obtained at discrete points in our
square domain.
We start with the mesh and describe the profile within the domain at discrete intervals,

The spacing between the two discrete points in the x & y direction is dx and dy respectively. Now we try to formulate the
problem in such a way that we can set it up as a computer program,

TRANSIENT STATE HEAT CONDUCTION:

In heat conduction under unsteady state conditions, the temperature of a body at any point varies with time as well as position
in one- dimensional and multidimensional systems. This type of heat conduction is also known as transient heat conduction.
The 2D transient heat conduction equation is given by:

In Method of Characteristics, we learned that when having a time derivative in the governing equation we make use of the time marching solution. This solution is implemented by making use of numerical schemes at each and every discrete points on the
square domain

Time interval,
n is for explicit schemes
n+1 is for implicit schemes
The final equations we end up before developing it into a computer program are as follows:
For the Explicit method,

For the Implicit method,

CODE:

Steady state analysis:

clear all
close all
clc
% Steady state analysis
n = 50; % number of grid points
x = linspace(0,1,n);
dx = x(2)-x(1);
dy = dx;
y = linspace(0,1,n);
T = 298*ones(n,n); % initial temperature condition
T(:,1) = 400; % K ,temperature left boundary condition
T(end,:) = 600; % K ,temperature top boundary condition
T(:,end) = 800; % K ,temperature right boundary condition
T(1,:) = 900; % K ,temperature bottom boundary condition
tol = 1e-4; % absolute error criteria
Told1 = T; Told2 = T; Told3 = T;
T1 = T; T2 = T; T3 = T;
% method: 1.Jacobi Method, 2.Gauss Seidel Method, 3.Successive Over
% Relaxation Method (SOR)
for method = 1:3
% jacobi method
if method == 1

jac_error = 9e9; % initial value of error
jacobi_iter = 0;
tic
while jac_error>tol
jacobi_iter = jacobi_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T1(i,j) = (Told1(i-1,j) + Told1(i+1,j) + Told1(i,j-1) + Told1(i,j+1))/4;
end
end
jac_error = max(max(abs(T1-Told1)));
Told1 = T1;

end
analysis_time1 = toc
figure(1)
contourf(x,y,T1,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Steady State Analysis (Jacobi Method) nTemperature Distribution over a plate nNumber of iterations = %d nAnalysis time = %f seconds',jacobi_iter,analysis_time1))
pause(0.02)
% gauss seidel method
elseif method == 2
gauss_error = 9e9; % initial value of error
gauss_iter = 1;
tic
while gauss_error>tol
gauss_iter = gauss_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T2(i,j) = (T2(i-1,j) + Told2(i+1,j) + T2(i,j-1) + Told2(i,j+1))/4;
end

end
gauss_error = max(max(abs(T2-Told2)));
Told2 = T2;
end
analysis_time2 = toc
figure(2)
contourf(x,y,T2,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Steady State Analysis (Gauss Seidel Method) nTemperature Distribution over a plate Number of iterations = %d nAnalysis time = %f seconds',gauss_iter,analysis_time2))
pause(0.02)

% SOR method
w = 2/(1+sin(pi*dx)); % relaxation parameter
elseif method == 3
sor_error = 9e9; % initial value of error
sor_iter = 0;
tic
while sor_error>tol
sor_iter = sor_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T3(i,j) = (1-w)*T3(i,j) + w*(T3(i-1,j) + Told3(i+1,j) + T3(i,j-1) + Told3(i,j+1))/4;
end
end
sor_error = max(max(abs(Told3-T3)));
Told3 = T3;
end
analysis_time3 = toc
figure(3)
contourf(x,y,T3,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Steady State Analysis (SOR Method) nTemperature Distribution over a plate nNumber of iterations = %d nAnalysis time = %f seconds',sor_iter,analysis_time3))
pause(0.03)
end
end

output:

Transient State Analysis:

clear all
close all
clc
% transient heat conduction
n = 50; % number of grid points
x = linspace(0,1,n);
dx = x(2) - x(1);
dy = dx;
y = linspace(0,1,n);
total_time = 10;
dt = 0.01;
total_time_steps = total_time/dt;
T = 298*ones(n,n); % initial temperature condition
T(:,1) = 400; % K ,temperature left boundary condition
T(end,:) = 600; % K ,temperature top boundary condition
T(:,end) = 800; % K ,temperature right boundary condition
T(1,:) = 900; % K ,temperature bottom boundary condition
% defining the corner temperature for smoother graph
T(1,1) = 500; % k top left boundary condition
T(n,1) = 650; % k bottom left boundary condition
T(1,n) = 700; % k top right boundary condition
T(n,n) = 850; % k bottom right boundary condition
tol = 1e-4;
Alpha = 0.002; % thermal diffusivity
k = Alpha*dt/(dx^2);
CFL = Alpha*dt/(dx^2);
Told1 = T; Told2 = T; Told3 = T; Told4 = T;
T1 = T; T2 = T; T3 = T; T4 = T;
T_time1 = T; T_time2 = T; T_time3 = T; T_time4 = T; % for time loop
% 1.transient explicit, 2.transient implicit jacobi, 3.transient implicit gauss seidel, 4.transient
% implicit SOR
for method = 1:4

% 1.transient explicit
if method ==1
cycle = 0;
tic
for t = 1:total_time_steps
for i = 2:(n-1)
for j = 2:(n-1)
T1(i,j) = Told1(i,j) + k*(Told1(i-1,j) + Told1(i+1,j) + Told1(i,j-1) + Told1(i,j+1)-4*Told1(i,j));
end
end
Told1 = T1;
cycle = cycle + 1;
end
convergence_time1 = toc;
figure(1)
contourf(x,y,T1,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Transient State Analysis (Explicit Method) nTemperature Distribution over a plate nNumber of cycles = %d nTime for convergence = %f seconds',cycle,convergence_time1))
pause(0.02)

% 2.transient implicit jacobi
elseif method == 2
tic
for t = 1:total_time_steps
jacobi_error = 9e9;
jacobi_iter = 0;
while (jacobi_error>tol)
jacobi_iter = jacobi_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T2(i,j) = T_time2(i,j)/(1+4*k) + (k/(1+4*k))*(Told2(i-1,j) + Told2(i+1,j) + Told2(i,j) + Told2(i,j+1));
end
end
jacobi_error = max(max(abs(T2-Told2)));
Told2 = T2;
end
T_time2 = T2;
end
convergence_time2 = toc;
figure(2)
contourf(x,y,T2,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Transient State Analysis (Implicit Jacobi Method) nTemperature Distribution over a plate nNumber of iterations = %d nTime for convergence = %f seconds',jacobi_iter,convergence_time2))
pause(0.02)

% 3.transient implicit gauss seidel
elseif method == 3
tic
for t = 1:total_time_steps
GS_error = 9e9;
GS_iter = 0;
while(GS_error>tol)
GS_iter = GS_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T3(i,j) = T_time3(i,j)/(1+4*k) + (k/(1+4*k))*(T3(i-1,j) + Told3(i+1,j) + T3(i,j-1) + Told3(i,j+1));
end
end
GS_error = max(max(abs(Told3-T3)));
Told3 = T3;
end
T_time3 = T3;
end
convergence_time3 = toc;
figure(3)
contourf(x,y,T3,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Transient State Analysis (Implicit Gauss Seidel Method) nTemperature Distribution over a plate nNumber of iterations = %d nTime for convergence = %f seconds',GS_iter,convergence_time3))
pause(0.02)

% 4.transient implicit SOR
elseif method == 4

w = 2/(1+sin(pi*dx)); % relaxation parameter
tic
for t = 1:total_time_steps
SOR_error = 9e9;
SOR_iter = 0;
while (SOR_error>tol)
SOR_iter = SOR_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T4(i,j) =(1-w)*Told4(i,j) + w*(T_time4(i,j))/(1+4*k) + w*(k/(1+4*k))*(T4(i-1,j) + Told4(i+1,j) + T4(i,j-1) + Told4(i,j+1));
end
end
SOR_error = max(max(abs(Told4-T4)));
Told4 = T4;
end
T_time4 = T4;
end
convergence_time4 = toc;
figure(4)
contourf(x,y,T4,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Transient State Analysis (Implicit SOR Method) nTemperature Distribution over a plate nNumber of iterations = %d nTime for convergence = %f seconds',SOR_iter,convergence_time4))
pause(0.02)
end

end

Output:

CONCLUSION:

we have implement the following methods for solving implicit equations like Jacobi, Gauss-seidel and Successive over-relaxation for the 2D heat conduction equation of steady state & Transient State Analysis.

Implicit method can be used to problems when converged results are needed in less computation time as they are
unconditionally stable. With increases in time step size the convergence rate is increased which leads to steady state with less no of iterations.
The explicit method results in time-accurate solutions as it is governed by stability criteria when not satisfied it leads to unstable solution and requires more iterations to attain a steady-state.