Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
AIM
To solve the 2D heat conduction equation under steady and transient conditions using the explicit and implicit methods.
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.
- Jacobi
- Gauss-seidel
- Successive over-relaxation
Your absolute error criteria is 1e-4
- You will write code for both implicit and explicit schemes (This is for the transient analysis only)
- Compare the convergence rate [Number of iterations] of the steady state and transient [explicit and implicit] simulations and justify the numerics.
INTRODUCTION
Heat conduction (or thermal conduction) is the transmission of heat from one object to another one that has a different temperature, when they are in direct contact with each other. Is project we are going to solve the 2D heat conduction equation under steady and transient conditions using the explicit and implicit methods. The heat equation is a second order 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 place where it is higher towards places where it is lower. It is a special case of the diffusion equation.
In physics and engineering contexts, especially in the context of diffusion through a medium, it is more common to fix a Cartesian coordinate system and then to consider the specific case of a function T(x, y, z, t) of three spatial variables (x, y, z) and time variable t. One then says that u is a solution of the heat equation.
`(∂^2 T)/(∂t^2 )=α.((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )+(∂^2 T)/(∂z^2 ))`
THEORY
Discretization
Discretization methods are used to divide the region of a continuous function i.e., the real solution to a system of differential equations into a discrete function, where the solution values are defined at each point in space and time. Discretization simply refers to the spacing between each point in solution space and then replace derivatives in the governing equation with difference quotients.
Finite Difference Method
In numerical analysis, finite-difference methods (FDM) are a class of numerical techniques for solving differential equations by approximating derivatives with finite differences. Based on finite differences, it gives the numerical solutions for the values at discrete points in the domain which are known as grid points.
Finite differences method are used to compute approximate derivatives. This is used in numerical analysis in numerical ordinary differential equations and numerical partial differential equations, which gives the numerical solution of ordinary and partial differential equations respectively. The resulting methods are called finite-difference methods. Finite difference representations of derivatives are derived from Taylor series expansions.
Point Iterative Technique
An iterative method uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the n-th approximation is derived from the previous ones. A specific implementation of an iterative method, including the termination criteria, is an algorithm of the iterative method. An iterative method is called convergent if the corresponding sequence converges for given initial approximations.
Steady State
The meaning of Steady-state is that the properties measured in any field of space do not vary with respect to time. In general, even after some time(t) the properties in a considered point of interest do not vary and will remain constant throughout the process. This is the basic concept of a steady-state system. So when a system is said to be in a Steady-state condition, then, it is providing a better case to fix the assumption values which makes the calculation better and is also majorly used in Engineering problems and when we see this concept applied in the domain of heat transfer, our region of interest would be in finding the heat transfer across the field to find out its temperature profile.
Steady-state conduction is the form of conduction that happens when the temperature difference driving the conduction is constant, so that after an equilibration time, the spatial distribution of temperatures in the conducting object does not change any further. 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).
The Steady State form of the 2-D Heat Conduction Equation is given by,
`(∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )=0` or `∇^2 T=0` (Laplacian Equation)
Transient State (Unsteady State)
Heat transfer is the phenomenon that causes the flow of heat energy from the elevated temperature source to the lower temperature side. The modes of heat transfer also can be classified as Conduction, Convection, and Radiation, and here we deal only with conduction. Any system when involving a change in energy within itself due to temperature difference starts as an Unsteady state condition wherein the change in temperature occurs with respect to time. This is because during the start of the process there would be a huge temperature difference between the system and the surroundings due to which it adapts an unsteady or transient condition of heat transfer initially, but later on, when the inlet energy flow and the energy going out (which causes the transient condition) becomes equal, there will be only steady-state heat transfer prevailing inside the system till the end until achieving a Steady-state of the temperature profile, provided there is no external influence of any properties over the system that causes a change in the temperature profile.
During any period in which temperatures changes in time at any place within an object, the mode of thermal energy flow is termed transient conduction. Another term is "non-steady-state" conduction, referring to the time-dependence of temperature fields in an object.
The Transient form of the 2-D Heat Conduction Equation is given by,
`(∂T)/(∂t)+α((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))=0` or `(∂T)/(∂t)+α(∇T)=0`
Where, ∇T = Laplacian of T.
Explicit Method
When a direct computation is made for an unknown quantity in terms of known quantity. This computation is called Explicit. This method requires stability criteria which when satisfied leads to stable solution. Because of the stability criteria small time steps has to be taken to obtain stable solutions. With decrease in time step the number of computation increases and hence the computational time required for the solution to converge is less and accuracy of the solution is high.
Implicit Method
When a computation is made for dependent variable are defined by set of equations and a matrix is involved to find the unknown variable. This computation is called Implicit method. This method is unconditionally stable as a result, no stability criteria is required to obtain a stable results and the computational time required for the solution to converge is high and accuracy of the solution is less relative to explicit method because this method uses simultaneous equations and uses guess values as the initial values to find the approximated values. If the guess values are closer to the desired value the computation time is less. As there is no stability criteria requirement larger time steps can be used.
ITERATIVE SOLVERS
In CFD, we may have a set of linear equations which we can solve by two no. of ways (i) Direct method and (ii) Iterative method.
(i) Direct Method
The direct method uses both Gauss-Seidel and Jacobian methods where it can solve a (n x n) matrix with a fixed number of steps. But, this method is inefficient for some engineering problems, since some may have millions of cells forming a huge matrix to solve. Hence for the solving of such kind of huge matrix, the system must require more space. Therefore, the usage of the Direct Method is not a feasible way to get the linear equations solved in CFD.
Thus we go for an alter for solving such equations, that is by using Iterative Methods.
(ii) Iterative Method
The iterative methods are a better and efficient aliter for the Direct method. Here we will be solving the equations in a way that takes the previously calculated values for a particular grid point at its corresponding time step(for transient condition only) or from its previous time step to calculate the function at that specific point with the help of the reframed governing equations that follows the criterion of different iterative solvers, like,
- a) Jacobian Iterative Solver
- b) Gauss-Seidel Iterative Solver
- c) Successive Over Relaxation Iterative Solver
Of which all the above-mentioned iterative solvers are of Implicit forms. Apart from this, for the Unsteady-state condition, we also do have an explicit form/method to solve the problem statement.
The following flow-chart represents the follow-up of the procedure, how an Iterative solver functions.
- a) Jacobian Iterative Solver
In numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges.
It is one of the types of Iterative solvers, which solves for the diagonal elements of the matrix formed by the linear equations, and the values are approximated in each of the iterations and are just inserted again into the equation until it converges. It generally forms its right-hand side with known values i.e., taking the values from the previous iterations only.
- b) Gauss-Seidel Iterative Solver
Gauss-Seidel method is an iterative method to solve a set of linear equations and is very much similar to Jacobi’s method. This method is also known as the Liebmann method or the method of successive displacement. The name successive displacement is because the second unknown is determined from the first unknown in the current iteration, and the third unknown is determined from the first and second unknowns.
It is also similar to the Jacobian method, and also strictly requires the diagonal matrix for it to perfectly converge. The only difference is that it also forms the right-hand side with values that were just calculated from that corresponding iteration step.
- c) Successive Over Relaxation Iterative Solver
The method of successive over-relaxation (SOR) is a variant of the Gauss-Siedel method for solving a linear system of equations, resulting in faster convergence. A similar method can be used for any slowly converging iterative process
This method is used over the Jacobian (or) Gauss-Seidel methods to enhance their results even more quickly with an additional component called Scaling Factor(alpha).
TWO - DIMENSIONAL STEADY STATE HEAT CONDUCTION EQUATION
The steady state Two -Dimensional Heat Conduction Equation,
`(∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )=0`
`∇^2 T=0` (Laplacian Equation)
This is the second order partial differential equation. Laplace's equation are in the simplest form of elliptic partial differential equation.
To discretize this second order partial differential equation in space along (x & y axis) to use the Central Difference scheme for space derivative. The partial derivative can be approximated using the derivative definition will be applied for space.
`(∂^2 T)/(∂x^2 )=(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/dx^2`
`(∂^2 T)/(∂y^2 )=(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/dy^2`
Substituting above Equations in steady state Two -Dimensional Heat Conduction Equation,
`(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/dx^2 +(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/dy^2 =0`
Rearranging the Equation,
`(2T_(i,j))/dx^2 +(2T_(i,j))/dy^2 =(T_(i+1,j)+T_(i-1,j))/dx^2 +(T_(i,j+1)+T_(i,j-1))/dy^2`
Simplifying the equation by considering,
T(i,j) = TP
T(i+1,j) = TR
T(i-1,j) = TL
T(i,j+1) = TT
T(i,j-1) = TB
`T_p [2((dx^2 )+(dy^2 ))/((dx^2 )⋅(dy^2 ) )]=(T_R+T_L)/dx^2 +(T_T+T_B)/dy^2`
Considering `K=(2((dx^2 )+(dy^2 )))/((dx^2)⋅(dy^2))`
Rearranging the Equations,
`T_(i,j)=(1/k) (T_(i+1,j)+T_(i-1,j))/dx^2 +(1/k) (T_(i,j+1)+T_(i,j-1))/dy^2`
or
`T_(i,j)=(1/k)[(T_(i+1,j)+T_(i-1,j))/dx^2 +(T_(i,j+1)+T_(i,j-1))/dy^2 ]`
Jacobian Iterative Solver
`T(i,j)=((T(i-1,j)+T(i+1,j))/(k⋅dx^2 )+(T(i,j-1)+T(i,j-1))/(k⋅dy^2 ))^(old)`
Gauss-Seidel Iterative Solver
`T(i,j)=((T(i-1,j)^(old)+T(i+1,j))/(k⋅dx^2 )+(T(i,j-1)^(old)+T(i,j-1))/(k⋅dy^2 ))`
Successive Over Relaxation Iterative Solver
`T(i,j)=(1-ω)T(i,j)^(old)+ω((T(i-1,j)^(old)+T(i+1,j))/(k⋅dx^2 )+(T(i,j-1)^(old)+T(i,j-1))/(k⋅dy^2 ))`
TRANSIENT STATE ANALYSIS OF 2-D HEAT CONDUCTION EQUATION
TRANSIENT STATE
The Two -Dimensional Transient Heat Conduction equation is,
`(∂T)/(∂t)=α⋅[(∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )]`
This is the second order partial differential equation. Laplace's equation are in the simplest form of elliptic partial differential equation and time dependent partial differential equation.
To discretize this equation in space and time to use the Forward Difference scheme for time derivative and Central Difference scheme for space derivative. The partial derivative can be approximated using the derivative definition will be applied both for space and time;
Here the terms represents,
`(∂T)/(∂t)=Time Derivative term`
`[(∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )]=Space Derivative Term`
α = Thermal Diffusivity
Time derivative
Applying Forward difference scheme for the time derivative;
`(∂^2 T)/(∂x^2 )=(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/dx^2`
`(∂^2 T)/(∂y^2 )=(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/dy^2`
Substituting the above equations in Two -Dimensional Transient Heat Conduction equation
`(T_(i,j)^(n+1)-T_(i,j)^n)/(△t)=α⋅[(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/dx^2 +(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/dy^2 ]^n`
Rearranging the Equation;
`T_(i,j)^(n+1)=T_(i,j)^n+α⋅((Δt)/dx^2 )⋅(T_(i+1,j)-2T_(i,j)+T_(i-1,j) )+α⋅((Δt)/dy^2 )⋅(T_(i,j+1)-2T_(i,j)+T_(i,j-1) ) ]^n`
Considering,
`K_1=α⋅((Δt)/dx^2 )`
`K_2=α⋅((Δt)/dy^2 )`
`T_(i,j)^(n+1)=T_(i,j)^n+[(K_1 )⋅(T_(i+1,j)-2T_(i,j)+T_(i-1,j) )+(K_2 )⋅(T_(i,j+1)-2T_(i,j)+T_(i,j-1) )]^n`
where n+1 and n are two consecutive steps in time and i,j are neighbouring points of the discretized x an Y - coordinate (grid points). The only unknown in the equation written in this way is `T_(i,j)^(n+1)` .To solve for this unknown to get an equation that allows us to advance in time as Time Marching Procedure Method. This equation is used in program to solve the 2-D Heat conduction equation numerically.
Explicit And Implicit Method
- The explicit method and implicit method are numerical analysis methods used to solve a time-dependent differential equation.
- The explicit method calculates the system status at a future time from the currently known system status.
(ie;), if know the state at n, can calculate the state at n+1.
- In explicit method the stability is so low that need to use a step size small enough to prevent divergence.
- The implicit method calculates the system status at a future time from the system statuses at present and future times.
- The implicit method has high stability and converges if setted with proper parameters.
- As the implicit method can use a sufficiently large step size, it is suitable for solving equations that involve a long time.
- Also, in non-linear equations such as contact, it is difficult to predict a future from the past state.
Explicit Equation
`T_(i,j)^(n+1)=T_(i,j)^n+[(K_1 )⋅(T_(i+1,j)-2T_(i,j)+T_(i-1,j) )+(K_2 )⋅(T_(i,j+1)-2T_(i,j)+T_(i,j-1) )]^n`
or
`T_p^(n+1)=T_p^n+[(K_1 )⋅(T_R-2T_P+T_L )+(K_2 )⋅(T_T-2T_P+T_B )]^n`
Implicit equation
`T_(i,j)^(n+1)=T_(i,j)^n+[(K_1 )⋅(T_(i+1,j)-2T_(i,j)+T_(i-1,j) )+(K_2 )⋅(T_(i,j+1)-2T_(i,j)+T_(i,j-1) )]^(n+1)`
or
`T_(i,j)^(n+1)=(T_(i,j)^n+(K_1 )⋅(T_(i+1,j)+T_(i-1,j) )+(K_2 )⋅(T_(i,j+1)+T_(i,j-1) )^(n+1))/(1+(2⋅(K_1+K_2 )) )`
`T_P^(n+1)=(T_P^n+(K_1 )⋅(T_R+T_L )+(K_2 )⋅(T_T+T_B )^(n+1))/(1+(2⋅(K_1+K_2 )) )`
Courant Number
The Courant-Friedrichs-Lewy (CFL) condition is absolutely essential in solving the partial differential equations. This non dimensionless number arises while solving the time marching problems. It expresses that the distance that any information travels during the time step length within the mesh must be lower than the distance between mesh elements.
In other words, information from a given cell or mesh element must propagate only to its immediate neighbors. The derivation of the CFL condition leads to the formula for Courant number and is given by,
`CFL=α⋅(Δt)/(Δx^2 )≤0.25;α⋅(Δt)/(Δy^2 )≤0.25`
Where,
CFL = Courant number, u
α = Thermal Diffusivity
Δt = Time step size
Δx = The length between mesh elements.
When the Courant number exceeds the value of 0.5 for Heat conduction equation the solution will be unstable and blew off, instabilities are amplified throughout the domain and may cause divergence of the simulation.
`CFL=α⋅(Δt)[1/(Δx^2 )+1/(Δy^2 )]≤0.5`
Courant number can be controlled by altering the time or space interval. When the space interval is reduced or time interval is increased the courant number increases. Since courant number governs the stability of a problem in an explicit method it is very important to choose the intervals precisely.
SOLVING AND MODELLING APPROACH
STEP 1: Defining the input parameters for the domain.
STEP 2: Defining the distance to the number of nodes provided.
STEP 3: Assigning the Initial Temperature to the respective number of nodes .
STEP 4: Applying Boundary conditions to the Left, Right, Top and Bottom of the domain and assigning to the corner nodes.
STEP 6: Updating the Initial and old Temperature values.
STEP7: Creating a function to Compute the Iterative methods.
STEP 8: Using the iterative method (Jacobi, Gauss seidal and SOR method)to determine the Temperature distribution at steady state for the respective number of nodes.
STEP 9: Updating the current temperature Value.
STEP 10: Calling a function to detrmine the Temperature distribution @ steady state.
STEP 12: Using contourf, the temperature distribution are potted for x and y Direction.
MATLAB CODING
STEADY STATE HEAT CONDUCTION
Function Code:
% Creating a function to solve Iterative Methods,
function [T,iteration]=Iterative_solver(tol,error,It_s,nx,ny,T_old,T,dx,dy,k)
iteration=0;
switch It_s % Switching to Iterative Methods
%% Solving Using Jacobi Method,
case 1
while(error>tol) % Convergence loop
for i=2:nx-1 % Spacial loop on x
for j=2:ny-1 % spacial loop on y
H=(T_old(i-1,j)+T_old(i+1,j))/(dx^2);
V=(T_old(i,j-1)+T_old(i,j+1))/(dy^2);
T(i,j)=k*(H+V);
end
end
error=max(max(abs(T-T_old)));
T_old=T; % Updating Temperature
iteration=iteration+1;
end
%% Solving Using Gauss Seidal method,
case 2
while(error>tol) % Convergence loop
for i=2:nx-1 % Spacial loop on x
for j=2:ny-1 % spacial loop on y
H=(T(i-1,j)+T_old(i+1,j))/(dx^2);
V=(T(i,j-1)+T_old(i,j+1))/(dy^2);
T(i,j)=k*(H+V);
end
end
error=max(max(abs(T-T_old)));
T_old=T; % Updating Temperature
iteration=iteration+1;
end
%% Solving Using Successive Over Relaxation Method (SOR),
case 3
omega=1.5;
while(error>tol) % Convergence loop
for i=2:nx-1 % Spacial loop on x
for j=2:ny-1 % spacial loop on y
H=(T(i-1,j)+T_old(i+1,j))/(dx^2);
V=(T(i,j-1)+T_old(i,j+1))/(dy^2);
T(i,j)=k*(H+V);
T(i,j)=T_old(i,j)*(1-omega)+(omega*T(i,j));
end
end
error=max(max(abs(T-T_old)));
T_old=T; % Updating Temperature
iteration=iteration+1;
end
end
end
Main Code:
%% Solving 2-D Temperature Heat Conduction Equation,
% Steady State Heat Conduction Equation,
clear all;close all;clc
%% Inputting the parameters for Domain,
l=1; % Length of the domain
nx=40; % Number of Nodes on X -Axis
ny=nx; % Number of Nodes on y -Axis
% Defining the length of domain repect to number of nodes,
x=linspace(0,l,nx);
y=linspace(0,l,ny);
% Finding the Grid Size along x and y axis,
dx=x(2)-x(1); % Along X-axis
dy=y(2)-y(1); % Along Y-axis
%% Defining the Temperature,
T_i=250; % Initial condition for Temperature
T=T_i*ones(nx,ny);
% Applying Boundary Conditions,
T_L=400;
T_T=600;
T_R=800;
T_B=900;
% Assigning Boundary Conditions on Nodes,
T(1,:)=600; % Temperature on Top
T(end,:)=900; % Temperature on Bottom
T(:,1)=400; % Temperature on Left
T(:,end)=800; % Temperature on Right
% Assigning Average Temperature on Corner Grids,
T(1,1)=(T_L+T_T)/2;
T(1,ny)=(T_R+T_T)/2;
T(nx,1)=(T_L+T_B)/2;
T(nx,ny)=(T_B+T_R)/2;
% Updating the Temperature,
T_old=T;
% Defining Error & Tolerance;
tol=1e-4; % Tolerence
error=1e5; % Error
It_s=1:3; % For Different Iterative Schemes (Jcobo - 1, Gauss Seidal- 2, SOR - 3)
% Calculation
k=(dx^2*dy^2)/(2*(dx^2+dy^2));
%% Calling a function to solve Iterative Solvers,
[T_J,iteration_jacobi]=Iterative_solver(tol,error,It_s(1),nx,ny,T_old,T,dx,dy,k); % Calling Jacobi Method
[T_GS,iteration_GS]=Iterative_solver(tol,error,It_s(2),nx,ny,T_old,T,dx,dy,k); % Calling Gauss Seidal Method
[T_SOR,iteration_SOR]=Iterative_solver(tol,error,It_s(3),nx,ny,T_old,T,dx,dy,k); % Calling SOR Method
%% Plotting,
% Plotting for Temperature Distribution by Jacobi Method
figure("Name",'STEADY STATE ANALYSIS')
contourf(x,y,T_J,'ShowText','on','LineWidth',0.5);
colorbar
colormap(jet)
set(gca,'YDir','reverse')
text_title=({'2-DSteady State Heat Conduction (Jacobi Method)';['No.of.Grid Points =',num2str(nx)];['Total Iterations =',num2str(iteration_jacobi)]});
title(text_title)
xlabel('X-Direction')
ylabel('Y-Direction')
% Plotting for Temperature Distribution by Gauss Seidal Method
figure("Name",'STEADY STATE ANALYSIS')
contourf(x,y,T_GS,'ShowText','on','LineWidth',0.5);
colorbar
colormap(jet)
set(gca,'YDir','reverse')
text_title=({'2-D Steady State Heat Conduction (Gauss Seidal Method)';['No.of.Grid Points =',num2str(nx)];['Total Iterations =',num2str(iteration_GS)]});
title(text_title)
xlabel('X-Direction')
ylabel('Y-Direction')
% Plotting for Temperature Distribution by Successive Over Relaxation (SOR) Method
figure("Name",'STEADY STATE ANALYSIS')
contourf(x,y,T_SOR,'ShowText','on','LineWidth',0.5);
colorbar
colormap(jet)
set(gca,'YDir','reverse')
text_title=({'2-D Steady State Heat Conduction (SOR Method)';['No.of.Grid Points =',num2str(nx)];['Total Iterations =',num2str(iteration_SOR)]});
title(text_title)
xlabel('X-Direction')
ylabel('Y-Direction')
Code Explanation
- Inputting the Domain parameters for length.
- Assigning the number of nodes on x and y axis as (nx and ny) respectively, based on the number of nodes the array for the x-axis and y-axis are assigned to the length of the domain using linspace command, equal intervals of arrays are created over the domain length.
- Determining the distance between the two nodes dx and dy , by differencing the second node to the first node.
- Defining the initial temperature values to the number of respective nodes using ones command. Ones command creating the value of one to the respective number of nodes.
- Applying Boundary conditions on top, bottom, left and right side of the domain using the respective commands and the corner node temperature are assigned by averaging the temperature on the nearby nodes.
- Updating the initial and old temperature values.
- Defining the error and tolerance criteria as 1e-4.
- Creating a function to solve iterative methods (Jacobi, Gauss seidal and SOR method) by passing the inputs as, tolerance, error, nx,ny,Temperature values and calculated terms.
- In the function the iterative methods are done using the convergence and nodel loops.
- The Switch case command is used to iterate for different methods based on given inputs respectively.
- The convergence loop- while loop is used to perform the error and tolerance which exceeds the numeric solutions are close to the exact solutions.
- For 2-D two for loops are assigned to perform nodel iterations on x and y axis.
- The error are determined for the iterated temperature values and the Temperatures values are updated.
- Using Iteration counter the number of iterations are counted for all the three iterative methods.
- By calling a function the temperature distribution for the different iterative methods are computed and plotted accordingly.
- Plotting the Temperature distribution profiles for the respective number of time step's dt for the domain along x and y axis.
- Using Contourf command the temperature distribution for the respective number of nodal points are plotted.
TRANSIENT STATE HEAT CONDUCTION
EXPLICIT:
%% Solving 2-D Temperature Heat Conduction Equation,
% Transient Heat Conduction Equation, By Using Explicit Method,
clear all;close all;clc
% Inputting the parameters for Domain,
l=1; % Length of the domain
nx=40; % Number of Nodes on X -Axis
ny=nx; % Number of Nodes on y -Axis
% Defining the length of domain repect to number of nodes,
x=linspace(0,l,nx); % Along X-Direction
y=linspace(0,l,ny); % Along Y- Direction
% Finding the Grid Size along x and y axis,
dx=x(2)-x(1);
dy=y(2)-y(1);
% Defining the Initial Temperature,
T_i=250; % Initial condition for Temperature
T=T_i*ones(nx,ny); % Assigning Temperature on nodes
% Applying Boundary Conditions,
T_L=400; % Temperature on the Left
T_T=600; % Temperature on the TOP
T_R=800; % Temperature on the RIGHT
T_B=900; % Temperature on the BOTTOM
% Assigning Boundary Conditions on Nodes,
T(1,:)=600; % Temperature on Top
T(end,:)=900; % Temperature on Bottom
T(:,1)=400; % Temperature on Left
T(:,end)=800; % Temperature on Right
% Assigning Average Temperature on Corner Grids,
T(1,1)=(T_L+T_T)/2;
T(1,ny)=(T_R+T_T)/2;
T(nx,1)=(T_L+T_B)/2;
T(nx,ny)=(T_B+T_R)/2;
% Updating the current Temperature,
T_old=T;
% Defining the Time Steps
t=0.15; % Time for Temperature Distribution
dt=1e-4; % Time steps
nt=t/dt; % Number of Time Steps
alpha=1.4; % Thermal Diffusivity
% Calculations,
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
CFL=k1+k2;
%% TIME MARCHING PROCEDURE,
Iteration=1;
for k=1:nt % Time loop
for i=2:nx-1 % Space loop
for j=2:ny-1 % Space loop
H=(T_old(i+1,j)-2*T(i,j)+T_old(i-1,j));
V=(T_old(i,j+1)-2*T(i,j)+T_old(i,j-1));
T(i,j)=T_old(i,j)+k1*(H)+k2*(V);
end
end
T_old=T;
Iteration=Iteration+1;
end
%% PLOTTING FOR EXPLICIT METHOD
figure('Name','TRANSIENT STATE ANALYSIS - EXPLICIT')
contourf(x,y,T,'ShowText','on','LineWidth',0.5);
colorbar
colormap(jet)
set(gca,'YDir','reverse')
title('2-D Transient State Heat Conduction (Explicit Method)')
text=sprintf(['Total Iterations =%d; ','CFL=%g\n ','Grid Points=%g; ','dt=%g; ','dx=%g'],...
Iteration,CFL,nx,dt,dx);
subtitle(text)
xlabel('X-Direction')
ylabel('Y-Direction')
Code Explanation:
- The Explicit-Scheme tells us that in order to find out the value of the node at the next time step we need the value of that node from its previous matrix and the value of its neighbouring nodes.
- The entire matrix is already known to us here and we just have to increment the value of our time.
- We don't require any iterative solvers here as we can reach convergence just by going up in time.
- The only thing we should be concerned about here is how long we need to increase our value of time.
- We can estimate that by using a convergence variable which will enable us to subtract the temperature matrix at the time, t = 0 and t = 1, and so on so that we can compare the value we got after subtracting it with one another and compare it with the value of convergence to determine whether or not it has reached convergence.
- The only iteration present here is for time-iteration.
IMPLICIT:
%% Solving 2-Dimensional Temperature Heat Conduction Equation,
% Transient Heat Conduction Equation, By Using Implicit Method,
clear all;close all;clc
%% Inputting the Domain Parmeters,
l=1; % Length of the domain
nx=40; % Number of Nodes on X -Axis
ny=nx; % Number of Nodes on y -Axis
% Defining the length of domain repect to number of nodes,
x=linspace(0,l,nx);
y=linspace(0,l,ny);
% Finding the Grid Size along x and y axis,
dx=x(2)-x(1); % On X-axis
dy=y(2)-y(1); % on Y-axis
%% Defining the Initial Temperature,
T_i=250; % Initial condition for Temperature
T=T_i*ones(nx,ny); % Assigning Temperature on nodes
% Applying Boundary Conditions,
T_L=400; % Temperature on the Left
T_T=600; % Temperature on the TOP
T_R=800; % Temperature on the RIGHT
T_B=900; % Temperature on the BOTTOM
% Assigning Boundary Conditions on Nodes,
T(1,:)=600; % Temperature on Top
T(end,:)=900; % Temperature on Bottom
T(:,1)=400; % Temperature on Left
T(:,end)=800; % Temperature on Right
% Assigning Average Temperature on Corner Grids,
T(1,1)=(T_L+T_T)/2;
T(1,ny)=(T_R+T_T)/2;
T(nx,1)=(T_L+T_B)/2;
T(nx,ny)=(T_B+T_R)/2;
% Updating the Temperatures,
T_old=T;
T_prev=T;
% Defining the Time Steps
t=0.2; % Time for Temperature Distribution
dt=1e-4; % Time steps
nt=t/dt; % Number of Time Steps
alpha=1.4; % Thermal Diffusivity
% Calculations,
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
CFL=k1+k2;
Term1=1/(1+(2*(k1+k2)));
Term2=k1*Term1;
Term3=k2*Term1;
% Defining the Tolerence,
tol=1e-4; % Tolerence Creteria
omega=1.2; % Relaxation Factor
Iteration=1; % Initial Iteration
%% Input User Selection Prompt,
prompt=(['Select Iterative Method\n (1).JACOBI ITERATIVE METHOD\n ' ...
'(2).GAUSS SEIDAL ITERATIVE METHOD\n ' ...
'(3).SOR ITERATIVE METHOD\n']);
Iterative_Method=input(prompt);
%% Solving 2-D Using Iterative Method,
% Time loop,
tic
for k=1:nt
error=1e5; % Initialize Error on each iteration time loop
if Iterative_Method==1 %% Jacobi Iterative Method
% Convergence loop
while error>tol
% Nodal loop
for i=2:nx-1
for j=2:ny-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_prev(i,j)*Term1)+(Term2*H)+(Term3*V));
end
end
error=max(max(abs(T-T_old)));
% Updating temperatures,
T_old=T;
Iteration=Iteration+1;
end
end
%% Gauss Seidal Iterative Method
if Iterative_Method==2
% Convergence loop
while error>tol
% Nodal loop
for i=2:nx-1
for j=2:ny-1
H=(T_old(i+1,j)+T(i-1,j));
V=(T_old(i,j+1)+T(i,j-1));
T(i,j)=((T_prev(i,j)*Term1)+(Term2*H)+(Term3*V));
end
end
error=max(max(abs(T-T_old)));
% Updating temperatures,
T_old=T;
Iteration=Iteration+1;
end
end
%% Successive Over Relaxation (SOR) Iterative Method
if Iterative_Method==3
% Convergence loop
while error>tol
% Nodal loop
for i=2:nx-1
for j=2:ny-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_prev(i,j)*Term1)+(Term2*H)+(Term3*V));
T(i,j)=(T_old(i,j)*(1-omega))+omega*T(i,j);
end
end
error=max(max(abs(T-T_old)));
% Updating temperatures,
T_old=T;
Iteration=Iteration+1;
end
end
T_prev=T;
sim_time_Jacobi=toc; % Simulation time for Jacobi Iteration
sim_time_Gs=toc; % Simulation time for Gauss seidal Iteration
sim_time_SOR=toc; % Simulation time for SOR Iteration
%% Plotting
figure(1)
contourf(x,y,T,'ShowText','on','LineWidth',0.5)
colorbar
colormap('jet')
set(gca,'YDir','reverse')
if Iterative_Method==1
title('2-D Transient State Heat Conduction (Implicit - Jacobi Method)')
text=sprintf(['Total Iterations =%d; ','CFL=%g\n ','Grid Points=%g; ','dt=%g; ','dx=%g\n','Simulation Time=%g'],...
Iteration,CFL,nx,dt,dx,sim_time_Jacobi);
subtitle(text)
elseif Iterative_Method==2
title('2-D Transient State Heat Conduction (Implicit - Gauss Seidal Method)')
text=sprintf(['Total Iterations =%d; ','CFL=%g\n ','Grid Points=%g; ','dt=%g; ','dx=%g\n','Simulation Time=%g'],...
Iteration,CFL,nx,dt,dx,sim_time_Gs);
subtitle(text)
elseif Iterative_Method==3
title('2-D Transient State Heat Conduction (Implicit - SOR Method)')
text=sprintf(['Total Iterations =%d; ','CFL=%g\n ','Grid Points=%g; ','dt=%g; ','dx=%g\n','Simulation Time=%g'],...
Iteration,CFL,nx,dt,dx,sim_time_SOR);
subtitle(text)
end
xlabel('X-Direction')
ylabel('Y-Direction')
end
Code Explanation:
- Inputting the Domain parameters for length.
- Assigning the number of nodes on x and y axis as (nx and ny) respectively, based on the number of nodes the array for the x-axis and y-axis are assigned to the length of the domain using linspace command, equal intervals of arrays are created over the domain length.
- Determining the distance between the two nodes dx and dy , by differencing the second node to the first node.
- Defining the initial temperature values to the number of respective nodes using ones command. Ones command creating the value of one to the respective number of nodes.
- Applying Boundary conditions on top, bottom, left and right side of the domain using the respective commands and the corner node temperature are assigned by averaging the temperature on the nearby nodes.
- Updating the initial and old temperature values.
- Defining the error and tolerance criteria as 1e-4.
- Defining the time steps to solve the temperature distribution from transient state to steady state for convergence of the solution.
- Assigning the thermal diffusivity value and relaxation factor to solve the iterative methods (Jacobi, Gauss seidal and Successive over relaxation Methods) numerically.
- Based on the selection the respective iterative methods are selected and computed accordingly.
- Using Time marching method the Explicit method for 2-D heat conduction equation are analysed numerically and updating the temperature values.
- Determining the Number of time steps (nt) by using the simulation run time for the respective time steps.
- Using "for" loop the time step (dt) are iterated for varying number of time steps.
- Using the Forward Differencing approximations the time derivatives are derived.
- Using the Central Differencing approximations the space derivatives are derived along x and y axis.
- By combing the time and space derivatives the respective the future time step for the respective temperature distribution for the domain's are determined.
- By using the "for" loop the time loop are iterated with to the respective space loop for the respective nodes.
- After the space loop iterated, updating the Temperatures for the respective time and space loop.
- Computing the simulation time for each iteration on time steps using "Tic Toc" command tic works with the toc function to measure elapsed time.
- The tic function records the current time, and the toc function uses the recorded value to calculate the elapsed time.
- Plotting the Temperature distribution profiles for the respective number of time step's (dt's) to the number of iterations.
- Using Contourf command the temperature distribution for the respective number of nodal points are plotted along x and y axis.
RESULT
STEADY STATE HEAT CONDUCTION
- The plotting describes, Temperature distribution along x and y axis for the steady state analysis of 2-D heat conduction equation.
- From the plotting the temperature at the bottom of domain as the higher temperature can be distributed over the domain.
- The Jacobi Iterative method are performed for steady state analysis and for 3084 iterations are required to get the convergence of steady state for the domain.
- The plotting describes, Temperature distribution along x and y axis for the steady state analysis of 2-D heat conduction equation.
- From the plotting the temperature at the bottom of domain as the higher temperature can be distributed over the x and y co-ordinates of a domain.
- The Gauss seidal Iterative method are performed for steady state analysis and for 1654 iterations are required to get the convergence of steady state for the domain.
- Compared to Jacobi method, the gauss seidal method converge faster.
- The plotting describes, Temperature distribution along x and y axis for the steady state analysis of 2-D heat conduction equation.
- From the plotting the temperature at the bottom of domain as the higher temperature can be distributed over the x and y co-ordinates of a domain.
- The Successive over relaxation iterative method are performed for steady state analysis and for 606 iterations are required to get the convergence of steady state for the domain.
- Compared to Jacobi method and gauss seidal method the SOR iterative method converge faster.
TRANSIENT STATE
EXPLICIT
- The plotting describes, Temperature distribution along x and y axis for the Transient state analysis of 2-D heat conduction equation.
- From the plotting the temperature at the bottom of domain as the higher temperature can be distributed over the x and y co-ordinates of a domain.
- The Explicit method are performed for Transient state analysis and for 1500 iterations are performed to get the convergence of steady state for the domain @ time step of 1e-4 @ dx= 0.02564.
- The CFL criteria for the Explicit method plays crucial role in the stability of the solution.
- If the CFL values exceeds 0.5 the solution will be unstable based on the respective time step.
IMPLICIT
- The plotting describes, Temperature distribution along x and y axis for the Transient state analysis of 2-D heat conduction equation implicitly.
- From the plotting the temperature at the bottom of domain as the higher temperature can be distributed over the domain.
- The Jacobi Iterative method are performed for transient state analysis and for 19732 number of iterations are required to get the convergence of steady state for the domain.
- The convergence loop are performed under the time loop "nt" where the older value of temperature and previous time step values of temperature are updated accordingly.
- The plotting describes, Temperature distribution along x and y axis for the Transient state analysis of 2-D heat conduction equation implicitly.
- From the plotting the temperature at the bottom of domain as the higher temperature can be distributed over the domain.
- The Gauss Seidal Iterative method are performed for transient state analysis and for 14223 number of iterations are required to get the convergence of steady state for the domain.
- The convergence loop are performed under the time loop "nt" where the older value of temperature and previous time step values of temperature are updated accordingly.
- Compared to Jacobi Iterative method the Gauss Seidal method converge faster.
- The plotting describes, Temperature distribution along x and y axis for the Transient state analysis of 2-D heat conduction equation implicitly.
- From the plotting the temperature at the bottom of domain as the higher temperature can be distributed over the domain.
- The Successive Over Relaxation Iterative method are performed for transient state analysis and for 9477 number of iterations are required to get the convergence of steady state for the domain.
- The convergence loop are performed under the time loop "nt" where the older value of temperature and previous time step values of temperature are updated accordingly.
- Compared to Jacobi Iterative method and Gauss Seidal method, the SOR Iterative method converge faster based on the relaxation factor.
- The relaxation factor are considered as a certain range to get the solution converge faster then gauss seidal method.
CONCLUSION
- The 2-Dimensional Steady state heat conduction equation are analysed numerically by using iterative methods (Jacobi, Gauss Seidel and Successive over Relaxation Methods).
- The 2-Dimensional Transient state heat conduction equation are analysed numerically by Explicit and Implicit methods.
- The Thermal diffusivity value are selected as α=1.4 to diffuse the temperature distribution at transient heat conduction.
- The Relaxation factor ω are selected as in the value of 1.2, if the relaxation factor overturned in SOR method the iterations are overshoot then the targeted value and taking longer time of iterations and computation time then Gauss Seidel method.
- The Transient and Steady state analysis of two Dimensional eat conduction equation are analysed and the temperature distribution along the domain for the defined boundary conditions are performed.
- In steady state analysis the Number of iterations to converge for steady state for faster in SOR iterative method then Gauss Seidel Iterative method then Jacobi method.
- In Transient state Explicit analysis the Equations are directly solved on the time loop to get the desired temperature distribution along x and y axis.
- In Transient state Implicit analysis the system of equations are solved Using Iterative methods to get the convergence of the steady state solutions.