Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
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.
In order to solve for dependent variable iterative solvers are used. Iterative solvers used are
1. Jacobi iterative solver
2. Gauss seidal iterative solver
3. Successive over relaxation iterative solver.
2D Heat conduction Equation
∂T∂t=α⋅(∂2T∂x2+∂2T∂y2)∂T∂t=α⋅(∂2T∂x2+∂2T∂y2)
The time derivative term ∂T∂x∂T∂x is hyperbolic in nature which can be discretized using forward differencing scheme and the spacial derivatives are elliptical in nature which can be discretized using central differencing scheme.
Discretization of Time derivative using forward difference scheme since the time derivative is varying in time direction which can be solved by marching method.
Discretization of spacial derivatives using central differencing scheme since the spacial derivatives are elliptical in nature, which uses information from neighbouring points to obtain a solution. Below are the discretization for x and y spacial derivatives.
From the PDE heat conduction equation,
where n is current time step and n+1 is next time step.
i,i+1 are the neighbouring points in x direction, j,j+1 are the neighbouring points in y direction.
△x△x is grid spacing in x direction. △y△y is grid spacing in y direction.
Iterative Solver:
1. Jacobi iterative solver:
In this method the temperature at n+1th iteration is found using old values of the temperature defined at nth iterations. Using the approximation obtained iterations procedure it continued untill the desired accuracy of the solution is obtained.
2. Gauss Seidal iterative solver:
In this type iterative solver, temperature at n+1 iterations uses values of n+1 iterations that have already been computed. By this way old approximated values are overwritten with new approximate values for the iterations performed. Unlike jacobi, Gauss seidal require one storage vector which can be overwritten with each computation. This helps for the solution to coverge faster than the jacobi method.
3. Successive over relaxation iterative solver.
In this type of iterative solver relaxation factor is used which helps in increasing the convergence rate. It goes by the formula
Tn+1=Tn⋅(1−ω)+ω⋅T(GS)orT(jacobi)Tn+1=Tn⋅(1-ω)+ω⋅T(GS)orT(jacobi)
where ωω is the relaxation factor . Its range is from 0<ωω<2
when ωω is between 0-1 it is called under relaxed state.
This condition is used when the temperature gradient is high and due to this condition the time step does not over shoot and helps to converge at a faster rate.
when ωω is between 1-2 it is called over relaxed state.
This condition is used when the temperature gradient is low and due to this condition larger time steps are used to obtain the steady state i.e., for the values to converge at faster rate.
TnTn is the value of temperature for previous iterations
T(GS) or T(jacobi) is the temperature of value for n+1th iteration using gauss seidal or jacobi iterative solver.
Steady State 2D Heat Conduction Equation :
Steady state is where the fluid properties does not change with respect to time.
From the 2D heat conduction equation since the system is in steady state time derivative is zero and diffusivity is constant.
Hence we the equation for steady state as:-
Discretization for Steady State:
Solving for T(i,j), we get
, 
Solving the 2D Steady Heat Conduction equation implicitly.
Formula for Jacobi iterative solver:-
Formula for Gauss Seidal Iterative Solver:-
(Tp)n+1=K1(TL+Tr)n+1+K2(TT+TB)n+1(Tp)n+1=K1(TL+Tr)n+1+K2(TT+TB)n+1
Formula for Successive over relaxation iterative solver:-
Transient state:
Unsteady state is where the fluid properties change with respect to time.
The formula for 2D transient heat conduction equation is :
∂T∂t=α⋅(∂2T∂x2+∂2T∂y2)∂T∂t=α⋅(∂2T∂x2+∂2T∂y2)
Discretization for Unsteady State heat conduction equation:-
Solving 2D transient heat conduction equation explicitly we get,
Tn+1P=(1−2⋅K1−2⋅K2)⋅TnP+K1(TL+TR)n+K2(TT+TB)nTPn+1=(1-2⋅K1-2⋅K2)⋅TPn+K1(TL+TR)n+K2(TT+TB)n
K1=α⋅△t△x2K1=α⋅△t△x2
K2=α⋅△t△y2K2=α⋅△t△y2
Stability Criteria:
K1+K2≤0.5K1+K2≤0.5 Stability criteria for 2D heat conduction equation when solving explicitly.
Solving for 2D Transient Heat Conduction equation implicitly we get,
Tn+1P=1K⋅[TnP+K1(TL+TR)n+1+K2(TT+TB)n+1]TPn+1=1K⋅[TPn+K1(TL+TR)n+1+K2(TT+TB)n+1]
K=1+2⋅K1+2⋅K2K=1+2⋅K1+2⋅K2
K1=α⋅△t△x2K1=α⋅△t△x2
K2=α⋅△t△y2K2=α⋅△t△y2
here n+1 is the next time step and n is the preceeding time step.
Formula for jacobi iterative solver:-
Formula for Gauss seidal
Formula for Successive over relaxation
Boundary Conditions:
Boundary conditions are defined at the boundaries of independent variable of the function. In this problem mesh size of 10x10 is considered for which the boundary conditions are defined at the extreme rows and columns of the mesh.
For the mesh size 10x10 shown above the corner boxes are average of the two boundaries intersecting at that point. The values other than the boundaries are taken as guess values which is in this case in 300K (initial condition for transient heat conduction equation).
Main Script:
%Code to study the behaviour of steady state and transient state analysis
%using explicit and implicit method.
clear all
close all
clc
%Defining spacial domain
nx = 10;ny=10;
%Specifying length of the domain.
L=1;
%defining space domain
x = linspace(0,L,nx);
y = linspace(0,L,ny);
%grid size
dx = x(2)-x(1);
dy = y(2)-y(1);
%defining temperature
t = 300*ones (nx,ny);
% [t] = 300*meshgrid(x,y);
t(1:end,1)=400;%<-- left boundary
t(end,1:end)=900;%<-- bottom boundary
t(1,1:end)=600;%<-- top boundary
t(1:end,end)=800;%<-- right boundary
% Edge temperature
t(1,1)= (t(1,2)+t(2,1))/2;
t(end,1) = (t(end,2)+t((end-1),1))/2;
t(1,end)= (t (1,(end-1))+t(2,end))/2;
t(end,end)=(t((end-1),(end))+ t(end,(end-1)))/2;
%tolerance
tol=1e-4;
%Making a copy of the initial temperature matrix.
told=t;
t_pre=t;
%Coefficient obtained from the derivation for steady state equation
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
%diffusivity co-efficient in the 2D heat conduction equation
alpha = 1.4;
% Calculation of time step
dt = (0.5)/((alpha/dx^2)+(alpha/dy^2));
% Co-efficient multipied for x diffuse term
k1 = dt*alpha/dx^2;
% Co-efficient multipied for y diffuse term
k2 = dt*alpha/dy^2;
%relaxation factor for successive over relaxation facotr for steady state
%heat conduction.
a = 1.5;
%relaxation factor for Successive over relaxation factor for transiet heat
%conduction
a1 = 1.15;
%Function to calculate the temperature of Steady state using jacobi iterative solver
[error_j,iter1] = jacobi(nx,ny,t,told,k,dx,dy,tol);
%Function to calculate the temperature of Steady state using Gauss Seidal iterative solver
[error_gs,iter2] = gs(nx,ny,t,told,k,dx,dy,tol);
%Function to calculate the temperature Steady state using Successive over relaxation iterative solver
[error_sor_gs,iter3] = sor_gs(nx,ny,t,told,k,dx,dy,tol,a);
%Function to calculate the temperature of Transient heat conduction explicitly
[error_dt,iter_dt]=transient_heat_conduction(nx,ny,k1,k2,told,t,tol);
%Function to calculate the temperature of Unsteady state using jacobi iterative solver
[error_jac,iter_jac]=transient_heat_conduction_jacobi(nx,ny,k1,k2,told,t,tol,t_pre);
%Function to calculate the temperature of Unsteady state using Gauss Seidal iterative solver
[error_gs_imp,iter_gs_imp]=transient_heat_conduction_gs(nx,ny,k1,k2,told,t,tol,t_pre);
%Function to calculate the temperature of Unsteady state using Successive over relaxation iterative solver
[error_sor_gs_imp,iter_sor_gs_imp]=transient_heat_conduction_sor_gs(nx,ny,k1,k2,told,t,tol,t_pre,a1);Explanation of the Main Script:
1. Specifying the number of grid points in the spacial domain x and y direction. In this case is 10 grid points in x and y direction.
2. Specifying the length of the domain in x and y space domain.
3. Defining the spacial domain for the length and grid points specified.
4. Calculating the grid spacing in the x and y space domain.
5. Defining temperature with the guess values 300K for the number of the grid points defined.
6. Applying the Boundary Conditions as mentioned in the problem.
7. The accuracy of the solution is measured with computed approximated solution and compared with the tolerance value so that the converged solution is in the desired accuracy.
8. A copy of the initial temperature matrix is made in told and t_pre which is used while calculating the temperature values of n+1 time steps explicitly and implicitly.
9. Calculating KK,K1andK2K1andK2. Here αα is taken 1.4.
10. Function is called to calculate the temperature values at n+1 time step for both steady and unsteady state using different iterative solvers and explicit method is also used for calculating temperature for unsteady state heat conduction equation.
Explanation of Functions defined in the main script:
In the functions dummy error is assigned to error variable so that the while loop condition would satisfy and iterations are calculated.
The error is calculated and updated for each iteration. The error calculated is coded to take only the absolute values since negative values will not statisfy the while loop condition and iterations will stopped and the results obtained will not be a converged results.
Number of iterations taken for the solution to converge is stored in the iter variable. This is done by increamenting the iter variable by 1 for each iteration.
As the iteration proceeds to calculate temperature at n+1 time step temperature at nth time step is required and these temperature values should be updated for each iteration. If the temperature values are not updated for each iteration the solution will not converge and while loop will run infinite times.
Contourf plot is used to check if the solutions are converged or not i.e., the plot does not change as the approximated temperature value changes in the order of tolerance value defined.
told stores the temperature values of the preceeding time step.
t stores the temperature values of the n+1 time step.
This logic is used in the below functions to calculate the temperature of Steady state Heat Conduction equation for the desired accuracy.
Functions are as follows:
Function to calculate the temperature of Steady state using jacobi iterative solver
%Function to calculate the temperature of Steady state using jacobi iterative solver
function [error_j,iter1] = jacobi(nx,ny,t,told,k,dx,dy,tol)
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter1 = 0;
% Dummy error value
error_j = 9e9;
%convergence loop
while (error_j>tol)
%Iteration loop
for i = 2:(nx-1)
for j = 2:(ny-1)
%Formula to calculate the temperature at each grid
%point
t(i,j)= (1/k)*(((told((i-1),j)+told((i+1),j))/dx^2)+((told(i,(j+1))+told(i,(j-1)))/dy^2));
end
end
%error for each iteration
error1 = abs(told-t);
error_j = max(error1,[],'all');
%updating the temperature values for the next iteration.
told = t;
%incrementing iter variable by 1 to display the no. of
%iterations taken for the solution to converge.
iter1 = iter1+1;
end
% Plot representing the converged results.
figure(1)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h);
axis ij;
title('number of iterations for Jacobi method for steady state heat conduction',iter1);
end
This functions uses jacobi iterative solver, temperature at n+1 time step is calculated using temperature at nth time step.
Function to calculate the temperature of Steady state using Gauss Seidal iterative solver
%Function to calculate the temperature of Steady state using Gauss Seidal iterative solver
function [error_gs,iter2] = gs(nx,ny,t,told,k,dx,dy,tol)
% Dummy error value
error_gs = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter2 = 0;
%convergence loop
while (error_gs>tol)
%Iteration loop
for i = 2:(nx-1)
for j = 2:(ny-1)
%Formula to calculate the temperature at each grid
%point
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));
end
end
%error for each iteration
error1 = abs(told-t);
error_gs = max(error1,[],'all');
%updating the temperature values for the next iteration.
told = t;
%incrementing iter variable by 1 to display the no. of
%iterations taken for the solution to converge.
iter2 = iter2+1;
end
% Plot representing the converged results.
figure(2)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h);
axis ij;
title('Number of iterations for Gauss seidal method steady state heat conduction',iter2)
end
This function uses Gauss seidal iterative solver in which temperature at n+1 time step is calculated using n+1 approximated temperature in the previous iterations.
Function to calculate the temperature of Steady state using Successive over relaxation iterative solver
%Function to calculate the temperature of Steady state using Successive over relaxation iterative solver
function [error_sor_gs,iter3] = sor_gs(nx,ny,t,told,k,dx,dy,tol,a)
%Dummy error value
error_sor_gs = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter3 = 0;
%Convergence loop
while (error_sor_gs>tol)
%iteration loop
for i = 2:(nx-1)
for j = 2:(ny-1)
%Formula to calculate the temperature at each grid
%point
t(i,j)= (1/k)*(((t((i-1),j)+told((i+1),j))/dx^2)+((told(i,(j+1))+t(i,(j-1)))/dy^2));
%Successive over relaxation Formula to calculate the
%temperature at each grid point.
t(i,j) = (told(i,j))*(1-a)+a*(t(i,j));
end
end
%error for each iteration
error1 = abs(told-t);
error_sor_gs = max(error1,[],'all');
%updating the temperature values for the next iteration.
told = t;
%incrementing iter variable by 1 to display the no. of
%iterations taken for the solution to converge.
iter3 = iter3+1;
end
% Plot representing the converged results.
figure(3)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h);
axis ij;
title('number of iteration for Successive over relaxation method for steady state heat conduction', iter3)
end
This function uses Successive over relaxation to calculate the temperature at n+1 time step. It uses temperature values calculated using Gauss seidal method and relaxation factor which helps in increasing the convergence rate.
Functions mentioned below uses Unsteady state head conduction equation formula to calculate the temperature values at n+1 time step.
Function to calculate the temperature of Unsteady state using explicit method
%Function to calculate the temperature of Unsteady state using explicit
%method
function [error_dt,iter_dt]=transient_heat_conduction(nx,ny,k1,k2,told,t,tol)
%dummy error
error_dt = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter_dt =0;
%convergence loop
while (error_dt>tol)
%Iteration loop
for i = 2:nx-1
for j = 2:ny-1
H = told(i-1,j)-2*told(i,j)+told(i+1,j);%<-- This line represent temperature values in x direction
V = told(i,j-1)-2*told(i,j)+told(i,j+1);%<-- This line represent temperature values in y direction
%Formula to calculate the temperature at each grid
%point
t(i,j)= told(i,j)+k1*H+k2*V;
end
end
%error for each iteration
error1 = abs(told-t);
error_dt = max(error1,[],'all');
%updating the temperature values for the next iteration.
told = t;
%incrementing iter variable by 1 to display the no. of
%iterations taken for the solution to converge.
iter_dt = iter_dt+1;
end
%Plot representing the converged results.
figure(4)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h);
axis ij
title('Number of iterations using transient heat conduction explicitly',iter_dt);
end
Function to calculate the temperature in Unsteady state using jacobi iterative solver
As we are now solving for the unsteady state heat conduction equation, steady state is achieved after a certain number of time steps.
Hence 1st while loop is used to check for the steady state of the solution by calculating the error between temperature between successive time steps and comparing it with the tolerance values. At steady state the temperature values do not change with increases in number of iterations and error value reduces below the tolerance value . With this the converge point is found out.
2nd while loop functions in the similar way as 1st while loop. But in 2nd while loop it finds the converged value of temperature for the particular time step.
This logic is used in the below functions which are used to solve 2D heat conduction equation implicitly.
%Function to calculate the temperature of Unsteady state using jacobi
%iterative solver
function [error_jac,iter_jac]=transient_heat_conduction_jacobi(nx,ny,k1,k2,told,t,tol,t_pre)
% dummy error
error_jac=9e9;
error_f = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter_jac=0;
%While loop to stop iterations when the temperature values converge
%whle computing for different time steps.
while (error_f>tol)
%Convergence loop each time step
while (error_jac>tol)
%Iteration loop
for i = 2:nx-1
for j = 2:ny-1
%Co-efficient obtained while deriving for transient heat conduction equation
term1= (1+(2*k1)+(2*k2))^-1;
%k1 and k2 represent the cfl value for x and y
%domain repectively.
H = told(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
V = told(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
%Formula to calculate the temperature at each grid
%point
t(i,j)= (t_pre(i,j)*term1)+(term1*k1*H)+(term1*k2*V);
end
end
%error for each iteration
error1 = abs(told-t);
error_jac = max(error1,[],'all');
%updating the temperature values for the next iteration.
told = t;
%incrementing iter variable by 1 to display the no. of
%iterations taken for the solution to converge.
iter_jac = iter_jac+1;
end
% error defined to obtain desired accuracy for the convergenc rate.
error_f = max(max(abs(t_pre-t)));
t_pre = t;
%Reseting the error value for the covergence loop to compute for the
%next time step.
error_jac = 9e9;
end
%Plot representing the converged results.
figure(5)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h)
axis ij
title('Number of iterations for implicit Jacobi',iter_jac);
end
Function to calculate the Temperature in Unsteady state using Gauss seidal iterative solver.
%Function to calculate the temperature of Unsteady state using gauss seidal
%iterative solver
function [error_gs_imp,iter_gs_imp]=transient_heat_conduction_gs(nx,ny,k1,k2,told,t,tol,t_pre)
%Definition of dummy error
error_f = 9e9;
error_gs_imp=9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter_gs_imp=0;
%While loop to stop iterations when the temperature values converge
%whle computing for different time steps.
while(error_f>tol)
%Convergence loop each time step
while (error_gs_imp>tol)
%Iteration loop
for i = 2:nx-1
for j = 2:ny-1
%Co-efficient obtained while deriving for transient heat conduction equation
term1= (1+(2*k1)+(2*k2))^-1;
%k1 and k2 represent the cfl value for x and y
%domain repectively.
H = t(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
V = t(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
%Formula to calculate the temperature at each grid
%point
t(i,j)= (t_pre(i,j)*term1)+(term1*k1*H)+(term1*k2*V);
end
end
%error for each iteration
error1 = abs(told-t);
error_gs_imp = max(error1,[],'all');
%updating the temperature values for the next iteration.
told = t;
%incrementing iter variable by 1 to display the no. of
%iterations taken for the solution to converge.
iter_gs_imp = iter_gs_imp+1;
end
% error defined to obtain desired accuracy for the convergenc rate.
error_f = max(max(abs(t_pre-t)));
%Updating the previous temperature values to computed values.
t_pre = t;
%Reseting the error value for the covergence loop to compute for the
%next time step.
error_gs_imp = 9e9;
end
%Plot representing the converged results.
figure(6)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h)
axis ij
title('Number of iterations for implicit Gauss seidal',iter_gs_imp);
end
Function to calculate the temperature in Unsteady state using successive over relaxation iterative solver
%Function to calculate the temperature of Unsteady state using successive
%over relaxation iterative solver
function [error_sor_gs_imp,iter_sor_gs_imp]=transient_heat_conduction_sor_gs(nx,ny,k1,k2,told,t,tol,t_pre,a1)
%Definition of dummy error
error_f = 9e9;
error_sor_gs_imp = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter_sor_gs_imp =0;
%While loop to stop iterations when the temperature values converge
%whle computing for different time steps.
while (error_f>tol)
%Convergence loop each time step
while (error_sor_gs_imp>tol)
%Iteration loop
for i = 2:nx-1
for j = 2:ny-1
%Co-efficient obtained while deriving for transient heat conduction equation
term1= (1+(2*k1)+(2*k2))^-1;
H = t(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
V = t(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
%Formula to calculate the temperature at each grid
%point
t(i,j)= (t_pre(i,j)*term1)+(term1*k1*H)+(term1*k2*V);
%Successive over relaxation Formula to calculate the
%temperature at each grid point.
t(i,j)= (told(i,j))*(1-a1)+ a1*t(i,j);
end
end
%error for each iteration
error1 = abs(told-t);
error_sor_gs_imp = max(error1,[],'all');
%updating the temperature values for the next iteration.
told = t;
%incrementing iter variable by 1 to display the no. of
%iterations taken for the solution to converge.
iter_sor_gs_imp = iter_sor_gs_imp+1;
end
%error defined to obtain desired accuracy for the convergenc rate.
error_f = max(max(abs(t_pre-t)));
%Updating the previous temperature values to computed values.
t_pre = t;
%Reseting the error value for the covergence loop to compute for the
%next time step.
error_sor_gs_imp = 9e9;
end
%Plot representing the converged results.
figure(7)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h)
axis ij
title('Number of iteration for implicit SOR',iter_sor_gs_imp);
end
Output Graph:
The above plots shows number of iterations required by the iterative solver to obtain the desired accurate solution.
Inference from the Plots:
Successive over relaxation iterative solver convergence rate is higher than the jacobi and gauss seidal method for both steady and unsteady state heat conduction equation. The number of iterations required by the successive over relaxation method is less compared to jacobi and gauss seidal method.
Below Table gives the number of iterations required by iterative solvers to converge.
| Steady State Heat Conduction Equation | |
| Iterative Solver | No. of iterations |
| Jacobi | 207 |
| Gauss Seidal | 111 |
| Successive over Relaxation | 28 |
| Unsteady State Heat Conduction Equation. | |
| Iterative Solver | No. of iterations |
| Jacobi | 1982 |
| Gauss Seidal | 1458 |
| Successive over Relaxation | 1148 |
Main Window
