Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

Aim :- Steady state analysis & Transient State Analysis.

Objective :- Solve the 2D heat conduction equation by using the point iterative techniques.

The Boundary conditions,

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

Theory :-

Solving 2D heat conduction equation of steady state :

`(delT)/(delt)-alpha((del^2T)/(delx^2)+(del^2T)/(dely^2)) = 0`

since, the temperature doesn't change with time the above equation reduces to,

`((del^2T)/(delx^2)+(del^2T)/(dely^2)) = 0`

Steps,

Diecretization,

`(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(delx^2)+(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(dely^2)= 0`

As we don't know value sat some nodes like T( i+1) and T(j+1) we take guess values and we get a set of coupled equations and are to be solved at once, as when there are guess values we use iterative solvers and those values converge under solver and we get a converged solution.we will say that our solution is converged when new iteration values and old values are same.

Here we use three type of iterative solvers ,

1) Jacobian Iterative Solver :

In jacobi method , each iteration is done using a set of previous values, even if new values are available in the middle of an iteration.Like when you are solving for n+1 iteration and for i th node you will use values of previous(n) iteration even if you have the updated values of the nodes in current iteration.so this method is called method of simultaneous displacements.

The 2D steady state heat conduction equation for this method is,

`T_(i,j)^(n+1) = 1/k((T_(i-1,j)^n +T_(i+1,j)^n)/(Deltax^2)+(T_(i,j-1)^n +T_(i,j+1)^n)/(Deltay^2))`

2) Gauss-Seidel iterative solver :

In Gauss-Seidel, as soon as you find a new value in the current iteration you will use it for finding other values in the current iteration. Like when you are solving for n+1 iteration and for i th node you will use all the available values upto i-1 node of current iteration(n+1) and use i+1 node values and above values of previous iteration(n).SO we get more converged values.so this method is called method of spontaneous displacements.

The 2D steady state heat conduction equation for this method is ,

`T_(i,j)^(n+1) = 1/k((T_(i-1,j)^n +T_(i+1,j)^n)/(Deltax^2)+(T_(i,j-1)^n +T_(i,j+1)^n)/(Deltay^2))`

3)Successive Over-Relaxation( SOR) method :

It is a variant of the Gauss-Seidel method for solving a linear system of equations, resulting in faster convergence.You can also use jacobi method .Here an extra term called relaxation factor (alpha)is multiplied to gauss seidel generated term to scale the convergence rate and solve the equation faster.To calculate new value we take old value and take a parameter omega and take rate of change.

The 2D steady state heat conduction equation for this method is

for scaling jacobian method,where Tj is the result of jacobian method.

`T_(i,j)^(n+1) = T_(i,j)^n(1-alpha)+alpha.(Tj)`

and for scaling Gauss-Seidel method,where Tgs is the result of Gauss-Seidel method.

`T_(i,j)^(n+1) = T_(i,j)^n(1-alpha)+alpha.(Tgs)`

when alpha = 0 we get new value is equal to old value which gives wrong result.

when alpha goes from 0 to 1 we are switching from no update to an update using gauss siedal or jacobi method.

when alpha >1 then we are taking the update from gauss siedal method and scaling up the convergence rate. The value we give is called over relaxition factor.

when alpha <1 ,The value we give is called under relaxition factor.This is used when the solution given by gauss siedel is over relaxed we give omega value less than 1 to get the converged solution.

MATLAB Code :-

% Solving the steady state 2D heat conduction problem using iterative solvers over a square domain.
clear all
close all
clc
%intial Conditions
L=1;
nx=20;
ny=nx;
tolerance=1e-4;
% TEMPERATURE matrix

error=999;
solver_number= input('enter_the_solver_number: ')

if solver_number==1
%function for steady_state_jacobian_iteration_method
steady_jacobian_iteration_method(L,nx,ny,error,tolerance)
end

if solver_number==2
% function for steady_state_gauss_seidel_iteration_method
steady_gauss_seidel_iteration_method(L,nx,ny,error,tolerance)
end

if solver_number==3
% function for steady_state_SOR_method applied on gauss seidel method
steady_gauss_SOR_method(L,nx,ny,error,tolerance)
end

Function for Jacobi interation method,

% Function to solve 2D steady state heat conduction equation using jacobi iterative solver.


function steady_jacobian_iteration_method(L,nx,ny,error,tolerance)

tic
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=L/(nx-1);
dy=L/(ny-1);
% TEMPERATURE matrix
T= 300*ones(nx,ny);

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T(:,1)=T_L;
T(:,end)=T_R;
T(1,:)=T_T;
T(end,:)=T_B;
%edge refining

T(1,1) = (T_L+T_T)/2;
T(1,end) = (T_T+T_R)/2;
T(end,1) = (T_L+T_B)/2;
T(end,end) = (T_B+T_R)/2;
k=(2*(dx^2+dy^2))/(dx^2*dy^2);
Told=T;
j_iteration_number=1;
%convergence loop
while (error>tolerance)

%nodal loop
for i= 2: (nx-1)
for j= 2:(ny-1)

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=max(max(abs(Told-T))); 
%updating old values
Told=T;
j_iteration_number= j_iteration_number+1;
%creating a plot
figure(1)
[C,h]=contourf(x,fliplr(y),T);
clabel(C,h);
colormap(jet)
colorbar;
time_required=toc
title_text = sprintf('solving 2D steady state heat equation using jacobi iterative solver \n j_iteration number=%d ; computation time: %0.3f s', j_iteration_number,time_required),
title(title_text);
xlabel('X axis')
ylabel('Y axis')
end
end

Function for Gauss Seidel interation method,

% Function to solve 2D steady state heat conduction equation using gauss seidel iterative solver.

function steady_gauss_seidel_iteration_method(L,nx,ny,error,tolerance)

tic
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=L/(nx-1);
dy=L/(ny-1);
T= 300*ones(nx,ny);

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T(:,1)=T_L;
T(:,end)=T_R;
T(1,:)=T_T;
T(end,:)=T_B;
%edge refining

T(1,1) = (T_L+T_T)/2;
T(1,end) = (T_T+T_R)/2;
T(end,1) = (T_L+T_B)/2;
T(end,end) = (T_B+T_R)/2;
k=(2*(dx^2+dy^2))/(dx^2*dy^2);
Told=T;

gs_iteration_number=1;
while (error>tolerance)

for i= 2: (nx-1)
for j= 2:(ny-1)

T(i,j)=(1/k)*((T(i-1,j)+Told(i+1,j))/dx^2 + (T(i,j-1)+Told(i,j+1))/(dy^2));

end
end
error=max(max(abs(Told-T)));
Told=T;
gs_iteration_number= gs_iteration_number+1;
figure(1)
[C,h]=contourf(x,fliplr(y),T);
clabel(C,h)
colormap(jet)
colorbar;
time_required=toc;
title_text = sprintf('solving 2D steady state heat equation using gauss seidel iterative solver \n iteration number=%d ; computation time: %0.1f s', gs_iteration_number,time_required)
title(title_text);
xlabel('X axis')
ylabel('Y axis')
end
end

Function for Succesive Over Relaxation interation method,

% Function to solve 2D steady state heat conduction equation using SOR iterative solver.

function steady_gauss_SOR_method(L,nx,ny,error,tolerance)

tic
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=L/(nx-1);
dy=L/(ny-1);
T= 300*ones(nx,ny);

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T(:,1)=T_L;
T(:,end)=T_R;
T(1,:)=T_T;
T(end,:)=T_B;
%edge refining

T(1,1) = (T_L+T_T)/2;
T(1,end) = (T_T+T_R)/2;
T(end,1) = (T_L+T_B)/2;
T(end,end) = (T_B+T_R)/2;
k=(2*(dx^2+dy^2))/(dx^2*dy^2);
Told=T;
alpha=1.72;
sor_iteration_number=1;
while (error>tolerance)

for i= 2: (nx-1)
for j= 2:(ny-1)

Tgs(i,j)=(1/k)*((T(i-1,j)+Told(i+1,j))/dx^2 + (T(i,j-1)+Told(i,j+1))/(dy^2));
T(i,j)=(1-alpha)*(Told(i,j))+alpha*(Tgs(i,j));
end
end
error=max(max(abs(Told-T)));
Told=T;
sor_iteration_number= sor_iteration_number+1;
figure(1)
[C,h]=contourf(x,fliplr(y),T);
clabel(C,h)
colormap(jet)
colorbar;
time_required=toc;
title_text = sprintf('solving 2D steady state heat equation using SOR iterative solver \niteration number=%d ; computation time: %0.1f s', sor_iteration_number,time_required);
title(title_text);
xlabel('X axis')
ylabel('Y axis')
end
end

Output Of Each Solver,

Jacobi Method

Gauss Seidel Method

SOR Method

In the above three plot results you can see that the number of iterations and the computation time for jacobi iteration method is high i.e,826 iterations and 111.894 seconds time and then it is for gauss seidel method with 442 iterations and 56.3 seconds time and then for SOR.On addition of over relaxation factor to the gauss seidel method i.e, SOR method you can see that the solution is calculated in just 59 iterations and 8.6 seconds time.The SOR method depends on relaxation factor and the solution is optimum with less number of iterations and computational time when the alpha value is 1.72 and the value above or below this the number iterations gets increased and when the value is above 1.99 the solution get blow up.

Solving 2D heat conduction equation of transient state :

The equation for 2D heat conduction equation of transient state is,

`(delT)/(delt)-alpha((del^2T)/(delx^2)+(del^2T)/(dely^2)) = 0`

where alpha is the thermal diffusivity.

Then discretize the equation by using forward difference scheme for time and central difference for space.

1). If we consider RHS terms of n th time step while calculating for (n+1) time step then we call it as explicit form.

`(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat) = alpha((T_(i-1,j)^n +T_(i+1,j)^n)/(Deltax^2)+(T_(i,j-1)^n +T_(i,j+1)^n)/(Deltay^2))`

As it is of explict form it easy to solve these equations as we use previous time step values which are already known . we solve these equations directly by keeping it in the time loop and space loop.

MATLAB Code :-

% Solving the transient state 2D heat conduction problem using explicit method over a square domain.

clear all 
close all
clc
tic
%intial Conditions
L=1;
nx=20;
ny=nx;
% TEMPERATURE matrix
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=L/(nx-1);
dy=L/(ny-1);
% TEMPERATURE matrix
T= 300*ones(nx,ny);

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T(:,1)=T_L;
T(:,end)=T_R;
T(1,:)=T_T;
T(end,:)=T_B;
%edge refining 

T(1,1) = (T_L+T_T)/2;
T(1,end) = (T_T+T_R)/2;
T(end,1) = (T_L+T_B)/2;
T(end,end) = (T_B+T_R)/2;
Tprev=T;
dt=1e-4;
t=0.15;
nt=t/dt;
alpha=1;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
 
for k=1:nt
    for i=2:nx-1
    for j=2:ny-1
   
    T(i,j)=Tprev(i,j)*(1-2*k1-2*k2)+k1*(Tprev(i-1,j)+Tprev(i+1,j))+k2*(Tprev(i,j-1)+Tprev(i,j+1));
    end
    end
    Tprev=T;
       figure(1)
[C,h]=contourf(x,fliplr(y),T);
clabel(C,h)
colormap(jet)
colorbar;
time_of_simulation=toc;
title_text = sprintf('solving 2D transient state heat equation explicitly \n time of simulation=%0.3f',time_of_simulation);
title(title_text)
xlabel('X axis')
ylabel('Y axis')

end

Output :-

Here we get the solution when it satisfy the CFL number condition.

2). when we consider RHS terms of (n+1)th time step while calculating for (n+1) time step then we call it as implicit form.

`(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat) = alpha((T_(i-1,j)^(n+1) -2T_(i,j)^(n+1)+T_(i+1,j)^(n+1))/(Deltax^2)+(T_(i,j-1)^n -2T_(i,j)^(n+1)+T_(i,j+1)^(n+1))/(Deltay^2))`

now to solve these equations as we get set of coupled equations using direct solvers it is difficult, we use iterative solvers.

MATLAB Code :-

% Solving the transient state 2D heat conduction problem using iterative solvers over a square domain.

clear all 
close all
clc
tic
%intial Conditions
L=1;
nx=20;
ny=nx;
tolerance=1e-4;

solver_number= input('enter_the_solver_number: ')

if solver_number==1
    %function for transient_state_jacobian_iteration_method
transient_jacobian_iteration_method(L,nx,ny,tolerance)
end

if solver_number==2
%    function for transient_state_gauss_seidel_iteration_method
transient_gauss_seidel_iteration_method(L,nx,ny,tolerance)
end

if solver_number==3
%   function for  transient SOR_iteration_method  applied on gauss seidel method
     transient_gauss_SOR_iteration_method(L,nx,ny,tolerance)
end

Function for Transient State Jacobi method,

% Function to solve 2D transient state heat conduction equation using jacobi iterative solver

function transient_jacobian_iteration_method(L,nx,ny,tolerance)
tic
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=L/(nx-1);
dy=L/(ny-1);
% TEMPERATURE matrix
T= 300*ones(nx,ny);

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T(:,1)=T_L;
T(:,end)=T_R;
T(1,:)=T_T;
T(end,:)=T_B;
%edge refining 

T(1,1) = (T_L+T_T)/2;
T(1,end) = (T_T+T_R)/2;
T(end,1) = (T_L+T_B)/2;
T(end,end) = (T_B+T_R)/2;
Tprev=T;
Told =T;
dt=1e-4;
t=0.15;
nt=t/dt;
alpha=1;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
 j_iteration_number=1
%time loop
for k=1:nt
    error =999
% convergence loop
    while(error>tolerance)
% nodal loop
    for i=2:nx-1
    for j=2:ny-1
        
         T(i,j)=(Tprev(i,j)+k1*(Told(i-1,j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1)))/(1+2*k1+2*k2);
         
    end
    end
    error=max(max(abs(T-Told)));
%updating old values
    Told=T;
    j_iteration_number= j_iteration_number+1
    end
%updating previous values
    Tprev=T;
    figure(1)
[C,h]=contourf(x,fliplr(y),T);
clabel(C,h);
colormap(jet)
colorbar;
time_required=toc
title_text = sprintf('solving 2D transient state heat equation using jacobi iterative solver \n j_iteration number=%d ; computation time: %0.3f s', j_iteration_number,time_required),
title(title_text);
xlabel('X axis')
ylabel('Y axis')
end
end

Function for Transient State Gauss Seidel Method

% Function to solve 2D transient state heat conduction equation using gauss seidel iterative solver
function transient_gauss_seidel_iteration_method(L,nx,ny,tolerance)
tic
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=L/(nx-1);
dy=L/(ny-1);
% TEMPERATURE matrix
T= 300*ones(nx,ny);

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T(:,1)=T_L;
T(:,end)=T_R;
T(1,:)=T_T;
T(end,:)=T_B;
%edge refining 

T(1,1) = (T_L+T_T)/2;
T(1,end) = (T_T+T_R)/2;
T(end,1) = (T_L+T_B)/2;
T(end,end) = (T_B+T_R)/2;
Tprev=T;
Told =T;
dt=1e-4;
t=0.25;
nt=t/dt;
alpha=1;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
 gs_iteration_number=1
%  timeloop
for k=1:nt 
    error =999
%     convergence loop
    while(error>tolerance)
%         nodal loop
    for i=2:nx-1
    for j=2:ny-1
        
         T(i,j)=(Tprev(i,j)+k1*(T(i-1,j)+Told(i+1,j))+k2*(T(i,j-1)+Told(i,j+1)))/(1+2*k1+2*k2);
         
    end
    end
    error=max(max(abs(T-Told)));
%     upgating old values
    Told=T;
    gs_iteration_number= gs_iteration_number+1;
    end
%     updating previous values
    Tprev=T;
    figure(1)
%     creating the plot
[C,h]=contourf(x,fliplr(y),T);
clabel(C,h);
colormap(jet)
colorbar;
time_required=toc
title_text = sprintf('solving 2D transient state heat equation using gauss seidel solver \n gs_iteration number=%d ; computation time: %0.3f s', gs_iteration_number,time_required),
title(title_text);
xlabel('X axis')
ylabel('Y axis')
end
end

Function for Transient State SOR Method

% Function to solve 2D transient state heat conduction equation using SOR solver
function transient_gauss_SOR_iteration_method(L,nx,ny,tolerance)
tic
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=L/(nx-1);
dy=L/(ny-1);
% TEMPERATURE matrix
T= 300*ones(nx,ny);

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T(:,1)=T_L;
T(:,end)=T_R;
T(1,:)=T_T;
T(end,:)=T_B;
%edge refining 

T(1,1) = (T_L+T_T)/2;
T(1,end) = (T_T+T_R)/2;
T(end,1) = (T_L+T_B)/2;
T(end,end) = (T_B+T_R)/2;
Tprev=T;
Told =T;
dt=1e-4;
t=0.15;
nt=t/dt;
alpha=1;
% Tgs=ones(nx,ny)
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
 SOR_iteration_number=1
omega=1.06
for k=1:nt
    error =999;
    while(error>tolerance)
    for i=2:nx-1
    for j=2:ny-1
        
         Tgs(i,j)=(Tprev(i,j)+k1*(T(i-1,j)+Told(i+1,j))+k2*(T(i,j-1)+Told(i,j+1)))/(1+2*k1+2*k2);
         T(i,j)=(1-omega)*Told(i,j)+omega*Tgs(i,j);
         
    end
    end
    error=max(max(abs(T-Told)));
    Told=T;
    SOR_iteration_number= SOR_iteration_number+1;
    end
    Tprev=T;
   figure(1)
[C,h]=contourf(x,fliplr(y),T);
clabel(C,h);
colormap(jet)
colorbar;
time_required=toc;
title_text = sprintf('solving 2D transient state heat equation using SOR iterative solver \n SOR_iteration number=%d ; computation time: %0.3f s', SOR_iteration_number,time_required);
title(title_text);
xlabel('X axis')
ylabel('Y axis')
end
end

Output :-

Jacobi Method

Gauss Seidel Method

SOR Method

From the above plots you can conclude that explicit method takes less time to execute but when it comes to stability it mainly depend on CFL condition if this condition is not satisfied the solution gets blow up.Where as in the implicit methods jacobian method gives the solution in 7970 iterations where as gauss method takes 10115 iterations but using SOR method you get the solution in 5098 iterations. SOR method is is faster than other two.Here in SOR method if you use alpha value more or less than 1.06 you can see the iteration count increases when it crosses 1.9 then you see that the solution got blew up.so its value should be given less than 2.you can also observe that for any values of alpha dx and dt the solution will be stable, like explicit it is not dependent of CFL condition.