Solving Steady and Unsteady 2D Heat Conduction Problem

Aim : To solve the 2D Heat Conduction equation using a Steady state solver and Transient solver using Iterative techniques like,

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Absolute error criteria is 1e-4

Assume that the domain is a unit square.

Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]

Boundary conditions for steady and transient case

1. Left:400K

2. Top:600K

3. Right:800K

4. Bottom:900K

The 2D Heat Conduction equation is given below:

where, first term is time derivative, second term is spatial derivative and `α`is thermal diffusivity.

Derivations for Steady State :-

A State of Fluid in which there is no change in the fluid properties with respect to time is known as Steady State.

So that,

We have to solve 2D Steady Heat Conduction equation using three iterative methods as following:

Jacobi Method

Gauss-Seidel Method

Successive Over Relaxation Method (SOR Method)

We have to discretize the steady heat conduction equation using Central Differencing Scheme.

After shuffling, we obtained,

Assume,

So that,

Steady State 2D Heat Conduction equation using Jacobi Mathod is given below:

In Jacobi Method, each iteration is done using a set of previous values, even if new values become available in the middle of an iteration.

Steady State 2D Heat Conduction equation using Gauss-Seidel is given below:

In Gauss-Seidel Method, as soon as a new iterate for a particular component is obtained, it will take it for all subsequent computation.

Steady State 2D Heat Conduction equation using Successive Over Relaxation is given below:

Where, `λ`is the Over Relaxation Factor.

Main MATLAB code for Steady State Analysis :

clear all
close all
clc

%Unit Square Length of Domain
L=1;
%Number of Grid Points along X-direction and Y-direction
nx =10;
ny = 10;
%Discretized the Domain in X-direction and Y- direction
x = linspace(0,L,nx);
dx = x(2)-x(1);       %OR L/(nx-1);
y = linspace(0,L,ny);
dy = y(2)-y(1);       %OR L/(ny-1);
%Error
error = 9e9;
%Tolerance
tol = 1e-4;

k = 2*((1/dx^2)+(1/dy^2));

%Initilized Temperature Array
T = ones(nx,ny);

%Defining Boundary Conditions
TL = 400;         %Left
TR = 800;         %Right
TT = 600;         %Top  
TB = 900;         %Bottom
%Applying Boundary Conditions
T(:,1)= TL;       %Left
T(:,nx)= TR;      %Right
T(1,:)= TT;       %Top
T(ny,:)= TB;      %Bottom
%Temperature at the Corners
T(1,1)=(TL+TT)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)= (TT+TR)/2;
T(nx,1)= (TL+TB)/2;

%Steady State Analysis
for iter_method= [1 2 3]
    if iter_method ==1
        Time1 = Steady_Jacobi(nx,ny,x,y,dx,dy,T,tol,k);
    elseif iter_method ==2
        Time2 = Steady_GaussSeidel(nx,ny,x,y,dx,dy,T,tol,k);
    elseif iter_method ==3
        Time3 = Steady_SOR(nx,ny,x,y,dx,dy,T,tol,k);  
    end
end

MATLAB Function for Jacobi Method :

function [out] = Steady_Jacobi(nx,ny,x,y,dx,dy,T,tol,k)

meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
Jacobi_Iterations =1;
error=9e9;

tic
%Converge Loop
while error>tol
    %Nodal Loop
    for i=2:nx-1
        for j=2:ny-1
            term1 = (T_old(i-1,j)+T_old(i+1,j))/(k*dx^2);
            term2 = (T_old(i,j-1)+T_old(i,j+1))/(k*dy^2);
            T(i,j)= term1+term2;
        end
    end 
    error = max(max(abs(T_old-T)));
    %Update Previous Timestep Value
    T_old = T;
    %Update Iterations
    Jacobi_Iterations = Jacobi_Iterations+1;
end
iter_time = toc;

figure(1)
[a,b]=contourf(x,y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Jacobi Method');
text2 = sprintf('Number of Iterations = %d ',Jacobi_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Steady State Analysis';text1;text2;text3;text4});

out= iter_time;
end

MATLAB Function for Gauss Seidel Method :

function [out] = Steady_GaussSeidel(nx,ny,x,y,dx,dy,T,tol,k)

meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
GS_Iterations =1;
error=9e9;

tic
%Converge Loop
while error>tol
    %Nodal Loop
    for i=2:nx-1
        for j=2:ny-1
            term1 = (T(i-1,j)+T_old(i+1,j))/(k*dx^2);
            term2 = (T(i,j-1)+T_old(i,j+1))/(k*dy^2);
            T(i,j)= term1+term2;
        end
    end 
    error = max(max(abs(T_old-T)));
    %Update Previous Timestep Value
    T_old = T;
    %Update Iterations
    GS_Iterations = GS_Iterations+1;
end
iter_time = toc;

figure(2)
[a,b]=contourf(x,y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Gauss Seidel Method');
text2 = sprintf('Number of Iterations = %d ',GS_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Steady State Analysis';text1;text2;text3;text4});

out= iter_time;
end

MATLAB Function for Successive Over Relaxation Method :

function [out] = Steady_SOR(nx,ny,x,y,dx,dy,T,tol,k)

meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
SOR_Iterations =1;
error=9e9;
%Over Relaxation Factor
Lambda = 1.49;

tic
%Converge Loop
while error>tol
    %Nodal Loop
    for i=2:nx-1
        for j=2:ny-1
            term1 = (T(i-1,j)+T_old(i+1,j))/(k*dx^2);
            term2 = (T(i,j-1)+T_old(i,j+1))/(k*dy^2);
            T(i,j)= (1-Lambda)*T_old(i,j)+(term1+term2)*Lambda;
        end
    end 
    error = max(max(abs(T_old-T)));
    %Update Previous Timestep Value
    T_old = T;
    %Update Iterations
    SOR_Iterations = SOR_Iterations+1;
end
iter_time = toc;

figure(3)
[a,b]=contourf(x,y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Successive Over Relaxation Method');
text2 = sprintf('Number of Iterations = %d ',SOR_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Steady State Analysis';text1;text2;text3;text4})'

out= iter_time;
end

Code Explanation :

• Defined different parameters in the main MATLAB code.
• Initialized the Temperature array.
• Applied boundary conditions at edges and corners.
• Called the different defined functions in the main program.
• In functions, "T" values are assigned to "T_old" for the calculation.
• Initialized the iterations.
• Created a Time loop, Converge loop, Nodal loop in such a way that it can calculate the temperatue distribution in the given domain for the given timesteps.
• Appropriate equations are used for that in context of the type of analysis.
• "tic" & "toc" is used to get the simulation time.
• Plotted a calculated data to get the worthful interpretation.
• Matlab command like "contourf","set(gca,'Ydir','reverse')","colormap(jet)","colorbar","clabel","xlabel", "ylabel","title" are used for that.

Result of Steady State Analysis :

From the above plots, we can see that number of iteration in the jacobi method is higher than the gauss seidel method and SOR method. This is same for the simulation time and the error. Therefore we can say that SOR method for Steady State analysis is better than gauss seidel and jacobi method as it is faster and accurate.

Derivations for Transient State :-

A State of Fluid in which there is a change in the fluid properties with respect to time is known as Unsteady State or Transient State.

Transient 2D Heat Conduction equation is given below:

Discretizing the time derivative using Forward Differencing Scheme and spatial derivative using Central Differencing Scheme.

Furthermore Simplifying,

Assume,

We obtained,

After arranging,

Implicit Form

Implicit functions are given in terms of both dependent and independent variables.

Here, `T_P^n` does not change with every iteration while, (n+1) terms change.

Implicit form of Transient Heat Conduction equation for Jacobi Method :

Implicit form of Transient Heat Conduction equation for Gauss Seidel Method :

Implicit form of Transient Heat Conduction equation for SOR Method :

Explicit Form

Explicit functions are given in terms of independent variables.

Explicit form of Transient Heat Conduction equation for Jacobi Method :

Explicit form of Transient Heat Conduction equation for Gauss Seidel Method :

Explicit form of Transient Heat Conduction equation for SOR Method :

Main MATLAB code for Transient State Analysis :

clear all
close all
clc

%Unit Square Length of Domain
L=1;
%Number of Grid Points along X-direction and Y-direction
nx =10;
ny = 10;
%Discretized the Domain in X-direction and Y- direction
x = linspace(0,L,nx);
dx = x(2)-x(1);       %OR L/(nx-1);
y = linspace(0,L,ny);
dy = y(2)-y(1);       %OR L/(ny-1);
t =10;
dt =1e-2; 
nt=t/dt;
%Error
error = 9e9;
%Tolerance
tol = 1e-4;
%Thermal Diffusivity
Alpha = 0.3;
%Assumption for Unsteady Constant (CFL number)
k1 = Alpha*(dt/dx^2);
k2 = Alpha*(dt/dy^2);

%Initilized Temperature Array
T = ones(nx,ny);

%Defining Boundary Conditions
TL = 400;         %Left
TR = 800;         %Right
TT = 600;         %Top  
TB = 900;         %Bottom
%Applying Boundary Conditions
T(:,1)= TL;       %Left
T(:,nx)= TR;      %Right
T(1,:)= TT;       %Top
T(ny,:)= TB;      %Bottom
%Temperature at the Corners
T(1,1)=(TL+TT)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)= (TT+TR)/2;
T(nx,1)= (TL+TB)/2;



%Transient State Analysis
%Implicit Approach
for iter_method= [1 2 3]
    if iter_method ==1
        Time1 = Transient_Implicit_Jacobi(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2);
    elseif iter_method ==2
        Time2 = Transient_Implicit_GaussSeidel(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2);
    elseif iter_method ==3
        Time3 = Transient_Implicit_SOR(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2);  
    end
end

%Explicit Approach
for iter_method= [1 2 3]
    if iter_method ==1
        Time4 = Transient_Explicit_Jacobi(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2);
    elseif iter_method ==2
        Time5 = Transient_Explicit_GaussSeidel(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2);
    elseif iter_method ==3
        Time6 = Transient_Explicit_SOR(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2);  
    end
end

Implicit Approach :

MATLAB Function for Jacobi Method :

function [out] = Transient_Implicit_Jacobi(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2)

[X,Y]=meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
T_prev_dt = T;

nt = t/dt;
I_Jacobi_Iterations =1;

tic
%Time Loop
for m = 1:nt
    error=9e9;
    %Converge Loop
    while error>tol
         %Nodal Loop
         for i=2:nx-1
             for j=2:ny-1
                 term1 = 1/(1+(2*k1)+(2*k2));
                 term2 = k1*term1;
                 term3 = k2*term1;
                 H = (T_old(i-1,j)+T_old(i+1,j));
                 V = (T_old(i,j-1)+T_old(i,j+1));
                 T(i,j)= (T_prev_dt(i,j)*term1)+ H*term2 + V*term3;
             end
         end 
         error = max(max(abs(T_old-T)));
         %Update old Value
         T_old = T;
         %Update Iterations
         I_Jacobi_Iterations = I_Jacobi_Iterations+1;
    end
    %Update Previous Timestep Value
    T_prev_dt = T;
end
iter_time = toc;

figure(1)
[a,b]=contourf(X,Y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Implicit Jacobi Method');
text2 = sprintf('Number of Iterations = %d ',I_Jacobi_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Transient State Analysis';text1;text2;text3;text4});

out= iter_time;
end

MATLAB Function for Gauss-Seidel Method :

function [out] = Transient_Implicit_GaussSeidel(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2)

[X,Y]=meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
T_prev_dt = T;

nt= t/dt;
I_GS_Iterations =1;

tic
%Time Loop
for k=1:nt
    error=9e9;
    %Converge Loop
    while error>tol
         %Nodal Loop
         for i=2:nx-1
             for j=2:ny-1
                 term1 = 1/(1+(2*k1)+(2*k2));
                 term2 = k1*term1;
                 term3 = k2*term1;
                 H = (T(i-1,j)+T_old(i+1,j));
                 V = (T(i,j-1)+T_old(i,j+1));
                 T(i,j)= (T_prev_dt(i,j)*term1)+ H*term2 + V*term3;
             end
         end 
         error = max(max(abs(T_old-T)));
         %Update old Value
         T_old = T;
         %Update Iterations
         I_GS_Iterations = I_GS_Iterations + 1;
    end
    %Update Previous Timestep Value
    T_prev_dt = T;
end
iter_time = toc;

figure(2)
[a,b]=contourf(X,Y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Implicit Gauss Seidel Method');
text2 = sprintf('Number of Iterations = %d ',I_GS_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Transient State Analysis';text1;text2;text3;text4});

out= iter_time;
end

MATLAB Function for Successive Over Relaxation Method :

function [out] = Transient_Implicit_SOR(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2)

[X,Y]=meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
T_prev_dt = T;

Lambda = 1.2;
nt= t/dt;
I_SOR_Iterations =1;


tic
%Time Loop
for k=1:nt
    error=9e9;
    %Converge Loop
    while (error>tol)
         %Nodal Loop
          for i=2:nx-1
             for j=2:ny-1
                 term1 = (1-Lambda)*T_old(i,j);
                 term2 = Lambda/(1+(2*k1)+(2*k2));     
                 H = (T(i-1,j)+T_old(i+1,j));
                 V = (T(i,j-1)+T_old(i,j+1));
                 term3 = T_prev_dt(i,j)+ (k1*H) + (k2*V); 
                 T(i,j)= term1 + (term2*term3);
             end
          end 
          error = max(max(abs(T_old-T)));
          %Update old Value
          T_old = T;
          %Update Iterations
          I_SOR_Iterations = I_SOR_Iterations + 1;
    end
    %Update Previous Timestep Value
    T_prev_dt = T;
end
iter_time = toc;

figure(3)
[a,b]=contourf(X,Y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Implicit SOR Method');
text2 = sprintf('Number of Iterations = %d ',I_SOR_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Transient State Analysis';text1;text2;text3;text4});

out= iter_time;
end

Result of Transient State Analysis with Implicit Approach :

Explicit Approach :

MATLAB Function for Jacobi Method :

function [out] = Transient_Explicit_Jacobi(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2)

[X,Y]= meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
T_prev_dt = T;

nt= t/dt;
E_Jacobi_Iterations =1;

tic
%Time Loop
for k=1:nt
    error=9e9;
    %Converge Loop
    while error>tol
         %Nodal Loop
         for i=2:nx-1
             for j=2:ny-1
                 term1 = (1-(2*k1)-(2*k2));
                 term2 = k1;
                 term3 = k2;
                 H = (T_old(i-1,j)+T_old(i+1,j));
                 V = (T_old(i,j-1)+T_old(i,j+1));
                 T(i,j)= (T_prev_dt(i,j)*term1)+ H*term2 + V*term3;
             end
         end 
         error = max(max(abs(T_old-T)));
         %Update old Value
         T_old = T;
         %Update Iterations
         E_Jacobi_Iterations = E_Jacobi_Iterations+1;
    end
    %Update Previous Timestep Value
    T_prev_dt = T;
end
iter_time = toc;

figure(4)
[a,b]=contourf(X,Y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Explicit Jacobi Method');
text2 = sprintf('Number of Iterations = %d ',E_Jacobi_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Transient State Analysis';text1;text2;text3;text4});

out= iter_time;
end

MATLAB Function for Gauss-Seidel Method :

function [out] = Transient_Explicit_GaussSeidel(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2)

[X,Y]=meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
T_prev_dt = T;

nt= t/dt;
E_GS_Iterations =1;

tic
%Time Loop
for k=1:nt
    error=9e9;
    %Converge Loop
    while error>tol
         %Nodal Loop
         for i=2:nx-1
             for j=2:ny-1
                 term1 = (1-(2*k1)-(2*k2));
                 term2 = k1;
                 term3 = k2;
                 H = (T(i-1,j)+T_old(i+1,j));
                 V = (T(i,j-1)+T_old(i,j+1));
                 T(i,j)= (T_prev_dt(i,j)*term1)+ H*term2 + V*term3;
             end
         end 
         error = max(max(abs(T_old-T)));
         %Update old Value
         T_old = T;
         %Update Iterations
         E_GS_Iterations = E_GS_Iterations+1;
    end
    %Update Previous Timestep Value
    T_prev_dt = T;
end
iter_time = toc;

figure(5)
[a,b]=contourf(X,Y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Explicit GaussSeidel Method');
text2 = sprintf('Number of Iterations = %d ',E_GS_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Transient State Analysis';text1;text2;text3;text4});

out= iter_time;
end

MATLAB Function for Successive Over Relaxation Method :

function [out] = Transient_Explicit_SOR(nx,ny,x,y,dx,dy,T,t,dt,tol,k1,k2)

[X,Y]=meshgrid(x,y);
%Assigning Temperature Array
T_old = T;
T_prev_dt = T;

Lambda = 1.3;
nt= t/dt;
E_SOR_Iterations =1;

tic
%Time Loop
for k=1:nt
    error=9e9;
    %Converge Loop
    while error>tol
         %Nodal Loop
         for i=2:nx-1
             for j=2:ny-1
                 term1 = (1-Lambda)*T_old(i,j);
                 term2 = Lambda*(T_prev_dt(i,j)*(1-2*k1-2*k2));
                 H = (T(i-1,j)+T_old(i+1,j));
                 V = (T(i,j-1)+T_old(i,j+1));
                 %term3 = (Lambda*k1*H) + (Lambda*k2*V)
                 T(i,j)= term1 + term2 + (Lambda*k1*H) + (Lambda*k2*V) ;
             end
         end 
         error = max(max(abs(T_old-T)));
         %Update Previous old Value
         T_old = T;
         %Update Iterations
         E_SOR_Iterations = E_SOR_Iterations+1;
    end
    %Update Previous Timestep Value
    T_prev_dt = T;
end
iter_time = toc;

figure(6)
[a,b]=contourf(X,Y,T);
set(gca,'Ydir','reverse')
colormap(jet);
colorbar;
clabel(a,b);
xlabel('Domain Length in X-direction');
ylabel('Domain Length in Y-direction');
text1 = sprintf('Explicit SOR Method');
text2 = sprintf('Number of Iterations = %d ',E_SOR_Iterations);
text3 = sprintf('Simulation Time = %f s',iter_time);
text4 = sprintf('Final Error = %f',error);
title({'Transient State Analysis';text1;text2;text3;text4});

out= iter_time;
end

Result of Transient State Analysis with Explicit Approach :

Conclusion :

As we can see, Successive Over Relaxation method has the best result in terms of accuracy, simulation time and number of iterations.

Explicit approach is much faster compared to Implicit approach with less amount of simulation time and number of iterations

Accuracy of Implicit approach is much higher compared to Explicit approach.