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
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 Seidel & 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 a function
u = (x, y, z, t)
of three spatial variables (x, y, z) (see Cartesian coordinate system) and the time variable t, the heat equation is
Where, 
is a real coefficient called the diffusivity of the medium.
Now by considering the assumption that the above equation is solved for a thin plate, dimensions in the Z direction are neglected, thereby the above equation is converted into the following form which is as follows:
The above equation can be calculated by considering two approaches:
1.Steady State equation:
Under these conditions, the rate of change of temperature remains 0. Therefore in this condition time, the derivative term is equated to zero.
Steady-state 2D heat conduction equation:
2.Transient State equation:
Under these conditions, both terms are considered therefore the equation is as follows:
Above two equations can be solved in 2 ways by using 3 iterative methods such as Jacobi, Gauss-Seidel and SOR.
1.Explicit Method
2.Implicit Method
Here in this project Transient state equation is solved by using both ways.
Time Loop Logic:
Boundary Condition Map:
Here,
TL = Left boundary
TT = Top boundary
TP= Middle/Point boundary
TB= Bottom boundary
TR= Right boundary
The assumption for solving the problem:
- Assume that the domain is a unit square.
- Assume nx = ny=10(for steady ) and [Number of points along the x-direction is equal to the number of points along the y direction]
- Time step nt =1400
- Boundary conditions for steady and transient case
- Left:400K
- Top:600K
- Right:800K
- Bottom:900K
- Initial temperature = 300K
- Tolerance=1e-4
- Starting error=9e9
- Thermal diffusivity(Alpha)=1.1 m^2/sec
- Domain length=1m
- CFL(for stability)= Alpha*(dt/(dx^2))
- To avoid lengthy execution of program omega assumed value 1.2
- Transient/Unsteady tolerance dt = 1e-3;
IMPLICIT Scheme Solving the Steady-state equation:
Considering Two dimensional Steady State Equation
here discretising the equation by applying the central differencing scheme.
Here the value of k is,
Solving the IMPLICIT steady state equation by Jacobi Method
MATLAB code:
%Solving the IMPLICIT Steady state 2D heat conduction (Jacobi Method)
close all
clear all
clc
%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 = 9e9;
tolerance = 1e-4;
k=2*((dx^2+dy^2)/(dx^2*dy^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 original values to T
T_intial = T;
T_old = T;
%Calculation of 2D steady heat conduction EQUATION by JACOBI method
iterative_solver = 1;
if iterative_solver == 1
jacobi_iteration = 1;
while(error > tolerance)
for i = 2:nx - 1
for j = 2:ny - 1
term1 = (T_old(i-1,j)+T_old(i+1,j))/(k*(dx^2));
term2 = (T_old(i,j-1)+T_old(i,j+1))/(k*(dy^2));
T(i,j) = (term1+term2);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
jacobi_iteration = jacobi_iteration + 1;
end
end
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03);
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Jacobi Iterations(IMPLICIT) = %d', jacobi_iteration));
Result
Solving the IMPLICIT Steady state 2D heat conduction (Gauss Seidel Method)
MATLAB code
%Solving the IMPLICIT Steady state 2D heat conduction (Gauss Seidel Method)
close all
clear all
clc
%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 = 9e9;
tolerance = 1e-4;
k=2*((dx^2+dy^2)/(dx^2*dy^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 original values to T
T_intial = T;
T_old = T;
%Calculation of 2D steady heat conduction EQUATION by Gauss-Seidel method
iterative_solver = 2;
if iterative_solver == 2
GS_iteration = 1;
while(error > tolerance)
for i = 2:nx - 1
for j = 2:ny - 1
term1 = (T(i-1,j)+T_old(i+1,j))/(k*(dx^2));
term2 = (T(i,j-1)+T_old(i,j+1))/(k*(dy^2));
T(i,j) = (term1+term2);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
GS_iteration = GS_iteration + 1;
end
end
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03);
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Gauss-Seidel Iterations(IMPLICIT) = %d', GS_iteration));
Result
Solving the IMPLICIT Steady state 2D heat conduction SOR(Successive over-relaxation Method)
MATLAB code
%Solving the IMPLICIT Steady state 2D heat conduction SOR(Successive Over Relaxation Method)
close all
clear all
clc
%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 = 9e9;
tolerance = 1e-4;
k=2*((dx^2+dy^2)/(dx^2*dy^2));
omega = 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 original values to T
T_intial = T;
T_old = T;
%Calculation of 2D steady heat conduction EQUATION by Successive Over Relaxation method
iterative_solver = 3;
if iterative_solver == 3
SOR_iteration = 1;
while(error > tolerance)
for i = 2:nx - 1
for j = 2:ny - 1
term1 = (T(i-1,j)+T_old(i+1,j))/(k*(dx^2));
term2 = (T(i,j-1)+T_old(i,j+1))/(k*(dy^2));
T(i,j) = omega*(term1+term2)+ (1 - omega)*T_old(i,j);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
SOR_iteration = SOR_iteration + 1;
end
end
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03);
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Successive Over-Relaxation Iterations(IMPLICIT) = %d', SOR_iteration));
Result
IMPLICIT Scheme Unsteady/Transient-state Condition Program:
Unsteady-state Equation discretisation:
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.
In Jacobi, the unknowns on the right-hand side are assumed to have an initial value, i.e. T at \'n+1\' at all nodes is assumed. Using the equation, then the LHS is calculated. On the next value of i,j similar method is applied to calculate LHS. Thus in the first iteration, value of T for time \'n+1\' is calculated at all nodes. The difference between T at \'n+1\' which is obtained and T at \'n+1\' which was assumed previously is taken. If it is greater than a specified tolerance value, iterations are carried out to calculate T at \'n+1\'. But in the second iteration, T at \'n+1\' is assumed as the value of T which was obtained in the previous iteration.
Once the value of T is found out iteratively whose error is less than the tolerance value, the value of T is calculated at the next time step \'n+2\' using the same iterative method.
Equation for Jacobi method are:
term1 = (1+(2*k1)+(2*k2))^(-1);
term2= k1*term1;
term3= k2*term1;
H=(Told(i-1,j)+Told(i+1,j));
V=(Told(i,j-1)+Told(i,j+1));
T(i,j) = (T_intial(i,j)*term1) + (H*term2) + (V*term3);
MATLAB code
%Solving the IMPLICIT Unsteady state 2D heat conduction by Jacobi Method
close all
clear all
clc
%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;
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 original values to T
T_old = T;
T_intial = T;
%Calculation of 2D Unsteady heat conduction EQUATION by Jacobi 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
term1 = (1 + (2*k1) + (2*k2))^(-1);
term2 = k1*term1;
term3 = k2*term1;
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)*term1) + (H*term2) + (V*term3);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
Jacobi_iteration = Jacobi_iteration + 1;
end
T_intial = T;
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03)
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Unsteady Jacobi Iterations(IMPLICIT) = %d', Jacobi_iteration));
endResult
Solving the IMPLICIT Unsteady state 2D heat conduction by Gauss Seidel Method
The gauss Seidel method is similar to the Jacobi method, but the only difference is- For the first iteration when the value of T at a particular node is calculated, this updated value at that node is used to calculate the value at the next node. But, the error term is calculated similar to Jacobi method, i.e. T at \'n+1\' calculated minus the T at \'n+1\' previously assumed.
Once the value of T is found out iteratively whose error is less than the tolerance value, the value of T is calculated at the next time step \'n+2\' using the same iterative method.
equations for gauss Seidel are
term1 = (1+(2*k1)+(2*k2))^(-1);
term2= k1*term1;
term3= k2*term1;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j) = (T_intial(i,j)*term1) + (H*term2) + (V*term3);
MATLAB code
%Solving the IMPLICIT Unsteady state 2D heat conduction by Gauss Seidel Method
close all
clear all
clc
%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;
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 original values to T
T_old = T;
T_intial = T;
%Calculation of 2D Unsteady heat conduction EQUATION by Gauss Seidel method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
GS_iteration = 1;
for k = 1:nt
error = 9e9;
while(error > tolerance)
for i = 2:nx - 1
for j = 2:ny - 1
term1 = (1 + (2*k1) + (2*k2))^(-1);
term2 = k1*term1;
term3 = k2*term1;
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)*term1) + (H*term2) + (V*term3);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
GS_iteration = GS_iteration + 1;
end
T_intial = T;
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03)
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Unsteady Gauss-Seidel Iterations(IMPLICIT) = %d', GS_iteration));
endResult
Solving the IMPLICIT Unsteady state 2D heat conduction by Successive over-relaxation Method(SOR)
SOR:
This method can be applied either by using Jacobi or Gauss Seidel but the only difference is consideration of factor called omega(scaling factor) which in turn makes the process faster by which the number of iterations required to solve will be decreased
The optimum value of omega is calculated by using the formula:omega=2/(1+sin(pi*dx))
To avoid lengthy execution of program omega assumed value 1.2
Equations for SOR:
term1 = (1+(2*k1)+(2*k2))^(-1);
term2= k1*term1;
term3= k2*term1;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j)= omega*((T_intial(i,j)*term1)+(H*term2)+(V*term3))+(1-omega)*T_old(i,j);
MATLAB code
%Solving the IMPLICIT Unsteady state 2D heat conduction by Successive Over Relaxation Method(SOR)
close all
clear all
clc
%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;
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.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 original values to T
T_old = T;
T_intial = T;
T_end = T;
%Calculation of 2D Unsteady heat conduction EQUATION by Successive Over Relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
SOR_iteration = 1;
for k = 1:nt
error = 9e9;
while(error > tolerance)
for i = 2:nx - 1
for j = 2:ny - 1
term1 = (1+(2*k1)+(2*k2))^(-1);
term2 = k1*term1;
term3 = k2*term1;
H = (T(i-1,j)) + (T_old(i+1,j));
V = (T(i,j+1)) + (T_old(i,j-1));
T(i,j)= omega*((T_intial(i,j)*term1)+(H*term2)+(V*term3))+(1-omega)*T_old(i,j);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
SOR_iteration = SOR_iteration + 1;
end
T_intial = T;
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03)
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Unsteady SOR Iterations(IMPLICIT) = %d', SOR_iteration));
end
Result
EXPLICIT Scheme - Unsteady/Transient state Condition Program
Unsteady-state Equation discretisation(EXPLICIT):
Solving the above equation Explicitly:
Explicit Approach is used in numerical analysis for obtaining numerical approximations to the solution of time-dependent ordinary or partial differential equations as in required in computer simulations of a physical process. Explicit methods calculate the state of a system at a later time from the state of the system at the current time.
In this method, each differential equation approach contains only one unknown and can be solved explicitly for this one unknown using this method.
Here in this equation time derivative terms are calculated using the first order differential equations and whereas the term 2 has second-order derivative terms and elliptical in nature so information needs to be taken from all the points equally, So, therefore, a central differencing scheme is applied for the discretising the term2.
The discretised equation for 2D heat conduction equation for the explicit form is as follows:
Solving the EXPLICIT Unsteady state (Transient state Condition)2D heat conduction Jacobi method.
MATLAB code
%Solving the EXPLICIT Unsteady state (Transient state Condition)2D heat conduction Jacobi method.
close all
clear all
clc
%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;
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 original values to T
T_intial = T;
T_old = T;
%Calculation of 2D Unsteady heat conduction EQUATION by Jacobi method
iterative_solver = 1;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
if iterative_solver == 1
Explicit_Unsteady_iteration = 1;
for k = 1:nt
error = 1;
while(error > tolerance)
for i = 2:nx - 1
for j = 2:ny - 1
term1 = T_old(i,j);
term2 = k1*(T_old(i+1,j) - 2*T_old(i,j) + (T_old(i-1,j)));
term3 = k2*(T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1));
T(i,j) = term1 + term2 + term3;
end
end
error = max(max(abs(T_old - T)));
T_old = T;
Explicit_Unsteady_iteration = Explicit_Unsteady_iteration + 1;
end
end
end
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03);
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Unsteady Iterations(EXPLICIT) Jacobi Method = %d', Explicit_Unsteady_iteration));
Result
Solving the EXPLICIT Unsteady state (Transient state Condition)2D heat conduction Gauss Seidel method.
MATLAB code
%Solving the EXPLICIT Unsteady state (Transient state Condition)2D heat conduction Gauss Seidel method.
close all
clear all
clc
%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;
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.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 original values to T
T_intial = T;
T_old = T;
%Calculation of 2D Unsteady heat conduction EQUATION by Gauss Seidel method
iterative_solver = 1;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
if iterative_solver == 1
Explicit_Unsteady_iteration = 1;
for k = 1:nt
error = 1;
while(error > tolerance)
for i = 2:nx - 1
for j = 2:ny - 1
term1 = T_old(i,j);
term2 = k1*(T(i+1,j) - 2*T_old(i,j) + (T(i-1,j)));
term3 = k2*(T(i,j+1) - 2*T_old(i,j) + T(i,j-1));
T(i,j) = term1 + term2 + term3;
end
end
error = max(max(abs(T_old - T)));
T_old = T;
Explicit_Unsteady_iteration = Explicit_Unsteady_iteration + 1;
end
end
end
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03);
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Unsteady Iterations(EXPLICIT) GS Method = %d', Explicit_Unsteady_iteration));
Result
Solving the EXPLICIT Unsteady state (Transient state Condition)2D heat conduction SOR method.
MATLAB code
%Solving the EXPLICIT Unsteady state (Transient state Condition)2D heat conduction SOR method.
close all
clear all
clc
%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;
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.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 original values to T
T_intial = T;
T_old = T;
%Calculation of 2D Unsteady heat conduction EQUATION by SOR method
iterative_solver = 1;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
if iterative_solver == 1
Explicit_Unsteady_iteration = 1;
for k = 1:nt
error = 1;
while(error > tolerance)
for i = 2:nx - 1
for j = 2:ny - 1
term1 = T_old(i,j);
term2 = k1*(T(i+1,j) - 2*T_old(i,j) + (T(i-1,j)));
term3 = k2*(T(i,j+1) - 2*T_old(i,j) + T(i,j-1));
T(i,j) = omega*(term1 + term2 + term3)+(1-omega)*T_old(i,j);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
Explicit_Unsteady_iteration = Explicit_Unsteady_iteration + 1;
end
end
end
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
pause(0.03);
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Unsteady Iterations(EXPLICIT) SOR Method = %d', Explicit_Unsteady_iteration));
Result
Conclusion
- 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 converging solutions in both the schemes (EXPLICIT & IMPLICIT) compared to Jacobi & Gauss Seidel.
- In Steady state equation SOR methods yields the results very quickly when compare to other two methods.whereas Jacobi takes more number of iterations and more time to solve the unknown.
- So the SOR method is most efficient when compared to the other two methods to solve the steady-state equation.
- For transient equations, the time required to solve is more when compared to steady-state equations.So, therefore it is recommended to solve the 2d heat conduction equation with steady state conditions.
- In transient state equation, explicit methods yield the results quickly when compared to the implicit method, but is recommended under the situation when CFL is less than 0.25.
- When transient state equation is solved implicitly Jacobi takes less time when compare to the other two methods whereas SOR is more stable compared to Gauss Seide