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

Objective:

To solve the 2D heat conduction equation in Steady-state & Transient State by using the point iterative techniques Jacobi, Gauss-seidel & Successive over-relaxation with absolute error criteria is 1e-4.

2D Transient Heat Conduction Equation

Steady-State Analysis:

Deriving the above equation w r t to steady-state, Temperature at point (i,j) can be calculated by the following equation

The steady-state solution is calculated by using the stencil containing points located at Left, Right, Top & Bottom of the preferred point. At what point of calculation, the adjacent values should be taken into consideration for calculation depends on the solver chosen

The program is written to solve this equation using Iterative Solvers and the results are discussed.

This program asks for which iterative method to be used for solving the equation. Please provide the input as listed below

1. Jacobi Method
2. GS Method
3. SOR Method - Requires additional input of w

Main Program to solve the Steady State Equation:

clear all
close all
clc

% Getting Solver Info
fprintf('Steady State Method   ')
fprintf('Enter 1 for Jacobi Method  ')
fprintf('Enter 2 for GS Method  ')
fprintf('Enter 3 for SOR Method  ')
prompt = ('enter iteration method:')
iterative_solver = input(prompt);


% Definig length
L = 1;

% number of grid points
nx = 10;
ny = 10;

% Defining No of Elements
x = linspace(0,L,nx);
y = linspace(0,L,ny);
  
% Mesh Spacing
dx = x(2)-x(1);
dy = y(2)-y(1);

% Initialization
% Assuming the initial values at all points as 300 K
T = 300 * ones(nx,ny);

% Boundary Condition - given value
T_T = 600;
T_L = 400;
T_B = 900;
T_R = 800;

% Assigning BC to Corner Points
T(1,1)   = (T_B + T_L)/2 ;
T(nx,1)  = (T_T + T_L)/2 ;
T(1,ny)  = (T_B + T_R)/2 ;
T(nx,ny)  = (T_T + T_R)/2 ;

% Assigninn BC to line points
T(1,2:nx-1)  = T_B;
T(nx,2:ny-1) = T_T;
T(2:ny-1,1)  = T_L;
T(2:ny-1,nx) = T_R;

% Copying of BC's)
Told= T;

% Calculating the value of k
k = 2 * (dx^2 + dy^2) / (dx^2 * dy^2);

% Convergence Criteria
error = 9e9;
tol = 1e-4;

% Jacobi Method
if iterative_solver == 1
    jacobi_iter = 1;
    while(error > tol)
        
        for i = 2:nx-1
        for j = 2:ny-1
            term1(i,j) = (Told(i-1,j) + Told(i+1,j))/dx^2;
            term2(i,j) = (Told(i,j+1) + Told(i,j-1))/dy^2;
            T(i,j) = (1/k) * (term1(i,j) + term2(i,j));
        end
        end
        error = max(max(abs(Told-T)))
        Told = T;
        jacobi_iter = jacobi_iter + 1;
        plot(x,T(3,:)) % for viewing the convergence graph
        xlabel('X axis')
        ylabel('U')
        title_text = sprintf('Iteration number (Jacobi Method) = %d',jacobi_iter);
        title(title_text)
        pause (0.03)
    end

figure(2)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
colorbar;
title_text = sprintf('Iteration number(Jacobi Method) = %d',jacobi_iter);
title(title_text)

end

% GS Method
if iterative_solver == 2
    GS_iter = 1;
    while(error > tol)
        
        for i = 2:nx-1
        for j = 2:ny-1
            term1(i,j) = (T(i-1,j) + Told(i+1,j))/dx^2;
            term2(i,j) = (T(i,j-1) + Told(i,j+1))/dy^2;
            T(i,j) = (1/k) * (term1(i,j) + term2(i,j));
        end
        end
        error = max(max(abs(Told-T)))
        Told = T;
        GS_iter = GS_iter + 1;
        plot(x,T(3,:)) % for viewing the convergence graph
        xlabel('X axis')
        ylabel('U')
        title_text = sprintf('Iteration number (GS Method)= %d',GS_iter);
        title(title_text)
        pause (0.03)
    end

figure(2)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
colorbar;
title_text = sprintf('Iteration number (GS Method)= %d',GS_iter);
title(title_text)

end

% SOR Method
if iterative_solver == 3
    SOR_iter = 1;
    prompt = ('Enter Relaxation Factor, W:')
    w = input(prompt);
    while(error > tol)
        
        for i = 2:nx-1
        for j = 2:ny-1
            term1(i,j) = (T(i-1,j) + Told(i+1,j))/dx^2;
            term2(i,j) = (T(i,j-1) + Told(i,j+1))/dy^2;
            T(i,j) = ((1-w)*Told(i,j)) +(w * ((1/k) * (term1(i,j) + term2(i,j))));
        end
        end
        error = max(max(abs(Told-T)))
        Told = T;
        SOR_iter = SOR_iter + 1;
        plot(x,T(3,:)) % for viewing the convergence graph
        xlabel('X axis')
        ylabel('U')
        title_text = sprintf('Iteration number (SOR Method)= %d',SOR_iter);
        title(title_text)
        pause (0.03)
    end

figure(2)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
colorbar;
title_text = sprintf('Iteration number (SOR Method)= %d',SOR_iter);
title(title_text)

end

Code Explanation:

The main program starts with default commands on clearing all the stored values & the screen and closing the windows.

An input is prompted to the user to chose which iterative solver should be used to solve the equation.

The value of L (length) is assigned initially. Following this, the number of elements along x & y is provided with this the value of " dx & dy " is calculated.

Following this, the initial values are applied followed by the Boundary Conditions. Boundary conditions at corners are provided by taking the average of two adjacent sides.

The T value is copied to Told which is used inside the solver.

Constant, k is calculated at this step and the convergence criteria are specified.

At this point, an " if " condition is used to determine which solver should be used to solve. The method which is used to solve depends on the value of the variable " iterative_solver ". If the value is equal to one, the Jacobi method is used, if the value is two, the GS method is used & if the value is equal to three, the SOR method id used which requires additional value " w ".

Jacobi Method:

A while loop is used to make the solver run until the preferred error value is obtained. Here two " for " loops are used to consider the values in both dimensions (x & y). The value of T at any point is obtained by substituting the values obtained at the previous iteration on the steady-state equation which was derived. Outside the two " for " loops, Told values are updated & Error value is calculated. The " for "loop is repeated until the preferred error value is reached. A plot is made to see the proceedings of the calculations. This plot is done with suitable labels and the title is named after the current iteration number to preview the values after each iteration.

After the while loop is ended, Contour is plotted to check the results obtained from the solver. Results are discussed later in the part.

GS Method:

All the loops are the same as that of the Jacobi Solver. But while calculating for the values of T at any point, The values of T at the current iteration is used to solve.

Contour is plotted similarly as that of the Jacobi Method.

SOR Method:

This method uses a Relaxation factor and that factor is used to calculate the solution in addition to the GS method solver.

Results - Steady State Solver:

Jacobi Method:

The solution is converged at 208th iteration while using Jacobi Method

GS Method:

The solution is converged at 112th iteration while using GM Method

SOR Method:

The solution is converged at 75th iteration while using SOR Method with w=1.2

The solution is converged at 60th iteration while using SOR Method with w=1.3

Steady-State Solution Inference:

From the above results, we can conclude that all the method leads to the same solution but the iteration count differs for each method. The Jacobi method takes more iteration to solve, whereas the GS method converges quickly when compared with Jacobi but the SOR method converges quickly that both methods. If the relaxation increases, the SOR method converges more quickly but up to certain values of the factor, beyond that limit, an infinite loop occurs and no solution is obtained.

With respect to number of Iterations to solve => Jacobi > G S > SOR

Explicit Transient State Analysis:

The solution is calculated by the time marching method and calculated using the solution at the previous time step, then it is called the explicit method.

Deriving the Heat Conduction equation w r t to Explicit Transient state, Temperature at point (i,j) can be calculated the following equation

Program is written to solve this equation

Main Program to solve the Steady State Equation:

clear all
close all
clc

% Definig length
L = 1;

% number of grid points
nx = 10;
ny = 10;

% Defining No of Elements
x = linspace(0,L,nx);
y = linspace(0,L,ny);
  
% Definig the mesh
dx = x(2)-x(1);
dy = y(2)-y(1);

% Initialization
% Assuming the initial values at all points as 300 K
T = 300 * ones(nx,ny);

% Boundary Condition - given value
T_T = 600;
T_L = 400;
T_B = 900;
T_R = 800;

% Assigning BC to Corner Points
T(1,1)   = (T_B + T_L)/2 ;
T(nx,1)  = (T_T + T_L)/2 ;
T(1,ny)  = (T_B + T_R)/2 ;
T(nx,ny)  = (T_T + T_R)/2 ;

% Assigninn BC to line points
T(1,2:nx-1)  = T_B;
T(nx,2:ny-1) = T_T;
T(2:ny-1,1)  = T_L;
T(2:ny-1,nx) = T_R;

% Copying of BC's)
Told= T;

% Defining Time limits
dt = 2e-3; % timestep size in sec
end_time = 1;
nt = end_time/dt; % number of time steps

% defining Constant
alpha = 1.4;

% Calculating the value of k
CFLx = alpha * dt / (dx^2);
CFLy = alpha * dt / (dy^2);

% CFL Calculations
CFL = CFLx + CFLy

% Determining the stability of solution 
if CFL <= 0.5
    fprintf('The solution is stable')
else
    fprintf('The solution is unstable')
end

% Initial non solution data
error = 9e9;
iter = 1;

for ti = 2:nt
    for i = 2:nx-1
    for j = 2:ny-1
        T(i,j) = Told(i,j) + ((T(i+1,j) - (2 * T(i,j)) + T(i-1,j))* CFLx) + ((T(i,j+1) - (2 * T(i,j)) + T(i,j-1)) * CFLy);
    end
    end
    error = max(max(abs(Told-T)))
    iter = iter + 1;
    Told = T;
    plot(x,T(3,:)) % for viewing the convergence graph
    xlabel('X axis')
    ylabel('U')
    title_text = sprintf('Iteration number (Explicit Method) = %d',iter);
    title(title_text)
    pause(0.0003)
end
figure(2)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
colorbar;
title_text = sprintf('Temperature at %d (sec)',end_time);
title(title_text)

Code Explanation:

The main program starts with default commands on clearing all the stored values & the screen and closing the windows.

The value of L (length) is assigned initially. Following this, the number of elements along x & y is provided with this the value of " dx & dy " is calculated.

Following this, the initial values are applied followed by the Boundary Conditions. Boundary conditions at corners are provided by taking the average of two adjacent sides.

The T value is copied to Told which is used inside the solver.

Time step size and end limit are given as input and with these values, the number of time steps is calculated.

Constant, " alpha " is provided here and with this alpha value CFL value along x & y coordinates are calculated to determine the stability of the solution. Depending upon the value of CFL, a message is displayed stating the stability of the solution.

Error and iteration variables are introduced to the solver.

The first " for " loop of the solver contains time as the main factor with an increment of dt. Inside the time loop, two " for " loops are used to calculate the value of T at any point. This value of T is obtained by substituting the values obtained at the previous timestep on the transient state equation which was derived. Outside the " for " loop, Told values are updated & Error value is calculated. The time loop is repeated for the number of time steps given. A plot is made to see the proceedings of the calculations. This plot is done with suitable labels and the title is named after the current iteration number.

Results - Transient Explicit Solver:

This result is obtained as a result of the program.

While the calculation was in progress we can notice that after two hundred time steps, the value of error was almost stable whereas after four hundred time steps there is no change in the value of error which means the steady solution is obtained.

Implicit Transient State Analysis:

The solution is calculated by the time marching method and calculated using the solution at the same time step when multiplied by the variable alpha & dt, then it is called the Implicit method.

Deriving the Heat Conduction equation w r t to Implicit Transient state, Temperature at point (i,j) is calculated by the following equation

The program is written to solve this equation using Iterative Solvers and the results are discussed.

This program asks for which iterative method to be used for solving the equation. Please provide the input as listed below

1. Jacobi Method
2. GS Method
3. SOR Method - Requires additional input of w

Main Program to solve the Implicit Transient State Equation:

clear all
close all
clc

% Getting Solver Infor
fprintf('Steady State Method   ')
fprintf('Enter 1 for Jacobi Method  ')
fprintf('Enter 2 for GS Method  ')
fprintf('Enter 3 for SOR Method  ')
prompt = ('enter iteration method:')
iterative_solver = input(prompt);

% Definig length
L = 1;

% number of grid points
nx = 10;
ny = 10;

% Defining No of Elements
x = linspace(0,L,nx);
y = linspace(0,L,ny);
  
% Definig the mesh
dx = x(2)-x(1);
dy = y(2)-y(1);

% Initialization
% Assuming the initial values at all points as 300 K
T = 300 * ones(nx,ny);

% Boundary Condition - given value
T_T = 600;
T_L = 400;
T_B = 900;
T_R = 800;

% Assigning BC to Corner Points
T(1,1)   = (T_B + T_L)/2 ;
T(nx,1)  = (T_T + T_L)/2 ;
T(1,ny)  = (T_B + T_R)/2 ;
T(nx,ny)  = (T_T + T_R)/2 ;

% Assigninn BC to line points
T(1,2:nx-1)  = T_B;
T(nx,2:ny-1) = T_T;
T(2:ny-1,1)  = T_L;
T(2:ny-1,nx) = T_R;

% Copying of BC's)
Told = T;
Tpre = T;

% Defining Time limits
dt = 2e-3; % timestep size in sec
end_time = 1;
nt = end_time/dt; % number of time steps

% defining Constant
alpha = 1.4;

% Calculating the value of k
CFLx = alpha * dt / (dx^2);
CFLy = alpha * dt / (dy^2);

% CFL Calculations
CFL = CFLx + CFLy

% Determining the stability of solution 
if CFL <= 0.5
    fprintf('The solution is stable')
else
    fprintf('The solution is unstable')
end

% defining Coeff terms
term1 = 1/(1 + (2 * CFLx) + (2 * CFLy));
term2 = CFLx * term1;
term3 = CFLy * term1;

% Convergence Criteria
error = 9e9;
tol = 1e-4;

% Jacobi Method
if iterative_solver == 1
    Jacobi_iter = 1;
    for ti = 2:nt
        while(error > tol)
            for i = 2:nx-1
            for j = 2:ny-1
                T(i,j) = (term1 * Tpre(i,j)) + (term2 * (Told(i-1,j) + Told(i+1,j))) + (term3 * (Told(i,j+1) + Told(i,j-1)));
            end
            end
            error = max(max(abs(Told-T)))
            Told = T;
            plot(x,T(3,:)) % for viewing the convergence graph
            xlabel('X axis')
            ylabel('U')
            title_text = sprintf('Iteration number (Implicit Jacobi Method) = %d at time step %d' ,Jacobi_iter ,ti);
            title(title_text)
            pause(0.0003)
            Jacobi_iter = Jacobi_iter + 1;
        end
        Tpre = T;
        error = 9e9;
    end
    figure(2)
    contourf(x,y,T);
    [a b] = contourf(x,y,T);
    clabel(a,b);
    colorbar;
    title_text = sprintf('Iteration number (Implicit Jacobi Method) = %d at time step %d' ,Jacobi_iter ,ti);
    title(title_text)
end

% GS Method
if iterative_solver == 2
    GS_iter = 1;
    for ti = 2:nt
        while(error > tol)
            for i = 2:nx-1
            for j = 2:ny-1
                T(i,j) = (term1 * Tpre(i,j)) + (term2 * (T(i-1,j) + Told(i+1,j))) + (term3 * (Told(i,j+1) + T(i,j-1)));
            end
            end
            error = max(max(abs(Told-T)))
            Told = T;
            plot(x,T(3,:)) % for viewing the convergence graph
            xlabel('X axis')
            ylabel('U')
            title_text = sprintf('Iteration number (Implicit GS Method) = %d at time step %d' ,GS_iter ,ti);
            title(title_text)
            pause(0.0003)
            GS_iter = GS_iter + 1;
        end
        Tpre = T;
        error = 9e9;
    end
    GS_iter = GS_iter - 1
    figure(2)
    contourf(x,y,T);
    [a b] = contourf(x,y,T);
    clabel(a,b);
    colorbar;
    title_text = sprintf('Iteration number (Implicit GS Method) = %d at time step %d' ,GS_iter ,ti);
    title(title_text)
end


% SOR Method
if iterative_solver == 3
    SOR_iter = 1;
    prompt = ('Enter Relaxation Factor, W:')
    w = input(prompt);
    T_GS = T;
    for ti = 2:nt
        while(error > tol)
            for i = 2:nx-1
            for j = 2:ny-1
                T_GS(i,j) = (term1 * Tpre(i,j)) + (term2 * (T(i-1,j) + Told(i+1,j))) + (term3 * (Told(i,j+1) + T(i,j-1)));
            	T(i,j) = ((1-w) * Told(i,j)) + (w * T_GS(i,j));
            end
            end
            error = max(max(abs(Told-T)))
            Told = T;
            plot(x,T(3,:)) % for viewing the convergence graph
            xlabel('X axis')
            ylabel('U')
            title_text = sprintf('Iteration number (Implicit SOR Method) = %d at time step %d' ,SOR_iter ,ti);
            title(title_text)
            pause(0.0003)
            SOR_iter = SOR_iter + 1;
        end
        Tpre = T;
        error = 9e9;
    end
    SOR_iter = SOR_iter - 1
    figure(2)
    contourf(x,y,T);
    [a b] = contourf(x,y,T);
    clabel(a,b);
    colorbar;
    title_text = sprintf('Iteration number (Implicit SOR Method with w = %d) = %d at time step %d' ,w,SOR_iter ,ti);
    title(title_text)
end

Code Explanation:

The main program starts with default commands on clearing all the stored values & the screen and closing the windows.

An input is prompted to the user to chose which iterative solver should be used to solve the equation.

The value of L (length) is assigned initially. Following this, the number of elements along x & y is provided with this the value of " dx & dy " is calculated.

Following this, the initial values are applied followed by the Boundary Conditions. Boundary conditions at corners are provided by taking the average of two adjacent sides.

The T value is copied to Told & Tpre which is used inside the solver.

Time step size and end limit are given as input and with these values, the number of time steps is calculated.

Constant, " alpha " is provided here and with this alpha value CFL value along x & y coordinates are calculated to determine the stability of the solution. Depending upon the value of CFL, a message is displayed stating the stability of the solution.

Error variable is introduced to the solver & the tolerance value is specified.

At this point, an " if " condition is used to determine which solver should be used to solve. The method which is used to solve depends on the value of the variable " iterative_solver ". If the value is equal to one, the Jacobi method is used, if the value is two, the GS method is used & if the value is equal to three, the SOR method id used which requires additional value " w ".

Jacobi Method:

The first " for " loop of the solver contains time as the main factor with an increment of dt. Inside the time loop, a while loop is used to make the solver run until the preferred error value is obtained. Here two " for " loops are used to consider the values in both dimensions. The value of T at any point is obtained by substituting the values obtained at the previous iteration on the same time step of the transient-state equation which was derived. Outside the two " for " loop, Told values are updated & Error value is calculated. The " for "loop is repeated until the preferred error value is reached. A plot is made to see the proceedings of the calculations. This plot is done with suitable labels and the title is named after the current iteration number.

After the while loop is ended, the Tpre value is updated and the error value is reset so that the while loop is executed again with the increment in the value of time. This is repeated until the specified time value is reached.

After the execution of the time loop, Contour is plotted to check the results obtained from the solver. Results are discussed later in the part.

GS Method:

All the loops are the same as that of the Jacobi Solver. But while calculating for the values of T at any point, The values of T at the current iteration is used to solve.

Contour is plotted similarly as that of the Jacobi Method.

SOR Method:

This method uses a Relaxation factor and that factor is used to calculate the solution in addition to the GS method solver.

Results - Implicit Transient State Solver:

Jacobi Method:

The solution is obtained at 2317th iteration (after 500 Time Steps) while using the Jacobi Method

GS Method:

The solution is obtained at 1785th iteration (after 500 Time Steps) while using GM Method

SOR Method:

The solution is obtained at 1514th iteration (after 500 Time Steps) while using the SOR Method with w=1.1

The solution is obtained at 1485th iteration (after 500 Time Steps) while using the SOR Method with w=1.2

Implicit Transient-State Solution Inference:

From the above results, we can conclude that all the method leads to the same solution but the iteration count differs for each method. The Jacobi method takes more iteration to solve, whereas the GS method converges quickly when compared with Jacobi but the SOR method converges quickly that both methods. If the relaxation increases, the SOR method converges more quickly but up to certain values of the factor, beyond that limit, an infinite loop occurs and no solution is obtained.

With respect to number of Iterations to solve => Jacobi > G S > SOR