Solving the steady and unsteady 2D heat conduction problem

REPORT – SOLVING THE STEADY AND UNSTEADY 2D HEAT CONDUCTION PROBLEM

AIM

Steady-State Analysis and Transient State Analysis

OBJECTIVE

• Solve the 2D heat conduction equation by using the point iterative techniques.
• Implement the following methods for solving implicit equations;
  • Jacobi
  • Gauss – seidel
  • Successive over-relaxation
• Write the transient analysis code for both implicit and explicit schemes.

GIVEN

The Boundary conditions for the problem are as follows;

• Top Boundary = 600 k
• Bottom Boundary = 900k
• Left Boundary = 400 k
• Right Boundary = 800 k
• Absolute error criteria = 1e - 4

SOLUTION

Thermal Conduction

Heat is transferred by conduction when adjacent atoms or molecules collide, or as several electrons move backwards and forwards from atom to atom in a diagnosed way so as not to form a macroscopic electric current, or as photons collide and scatter. Conduction is the most significant means of heat transfer within a solid or between solid objects in thermal contact.

The 2D temperature heat equation is a partial differential equation.

`(dT)/dt – alpha* ( (partial^2T)/(partialx^2) + ( partialT^2)/(partialy^2)) = 0`

Here, alpha is a positive coefficient called the thermal diffusivity of the medium.

Steady-State Conduction

Steady-state conduction is the form of conduction that happens when the temperature difference driving the conduction is constant so that the spatial distribution of temperatures in the conducting object does not change any further. Thus, all partial derivatives of temperature concerning space may either be zero or have nonzero values, but all derivatives of temperature at any point concerning time are uniformly zero. In steady-state conduction, the amount of heat entering any region of an object is equal to the amount of heat coming out, i.e.

`( partial^2T)/(partialx^2) + (partialT^2)/(partialy^2) = 0`

Transient Conduction

During any period in which temperatures change in time at any place within an object, the mode of thermal energy flow is termed transient conduction. Another term is "non-steady-state" conduction, referring to the time-dependence of temperature fields in an object. Non-steady-state situations appear after an imposed change in temperature at a boundary of an object. The formula for the 2D transient heat conduction equation is:

`(partialT)/(partialt) = (alpha)*((partial^2T)/(partialx^2) + (partial^2T)/(partialy^2))`

EXPLICIT AND IMPLICIT METHODS are approaches used in numerical analysis for obtaining numerical approximations to the solutions of time-dependent ordinary and partial differential equations, as is required in computer simulations of physical processes.

Explicit Method

Implicit Method

· Explicit method calculates the state of a system at a later time from the state of the system at the current time. · Explicit methods require small time steps to keep the error in the result bounded. · With the decrease in time step the number of computations increases and hence the computational time required for the solution to converge is less and the accuracy of the solution is high. · This method is considered conditionally stable because if the time step is too large, there will be a significant numerical error.

· Implicit method finds a solution by solving an equation involving both the current state of the system and the later state. · In the implicit method, general equilibrium is achieved by iteration at first, followed by evaluation of local element variables. · This method considers simultaneous equations and uses guess values as the initial values to find the approximated values. If the guess values are closer to the desired value the computation time will be less. · Relative to the explicit method, the computational time required for the solution to converge is high and the accuracy of the solution is less. · For the implicit method, there is no stability criteria requirement so more significant time steps can be used. · Implicit methods require an extra computation since the matrix coefficients should recalculate with different time steps.

Iterative Method

In computational mathematics, an iterative method is a mathematical procedure that uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the n-th approximation is derived from the previous ones. A specific implementation of an iterative method, including the termination criteria, is an algorithm of the iterative method. An iterative method is called convergent if the corresponding sequence converges for given initial approximations.

Iterative Solver

1. Jacobi iterative solver

The Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges.

2. Gauss-Seidel iterative solver

The Gauss-seidel method, or the method of successive displacement, is an iterative method used to solve a system of linear equations similar to the Jacobi method. The temperature at n+1 iterations uses values of n+1 iterations that have already been computed. In this path, old approximated values are overwritten with new approximate values for the iterations performed. Unlike Jacobi, Gauss-Seidel requires one storage vector which can be overwritten with each computation. This helps the solution to converge faster than the Jacobi method.

3. Successive over-relaxation iterative solver

SOR is a variant of the Gauss-Seidel method for solving a linear system of equations, resulting in faster convergence. This exploration takes the form of a weighted average between the previous iterate and the computed Gauss-Seidel iterate successively for each component.

Formula: T_n+1 = T_n . (1-w) + w. T(GS)_n+1 or T(Jacobi)_n+1

Here, w is the relaxation factor.

Range: 0 < w < 2

• 0 < w < 1 (Under-Relaxed state)

This condition is used when the temperature gradient is high and due to this condition, the time step does not overshoot and helps to converge at a faster rate.

• 1 < w < 2 (Over-Relaxed state)

This condition is used when the temperature gradient is low and due to this condition more significant time steps are used to obtain the steady-state, i.e. for the values to converge at a faster rate.

Nomenclature

• T_P = T_{i, j}
• T_L = T_{i-1, j} (Left boundary)
• T_R = T_{i+1, j} (Right boundary)
• T_T = T_{i, j+1} (Top boundary)
• T_B = T_{i, j-1} (Bottom boundary)

Equations used

• 2D heat conduction equation:

`(dT)/dt – alpha* ( (partial^2T)/(partialx^2) + ( partialT^2)/(partialy^2)) = 0`

• In discrete form:

Here,

• n represents the current time step and n+1 represents the next time step.
• i, i+1 are the neighbouring points in the x-direction.
• j, j+1 are the neighbouring points in the y-direction.
• `/_x` is grid spacing in the x-direction.
• `/_y` is grid spacing in the y-direction.

The time derivative term is hyperbolic in nature which is to be discretised using a forward differencing scheme and the space derivatives are elliptical in nature which is to be discretised using a central differencing scheme.

• Steady State 2D Heat Conduction Equation

`( partial^2T)/(partialx^2) + (partialT^2)/(partialy^2) = 0`

• Steady State Discretization

(T_L – 2.T_P + T_R)/ `/_x^2` + (T_T – 2.T_P + T_B)/ `/_y^2` = 0;

• Therefore,

T_P = K1 (T_T + T_R) + K2 (T_T + T_B)

Where,

K1 = `(/_y^2)/(2.(/_x^2 + /_y^2))` and K2 = `(/_x^2)/(2.(/_x^2 + /_y^2))`

• Using iterative solvers for solving the conduction equation steady state:

Jacobi iterative solver

(T_P)ⁿ⁺¹ = K1 (T_L + T_R)ⁿ + K2 (T_T + T_B)ⁿ

Gauss-Seidel iterative solver

(T_P)ⁿ⁺¹ = K1 (T_L + T_R)ⁿ⁺¹ + K2 (T_T + T_B)ⁿ⁺¹

Successive over relaxation iterative solver

(T_P)ⁿ⁺¹ = (T_P)ⁿ.(1-w) + w.( K1(T_L +T_R)ⁿ⁺¹ + K2(T_T + T_B)ⁿ⁺¹)

• Transient State

`(partialT)/(partialt) = alpha*((partial^2T)/(partialx^2) + (partial^2T)/(partialy^2))`

• Unsteady-State Discretization

• Solving Explicitly

(T_P)ⁿ⁺¹ = (1 - 2.K₁ – 2.K₂).(T_P)ⁿ + K₁(T_L + T_R)ⁿ + K₂(T_T + T_B)ⁿ

Here,

K₁ = `alpha . (/_t)/(/_x^2)`

K₂ = `alpha . (/_t)/(/_y^2)`

• Stability Criteria

K₁ + K₂ `<=` 0.5

• Solving Implicitly

(1 - 2.K₁ – 2.K₂).(T_P)ⁿ⁺¹ = (T_P)ⁿ + K₁ (T_L + T_R)ⁿ⁺¹ + K₂ (T_T + T_B)ⁿ⁺¹

Here,

K₁ = `alpha . (/_t)/(/_x^2)`

K₂ = `alpha . (/_t)/(/_y^2)`

Jacobi iterative solver

(T_P)ⁿ⁺¹ = ((T_P)ⁿ + K₁ ((T_L)ⁿ + (T_R)ⁿ) + K₂ ((T_T)ⁿ + (T_B)ⁿ) / (1 - 2.K₁ – 2.K₂)

Gauss-Seidel iterative solver

(T_P)ⁿ⁺¹ = ((T_P)ⁿ + K₁ ((T_L)ⁿ⁺¹ + (T_R)ⁿ) + K₂ ((T_T)ⁿ + (T_B)ⁿ⁺¹) / (1 - 2.K₁ – 2.K₂)

Successive over-relaxation iterative solver

(T_P)ⁿ⁺¹ = (T_P)ⁿ . (1 - w) + w. [((T_P)ⁿ + K₁ ((T_L)ⁿ⁺¹ + (T_R)ⁿ) + K₂ ((T_T)ⁿ + (T_B)ⁿ⁺¹) / (1 - 2.K₁ – 2.K₂)]

ASSUMPTIONS

• Consider computation area as a square.
• nx = ny; i.e., the number of nodes in the x-direction is equal to the number of nodes in the y-direction.
• nx = ny = 10
• Apart from the boundary conditions, the inner nodes have an initial temperature value of 300K.

Boundary Conditions

• Top Boundary = 600 k
• Bottom Boundary = 900k
• Left Boundary = 400 k
• Right Boundary = 800 k

CODE

Code to study the behaviour of STEADY and UNSTEADY states for 2D heat conduction equation:

MAIN CODE

clear all
close all 
clc

% Defining domain
nx = 10;
ny = 10;
x = linspace(0,1,nx);
y = linspace(0,1,ny);

% Grid size
dx = x(2)-x(1);
dy = y(2)-y(1);

% Defining temperature
T = 300*ones(nx,ny);

% Boundary Conditions
T(1,1:end) = 600;
T(end,1:end) = 900;
T(1:end,1) = 400;
T(1:end,end) = 800;

% Temperature at edges
T(1,1) = (T(1,2)+T(2,1))/2;
T(end,1) = (T(end-1,1)+T(end,2))/2;
T(1,end) = (T(1,end-1)+T(2,end))/2;
T(end,end) = (T(end-1,end)+T(end,end-1))/2;

% Absolute error criteria
tol = 1e-4;

% Initial temperature
Told = T;

% Some constants
k = 2*(dx^2+dy^2)/(dx^2*dy^2);

% Relaxation Factor
w = 1.045 ;

% alpha
% Thermal diffusivity of the medium
A = 1.4;

% time interval 
dt = 1/(2*((A/dx^2)+(A/dy^2)));

% constants
k1 = dt*A/dx^2;
k2 = dt*A/dy^2;

T_initial = T;

% Relaxation factor for unsteady state
W = 1.15;

% Function for respective solution
[error_j,iter1] = jacobi_steady_function(nx,ny,T,Told,dx,dy,x,y,k,tol);
[error_g,iter2] = gauss_seidel_function(tol,nx,ny,T,Told,dx,dy,x,y,k);
[error_s,iter3] = successive_over_relaxation_function(tol,nx,ny,T,Told,dx,dy,x,y,k,w);
[error,iter] = unsteady_explicit_function(nx,ny,k1,k2,Told,T,tol);
[error_ja,iter_ja] = unsteady_jacobi_function(nx,ny,k1,k2,T,Told,tol,T_initial)
[error_gs,iter_gs] = unsteady_gauss_seidel_function(nx,ny,k1,k2,T,Told,tol,T_initial)
[error_sor,iter_sor] = unsteady_successive_over_relaxation_function(nx,ny,k1,k2,T,Told,tol,T_initial,W)

FUNCTION CODE

• Jacobi iteration method for steady-state

function [error_j,iter1] = jacobi_steady_function(nx,ny,T,Told,dx,dy,x,y,k,tol)

% dummy error value 
error_j = 9e9;

% Definig iteration variable to calculate no of iterations
iter1 = 0;

    % convergence loop 
    while (error_j > tol)

          % iteration loop
          for i = 2:nx-1

               for j = 2:ny-1

                       % Formula to calculate the temprature at each grid point
                       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 for each iteration 
          error_j = max(max(abs(Told-T)));

          % updating the temperature values for the next iteration
           Told = T;

           % Counting the no of iterations
            iter1 = iter1 + 1;
    
    end

    % Plotting
    figure(1)
    contourf(x,y,T);
    clabel(contourf(x,y,T));
    colorbar
    colormap("jet")
    set(gca,'ydir','reverse')
    title(sprintf('No. of Iterations = %d',iter1));

end

OUTPUT

• Gauss-Seidel iterative method for steady-state

function [error_g,iter2] = gauss_seidel_function(tol,nx,ny,T,Told,dx,dy,x,y,k)

% Definig iteration variable to calculate no of iterations
iter2 = 0;

% dummy error value 
error_g = 9e9;

% convergence loop 
while (error_g > tol)

    % iteration loop
    for i = 2:(nx-1)

        for j = 2:(ny-1)

            % Formula to calculate the temprature at each grid point
            T(i,j) = (1/k)*(((T((i-1),j)+T((i+1),j))/dx^2)+((T(i,(j-1))+T(i,(j+1)))/dy^2));

        end
        
    end

    % error for each iteration 
    error_g = max(max(abs(Told-T)));

    % updating the temperature values for the next iteration
    Told = T;

    % Counting the no of iterations
    iter2 = iter2 + 1;
   
end
    % Plotting
    figure(2)
    contourf(x,y,T);
    clabel(contourf(x,y,T));
    colorbar
    colormap("jet")
    set(gca,'ydir','reverse')
    title(sprintf('No. of Iterations = %d',iter2));

end

OUTPUT

• Successive over-relaxation method for steady-state

function [error_s,iter3] = successive_over_relaxation_function(tol,nx,ny,T,Told,dx,dy,x,y,k,w)

% Definig iteration variable to calculate no of iterations
iter3 = 0;

% dummy error value 
error_s = 9e9;

% convergence loop 
while (error_s > tol)

    % iteration loop
    for i = 2:(nx-1)

        for j = 2:(ny-1)

            % Formula to calculate the temprature at each grid point
            T(i,j) = Told(i,j)*(1-w)+w*((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 for each iteration 
    error_s = max(max(abs(Told-T)));

    % updating the temperature values for the next iteration
    Told = T;

    % Counting the no of iterations
    iter3 = iter3 + 1;
   
end
% Plotting
    figure(3)
    contourf(x,y,T);
    clabel(contourf(x,y,T));
    colorbar
    colormap("jet")
    set(gca,'ydir','reverse')
    title(sprintf('No. of Iterations = %d',iter3));

end

OUTPUT

• Solving for unsteady-state explicitly

function [error,iter] = unsteady_explicit_function(nx,ny,k1,k2,Told,T,tol)

% Dummy Error
error = 9e9;

% Iteration Value
iter = 0;

% Convergence loop
while (error > tol)

    for i = 2:nx-1

        for j = 2:ny-1

            X = Told(i-1,j) - 2*Told(i,j) + Told(i+1,j);
            Y = Told(i,j-1) - 2*Told(i,j) + Told(i,j+1);

            T(i,j) = Told(i,j) + k1*X + k2*Y;

        end

    end

    % Max Error for the iteration 
    error = max(max(abs(Told - T)));

    % Updating Temperature Value
    Told =T;

    % No. of iteration
    iter = iter + 1;

end

    % Plotting
    figure(4)
    [C,H] = contourf(T);
    clabel(C,H);
    axis ij
    title(sprintf('no. of iterations = %d', iter));
    
end

OUTPUT

• Unsteady state using Jacobi iteration

function [error_ja,iter_ja] = unsteady_jacobi_function(nx,ny,k1,k2,T,Told,tol,T_initial)

% Dummy Error
error_ja = 9e9;
error_t = 9e9;

% Iteration Value
iter_ja = 0;

% Time loop
while (error_t > tol)

    % Convergence loop
    while (error_ja > tol)
        
        % Space loop
        for i = 2:nx-1

            for j = 2:ny-1

              term1 = (1+(2*k1)+(2*k2))^-1;

              H = Told(i-1,j) + Told(i+1,j);
              V = Told(i,j-1) + Told(i,j+1);

              T(i,j) = (T_initial(i,j)*term1) + (term1*k1*H) + (term1*k2*V);

            end
            
        end

        % error for convergence loop 
        error_ja = max(abs(Told-T));

        % Temperature value
        Told = T;

        % Calculating no. of iteration
        iter_ja = iter_ja + 1;

    end

    % error for time loop
    error_t = max(max(abs(T_initial - T)));

    % Assigning value
    T_initial = T;

    
    error_ja = 9e9;

end

% Plotting
figure(5)
contourf(T);
clabel(contourf(T))
axis ij
title(sprintf('No. of iterations = %d',iter_ja))

end

OUTPUT

• Unsteady state using Gauss-Seidel iteration

function [error_gs,iter_gs] = unsteady_gauss_seidel_function(nx,ny,k1,k2,T,Told,tol,T_initial)

% Dummy Error
error_gs = 9e9;
error_t = 9e9;

% iteration value
iter_gs = 0;

% time loop
while (error_t > tol)

    % convergence loop
    while (error_gs > tol)
        
        % space loop
        for i = 2:nx-1

            for j = 2:ny-1

              term1 = (1+(2*k1)+(2*k2))^-1;

              H = T(i-1,j) + T(i+1,j);
              V = T(i,j-1) + T(i,j+1);

              T(i,j) = (T_initial(i,j)*term1) + (term1*k1*H) + (term1*k2*V);

            end
            
        end

        % max error
        error_gs = max(abs(Told-T));

        % temperature value
        Told = T;

        % no. of iterations
        iter_gs = iter_gs + 1;

    end

    % max error for time loop
    error_t = max(max(abs(T_initial - T)));

    T_initial = T;
    
    error_gs = 9e9;

end

% Plotting
figure(6)
contourf(T);
clabel(contourf(T))
axis ij
title(sprintf('No. of iterations = %d',iter_gs))

end

OUTPUT

• Unsteady state using Successive over relaxation iteration

function [error_sor,iter_sor] = unsteady_successive_over_relaxation_function(nx,ny,k1,k2,T,Told,tol,T_initial,W)

% Dummy Error
error_sor = 9e9;
error_t = 9e9;

% iteration
iter_sor = 0;

% time loop
while (error_t > tol)

    % convergence loop
    while (error_sor > tol)
        
        % space loop
        for i = 2:nx-1

            for j = 2:ny-1

              term1 = (1+(2*k1)+(2*k2))^-1;

              H = T(i-1,j) + T(i+1,j);
              V = T(i,j-1) + T(i,j+1);

              T(i,j) = (T_initial(i,j)*term1) + (term1*k1*H) + (term1*k2*V);

              T(i,j) = (Told(i,j))*(1-W) + W*T(i,j);

            end
            
        end

        % max error
        error_sor = max(abs(Told-T));

        % temperature
        Told = T;

        % no. of iterations
        iter_sor = iter_sor + 1;

    end

    error_t = max(max(abs(T_initial - T)));

    T_initial = T;
    
    error_sor = 9e9;

end

% Plotting
figure(7)
contourf(T);
clabel(contourf(T))
axis ij
title(sprintf('No. of iterations = %d',iter_sor))

end

OUTPUT