BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

KKasinisha K7 Aug 2022
MATLAB

Solving the steady and unsteady 2D heat conduction problem

Aim:

To solve Steady and Unsteady state 2D Heat Conduction Equation by using the point iterative techniques such as Jacobi, Gauss Seidel & Successive over-relaxation for both implicit and explicit schemes using MATLAB.

Given Data and assumptions for solving the problem:

1. Assume that the domain is a unit square.

2. Assume nx = ny=10(for steady ) and 41(for unsteady) [Number of points along the x-direction is equal to the number of points along the y direction]

3. Boundary conditions for steady and transient case

1. Left boundary:400K

2. Top boundary:600K

3. Right boundary:800K

4. Bottom boundary:900K

4.Tolerance=1e-4

5.Starting error=9e9

6.Thermal diffusivity(Alpha)=1.88e-5 m^2/sec

7.Doamin length=1m

8.CFL(for stability)=alpha*(dt/(dx^2))

THEORY:

The heat equation is a partial differential equation that describes how the distribution of some quantity (such as heat) evolves over time in a solid medium, as it spontaneously flows from places where it is higher towards places where it is lower.

For a function of three spatial variables (see Cartesian coordinate system) and the time variable , the heat equation is

where,

is a real coefficient called the diffusivity of the medium.

Types of Heat Conduction schemes

This transfer of Heat energy may occur under two conditions

  • Steady-state condition
  • Unsteady/Transient state condition

Unsteady state conditions are a precursor to steady state conditions. No system exists initially under steady state conditions. Some time must pass, after heat transfer is initiated, before the system reaches steady state. During that period of transition the system is under unsteady state conditions.

Time Loop Logic:

Boundary Condition Map:

Governing Equations

In order to solve the 2D Heat Conduction problem, the governing equations for respective conditions with respect to cartesian coordinates should be discretized to enter that in a code which is given under MAIN PROGRAM below.

MAIN PROGRAM:

Below there is a sequence of seven individual programs, starting with IMPLICIT Scheme of Steady & Unsteady-state conditions for three iterative methods Jacobi, Gauss-siedel & SOR followed by EXPLICIT Scheme for Unsteady/Transient Conditions.

IMPLICIT Scheme - Steady-state Condition Program:

Steady-state Equation discretization:

In steady-state conduction, the amount of heat entering any region of an object is equal to the amount of heat coming out (if this were not so, the temperature would be rising or falling, as thermal energy was tapped or trapped in a region).

`(∂ ^2 T)/ (∂ x^ 2) +( ∂^ 2 T)/( ∂ y^ 2) = 0`

By dioscretizing the equation above it becomes,

`T_p=1/k[(T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)]`

Hera the value of k is,

`K=2((Deltax^2+Deltay^2)/(Deltax^2Deltay^2))`

if we further simply this the equation for code becomes

`0.25*(T_L+T_R+T_T+T_B)`

Matlab code:

Jacobi Iteration

close all
clear all
clc

%Solving the Steady state 2D heat conduction (Jacobi Method)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1)  = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1)  = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_intial = T;
T_old = T;


%Calculation of 2D steady heat conduction EQUATION by JACOBI method
iterative_solver = 1;

if iterative_solver == 1
    jacobi_iteration = 1;
    
    while(error > tolerance)
        
        for i = 2:nx - 1
            
            for j = 2:ny - 1
                T(i,j) = 0.25*(T_old(i-1,j) + T_old(i+1,j) + T_old(i,j-1) + T_old(i,j+1));
            end
            
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        jacobi_iteration = jacobi_iteration + 1;  
        
    end
end

%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')

xlabel('X Axis')
ylabel('Y Axis')
title(sprintf('No. of Jacobi Iterations(IMPLICIT) = %d', jacobi_iteration));
pause(0.03);

Gauss-Siedel Iteration

close all
clear all
clc

%Solving the Steady state 2D heat conduction (Gauss-seidel Method)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1)  = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1)  = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_intial = T;
T_old = T;


%Calculation of 2D steady heat conduction EQUATION by Gauss-Seidel method
iterative_solver = 2;

if iterative_solver == 2
    Gauss_seidel_iteration = 1;
    
    while(error > tolerance)
        
        for i = 2:nx - 1
            
            for j = 2:ny - 1
                T(i,j) = 0.25*(T(i-1,j) + T_old(i+1,j) + T(i,j-1) + T_old(i,j+1));
            end
            
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        Gauss_seidel_iteration = Gauss_seidel_iteration + 1;  
        
    end
end

%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')

xlabel('X Axis')
ylabel('Y Axis')
title(sprintf('No. of Gauss-Seidel Iterations(IMPLICIT) = %d', Gauss_seidel_iteration));
pause(0.03);

Successive over relaxation

close all
clear all
clc

%Solving the Steady state 2D heat conduction (Successive over-relaxation Method)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
omega = 1.2;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1)  = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1)  = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_intial = T;
T_old = T;


%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
iterative_solver = 3;

if iterative_solver == 3
    Successive_Over_Relaxation_iteration = 1;
    
    while(error > tolerance)
        
        for i = 2:nx - 1
            
            for j = 2:ny - 1
                T(i,j) = omega*(0.25*(T(i-1,j) + T_old(i+1,j) + T(i,j-1) + T_old(i,j+1))) + (1 - omega)*T_old(i,j);
            end
            
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        Successive_Over_Relaxation_iteration = Successive_Over_Relaxation_iteration + 1;  
        
    end
end

%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')

xlabel('X Axis')
ylabel('Y Axis')
title(sprintf('No. of Successive Over-Relaxation Iterations(IMPLICIT) = %d', Successive_Over_Relaxation_iteration));
pause(0.03);

Output for IMPLICIT Scheme - Steady-state Condition

Jacobi Iteration

Gauss-Siedel Iteration

Successive over-relaxation

IMPLICIT Scheme Unsteady/Transient-state Condition Program:

Unsteady-state Equation discretization:

Non-steady-state situations appear after an imposed change in temperature at a boundary as well as inside of an object, as a result of a new source or sink of heat suddenly introduced within an object, causing temperatures near the source or sink to change in time.

`(∂ T)/( ∂ t) + α ( ∂^ 2 T)/( ∂ x^ 2) +( ∂^ 2 T)/( ∂ y^ 2 ) = 0`

if we discretize the equation it becomes,

`(T_P^(n+1)-T_p^n)/(Deltat)=alpha(T_L-2T_P+T_R)/(∂x^2)+(T_T-2T_P+T_B)/(0y^2)`

if we simplify it further,

`T_P^(n+1)=1/(1+k1+k2)(T_p^n+k1(T_L+T_R)^(n+1)+k2(T_T+T_B)^(n+1))`

`k1=alpha(Deltat)/(Deltax^2)k2=alpha(Deltat)/(Deltay^2)`

Alpha is diffusivity teram.

Jacobi Iteration

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Jacobi Method(IMPLICIT SCHEME)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
dt = 1e-3;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1)  = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1)  = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_old = T;
T_intial = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

Jacobi_iteration = 1;
for k = 1:nt
    error = 9e9; 
    
    while(error > tolerance)
        
        for i = 2:nx - 1
            
            for j = 2:ny - 1
                    
                Term_1 = 1/(1 + (2*k1) + (2*k2));
                Term_2 = k1*Term_1;
                Term_3 = k2*Term_1;
                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_intial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);
   
            end
            
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        Jacobi_iteration = Jacobi_iteration + 1;  
        
    end
    
    T_intial = T;
    
    %Plotting
    figure(1)
    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca, 'ydir', 'reverse')

    xlabel('X-Axis')
    ylabel('Y-Axis')
    title(sprintf('No. of Unsteady Jacobi Iterations(EXPLICIT) = %d', Jacobi_iteration));
    pause(0.03)
end

Gauss-Siedel Iteration

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Gauss Seidel Method(IMPLICIT SCHEME)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
dt = 1e-3;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1)  = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1)  = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_old = T;
T_intial = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

Gauss_Seidel_iteration = 1;
for k = 1:nt
    error = 9e9; 
    
    while(error > tolerance)
        
        for i = 2:nx - 1
            
            for j = 2:ny - 1
                    
                Term_1 = 1/(1 + (2*k1) + (2*k2));
                Term_2 = k1*Term_1;
                Term_3 = k2*Term_1;
                H = (T(i-1,j)) + (T_old(i+1,j));
                V = (T(i,j+1)) + (T_old(i,j-1));
                    
                T(i,j) = (T_intial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);
   
            end
            
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        Gauss_Seidel_iteration = Gauss_Seidel_iteration + 1;  
        
    end
    
    T_intial = T;
    
    %Plotting
    figure(1)
    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca, 'ydir', 'reverse')

    xlabel('X-Axis')
    ylabel('Y-Axis')
    title(sprintf('No. of Unsteady Gauss-Seidel Iterations(IMPLICIT) = %d', Gauss_Seidel_iteration));
    pause(0.03)
end

Successive over-relaxation

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Successive over-relaxation Method(IMPLICIT SCHEME)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
dt = 1e-3;
omega = 1.1;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1)  = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1)  = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_old = T;
T_intial = T;
T_end = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

Successive_over_relaxation_iteration = 1;
for k = 1:nt
    error = 9e9; 
    
    while(error > tolerance)
        
        for i = 2:nx - 1
            
            for j = 2:ny - 1
                    
                Term_1 = 1/(1 + (2*k1) + (2*k2));
                Term_2 = k1*Term_1;
                Term_3 = k2*Term_1;
                H = (T(i-1,j)) + (T_old(i+1,j));
                V = (T(i,j+1)) + (T_old(i,j-1));
                    
                T(i,j) = (T_intial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);
                T_end(i,j) = ((1-omega)*T_old(i,j)) + (T(i,j)*omega);
                
            end
            
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T_end;
        Successive_over_relaxation_iteration = Successive_over_relaxation_iteration + 1;  
        
    end
    
    T_intial = T_end;
    
    %Plotting
    figure(1)
    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca, 'ydir', 'reverse')

    xlabel('X-Axis')
    ylabel('Y-Axis')
    title(sprintf('No. of Unsteady SOR Iterations(IMPLICIT) = %d', Successive_over_relaxation_iteration));
    pause(0.03)
end

Output for IMPLICIT Scheme Unsteady/Transient-state Condition

Jacobi Iterations

Youtube video link:

https://youtu.be/uD0yAyo4Q3U

Gauss-Siedel Iteration

Youtube video link:

https://youtu.be/g6tfDddGLFk

Successive over-relaxation

Youtube video link:

https://youtu.be/B-C5RVvGS8I

EXPLICIT Scheme - Unsteady/Transient state Condition Program

Unsteady-state Equation discretization(EXPLICIT):

`(∂T)/(∂t)=alpha(∂^2T)/(∂x^2)+(∂^2T)/(∂y^2)`

if we discretize the equation it becomes,

`(T_p^(n+1)-T_p^n)/(Deltat)=alpha(T_L-2T_P+T_R)/(∂x^2)+(T_T-2T_P+T_B)/(∂y^2)`

if we simplify it further,

`T_P^(n+1)=T_P^n+k1(T_L-2T_P+T_R)^n+k2(T_T-2T_P+T_B)^n`

Explicit Scheme:

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction (EXPLICIT SCHEME)

%Input Parameters
%Number of grid points
nx = 10;
ny = nx;
nt = 1400;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

%Absolute error criteria > tolerance
error = 9e9;
tolerance = 1e-4;
dt = 1e-3;
omega = 1.1;

%Defining Boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx,ny);

T(2:ny - 1, 1)  = T_L;
T(2:ny - 1, nx) = T_R;
T(1, 2:nx - 1)  = T_T;
T(ny, 2:nx - 1) = T_B;

%Calculating Average temperature at corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

%Assigning orginal values to T
T_intial = T;
T_old = T;


%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
iterative_solver = 1;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

if iterative_solver == 1
    Steady_Explicit = 1;
    for k = 1:nt
        error = 1; 
    
        while(error > tolerance)
        
            for i = 2:nx - 1
            
                for j = 2:ny - 1
                    Term_1 = T_old(i,j);
                    Term_2 = k1*(T_old(i+1,j) - 2*T_old(i,j) + (T_old(i-1,j)));
                    Term_3 = k2*(T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1));
                    T(i,j) = Term_1 + Term_2 + Term_3;
                end
            
            end
        
            error = max(max(abs(T_old - T)));
            T_old = T;
            Steady_Explicit = Steady_Explicit + 1;  
        
        end
    end
end

%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')

xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('No. of Unsteady Iterations(EXPLICIT) = %d', Steady_Explicit));
pause(0.03);

Output for Explicit Scheme

Results:

IMPLICIT-STEADY STATE ANALYSIS:

JACOBI 208

GAUSS-SIEDEL 112

SUCCESSIVE OVER RELAXATION 46

IMPLICIT-TRANSIENT STATE ANALYSIS:

JACOBI 3880

GAUSS-SIEDEL 3612

SUCCESSIVE OVER RELAXATION 3290

EXPLICIT ANALYSIS: 1941

From the above outputs of various methods & comparing their No. of iterations, we can conclude that

  1. Jacobi > Gauss-seidel > Successive over-relaxation according to the No. of iterations computed for all methods & Schemes.
  2. When comparing Steady state and Unsteady state, it is proven that Steady state is the most successful one as it takes least number of iterations and least time to compute.
  3. SOR iteration method is the most effective iteration method as it took least No. of iterations for coverging solutions in both the schemes (EXPLICIT & IMPLICIT) compared to Jacobi & Gauss-siedel.
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.