Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
2-D Heat conduction equation:
`(∂T)/(∂t)= α((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))`
For 2D-steady state:
`((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))=0`
The above equation is of elliptical type so the given point is influenced and will influence the whole space with infinite speed. So, central difference scheme is used.
`(Tl-2Tp+Tr)/(∆x^2 )+ (Tt-2Tp+Tb)/(∆y^2 )=0`
`2Tp/(∆x^2 )+ 2Tp/(∆y^2 )= (Tl+Tr)/(∆x^2 )+ (Tt+Tb)/(∆y^2 )`
`Tp((2(∆x^2+ ∆y^2 ))/(∆x^2.∆y^2 ))= (Tl+Tr)/(∆x^2 )+ (Tt+Tb)/(∆y^2 )`
`Tp= (∆y^2)/(2(∆x^2+∆y^2 )) .(Tl+Tr)+ (∆x^2)/(2(∆x^2+∆y^2 )) .(Tt+Tb)`
Therefore,
`Ti,j=k1.(T_(i-1,j)+T_(i+1,j) )+k2.(T_(i,j-1)+T_(i,j+1) )`
`where, (∆y^2)/(2(∆x^2+∆y^2) ) =k1 and (∆x^2)/(2(∆x^2+∆y^2) ) =k2`
The above steady state equation can be solved using Jacobi, gauss-seidal and successive over relaxation method in this project:
1. Jacobi method:
`T_(i,j)=k1(T_(old(i-1,j)) +T_(old(i+1,j)) )+k2(T_(old(i,j-1)) +T_(old(i,j+1)) )`
2. Gauss-seidal method:
`T_(i,j)=k1(T_((i-1,j) )+T_((i+1,j) ) )+k2(T_((i,j-1) )+T_((i,j+1) ) )`
Here even for upper diagonal matrix old values are taken because those values are not updated when we are in the loop.
3. Successive over relaxation method:
`T_(i,j)= (1-w) T_(old(i,j)) +wT_(i,j)`
Where in the right hand side can be taken from either Jacobi method or Gauss-seidal method. Since gauss-seidal method has faster convergence compared to Jacobi in this project Gauss-seidal is taken.
W is called relaxation factor
w>1 --- over relaxation and w<1 under relaxation also w≠0.
Matlab code:
clear all
close all
clc
%2D steady state heat conduction
nx = 20;
ny = 20;
dx = 1/(nx-1);
dy = 1/(ny-1);
[x,y] = meshgrid(1:1:nx,1:1:ny); %creating a mesh grid
%entering the boundary conditions and forming a temperatue matrix
T = 300*ones(nx,ny);
T(1,1) = 500;
T(1,ny) = 700;
T(nx,1) = 650;
T(nx,ny) = 850;
T(2:nx-1,1) = 400;
T(2:nx-1,ny) = 800;
T(1,2:ny-1) = 600;
T(nx,2:ny-1) = 900;
method = input('press 1 for jacobi, 2 for gauss-seidal and 3 for succesive over relaxation: ');
T_old = T;
tol = 1e-4;
error = 9e9;
k1 = (dy^2)/(2*(dx^2 +dy^2));
k2 = (dx^2)/(2*(dx^2 +dy^2));
%solving using jacobi method
if method == 1
iter_jacobi = 1;
tic;
while error > tol
for i = 2 : nx-1
for j = 2 : ny-1
T(i,j) = k1*(T_old(i-1,j) + T_old(i+1,j)) + k2*(T_old(i,j-1) + T_old(i,j+1));
end
end
error = max((max(abs(T_old - T))));
T_old = T;
iter_jacobi = iter_jacobi +1; % records no of iterations
end
jacobi_time = toc; %records computation time
[c,h] = contourf(x,y,T); % generate filled contour plot specifying x,y coordinates
clabel(c,h);
colormap(jet);
colorbar;
set(gca,'YDIR','reverse');
xlabel('x axis (L = 1m)');
ylabel('y axis (L = 1m)');
title({'STEADY STATE SOLOUTION USING JACOBI METHOD ', ['No of Iterations = ',num2str(iter_jacobi)], ['Computation Time = ',num2str(jacobi_time)]},'color','b')
end
%solving using gauss-seidal method
if method == 2
iter_gs = 1;
tic;
while error > tol
for i = 2 : nx-1
for j = 2 : ny-1
T(i,j) = k1*(T(i-1,j) + T(i+1,j)) + k2*(T(i,j-1) + T(i,j+1));
end
end
error = max((max(abs(T_old - T))));
T_old = T;
iter_gs = iter_gs +1;
end
gauss_time = toc;
[c,h] = contourf(x,y,T);
clabel(c,h);
colormap(jet);
colorbar;
set(gca,'YDIR','reverse');
xlabel('x axis (L = 1m)');
ylabel('y axis (L = 1m)');
title({'STEADY STATE SOLUTION USING GAUSS-SEIDAL METHOD ', ['No of Iterations = ',num2str(iter_gs)], ['Computation Time = ',num2str(gauss_time)]},'color','b')
end
%solving using successive-over-relaxation method
%taking w = 1.5
if method == 3
iter_sor = 1;
w = 1.2;
tic;
while error > tol
for i = 2 : nx-1
for j = 2 : ny-1
T(i,j) = k1*(T(i-1,j) + T(i+1,j)) + k2*(T(i,j-1) + T(i,j+1));
T(i,j) = (1-w)*T_old(i,j) + w*T(i,j);
end
end
error = max((max(abs(T_old - T))));
T_old = T;
iter_sor = iter_sor +1;
end
sor_time = toc;
[c,h] = contourf(x,y,T);
clabel(c,h,'fontsize',10);
colormap(jet);
colorbar;
set(gca,'YDIR','reverse');
xlabel('x axis (L = 1m)');
ylabel('y axis (L = 1m)');
title({'STEADY STATE SOLOUTION USING SOR METHOD ', ['No of Iterations = ',num2str(iter_sor)], ['Computation Time = ',num2str(sor_time)]},'color','b')
end
Run the code in matlab and press 1 for jacobi method and press enter:- This solves steady state problem using jacobi method until error>tol iteration takes place and temperature matrix gets updated each time. Final temperature profile is shown as 2D-contour plot along with number of iterations and computational time.
Run the code in matlab and press 2 for Gauss-seidal method and press enter:-This solves steady state problem using Gauss-seidal method.
Run the code in matlab and press 3 for SOR method and press enter:- This solves steady problem using SOR method.
Code explaination:
Step 1 : To generate initial Temperature matrix including all the boundary conditions.
Step 2 : Updating the temperature matrix by running through the loop solving discretized equation each time for a selected method.
step 3 : plotting the temperature in 2D-contour plot along with Temperature values, Number of iterations and computational time.
Results:
For number of grid points 20X20 and using Jacobi method contour plot is found to be:
For number of grid points 20X20 and using Gauss-seidal method contour plot is found to be:
For number of grid points 20X20 and using successive over relaxation method with relaxation factor(w=1.2) the contour plot is found to be:
comparison of results:
|
Method |
Number of iterations |
Computation time(s) |
|
Jacobi |
826 |
0.034061 |
|
Gauss-seidal |
442 |
0.020432 |
|
Successive-over-relaxation |
304 |
0.018716 |
Conclusion:
From the above results it can be concluded that successive over relaxation method is faster compared to gauss-seidal method and Jacobi method.
Number of iterations for Jacobi > Gauss-seidal > successive over relaxation
Computational time for Jacobi > Gauss-seidal > successive over relaxation
With increase in number of grid points we get fine results and the corner lines in the contour plot converges to a point.
Note: when number of grid points are taken 100x100
The above plot has smooth curves and parallel lines at the corners are converged to a point, this analysis took long time.
Independent of material the above temperature variation is obtained at steady state for a given specific set of boundary conditions.
2D-Transient Heat conduction:
`(∂T)/(∂t)= α((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))`
The above equation on the left hand side is time marching type for which backward difference scheme is used and right hand side has elliptical type for which central difference scheme is used.
This can be solved either by using explicit scheme or by implicit scheme.
Matlab code:
%solving Transient 2D heat conduction equation
clear all
close all
clc
%grid size inputs
nx = 30;
ny = 30;
dx = 1/(nx-1);
dy = 1/(ny-1);
[x,y] = meshgrid(1:1:nx,1:1:ny);
%entering the boundary conditions and initial conditions for forming a temperatue matrix
T = 300*ones(nx,ny);
T(1,1) = 500;
T(1,ny) = 700;
T(nx,1) = 650;
T(nx,ny) = 850;
T(2:nx-1,1) = 400;
T(2:nx-1,ny) = 800;
T(1,2:ny-1) = 600;
T(nx,2:ny-1) = 900;
dx = 1/(nx-1);
dy = 1/(ny-1);
method = input('press 1 for explicit scheme and 2 for implicit scheme: ');
if method == 1
T_prev = T;
alpha = 1.4;
dt = 2.12e-4;
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);
CFl = k1+k2; % CFL must be less than 0.5
iter_exp = 0;
tic;
for t = 1:1700
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
T_prev = T;
iter_exp = iter_exp +1;
end
exp_time = toc;
[C,H] = contourf(x,y,T);
clabel(C,H,'fontsize',10);
colormap(jet);
colorbar;
set(gca,'YDIR','reverse');
xlabel('xaxis')
ylabel('yaxis')
title({'Transient 2d heat conduction Explicit scheme',['No of Iterations = ',num2str(iter_exp)],['Computational Time = ',num2str(exp_time)]},'color','b');
end
if method == 2
% implicit scheme
T_prev = T;
alpha = 1.4;
dt = 2.12e-4;
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);
cfl = k1 + k2;
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
tol = 1e-4;
T_old = T;
t_max = 1700;
procedure = input('press 1 for jacobi, 2 for gauss-seidal and 3 for SOR: ');
if procedure == 1
iter_jac = 0;
tic;
%jacobi
for t = 1 :t_max
error = 9e9;
while error > tol
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = term1*T_prev(i,j) + term2*(T_old(i-1,j) +T_old(i+1,j)) +term3*(T_old(i,j-1) + T_old(i,j+1));
end
end
error = max(max(abs(T_old - T)));
T_old = T;
e(i) = error;
iter_jac = iter_jac + 1;
end
T_prev = T;
end
jac_time = toc;
[c,h] = contourf(x,y,T);
clabel(c,h);
colormap(jet);
colorbar;
set(gca,'YDIR','reverse');
xlabel('x axis');
ylabel('y axis');
title({'Transient 2d heat conduction Implicit scheme using jacobi ',['No of Iterations = ',num2str(iter_jac)],['Computational Time = ',num2str(jac_time)]},'color','b');
end
if procedure == 2
%gauss-seidal
iter_gauss = 0;
tic;
for t = 1 : t_max
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
error = 9e9;
while error > tol
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = term1*T_prev(i,j) + term2*(T(i-1,j) +T(i+1,j)) +term3*(T(i,j-1) + T(i,j+1));
end
end
error = max(max(abs(T_old - T)));
T_old = T;
iter_gauss = iter_gauss + 1;
end
T_prev = T;
end
gauss_time = toc;
[c,h] = contourf(x,y,T);
clabel(c,h);
colormap(jet);
colorbar;
set(gca,'YDIR','reverse');
xlabel('x axis');
ylabel('y axis');
title({'Transient 2d heat conduction Implicit scheme using Gauss seidal ',['No of Iterations = ',num2str(iter_gauss)],['Computational Time = ',num2str(gauss_time)]},'color','b');
end
if procedure == 3
%successive over relaxation
w = 1.2;
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
iter_sor = 0;
tic;
for t = 1 : t_max
error = 9e9;
while error > tol
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (1-w)*T_old(i,j) + w*(term1*T_prev(i,j) + term2*(T(i-1,j) +T(i+1,j)) +term3*(T(i,j-1) + T(i,j+1)));
end
end
error = max(max(abs(T_old - T)));
T_old = T;
iter_sor = iter_sor + 1;
end
T_prev = T;
end
sor_time = toc;
[c,h] = contourf(x,y,T);
clabel(c,h);
colormap(jet);
colorbar;
set(gca,'YDIR','reverse');
xlabel('x axis');
ylabel('y axis');
title({'Transient 2d heat conduction Implicit scheme using SOR ',['No of Iterations = ',num2str(iter_sor)],['Computational Time = ',num2str(sor_time)]},'color','b');
end
end
Run the above code in matlab press 1 to solve the 2d transient heat conduction using explicit scheme and press enter
Run the above code in matlab again and press 2 to solve using implicit scheme and press enter, now select which 1 or 2 or 3 for solving implicit scheme using jacobi or gauss-seidal or SOR methods
Note: Code is written in this way to make it more interactive.
Code explaination:
step1: Generating initial temperature matrix using all the boundary conditions and initial assumptions and corner point are given average temperature values of edge temperatures.
step2: writing code for explicit scheme which takes initial temperature matrix and runs it in a loop until error>tol.
step3: plotting 2D-contour plot for explicit scheme along with number of iterations and computational time.
step4: Writing code for implicit scheme which has jacobi method, Gauss-seidal method and successive over relaxation method and repeat the same process.
Explicit scheme:
`(∂T)/(∂t)= α((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))`
on discretization
`(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 ))^(n)`
`T_(i,j)^(n+1)= T_(i,j)^n+(∝∆t)/((∆x)^2 ) (T_(i-1,j)-2T_(i,j)+T_(i+1,j) )^(n)+ (∝∆t)/((∆y)^2 ) (T_(i,j-1)-2T_(i,j)+T_(i,j+1) )^(n)`
`Let k1 = (∝∆t)/((∆x)^2 ) and k2 = (∝∆t)/((∆y)^2 ) ~~ CFL = k1 +k2`
`T_(i,j)^(n+1)=(1-2k1-2k2) T_(i,j)^n+k1(T_(i-1,j)+T_(i+1,j) )^(n)+k2(T_(i,j-1)+T_(i,j+1) )^(n)`
In case of explicit scheme CFL must be less than 0.5
Result:
On solving the Transient 2d heat conduction equation for CFL = 0.4992 for 1700 iteration result is found to be:
Time step = 2.12e-4
Alpha = 1.4
CFL = 0.4992
Number of iterations = 1700
Computation time = 0.05631 s
Implicit scheme:-
`(∂T)/(∂t)= α((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))`
on discretization
`(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 ))^(n+1)`
`Let k1 = (∝∆t)/((∆x)^2 ) and k2 = (∝∆t)/((∆y)^2 )`
`T_(i,j)^(n+1)=T_(i,j)^n+k1(T_(i-1,j)-2T_(i,j)+T_(i+1,j) )^(n+1)+k2(T_(i,j-1)-2T_(i,j)+T_(i,j+1) )^(n+1)`
`(1+2k1+2k2)T_(i,j)^(n+1)=T_(i,j)^n+k1(T_(i-1,j)+T_(i+1,j) )^(n)+k2(T_(i,j-1)+T_(i,j+1) )^(n)`
`T_(i,j)^(n+1)=1/((1+2k1+2k2) )( T_(i,j)^n+k1(T_(i-1,j)+T_(i+1,j) )^(n)+k2(T_(i,j-1)+T_(i,j+1) )^(n))`
`Let Term1 = 1/((1+2k1+2k2) )`
Term2 = k1 x Term1
Term3 = K2 x Term1
The above equation can be solved using iterative solvers like Jacobi, gauss-seidal and successive over relaxation.
Jacobi:
`T_(i,j)=Term1.T_(i,j)^prev+Term2.(T_(old(i-1,j))+T_(old(i+1,j)) )+Term3.(T_(old(i,j-1))+T_(old(i,j+1)) )`
Gauss-seidal:
`T_(i,j)=Term1.T_(i,j)^prev+Term2.(T_((i-1,j))+T_((i+1,j)) )+Term3.(T_((i,j-1))+T_((i,j+1)) )`
Successive over relaxation:
`T_(i,j)=(1-w) T_(i,j)^prev+w(Term1.T_(i,j)^prev+Term2.(T_((i-1,j) )+T_((i+1,j) ) )+Term3.(T_((i,j-1) )+T_((i,j+1) ) ))`
Results:
For Transient heat conduction implicit scheme using Jacobi method
Time step = 2.12e-4
Alpha = 1.4
CFL = 0.4992
Number of iterations: 14407
Computational time = 0.37344 s
For Transient heat conduction implicit scheme using Gauss-seidal method:
Time step = 2.12e-4
Alpha = 1.4
CFL = 0.4992
Number of iterations: 1407
Computational time = 0.30161 s
For Transient heat conduction implicit scheme using Successive over relaxation method:
Time step = 2.12e-4
Alpha = 1.4
CFL = 0.4992
Number of iterations = 7383
Computational time = 0.23456 s
Conclusions:
- Among implicit methods number of iterations and computational time for Jacobi method > Gauss-seidal method > Successive-over-relaxation method
- For same number of time iterations explicit method is much faster compared to all implicit methods with less number of iterations and less computational time.
- Thus successive over relaxation method is fastest among implicit methods when relaxation factor is chosen carefully.
- Implicit schemes are unconditionally stable and will also work for higher value of CFL.
- Explicit method is used to know the time accurate transient result when the desired end result might not be steady.
- Implicit method is suitable when steady state is desired end result but otherwise for large time steps there might be errors.