Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Question:
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
Objectives:
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.
Solution:
We have temperature conduction equation for transient heat conduction(3d),
`(∂T)/(∂t)+α((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )+(∂^2 T)/(∂z^2 ))=0`
But here we want to solve for 2d, so equation becomes;
`(∂T)/(∂t)+α((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))=0`
Now to solve this challenge we have 2 states or types transient/unsteady state and steady state;
Where for transient state equation is;
`(∂T)/(∂t)+α((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))=0`
Which we can also write in form of Laplacian `(∂T)/(∂t)+α∇(T)=0`
Where as for steady state the equation is;
`(∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )=0`
In Laplacian form it becomes;
But as solution we have to solve this equation using the implicit & explicit form, where implicit requires iterative solvers so we have to derive this equation likewise in those form and then go ahead for solution
Transient state:
So, let’s start using the transient/unsteady state first;
If we see above the equation of transient state there is dt which means time derivative and states that we need to solve this using the time marching method.
So the equation becomes;
`(T_(i,j)^(n+1)-T_(i-1,j)^n)/(∆t)+α((T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(∆x^2 )+(T_(i,j+1)-2T_((i,j) )+T_(i,j-1))/(∆y^2 ))=0`
Here we can see for first term we apply the forward differencing and central differencing to the spacing and also a elliptic term.
Where let’s see in diagram below what are these i and j are;
Fig 1. Elliptic term nodes
Or else we can also name those points as p, l, r, t, and b meaning point(I,j), left(i-1,j), right(i+1,j), top(I,j+1), bottom(I,j-1)
Now if we simplify the above term we get;
`T_(i,j)^(n+1)=T_(i,j)^n+α∆t((T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(∆x^2 )+(T_(i,j+1)-2T_((i,j) )+T_(i,j-1))/(∆y^2 ))`
Now if we introduce the new variables k1 and k2
Where, `k1=(α∆t)/(∆x^2 )` and `k2=(α∆t)/(∆y^2 )` this gives equation;
`T_(i,j)^(n+1)=T_(i,j)^n+k1(T_(i-1,j)-2T_(i,j)+T_(i+1,j) )+k2(T_(i,j-1)-2T_(i,j)+T_(i,J+1))`
The above term we get is for 2d transient conduction using explicit scheme
Now for implicit we have to simplify it more (let’s use p, l, r, t and b notations)
And as we want for implicit the term gets n+1 which represents position, so we took all term to LHS
`T_p^(n+1)-K1T_L^(n+1)+2k1T_P^(n+1)-k1T_r^(n+1)-k2T_t^(n+1)+2K2T_p^(n+1)-k2T_b^(n+1)=T_p^n`
Now if we take Tp common we got;
`T_p^(n+1) (1+2k1+2k2)-K1T_L^(n+1)-k1T_r^(n+1)-k2T_t^(n+1)-k2T_b^(n+1)=T_p^n`
Lets re arrange the term;
`T_p^(n+1) (1+2k1+2k2)=T_p^n+K1T_L^(n+1)+k1T_r^(n+1)+k2T_t^(n+1)+k2T_b^(n+1)`
Now keeping Tpn+1 term alone at LHS and rest of all at RHS
`T_p^(n+1)=(1/((1+2k1+2k2) ))[T_p^n+K1T_L^(n+1)+k1T_r^(n+1)+k2T_t^(n+1)+k2T_b^(n+1) ]`
The term above we get is 2d transient conduction using Implicit method
Steady state Derivation:
Here for steady state we have equation;
`(∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )=0`
Using central differencing we first discreate the equation;(here we use p, l, r, t, and b as node points)
`(T_l-2T_p+T_r)/(∆x^2 )+(T_t-2T_p+T_b)/(∆y^2 )=0`
After simplifying above equation, we get;
`(2T_p)/(∆x^2 )+(2T_p)/(∆y^2 )=((T_l+T_r)/(∆x^2 ))+((T_t+T_b)/(∆y^2 ))`
And if we take Tp as common from LHS equation becomes;
`T_p ((2(∆x^2+∆y^2 ))/(∆x^2 ∆y^2 ))=RHS`
Now we introduce new variable k;
`k=(2(∆x^2+∆y^2 ))/(∆x^2 ∆y^2 )`
And if we rearrange the equation we get;
`T_p=1/k (((T_l+T_r)/(∆x^2 ))+((T_t+T_b)/(∆y^2 )))`
Above equation is equation for steady state conduction
Iterative solvers:
As of now we derived the equations for both transient and steady state now lets se what are the iterative solvers which are require to solve the implicit methods.
There are three types of iterative solvers;
- Jacobi: In Jacobian iterative method the new value is not updated until the whole set of equation completes ongoing iteration using guess values or previous iteration value
- Gauss Seidel: In this method of iteration the value get from the ongoing iteration is used in that same running iteration, so number of fresh values it get from iteration it utilizes and rest of values it takes as guess values or from previous iteration
- Successive over relaxation: As its name suggest this method is dependent on relaxation factor(ω) its divergence and convergence all depends on this relaxation factor.
Lets see how our equation becomes for various iterative solvers;
2D Steady state heat conduction (Implicit):
Jacobi solver:
`T_(i,j)=1/k (((T_(i-1,j)^(old)+T_(i+1,j)^(old))/(∆x^2 ))+((T_(i,j-1)^(old)+T_(i,j+1)^(old))/(∆y^2 )))`……………………..(1)
Gauss seidel:
`T_(i,j)=1/k (((T_(i-1,j)+T_(i+1,j)^(old))/(∆x^2 ))+((T_(i,j-1)+T_(i,j+1)^(old))/(∆y^2 )))`………………………(2)
Successive over relaxation:
`T_(i,j)=(1-ω) T_(i,j)^(old)+ω(1/k (((T_(i-1,j)+T_(i+1,j)^(old))/(∆x^2 ))+((T_(i,j-1)+T_(i,j+1)^(old))/(∆y^2 ))))` ……(3)
Where Told is value from previous iteration or initial guess value and T is fresh value
2D Unsteady/transient state heat conduction (implicit):
Jacobi:
`T_(i,j)=(1/((1+2k1+2k2) ))[T_(i,j)^(old)+K1T_(i-1,j)^(old)+k1T_(i+1,j)^(old)+k2T_(i,j+1)^(old)+k2T_(i,j-1)^(old) ]` …(4)
Gauss Seidel:
`T_(i,j)=(1/((1+2k1+2k2) ))[T_(i,j)^(prev)+K1T_(i-1,j)+k1T_(i+1,j)^(old)+k2T_(i,j-1)+k2T_(i,j+1)^(old) ]` …..(5)
Successive over relaxation:
`T_(i,j)=(1-ω) T_(i,j)^(old)+ω[(1/((1+2k1+2k2) ))[T_(i,j)^(prev)+K1T_(i-1,j)+k1T_(i+1,j)^(old)+k2T_(i,j-1)+k2T_(i,j+1)^(old) ]]`….(6)
Code explanation:
- While creating code for various types and methods one thing is common which is its inputs like boundary conditions, tolerance, grid etc.
- So, we first define all those inputs in code where we take nx and ny as grids in x and y direction respectively (which should be equal according to given condition as we are considering the square plate)
- Next is tolerance where our error tolerance should be 1e-4 which means when error comes at this point solution is said to converged.
- Then we create the array at x and y direction which is having number of values equal to nx and ny respectivey. And like wise we create the dx and dy also using difference between to values of x and y respectively
- Now its time to define the boundary conditions which we have to create using the references given from question. which was our initial boundary condition now lets see how we created the condition;
- First of all we create the array of nx*ny with all ones in array.
- Next we multiply its whole first row by 900 (using T(1,:))
- Then similarly we multiply the last row by 600 (using T(end,:))
- Likewise we do it for left and right 800 and 400 respectively (using T(:,1) and T(:,end) respectively)
- And for corners we select those points and take average of two overlapped values meeting at that particular point for ex. We take T(1,1) which means 1st point of array where 400 and 900 column intersecting each other so we take 900+400/2 which gives us 650.
- As of now the initial boundary condition is defined, now it is a time where various methods and types have different-different equations and for these equation there are different for loops so till this step the code is same for all types and methods now let’s see state vise changes done in code
- Steady state 2D heat conduction:
Here we use implicit scheme which means there are iterative solvers. Where the format of “nested-for-loop” is same for all three solvers just inside equation varies as per solver (explanation going from inner loop to outer loop);
In steady state we first start with nested for loop in which first for loop defines the ‘i’ indices and second for loop defines the ‘j’ indices with nx and ny respectively where after defining those 2 indices we entered the equation there which we stated above during derivations which are eqn(1) for Jacobi, eqn(2) for Gauss seidel, eqn(3) for SOR.
The we end those both loops and create the convergence loop which is while loop which makes equation runs till the solution didn’t get converged in which we say while error is greater that tolerance value the for loop should be run, and also we update the values of Told here and calculates error value and also includes all plot commands in this loops.
- Transient state 2D heat conduction:
- Explicit: In explicit there are no iterative solvers just we have to use the normal transient equation which we derived for explicit scheme. And one change we have to o in the above steady state used for loop that is we have to add the new for loop which we can say time loop and then keep all those convergences and nested for loop inside that time loop (for loop is same how we expressed in steady state implicit)
- Implicit: Here all syntax for ‘nested for loop’ is similar we used in steady state implicit we just have to add time loop outside nested loop as this is transient state. And here we use eqn(4) for Jacobi iterative solver, eqn(5) for Gauss seidel iterative solver, and eqn(6) for SOR iterative solver.
- Contour creation is our final task where we create the contour of the heat conduction in convergence loop. We take it in convergence loop because we want to see how heat conduction going to take place there. So we use contour command to create the 2d contour.
Results:
So first we will check the results for the Steady state 2D heat conduction;
Fig 2. Steady state using Jacobi solver
Fig 3. Steady state using Gauss seidel solver
Fig 4. Steady state using SOR solver
Where from above contours we can see above (in title) the number of iterations required for converge the solution in each and every solver
|
Iterative Solver |
Number of iteration takes to converge |
|
Jacobi Iterative solver |
4590 |
|
Gauss seidel iterative solver |
2466 |
|
Successive over relaxation solver |
1709 |
Transient/unsteady state 2D heat conduction
In transient we use two schemes implicit and explicit scheme so we will check one by one;
Fig 5. Transient state using explicit scheme
Fig 6. Transient state using Jacobi iterative solver
Fig 7. Transient state using Gauss seidel iterative solver
Fig 8. Transient state using SOR iterative solver
Where from above contours we can see above (in title) the number of iterations required for converge the solution in each and every solver
|
Iterative Solver |
Number of iterations takes to converge |
|
Jacobi Iterative solver (Implicit) |
12524 |
|
Gauss seidel iterative solver (Implicit) |
8769 |
|
Successive over relaxation solver (Implicit) |
7639 |
|
Explicit scheme |
5757 |
If we see above all three solvers in implicit of transient as well as steady the number of iterations SOR takes is less than gauss seidel which is less than Jacobian solver but if we compare explicit with implicit method then explicit takes too less iterations to converge as compare to all the three iterative solver that we run in implicit
Comparing steady state to transient state, yes steady state converges faster and use less computation power because it is not time dependent and off course cannot be change with respect to time. This method uses mean values for computing output and for it is less computationally expensive and a lot faster than transient state. How ever transient simulation method calculates instantaneous values at each iteration and performing such simulation requires higher computational resources and off course more time much more data has to be produced.
Conclusion:
In the end we can conclude that the steady state solution converges faster than transient or unsteady state while in solvers we can say that the SOR converges faster than rest of two solvers and it is totally depends on relaxation factor, for faster convergence the factor should be higher while we also have to keep eye that should not get diverge.
Transient state implicit code:
clear all
close all
clc
% transient state heat conduction using Implicit scheme
% inputs
% number of grids in x and y direction
nx = 50;
ny = nx;
% tolerance value define
tol = 1e-4;
% creating array
x = linspace(0,1,nx);
y = linspace(0,1,ny);
% distance between grid
dx = 1/(nx-1);
dy = 1/(ny-1);
dt = 1;
% relaxation factor
alpha = 0.0001;
omega = 1.4;
% temperature grids define
T = 298*ones(nx,ny)
T(1,:) = 900
T(end,:) = 600
T(:,1) = 400
T(:,end) = 800
T(1,1) = (900+400)/2
T(1,end) = (900+800)/2
T(end,1) = (400+600)/2
T(end,end) = (800+600)/2
�lculation of temperature using explicit
k1 = (alpha*dt)/(dx^2)
k2 = (alpha*dt)/(dy^2)
CFL = alpha*dt*(1/dx^2+1/dy^2)
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = T;
method =input('Enter the number for the solver method: Jacobi=1 Gauss Seidal=2 SOR=3 Solver type: ')
%% jacobi method
tic;
if method == 1
ct = 1;
for nt = 1:1000
error = 9e9;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
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_prev_dt(i,j)*term1+term2*H+term3*V;
end
end
error = max(max(abs(Told-T)))
Told = T;
ct = ct+1;
end
T_prev_dt = T;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unstedy state heat conduction Method:implicit Method jacobi Number of iterations: %d(Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f', ct-1,nt*dt,time_counter,CFL));
xlabel('xx');
ylabel('yy');
end
end
%% Gauss seidal method
tic;
if method == 2
ct = 1;
for nt = 1:1000
error = 9e9;
while (error>tol)
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)+Told(i+1,j));
V = (T(i,j-1)+Told(i,j+1));
T(i,j) = T_prev_dt(i,j)*term1+term2*H+term3*V;
end
end
error = max(max(abs(Told-T)))
Told = T;
ct = ct+1;
end
T_prev_dt = T;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unstedy state heat conduction Method:implicit Method Gauss seidal Number of iterations: %d(Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f', ct-1,nt*dt,time_counter,CFL));
xlabel('xx');
ylabel('yy');
end
end
%% SOR successive over relsxation
tic;
if method == 3
ct = 1;
for nt = 1:1000
error = 9e9;
while (error>tol)
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)+Told(i+1,j));
V = (T(i,j-1)+Told(i,j+1));
T(i,j) = ((1-omega)*Told(i,j))+(omega*(T_prev_dt(i,j)*term1+term2*H+term3*V));
end
end
error = max(max(abs(Told-T)))
Told = T;
ct = ct+1;
end
T_prev_dt = T;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unstedy state heat conduction Method:implicit Method SOR Number of iterations: %d(Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f', ct-1,nt*dt,time_counter,CFL));
xlabel('xx');
ylabel('yy');
end
end
Transient state Explicit code:
close all
clear all
clc
%transient state 2D heat conduction using explicit scheme
%inputs
%number of grids in x and y direction
nx = 50;
ny = nx;
%olerance value define
tol = 1e-4;
%creating array
x = linspace(0,1,nx);
y = linspace(0,1,ny);
%distance between grid
dx = 1/(nx-1);
dy = 1/(ny-1);
dt = 1;
%relaxation factor
alpha = 0.0001;
%omega = 1.2;
%temperature grids define
T = 298*ones(nx,ny)
T(1,:) = 900
T(end,:) = 600
T(:,1) = 400
T(:,end) = 800
T(1,1) = (900+400)/2
T(1,end) = (900+800)/2
T(end,1) = (400+600)/2
T(end,end) = (800+600)/2
�lculation of temperature using explicit
k1 = (alpha*dt)/(dx^2)
k2 = (alpha*dt)/(dy^2)
�L calculation
CFL = alpha*dt*(1/dx^2+1/dy^2)
[xx,yy] = meshgrid(x,y);
Told = T;
ct = 1;
tic;
for nt = 1:1000
error = 9e9;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = Told(i-1,j)-2*Told(i,j)+Told(i+1,j);
term2 = Told(i,j-1)-2*Told(i,j)+Told(i,j+1);
T(i,j) = Told(i,j)+k1*term1+k2*term2;
end
end
error = max(max(abs(Told-T)))
Told = T;
ct = ct+1;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unstedy state heat conduction Method:Explicit Method Number of iterations: %d(Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f', ct-1,nt*dt,time_counter,CFL));
xlabel('xx');
ylabel('yy');
end
end
steady state code:
clear all
close all
clc
% steady state heat conduction using Implicit scheme
% inputs
% number of grids in x and y direction
nx = 50;
ny = nx;3
% tolerance value define
tol = 1e-4;
% creating array
x = linspace(0,1,nx);
y = linspace(0,1,ny);
% distance between grid
dx = 1/(nx-1);
dy = 1/(ny-1);
dt = 1;
% relaxation factor
alpha = 0.0001;
omega = 1.2;
% temperature grids define
T = 298*ones(nx,ny)
T(1,:) = 900
T(end,:) = 600
T(:,1) = 400
T(:,end) = 800
T(1,1) = (900+400)/2
T(1,end) = (900+800)/2
T(end,1) = (400+600)/2
T(end,end) = (800+600)/2
�lculation of temperature using explicit
k = 2*((dx^2+dy^2)/(dx^2*dy^2))
�L number
CFL = alpha*dt*(1/dx^2+1/dy^2)
% mesh grid
[xx,yy] = meshgrid(x,y);
Told = T;
method = input('Enter the number for the solver method: Jacobi=1 Gauss Seidal=2 SOR=3 Solver type: ')
%% jacobi method
tic;
if method == 1
ct = 1;
error = 9e9;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = ((Told(i-1,j)+Told(i+1,j))/(k*dx^2));
term2 = (Told(i,j-1)+Told(i,j+1))/(k*dy^2);
T(i,j) = term1+term2;
end
end
error = max(max(abs(Told-T)))
Told = T;
ct = ct+1;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D stedy state heat conduction Method:implicit Method jacobi Number of iterations: �omputation Time: %0.2f CFL Number: %0.3f', ct-1,time_counter,CFL));
xlabel('xx');
ylabel('yy');
end
end
%% Gauss seidal method
tic;
if method == 2
ct = 1;
error = 9e9;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = ((T(i-1,j)+Told(i+1,j))/(k*dx^2));
term2 = (T(i,j-1)+Told(i,j+1))/(k*dy^2);
T(i,j) = term1+term2;
end
end
error = max(max(abs(Told-T)))
Told = T;
ct = ct+1;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D stedy state heat conduction Method:implicit Method Gauss seidal Number of iterations: �omputation Time: %0.2f CFL Number: %0.3f', ct-1,time_counter,CFL));
xlabel('xx');
ylabel('yy');
end
end
%% SOR successive over relxation
tic;
if method == 3
ct = 1;
error = 9e9;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = ((T(i-1,j)+Told(i+1,j))/(k*dx^2));
term2 = (T(i,j-1)+Told(i,j+1))/(k*dy^2);
T(i,j) = (1-omega)*Told(i,j)+omega*(term1+term2);
end
end
error = max(max(abs(Told-T)))
Told = T;
ct = ct+1;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D stedy state heat conduction Method:implicit Method SOR Number of iterations: �omputation Time: %0.2f CFL Number: %0.3f', ct-1,time_counter,CFL));
xlabel('xx');
ylabel('yy');
end
end