Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
AIM:- To Solve the steady and unsteady 2D heat conduction problem using the Iterative method
Objective:-
- To perform Steady-state analysis & Transient State Analysis.
- To understand the different iterative techniques used to solve the steady and unsteady 2D heat conduction problem
Theory:-
Heat conduction (or thermal conduction) in general is the movement of heat from one object to another one that has different temperatures when they are touching each other. Or the movement of heat from one end to the other positions/points on an object.
The heat conduction equation is a second-order partial differential equation that describes the distribution of heat energy in a given body over time.
2D heat conduction
`(del^2T)/(del t^2) + alpha ((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
Steady-state 2D heat conduction- In a steady-state heat transfer, the heat flow rate remains constant at any point in time. Steady-state heat transfer can be described in terms of the temperature as the variable; the temperature of the system remains constant over time in steady-state heat transfer after thermal equilibrium is attained.
`(del^2T)/(delx^2) + (del^2T)/(dely^2) = 0`
We can discretize the above equation into an algebraic form so that it can be used for computation,
For discretization, we ssume
`T_p = T_(i,j)` - Current point for which calculating
`T_(R) = T_(i+1,j)` - point on the right of current point
`T_(L) = T_(i-1,j)` - point on the right of current point
`T_(T) = T_(i,j+1)` - Top point of current point
`T_(B) = T_(i,j-1)` - Bottom point of current point
Now the Equation can written as,
`implies (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`
`implies (T_(R)−2T_(P)+T_(L))/dx^2+(T_(T)−2T_(P)+T_(B))/dy^2=0`
On Rearranging,
`implies T_(p) = (dy^2)/(2.(dx^2+dy^2)).(T_(R)+T_(L)) + (dx^2)/(2.(dx^2+dy^2)).(T_(T)+T_(B))`
Considering,
`K_(X) = (dy^2)/(2.(dx^2+dy^2))` and
`K_(Y) = (dx^2)/(2.(dx^2+dy^2))`
Therefore,
`T_(P) = K_(X).(T_(R)+T_(L))+K_(Y).(T_(T)+T_(B))`
Transient state 2D heat conduction- Transient heat transfer can also be called unsteady-state heat transfer. This type of heat transfer only exists for a short duration of time. In transient heat transfers, the heat energy transferred through the medium is not constant. The heat flow rate keeps changing, and the cause of the varying rates of heat transfer can be either fluctuating temperature differences over the medium or changes in properties of the medium.
In the physical world, heat transfer starts as transient and then reaches a steady state until thermal equilibrium is reached.
`(del^2T)/(del t^2) + alpha ((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
Where, `alpha` is thermal diffusivity
We can discretize the above equation into an algebraic form so that it can be used for computation,
`implies (T_(i,j)^(n+1)-T_(i,j)^(n))/dt+ alpha. ((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`
`implies (T_(NP)-T_(OP))/dt + alpha.((T_(OR)-2T_(OP)+T_(OL))/dx^2 + (T_(OT)-2T_(OP)+T_(OB))/dy^2=0`
`implies T_(NP)= T_(OP) + alpha.dt((T_(OR)-2T_(OP)+T_(OL))/dx^2 + (T_(OT)-2T_(OP)+T_(OB))/dy^2)`
Here, current time step (i.e. at n) is showned by suffix of 'O' representing Old time, and for the next time step (i.e. n+1) is by suffix of 'N'
considering,
`K_(XT) = (alpha.dt)/dx^2` and
`K_(YT) = (alpha.dt)/dy^2`
Therefore,
Implicit,
`T_(NP)= T_(OP) + K_(XT).(T_(OR)-2T_(OP)+T_(OL)) + K_(YT).(T_(OT)-2T_(OP)+T_(OB)`
Explicit,
`T_(NP)= T_(OP) + K_(XT).(T_(NR)-2T_(NP)+T_(NL)) + K_(YT).(T_(NT)-2T_(NP)+T_(NB)`
Iterative Solver:
1. Jacobi method- In this method, the temperature at (n+1) iteration is found using old values of the temperature defined at nth iterations. Using the approximation obtained iterations procedure is continued until the desired accuracy of the solution is obtained.
2. Gauss-Seidel method- In this type of iterative solver, the temperature at n+1 iterations uses values of n+1 iterations that have already been computed. In this way, old approximated values are overwritten with new approximate values for the iterations performed. Unlike Jacobi, Gauss seidel requires one storage vector which can be overwritten with each computation.
3. Successive over-relaxation method- 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-w)+T(GS)
where w is between O-1 it is called under 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 a faster rate
Tn">Tn is the value of temperature for previous iterations.
T(GS) or T(Jacobi) is the temperature of value for (n+1) iteration using gauss seidel or Jacobi iterative solver.
THERMAL DIFFUSIVITY AND CFL CRITERIA -
The Thermal Diffusivity (`alpha`) is an important factor in the solution, as it is related to the CFL Criteria shown below,
CFL = `(alpha.dt) ((1)/(delx^2) + (1)/(dely^2)) = 0`
Thermal Diffusivity of material is an important factor and is interrelated with the dt, dx, and the CFL number. There will be an optimum CFL i.e. optimum value of dt, dx, and dy.
For the smaller value of the `alpha` the size of dt will be higher in order to maintain stability i.e. here CFL no., but it reduces the accuracy of the solution, and for the larger value of `alpha` the size of dt will be small for the optimum CFL, but it can increase the simulation time. Since the implicit method is unconditionally stable at any value of time step, still, the CFL condition can be used in order to determine the optimum 'dt' in order to maintain time accuracies.
MATLAB Program:- For Steady-State 2D Heat Conduction using Implicit Method-
close all
clear all
clc
% Steady State Heat Conduction || d^2T/dx^2 + d^2T/dy^2
nx = 10;
ny = 10;
T = 450*ones(nx,ny);
TL = 400;
TR = 800;
TT = 600;
TB = 900;
% Boundary Conditions
T(:,1) = TL;
T(:,ny) = TR;
T(1,:) = TT;
T(nx,:) = TB;
% Edge Refinement
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,ny) = (T(1,(ny-1))+T(2,ny))/2;
T(nx,1) = (T((nx-1),1)+T(nx,2))/2;
T(nx,ny) = (T(nx,(ny-1))+T((nx-1),ny))/2;
% Saving a Copy
Told = T;
T_initial = T;
error = 9e9;
tol = 1e-4;
% Iteration Solvers - Jacobi || GS || SOR
% 1- Jacobi Method
jacobi_iter = 0;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = 0.25*(Told(i+1,j) + Told(i-1,j) +Told(i,j+1) + Told(i,j-1));
end
end
error = max(max(abs(T-Told)));
Told = T;
jacobi_iter = jacobi_iter +1;
figure(1)
contourf(flipud(T),'ShowText','ON')
title(sprintf('Steady State Jacobi Iteration = %d',jacobi_iter));
end
colorbar
% Reinitialising
T = T_initial;
Told = T;
error = 9e9;
%2- Gauss Seidel Method
gs_iter = 0;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = 0.25*(Told(i+1,j) + T(i-1,j) +Told(i,j+1) + T(i,j-1));
end
end
error = max(max(abs(T-Told)));
Told = T;
gs_iter = gs_iter +1;
figure(2)
contourf(flipud(T),'ShowText','ON')
title(sprintf('Steady State GS Iteration= = %d',gs_iter));
pause(0.003)
end
colorbar
% Reinitialising
T = T_initial;
Told = T;
error = 9e9;
% 3- SOR
sor_iter = 0;
omega = 1.5;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
TGS(i,j) = 0.25*(Told(i+1,j) + T(i-1,j) + T(i,j-1) +Told(i,j+1));
T(i,j) = (1-omega)*Told(i,j) + omega*TGS(i,j);
end
end
error = max(max(abs(T-Told)));
Told = T;
sor_iter = sor_iter +1;
figure(3)
contourf(flipud(T),'ShowText','ON')
title(sprintf('Steady State SOR Iteration= = %d',sor_iter));
pause(0.003)
end
colorbar
Output-
- Steady-State Analysis using Jacobi iterative solver
- Steady-State Analysis via Gauss-Siedel iterative solver
- Steady-State Analysis using Successive Over Relaxation with Gauss-Seidel solver
MATLAB Program:- For Transient-State 2D Heat Conduction using Implicit Method
close all
clear all
clc
% Transient State Heat Conduction || dt^2/dT^2 + aplha(d^2T/dx^2 + d^2T/dy^2)
nx = 10;
ny = 10;
T = 450*ones(nx,ny);
TL = 400;
TR = 800;
TT = 600;
TB = 900;
% Boundary Conditions
T(:,1) = TL;
T(:,ny) = TR;
T(1,:) = TT;
T(nx,:) = TB;
% Edge Refinement
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,ny) = (T(1,(ny-1))+T(2,ny))/2;
T(nx,1) = (T((nx-1),1)+T(nx,2))/2;
T(nx,ny) = (T(nx,(ny-1))+T((nx-1),ny))/2;
% Inputs
% Time Step
dt = 1e-2;
% Grid Spacing
dx = 1e-2;
dy = 1e-2;
% Thermal Diffusibility
alpha = 1.5;
k1 = alpha*dt/dx^2
k2 = alpha*dt/dy^2
% Saving a Copy
Told = T;
T_initial = T;
T_prev = T;
error = 9e9;
conv = 9e9;
tol = 1e-4;
% Iteration Solvers - Jacobi || GS || SOR
% 1- Jacobi Method
jacobi_iter = 0;
jacobi_time_iter = 0;
% Time Loop
while (conv > tol)
% Space Loop
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (1+2*k1+2*k2)^-1 *(T_prev(i,j)+ k1*(Told(i+1,j)+ Told(i-1,j)) +k2*(Told(i,j+1) + Told(i,j-1)));
end
end
error = max(max(abs(T-Told)));
Told = T;
jacobi_iter = jacobi_iter +1;
figure(1)
contourf(flipud(T),'ShowText','ON')
title(sprintf('Transient Implicit Jacobi Iteration = %d \n Time Iteration = %d',jacobi_iter,jacobi_time_iter));
pause(0.003)
end
conv = max(max(abs(T-T_prev)));
T_prev = T;
jacobi_time_iter = jacobi_time_iter +1;
error = 9e9;
end
colorbar
% Reinitialising
T = T_initial;
Told = T;
T_prev = T;
error = 9e9;
conv = 9e9;
%2- Gauss Seidel Method
gs_iter = 0;
gs_time_iter =0;
% Time Loop
while (conv > tol)
% Space Loop
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (1+2*k1+2*k2)^-1 *(T_prev(i,j)+ k1*(Told(i+1,j)+ T(i-1,j)) +k2*(Told(i,j+1) + T(i,j-1)));
end
end
error = max(max(abs(T-Told)));
Told = T;
gs_iter = gs_iter +1;
figure(2)
contourf(flipud(T),'ShowText','ON')
title(sprintf('Transient Implicit Gauss Siedel Iteration = %d \n Time Iteration = %d',gs_iter,gs_time_iter));
pause(0.003)
end
conv = max(max(abs(T-T_prev)));
T_prev = T;
gs_time_iter = gs_time_iter +1;
error = 9e9;
end
colorbar
% Reinitialising
T = T_initial;
Told = T;
T_prev = T;
error = 9e9;
conv = 9e9;
% 3- SOR
sor_iter = 0;
sor_time_iter = 0;
omega = 1.5;
% Time Loop
while (conv > tol)
% Space Loop
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
TGS(i,j) = (1+2*k1+2*k2)^-1 *(T_prev(i,j)+ k1*(Told(i+1,j)+ T(i-1,j)) +k2*(Told(i,j+1) + T(i,j-1)));
T(i,j) = (1-omega)*Told(i,j) + omega*TGS(i,j);
end
end
error = max(max(abs(T-Told)));
Told = T;
sor_iter = sor_iter +1;
figure(3)
contourf(flipud(T),'ShowText','ON')
title(sprintf('Transient Implicit SOR Iteration = %d \n Time Iteration = %d',sor_iter,sor_time_iter));
pause(0.003)
end
conv = max(max(abs(T-T_prev)));
T_prev = T;
sor_time_iter = sor_time_iter +1;
error = 9e9;
end
colorbarOutput-
- Transient Analysis with implicit using Jacobi iterative solver
- Transient Analysis with implicit using Gauss-Siedel iterative solver
- Transient Analysis with implicit using SOR with Gauss-Seidel solver.
MATLAB Program:- For Transient-State 2D Heat Conduction using Explicit Method
close all
clear all
clc
% Transient State Heat Conduction || dt^2/dT^2 + aplha(d^2T/dx^2 + d^2T/dy^2)
nx = 10;
ny = 10;
T = 450*ones(nx,ny);
TL = 400;
TR = 800;
TT = 600;
TB = 900;
% Boundary Conditions
T(:,1) = TL;
T(:,ny) = TR;
T(1,:) = TT;
T(nx,:) = TB;
% Edge Refinement
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,ny) = (T(1,(ny-1))+T(2,ny))/2;
T(nx,1) = (T((nx-1),1)+T(nx,2))/2;
T(nx,ny) = (T(nx,(ny-1))+T((nx-1),ny))/2;
% Inputs
% Time Step
dt = 1e-2;
% Grid Spacing
dx = 1.3e-2;
dy = 1.3e-2;
% Thermal Diffusibility
alpha = 0.004;
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
cfl = k1 + k2
% Saving a Copy
T_prev = T;
conv = 9e9;
tol = 1e-4;
time_iter = 0;
% Time Loop
while (conv > tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (1-2*k1-2*k2)*T_prev(i,j) + k1*(T_prev(i+1,j)+ T_prev(i-1,j)) +k2*(T_prev(i,j+1) + T_prev(i,j-1));
end
end
conv = max(max(abs(T-T_prev)));
T_prev = T;
time_iter = time_iter +1;
figure(1)
contourf(flipud(T),'ShowText','ON')
title(sprintf('Transient Explicit Iteration = %d',time_iter));
pause(0.003)
end
colorbar
Output-
- Transient Analysis with Explicit method
Result/Observation:-.
- The CFL number should be maintained at less than 0.5 so that we can obtain converged results.
- The implicit method is more accurate than the explicit method even though the number of iterations is much lesser since the implicit method uses all the points of the domain which makes it more accurate
- The number of iterations for steady-state is much lesser than the unsteady-state. That is because to change the unsteady state to steady-state we have to time march much longer than the steady-state therefore the iteration numbers are more
- The observation we made after simulation-
| Solving Method | Iterations |
| Steady-State Analysis using Jacobi iterative solver | 198 |
| Steady-State Analysis via Gauss-Siedel iterative solver | 108 |
| Steady-State Analysis using Successive Over Relaxation with Gauss-Seidel solver | 27 |
| Transient Analysis with implicit using Jacobi iterative solver | 438 |
| Transient Analysis with implicit using Gauss-Siedel iterative solver | 250 |
| Transient Analysis with implicit using SOR with Gauss-Seidel solver. | 74 |
|
Transient Analysis with Explicit method |
209 |
- From the Graphs obtained for all the methods, we can compare the total iteration taken to converge the result,
For Steady State - `I_(SOR) < I_(GS) < I_(JACOBI)`
For Transient State(Implicit methods) - `I_(SOR) < I_(GS) < I_(JACOBI)`
- Hence, we see that in both Steady-state and transient analysis the SOR method is the most efficient and fastest, taking the least Iterations to converge the solution. But, it is also important to consider that the Steady-State Analysis does not include the Thermal Diffusivity and the time derivative terms therefore they are faster than the transient analysis.
- The explicit method doesn't require an Iterative solver it solves directly, it derives the value of the new time step by the values at the previous timestep. Therefore, the iterations or the no. of solving the stiffness matrix is limited to no. of time steps.