Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
OBJECTIVE:
To solve Steady and Unsteady state 2D Heat Conduction Equation by using the point iterative techniques such as Jacobi, Gauss Seidel & Successive over-relaxation for both implicit and explicit schemes using MATLAB.
Problem statement:
- Computation area is to be a square.
- Assume nx = ny [No. of points on x-direction is equal to no. of points in y-direction]
- We assume nx = ny = 10
- Boundary Conditions are
Left = 400k
Right = 600k
Top = 800k
Bottom = 900k - Apart from the boundary conditions, all other nodes are given an initial temperature of 300k.
- Various methods to be employed in solving the equation are given below:
- Various States: Steady-state & Unsteady-state Condition
- Various Schemes: Jacobi, Gauss-seidel & Successive over-relaxation
- Various Methods to be used in Unsteady-state: IMPLICIT & EXPLICIT
- Various States: Steady-state & Unsteady-state Condition
By using Point iterative method, solve the 2D Heat Conduction Equation and plot the Temperature Contour for all the three iterative methods namely Jacobi, Gauss-Siedel & SOR for Steady and Unsteady conditions under IMPLICIT scheme also additionally solving the EXPLICIT scheme for Unsteady Condition.
THEORY:
The heat equation is a partial differential equation that describes how the distribution of some quantity (such as heat) evolves over time in a solid medium, as it spontaneously flows from places where it is higher towards places where it is lower.
For the heat equation for Unsteady 2D Heat Conduction
`(delT)/(delt) = alpha*((del^2T)/(delx^2) + (del^2T)/(dely^2))`
For a steady state condition, the heat conduction equation is represented as :
`((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
These equation are partial differential equation, and hence the numerical discretization techniques can be used to solve these equations using the numerical analysis.
The boundary conditions for the problem are as given below :
Implicit form For Unsteady/Transient analysis:
If the of the second order derivative terms are computed at the current time step (n+1), then the equation is of the implicit form.
`(T_p^(n+1) - T_p^n)/(Deltat) = alpha*[(T_L^(n+1) - 2*T_P^(n+1) + T_R^(n+1))/(Deltax^2) + (T_T^(n+1) - 2*T_P^(n+1) + T_B^(n+1))/(Deltay^2)]`
Rearranging,
`T_p^(n+1) = T_p^n + alpha*(Deltat)*[(T_L^(n+1) - 2*T_P^(n+1) + T_R^(n+1))/(Deltax^2) + (T_T^(n+1) - 2*T_P^(n+1) + T_B^(n+1))/(Deltay^2)]`
Lets consider `(alpha*(Deltat)) / (Deltax^2) = k_1` and `(alpha*(Deltat)) / (Deltay^2) = k_2`
Subtituting these values in the equation and rearranging the equation :
`T_p^(n+1) =( T_p^n + k_1*(T_L^(n+1) + T_R^(n+1)) + k_2*(T_T^(n+1) + T_B^(n+1)))/(1+2*k_1+2*k_2)`
For implicit form of the equation the RHS is still unknown, therfore system of equations are converted to matrix form and then solved. In other ways the implicit form can also be solved using iterative solvers. The iterative solvers make guess for the unknown values and then compute the error value, and until the residual error comes below a specified limit, the guess values are modified using a transformation matrix and a constant. The residual error or otherwise called as the tolerance is specified by the user before the start of the iterations.
The various methods used by the iterative solvers are :
- Jacobi method :
Jacobi iterative method is considered as an iterative algorithm which is used for determining the solutions for the system of linear equations. In this method, an approximate value (guess) is used and then the entire equation is solved using this guess value. Until it converges, the process is iterated. With the Jacobi method, the values of `T_i^n` obtained in the `n^(th)` iteration remains unchanged until the entire `(n+1)^(th)` iteration has been calculated.
- Gauss Siedel method :
With the Gauss-Seidel method, we use the new values `T_i^(n+1)`as soon as they are known. For example,
once we have computed `T_1^(n+1)` from the first equation, its value is then used in the second equation to obtain the new `T_2^(n+1)` and so on.
In simpler words the difference between the Gauss–Seidel and Jacobi methods is that the Jacobi method uses the values obtained from the previous step while the Gauss–Seidel method always applies the latest updated values during the iterative procedures.
- Successive over relaxation (SOR) method :
This iterative method uses a value ωω (omega) to guess the value for next iteration. By employing the below mentioned formula the values are computed.
`T_i^(n+1) = T_i^n*(1-omega) + omega*(T_(GS))`
In this `T_(GS)` is the value computed using the gauss-siedel iteration method formula. In this manner the gauss siedel method is extended to the relaxation methods.
Transient Analysis using the implicit scheme
In this case, since the RHS of the equation are unknown, the iterative solvers are used. The three iterative methods indicated above are used in the following MATLAB program.
close all
clear all
clc
% input parameters
nx=10; % number of nodes in the x-direction
ny=10; % number of nodes in the y-direction
dt=0.0001 ; % timestep
alpha = 1.4; % thermal diffusivity
T = 300*ones(ny,nx); % initial guess values for temperature
% initial conditions
for i = 1:nx
T(1,i) = 600;
T(ny,i)=900;
end
for j=1:ny
T(j,1) = 400;
T(j,nx) = 800;
end
T(1,1)= (T(1,2)+ T(2,1))/2;
T(1,nx) = (T(1,nx-1)+T(2,nx))/2;
T(ny,1) = (T(ny-1,1) + T(ny,2))/2;
T(ny,nx) = (T(ny,nx-1)+T(ny-1,nx))/2;
%creating a 2D mesh
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
Told = T;
% selection of the types of solver
iteration = 3;
% tolerance for convergence of temperture
tol = 1e-4;
k1 = alpha*(dt/(dx^2));
k2 = alpha*(dt/(dy^2));
Tprev = T;
iter = 1 ;
%Joacobi method
if iteration ==1
for k=1:1500 % time loop
error = 9e9;
while(error > tol) % convergence loop
for i=2:ny-1
for j=2:nx-1
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
T(i,j) = (term1*Tprev(i,j)) + (term2*(Told(i+1,j)+Told(i-1,j))) + (term3*(Told(i,j-1)+Told(i,j+1)));
end
end
error = max(max(abs(Told-T))); % calculating maximum error value
Told = T; % updating the old temperature values
iter = iter +1;
end
Tprev = T;
end
end
%Gauss-siedel method
if iteration ==2
for k=1:1500 % time loop
error = 9e9;
while(error > tol) % convergence loop
for i=2:ny-1
for j=2:nx-1
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
T(i,j) = (term1*Tprev(i,j)) + (term2*(Told(i+1,j)+T(i-1,j))) + (term3*(T(i,j-1)+Told(i,j+1)));
end
end
error = max(max(abs(Told-T))); % calculating error values
Told = T; % updating the old temperature values
iter = iter +1;
end
Tprev = T;
end
end
%SOR method
if iteration ==3
omega = 1.013; % Value of omega
for k=1:1500 % time loop
error = 9e9;
while(error > tol) % convergence loop
for i=2:ny-1
for j=2:nx-1
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
T(i,j) = ((1-omega)*Told(i,j)) + omega*((term1*Tprev(i,j)) + (term2*(Told(i+1,j)+T(i-1,j))) + (term3*(T(i,j-1)+Told(i,j+1))));
end
end
error = max(max(abs(Told-T))); % calculating error value
Told = T; % updating the old temperature values
iter = iter +1;
end
Tprev = T;
end
end
% plotting the temperature distribution
figure(1)
colormap(jet);
A=contourf(x,y,T);
clabel(A);
set(gca,'YDir','reverse');
xlabel('X-Direction');
ylabel('Y-Direction');
grid on ;
title(' Transient analysis of 2D heat equation using implicit scheme');
B = sprintf('number of timestep = %d , Number of iteration = %d',k,iter);
subtitle(B);In the code, all the three iterative solvers are used. The parameter 'iteration' is used to select solver solver. iteration=1 is the iterative solver which uses the jacobi method. iteration=2 is the iterative solver which uses the gauss siedel method and hte iterative solver=3, uses the successive over relaxation method.
The formulas and the output temperature distribution for various methods is given below :
Jacobi method
Formula :
`T_p^(n+1) = (T(prev)_p^n + k_1*(T(old)_R^(n+1) + T(old)_L^(n+1)) + k_2*(T(old)_T^(n+1) + T(old)_B^(n+1)))/(1+2*k_1 + 2*k_2)`
The T(old) represent the guess values from the previous iteration. These values are updated after every iteration and not within an iteration. The T(prev) represents the temperature values from the previous timestep. It remains constant at a particular time step
Output Plot :
Gauss-siedel method
Formula :
`T_p^(n+1) = (T(prev)_p^n + k_1*(T(old)_R^(n+1) + T_L^(n+1)) + k_2*(T_T^(n+1) + T(old)_B^(n+1)))/(1+2*k_1 + 2*k_2)`
In this case the left point and the top point values are calculated in the iteration and guess value is updated within an iteration, Whereas in jacobi the old guess values are used and hence are not updated within an iteration.
Output Plot :
Successive over relaxation method
Formula :
`T_p^(n+1) = (1-omega)*T(old)_p^n + omega*((T(prev)_p^n + k_1*(T(old)_R^(n+1) + T(old)_L^(n+1)) + k_2*(T(old)_T^(n+1) + T(old)_B^(n+1)))/(1+2*k_1 + 2*k_2))`
The ω is the relaxation factor.
Output Plot :
| Name of The Method used | Number of Iterations |
| Jacobi | 5755 |
| Gauss-siedel | 5226 |
| Successive over relaxation | 4688 |
This data is based on running the iterative solver for 1500 timesteps, with each timestep of 0.0001 seconds. The trend in terms of the number of iteration for convergence is :
Jacobi method > Gauss-siedel method > SOR method.
Steady state analysis of 2D heat conduction problem
For a steady state condition, the heat conduction equation is represented as :
`((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
close all
clear all
clc
% input parameters
nx=10; % number of nodes in the x-direction
ny=10; % number of nodes in the y-direction
dt=1e-4 ; % timestep size
alpha = 1.4; % thermal diffusivity
T = 300*ones(nx,ny);
% initial conditions
for i = 1:nx
T(1,i) = 600;
T(ny,i)=900;
end
for j=1:ny
T(j,1) = 400;
T(j,nx) = 800;
end
T(1,1)= (T(1,2)+ T(2,1))/2;
T(1,nx) = (T(1,nx-1)+T(2,nx))/2;
T(ny,1) = (T(ny-1,1) + T(ny,2))/2;
T(ny,nx) = (T(ny,nx-1)+T(ny-1,nx))/2;
%creating a 2D mesh
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
Told = T;
% selection of the iterative solver method
iteration = 3;
% tolerance for convergence of temperture
tol = 1e-4;
error = 1e4;
k = ((2*(dy^2 + dx^2))/(dy^2*dx^2));
%Jacobi method
if iteration == 1
iter = 1 ;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term_1 = ((Told(i+1,j) + Told(i-1,j))/(k*dx^2));
term_2 = ((Told(i,j-1) + Told(i,j+1))/(k*dy^2));
T(i,j) = term_1 + term_2;
end
end
error = max(max(abs(Told-T)));
Told = T;
iter = iter +1;
end
end
%Gauss-siedel method
if iteration == 2
iter = 1 ;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term_1 = ((Told(i+1,j) + T(i-1,j))/(k*dx^2));
term_2 = ((T(i,j-1) + Told(i,j+1))/(k*dy^2));
T(i,j) = term_1 + term_2;
end
end
error = max(max(abs(Told-T)));
Told = T;
iter = iter +1;
end
end
%succesive over-relaxation method
omega=1.5;
if iteration == 3
iter = 1 ;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term_1 = ((Told(i+1,j) + T(i-1,j))/(k*dx^2));
term_2 = ((T(i,j-1) + Told(i,j+1))/(k*dy^2));
T(i,j) = ((Told(i,j)*(1-omega)) + omega*(term_1 + term_2));
end
end
error = max(max(abs(Told-T)));
Told = T;
iter = iter +1;
end
end
% plotting the temperature distribution
figure(1)
colormap(jet);
A=contourf(x,y,T);
clabel(A);
set(gca,'YDir','reverse');
xlabel('X-Direction');
ylabel('Y-Direction');
grid on ;
title('Steady State Analysis of 2D Heat conduction');
B = sprintf('Iterative method = %d , Number of iterations = %d , Error = %d',iteration,iter,error);
subtitle(B);This code uses three iterative solvers to solve the PDE. The parameter 'iteration' is used to select solver solver. iteration=1 is the iterative solver which uses the jacobi method. iteration=2 is the iterative solver which uses the gauss siedel method and hte iterative solver=3, uses the successive over relaxation method.
Since its a steady state the time loop is not found as in the case of the unsteady analysis code. Since its a steady state equation the formula also changes.
The formulas and the output plot when the three iterative solvers are used are as shown.
Jacobi method
Formula :
`T_p = 1/k*[((T(old)_R + T(old)_L)/(Deltax^2)) + ((T(old)_T + T(old)_B)/(Deltay^2))]`
Where `k = 2*((Deltax^2 + Deltay^2)/(Deltax^2*Deltay^2))`
Output Plot :
Gauss-Siedel method
Formula :
`T_p = 1/k*[((T(old)_R + T_L)/(Deltax^2)) + ((T_T + T(old)_B)/(Deltay^2))]`
Output Plot :
Successive over relaxation (SOR) method
Formula :
`T_p = (1-omega)*T(old)_p + omega*(1/k*[((T(old)_R + T_L)/(Deltax^2)) + ((T_T + T(old)_B)/(Deltay^2))])`
In this case, the relaxation factor(ω) is kept at 1.5. It is at this value that the number of iteration required to attain convergence is minimum.
Output Plot :
The number of iterations required to reach convergence is as follows :
| Name of The Method used | Number of Iterations |
| Jacobi | 208 |
| Gauss-siedel | 112 |
| Successive over relaxation | 29 |
It can be concluded that the number of iteration required is least in SOR method. The maximum number of iteration is seen in jacobi method. So, the trend in terms of the number of iteration for convergence is :
Jacobi method > Gauss-siedel method > SOR method.