2D Heat Conduction Problem - Midterm Project
We are solving 2d heat conduction equation in this project. We are doing it in Steady state and Transient state conditions. The iterative techniques include Gauss Siedal, Jacobi, SOR methods.
Steady State 2D Heat conduction
Main Equations
We know 2D heat conduction equation is given by
`alpha((del^2T)/(delx^2) + (del^2T)/(dely^2)) = (delT)/(delt)`
Descritising the above equation we get
`(del^2T)/(delx^2) = ((T_(i-1,j) - 2T_(i,j) + T_(i+1,j))/(Deltax^2) )`
`(del^2T)/(dely^2) = ((T_(i,j-1) - 2T_(i,j) + T_(i,j+1))/(Deltay^2) )`
Then the equation we get is
`(T_(i-1,j) - 2T_(i,j) + T_(i+1,j))/(Deltax^2) + (T_(i,j-1) - 2T_(i,j) + T_(i,j+1))/(Deltay^2) = 0`
Here we are assuming nx = nx and therefore `Deltax = Delta y`
therefore
`T_(i,j) =1/4 (T_(i-1,j) + T_(i+1,j) + T_(i,j-1) + T_(i,j+1))`
Main code **(Please run the codes in MATLAB)
clear all
close all
clc
%inputs
nx = 10;
ny = 10;
nx = ny;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
t = ones(nx,ny);
%Assigning BC
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
t(1,2:end-1) = 600;
t(end,2:end-1) = 900;
t(2:end-1,1) = 400;
t(2:end-1,end) = 800;
t(1,1) = ((T_top + T_left)/2);
t(end,1) = ((T_left + T_bottom)/2);
t(1,end) = ((T_top + T_right)/2);
t(end,end) = ((T_right + T_bottom)/2);
iterative_solvers = 3;
error = 9e9;
tol = 1e-4
%input solver number
%jacobi = 1
%gauss seidal = 2
%sor = 3
%steady state heat conduction
solver = 1
%Jacobi Method
if solver == 1
t = Steady_jacobi(x,y,t,tol,nx,ny,dx,dy);
%gauss seidel
else if solver == 2
t = Steady_gauss(x,y,t,tol,nx,ny,dx,dy);
%SOR
else if solver == 3
t = Steady_sor(x,y,t,tol,nx,ny,dx,dy);
end
end
endHere we assigned the input values in the above code. We formed a square domain and gave values to it. Later on we specified the boundary condition for the square domain.
Then we specified the number of iterative solver, error and tolerence. we assigned number 1, 2 and 3 for Jacobi, Gauss Seidal and SOR solver respectively. The function files for the iterative methods are given below it with their respective inputs from the main code.
Function file : Steady State Jacobi
%Steady state jacobi solver
function t = Steady_jacobi(x,y,t,tol,nx,ny,dx,dy);
[X,Y] = meshgrid(x,y);
jacobi_iteration = 0;
error = 9e9;
t_old = t;
tic
k = (2/(dx)^2) + (2/(dy)^2);
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
value_1 = (t_old(i-1,j) + t_old(i+1,j))/(k*dx^2);
value_2 = (t_old(i,j-1) + t_old(i,j+1))/(k*dy^2);
t(i,j) = value_1 + value_2;
end
end
error = max(max(abs(t_old - t)));
t_old = t;
jacobi_iteration = jacobi_iteration + 1 ;
end
time = toc
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temperature distribution \\n Jacobi method : %d iterations\', jacobi_iteration))
xlabel(\'X\')
ylabel(\'Y\')
clabel(C,\'FontSize\',10) % for temp values
set(gca, \'YDIR\',\'reverse\')
colormap(jet)
hold on;
%time to solve
display(\'Steady State 2D Heat Conduction(Jacobi method)\');
fprintf(\'Computational Time : %f s\\n\', time);
fprintf(\'Number of iterations : %f\', jacobi_iteration);
end
Here we form a function file, first we specify a meshgrid and give the iteration number and the error. Then we specify the old value t_old and specify the value of k. Now we form a while loop for which the error is greater than the tolerence. Then we specify the equation inside the for loops and later on update the t_old value and iterations value. This is then plotted in a countour graph giving the headings and other important datas.
Output
Computation time: 0.002710 sec
Number of iteration: 216.000
Function file : Steady Gauss Seidal
%Steady state H C Gauss Seidal Solver
function t = Steady_gauss(x,y,t,tol,nx,ny,dx,dy)
%creating mesh
[X,Y] = meshgrid(x,y);
Gauss_iteration = 0;
error = 9e9;
t_old = t;
tic
k = (2/(dx)^2) + (2/(dy)^2);
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
value_1 = (t(i-1,j) + t_old(i+1,j))/(k*dx^2);
value_2 = (t(i,j-1)+ t_old(i,j+1))/(k*dy^2);
t(i,j) = value_1 + value_2;
end
end
error = max(max(abs(t_old - t)));
t_old = t;
Gauss_iteration = Gauss_iteration + 1;
end
time = toc
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temperature distribution \\n Gauss Seidal method : %d iterations\', Gauss_iteration))
xlabel(\'X\')
ylabel(\'Y\')
clabel(C,\'FontSize\',10) % for temp values
set(gca, \'YDIR\',\'reverse\')
colormap(jet)
hold on;
%time to solve
display(\'Steady State 2D Heat Conduction(Gauss Seidal method)\');
fprintf(\'Computational Time : %f s\\n\', time);
fprintf(\'Number of iterations : %f\', Gauss_iteration);
endHere we form a function file, first we specify a meshgrid and give the iteration number and the error. Then we specify the old value t_old and specify the value of k. Now we form a while loop for which the error is greater than the tolerence. Then we specify the equation inside the for loops and later on update the t_old value and iterations value. This is then plotted in a countour graph giving the headings and other important datas.
Output
Computation time: 0.002383 sec
Number of iteration: 116.000
Function file : Steady State SOR
%Steady State Heat Conduction S.O.R Solver
function t = Steady_sor(x,y,t,tol,nx,ny,dx,dy)
[X,Y] = meshgrid(x,y);
SOR_iteration = 0;
error = 9e9;
t_old = t;
tic
k = (2/(dx)^2) + (2/(dy)^2);
%relaxation parameter as from the below website % w = omega
%https://www.sciencedirect.com/science/article/pii/S0893965908001523
w = 2/(1+sin(pi*dx));
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
value_1 = (t(i-1,j) + t_old(i+1,j))/(k*dx^2);
value_2 = (t(i,j-1) + t_old(i,j+1))/(k*dy^2);
t(i,j) = w * (value_1 + value_2) + (1-w)*t_old(i,j);
end
end
error = max(max(abs(t_old - t)));
t_old = t;
SOR_iteration = SOR_iteration + 1 ;
end
time = toc
figure(1)
C = contourf(X,Y,t);
title(sprintf(\'Temperature distribution \\n SOR method : %d iterations\', SOR_iteration))
xlabel(\'X\')
ylabel(\'Y\')
clabel(C,\'FontSize\',10) % for temp values
set(gca, \'YDIR\',\'reverse\')
colormap(jet)
hold on;
%time to solve
display(\'Steady State 2D Heat Conduction(SOR method)\');
fprintf(\'Computational Time : %f s\\n\', time);
fprintf(\'Number of iterations : %f\', SOR_iteration);
endHere we form a function file, first we specify a meshgrid and give the iteration number and the error. Then we specify the old value t_old and specify the value of k. Now we form a while loop for which the error is greater than the tolerence. Then we specify the equation inside the for loops and later on update the t_old value and iterations value. This is then plotted in a countour graph giving the headings and other important datas.
Output
Computation time: 0.002404 sec
Number of iteration: 30.000
Result : Steady State
From the above given outputs we can derive an inference that out of the three iterative methods, SOR gives us the result faster followed by Gauss Seidal and later followed by Jacobi method.
Transient State 2D Heat conduction
Main equations
We have 2D heat conduction equation for unsteady or transient condition as
`alpha((del^2T)/(delx^2) + (del^2T)/(dely^2)) = (delT)/(delt)`
Descritising the spacial terms by central differencing we get
`(del^2T)/(delx^2) = ((T_(i-1,j) - 2T_(i,j) + T_(i+1,j))/(Deltax^2) )`
`(del^2T)/(dely^2) = ((T_(i,j-1) - 2T_(i,j) + T_(i,j+1))/(Deltay^2) )`
Descritising time derivative by forward differencing we get
`(delT)/(delt) = ((T^(n+1))_(i,j) - (T^n)_(i,j))/(Deltat)`
and summing up we get,
`(T_(i-1,j) - 2T_(i,j) + T_(i+1,j))^(n+1)/(Deltax^2) + (T_(i,j-1) - 2T_(i,j) + T_(i,j+1))^(n+1)/(Deltay^2) = ((T^(n+1))_(i,j) - (T^n)_(i,j))/(Deltat)`
Now for implicit we have :
Here we defined 2 constants k1 and k2 such that
`k1 = alpha(Deltat)/(Deltax^2)`
`k2 = alpha(Deltat)/(Deltay^2)`
and we know
`Deltax = Deltay`
finally we get by substitution that
`(T^(n+1))_(i,j) = ((T^n)_(i,j) + (k1(T_(i-1,j) + T_(i+1,j)) + k2(T_(i,j-1) + T_(i,j+1)))^(n+1))/(1+2*k1+2*k2`
Now for Explicit we have :
`alpha((T_(i-1,j)-2*T_(i,j)+T_(i+1,j))^n/(Deltax^2) +(T_(i,j-1)-2*T_(i,j)+T_(i,j+1))^n/(Deltay^2)) = ((T^(n+1))_(i,j) - (T^n)_(i,j))/(Delta*t)`
Rearranging the above equation and substituting `lambda = (alpha*Deltat)/(Deltax^2`, we get:
`(T^(n+1))_(i,j) = lambda((T^n)_(i-1,j) + (T^n)_(i,j-1)) + (1-4*lambda)*(T^n)_(i,j)+lambda*((T^n)_(i+1,j)+(T^n)_(i,j+1)`
The greek letter \'lambda\' is called as the CFL number. This number is used to check the stability of solution
Main code **(Please run the codes in MATLAB)
clear all
close all
clc
%Transient method
%inputs
nx = 10;
ny = 10;
nx = ny;
%omega = thermal diffusivity \'w\'
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
T = 301*ones(nx,ny);% temp as 28 degree celcious(converted into kelvin gives 301k)
%Assigning BC
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
T(1,2:end-1) = 600;
T(end,2:end-1) = 900;
T(2:end-1,1) = 400;
T(2:end-1,end) = 800;
T(1,1) = ((T_top + T_left)/2);
T(end,1) = ((T_left + T_bottom)/2);
T(1,end) = ((T_top + T_right)/2);
T(end,end) = ((T_right + T_bottom)/2);
alpha = 1.5;
tol = 1e-4;
error = 9e9;
t_start = 0;
t_end = 1/(alpha);
dt_1 = 1e-2;%for implicit
dt_2 = 1e-3;%for explicit
t_1 = [t_start:dt_1:t_end];%for implicit
t_2 = [t_start:dt_2:t_end];%for explicit
nt_1 = length(t_1);%for implicit
nt_2 = length(t_2);%for explicit
alpha = 1.5;
k1 = (alpha*dt_1)/dx^2;%for implicit
k2 = (alpha*dt_1)/dy^2;%for implicit
k3 = (alpha*dt_2)/dx^2;%CFL for explicit
k4 = (alpha*dt_2)/dy^2;%CFL for explicit
k = 1/(1+2*k1+2*k2);
T_old = T;
T_prev = T_old;
%SOLVERS
%Explicit = 1
%implicit:
%jacobi = 2
%gauss seidal = 3
%sor = 4
%please enter the solver number
solver = 1;
%Explicit solver
if solver == 1
T = Explicit_transient(T,nt_2,x,y,nx,ny,error,k3,k4);
%implicit solvers
%Jacobi
else if solver == 2
T = Implicit_transient_jacobi(T,x,y,tol,error,nx,ny,k1,k2,nt_1);
%Gauss Seidal
else if solver == 3
T = Implicit_transient_gauss(T,x,y,tol,error,k1,k2,nx,ny,nt_1);
%SOR
else if solver == 4
T = Implicit_transient_SOR(T,x,y,tol,nx,ny,error,k1,k2,nt_1);
end
end
end
endHere we assigned the input values in the above code. We formed a square domain and gave values to it. Later on we specified the boundary condition for the square domain. After that we define the value of alpha, which is required for the calculation of k1, k2, k3 and k4. We specify the tolerence, error and t. The t begins with t_start and ends with t_end with a dt number of divisions or time steps. The dt_1 value is for implicit methods while dt_2 value is for explicit method. Similarly t_1,nt_1 and t_2,nt_2 are for implicit and explicit methods respectively. Then k1, k2, k3 and k4 are calculated from the above input and values of T is updated as T_old and T_prev.
Then we specified the number of iterative solver, we assigned number 1, 2, 3 and 4 for Explicit, Jacobi, Gauss Seidal and SOR solver respectively (out of which Jacobi, Gauss Seidal and SOR solver are coming under implicit scheme) The function files for the iterative methods are given below it with their respective inputs from the main code.
%Explicit solver
function T = Explicit_transient(T,nt_2,x,y,nx,ny,error,k3,k4)
%Creating mesh
[X,Y] = meshgrid(x,y);
%error = 9e9;
%Creating a copy of T
T_old = T;
T_prev = T;
% k3 and k4 are the CFL Numbers
tic
iteration = 1;
%Time loop
for k =1:nt_2;
error = 9e9;
%convergence loop
for i = 2:nx-1 ;%bc
for j = 2:ny-1;
value1 = T_prev(i,j);
value2 = k3*(T_prev(i-1,j) - 2*T_prev(i,j) + T_prev(i+1,j));
value3 = k4 * (T_prev(i,j-1) - 2*T_prev(i,j) + T_prev(i,j+1));
T(i,j) = value1 + value2 + value3;
end
end
error = max(max(abs(T_prev - T)));
T_prev = T;
iteration = iteration + 1;
end
time = toc;
%Plotting and Display Functions
figure(1)
[X,Y] = meshgrid(x,y);
C = contourf(X,Y,T);
set(gca, \'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'FontSize\',10); % for temp values
title(sprintf(\'Temperature distribution \\n Explicit method : %d iterations\',iteration));
xlabel(\'X\');
ylabel(\'Y\');
hold on
%time to solve
display(\'Transient 2D Heat Conduction(Explicit method)\');
fprintf(\'Computational Time : %fs \\n\', time);
fprintf(\'Number of iteration : %d \\n\',iteration);
fprintf( \'Error : %f \\n\', error);
end%Transient Implicit Solvers
%Jacobi Solver
function T = Implicit_transient_jacobi(T,x,y,tol,error,nx,ny,k1,k2,nt_1)
%copy of T for iteration
T_old = T;
T_prev = T_old;
value1 = (1+(2*k1)+(2*k2))^-1;
value2 = k1*value1;
value3 = k2*value1;
tic
jacobi_iteration = 1;
%time loop
for k = 1:nt_1;
error = 9e9;
%convergence loop
while(error>tol)
for i = 2:nx-1 %nodal loop
for j = 2:ny-1
EQ1 = T_old(i-1,j) + T_old(i+1,j);
EQ2 = T_old(i,j-1) + T_old(i,j+1);
T(i,j) = (T_prev(i,j)* value1) + (EQ1*value2) + (EQ2*value3);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
jacobi_iteration = jacobi_iteration + 1 ;
end
T_prev = T_old;
end
time = toc;
%Plotting and Display Functions
figure(1)
[X,Y] = meshgrid(x,y);
C = contourf(X,Y,T);
set(gca, \'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'FontSize\',10); % for temp values
title(sprintf(\'Temperature distribution \\n Jacobi method : %d iterations\', jacobi_iteration));
xlabel(\'X\');
ylabel(\'Y\');
hold on;
%time to solve
display(\'Transient 2D Heat Conduction(Jacobi method)\');
fprintf(\'Computational Time : %f s\\n\', time);
fprintf(\'Number of iterations : %d \\n\', jacobi_iteration);
fprintf(\'Error : %f \\n\', error);
endIn this function file, first we take input values from the main code file. Then we create a mesh. After that we create a copy of T, iteration number is counted and time loop is begun from 1 to nt_1. Then the convergence loop is begun specifying all the above mentioned equations. Also convergence loop is inside a while loop which ensures that error is greater than tolerence The convergence loop is ended and error ,T_old are given out. Then the iteration counter is incremented by one and time loop is also ended. The last time is notted and the plotting and display functions are specified.
Output
Number of iterations : 1367
Error : 0.000052
%Transient Implicit Solvers
%Gauss seidal Solver
function T = Implicit_transient_gauss(T,x,y,tol,error,k1,k2,nx,ny,nt_1)
value1 = (1+(2*k1)+(2*k2))^-1;
value2 = k1*value1;
value3 = k2*value1;
%copy of T for iteration
T_old = T;
T_prev = T_old;
tic
gauss_iteration = 1;
%time loop
for k = 1:nt_1;
error = 9e9;
%convergence loop
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
EQ1 = (T(i-1,j) + T(i+1,j));
EQ2 = (T(i,j-1) + T(i,j+1));
T(i,j) = (T_prev(i,j)*value1)+(EQ1*value2) + (EQ2*value3);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
gauss_iteration = gauss_iteration + 1 ;
end
T_prev = T_old;
end
time = toc;
%Plotting And Display Functions
figure(1)
[X,Y] = meshgrid(x,y);
C = contourf(X,Y,T);
set(gca, \'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'FontSize\',10); % for temp values
title(sprintf(\'Temperature distribution \\n Gauss Seidal method : %d iterations\', gauss_iteration));
xlabel(\'X\');
ylabel(\'Y\');
hold on
%time to solve
display(\'Transient 2D Heat Conduction(Gauss Seidal method)\');
fprintf(\'Computational Time : %fs \\n\', time);
fprintf(\'Number of iteration : %d \\n\',gauss_iteration);
fprintf( \'Error : %f \\n\', error);
endIn this function file, first we take input values from the main code file. Then we create a mesh. After that we create a copy of T, iteration number is counted and time loop is begun from 1 to nt_1. Then the convergence loop is begun specifying all the above mentioned equations. Also convergence loop is inside a while loop which ensures that error is greater than tolerence The convergence loop is ended and error ,T_old are given out. Then the iteration counter is incremented by one and time loop is also ended. The last time is notted and the plotting and display functions are specified.
Output
Computational Time : 0.004772s
Number of iteration : 850
Error : 0.000034
Function file : Transient Implicit SOR Solver
%Transient Implicit Solvers
%SOR Solver
function T = Implicit_transient_SOR(T,x,y,tol,nx,ny,error,k1,k2,nt_1)
%copy of T for iteration
T_old = T;
T_prev = T_old;
z = cos(pi/nx)+cos(pi/ny);
w = 1.2;%OMEGA
value1 = (1+(2*k1)+(2*k2))^-1;
value2 = k1 * value1;
value3 = k2 * value1;
tic
SOR_iteration = 1;
%time loop
for k = 1:nt_1;
error = 9e9;
%convergence loop
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
EQ1 = (T(i-1,j) + T(i+1,j));
EQ2 = (T(i,j-1) + T(i,j+1));
T(i,j) = w*((T_prev(i,j)*value1) + (EQ1*value2) + (EQ2*value3)) -((w-1)*T_prev(i,j));
end
end
error = max(max(abs(T_old - T)));
T_old = T;
end
SOR_iteration = SOR_iteration + 1;
T_prev = T_old;
end
time = toc;
%Plotting and Display Functions
figure(1)
[X,Y] = meshgrid(x,y);
C = contourf(X,Y,T);
set(gca, \'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'FontSize\',10); % for temp values
title(sprintf(\'Temperature distribution \\n SOR method : %d \\n iterations\',SOR_iteration));
xlabel(\'X\');
ylabel(\'Y\');
hold on;
%time to solve
display(\'Transient 2D Heat Conduction(SOR method)\');
fprintf(\'Computational Time : %f s\\n\', time);
fprintf(\'Number of iterations : %d \\n\',SOR_iteration);
fprintf(\'Error : %f\', error);
endIn this function file, first we take input values from the main code file. Then we create a mesh. After that we create a copy of T and we give a value for omega(1.2). Later iteration number is counted and time loop is begun from 1 to nt_1. Then the convergence loop is begun specifying all the above mentioned equations. Also convergence loop is inside a while loop which ensures that error is greater than tolerence The convergence loop is ended and error ,T_old are given out. Then the iteration counter is incremented by one and time loop is also ended. The last time is notted and the plotting and display functions are specified.
Output
Computational Time : 0.003776 s
Number of iterations : 68
Error : 0.000000
Result : Transient State
After the analysis, it was duly notted that the Explicit method was much faster (doing 668 iterations at 0.004574 sec).
Out of the three implicit method, SOR solver was notted to be faster than the other two(doing 68 iterations at 0.003757 sec). It is followed by Gauss Seidal solver and later on followed by Jacobi solver.
Comparing the two methods , implicit method is much more accurate as it uses many solvers rather than just one. Also it uses a wide range of values from the whole domain to converge to an answer.
**(Please run the codes in MATLAB only)