Solving the Steady and Unsteady 2D Heat Conduction Problem

Objective: -

In this project, we are going to solve the two dimensional Heat Conduction for steady and unsteady states.

For the Steady state, we have to solve it using three iterative approaches namely,

Jacobi Method

Gauss-Seidel Method and

Successive Over Realxation Method

For Transient state, we have to solve using the same iterative methods but in both Implicit and Explicit forms.

Solution: -

Steady and Unsteady Analyses are an integral part of solving and understanding Fluid flow and its properties. For the current project, we will be computing these anlyses for the 2 Dimensional Heat conduction problem.

Heat transfer is the transfer of thermal energy from a body, at a high temperature, to another at a lower temperature. This transfer of thermal energy may occur under steady or unsteady state conditions. In reality there is no 2 dimensional heat conduction, there is only the 3 Dimensional Heat Conduction. The 2-D Heat Conduction analysis is generally done to observe the flow and its properties in 2-D.

The 2-D Heat conduction Equation is given below:

`(deltaT)/(deltat) - alpha((delta^2T)/(deltax^2) + (delta^2T)/(deltay^2)) = 0`

where the first term, second term are the time derivative and the spatial derivatives respectively. `alpha` is the Thermal Diffusivity.

Derivations for Steady and Transient States:

Steady State:

A Fluid state in which there is no change in the fluid property with respect to time is known as a Steady State.

In the current case, the fluid property in question is temperature. So,

`(deltaT)/(deltat) = 0`

When the above equation is applied to the main Heat Conduction equation, we have

`(delta^2T)/(deltax^2) + (delta^2T)/(deltay^2) = 0`

assuming that the Thermal Conductivity `(alpha)` is constant.

According to the problem statement, we have to solve the Steady state heat conduction using three different iterative solvers namely,

Jacobi Method

Gauss-Seidel Method and

Successive Over Realxation Method

Upon discretising the above steady state equation using, Central Differencing Scheme (CDS), we obtain

`(T_L - 2*T_P + T_R)/(Deltax^2) + (T_T - 2*T_P + T_B)/(Deltay^2) = 0`

where

`T_P = T_(i,j)`

`T_L = T_(i-1,j)`

`T_R = T_(i+1,j)`

`T_T = T_(i,j+1)`

`T_B = T_(i,j-1)`

After Shuffling according to the common terms, we obtain

`(2*T_P)/(Deltax^2) + (2*T_P)/(Deltay^2) = (T_L+T_R)/(Deltax^2) + (T_T+T_B)/(Deltay^2)`

Assume `k = 2*[(1)/(Deltax^2) + (1)/(Deltay^2)]`, we have

`k*T_P = (T_L+T_R)/(Deltax^2) + (T_T+T_B)/(Deltay^2)`

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

Representing the above equation in an iterative form,

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

The above formulated equation is the Jacobi Method. Here, `n` represents the generation of the value that is being calculated as this is an iterative method. It is the same with all the other methods as well.

The Gauss-Seidel Method is similar to the Jacobi method with only one difference. In this method, we use the previously calculated value for the next generation of that value. The Steady State Heat Conduction equation using Gauss-Seidel method is shown below

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

The Successive Over-Relaxation (SOR) Method is a variant of the Gauss-seidel methd for solving a linear system of equations. The Steady State Heat Conduction equation using the SOR method is given below

`T_P^(n+1) = T_P^n*(1-lamda) + lambda*[(T_L+T_R)/(Deltax^2) + (T_T+T_B)/(Deltay^2)]^(n+1)`

where `lambda` is the over-relaxation factor

Transient State

A Fluid state in which there is a change in the fluid property with respect to time is known as an Unsteady or a Transient State. Transient 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 Transient state conditions.

The general representation of the Transient 2-D Heat conduction is given below,

`(deltaT)/(deltat) - alpha((delta^2T)/(deltax^2) + (delta^2T)/(deltay^2)) = 0`

Upon discretising the above equation using Forward Differencing Scheme (FDS) for the time derivative and Central Differencing Scheme (CDS) for the spatial derivatives, we obtain

`(T_P^(n+1)-T_P^n)/(Deltat) - alpha*((T_L - 2*T_P + T_R)/(Deltax^2) + (T_T - 2*T_P + T_B)/(Deltay^2)) = 0`

Upon simplifying, we get

`T_P^(n+1) = T_P^n + Deltat*alpha*((T_L - 2*T_P + T_R)/(Deltax^2) + (T_T - 2*T_P + T_B)/(Deltay^2))`

Assume `k_1 = (alpha*Deltat)/(Deltax^2)` and `k_2 = (alpha*Deltat)/(Deltay^2)`. We get,

`T_P^(n+1) = T_P^n + k_1*(T_L+2*T_P+T_R) + k_2*(T_T+2*T_P+T_B)`

As mentioned in the problem statement, we will be solving the Transient Heat conduction equation using the given iterative methods implementing them in both implicit and explicit forms.

Implicit Form

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

In the implicit form, the Transient equation is represented as

`T_P^(n+1) = T_P^n + k_1*(T_L^(n+1)+2*T_P^(n+1)+T_R^(n+1)) + k_2*(T_T^(n+1)+2*T_P^(n+1)+T_B^(n+1))`

After arranging all the `(n+1)` terms and simplifying the equation, we obtain

`T_P^(n+1) = T_P^n + (k_1*(T_L^(n+1)+T_R^(n+1)) + k_2*(T_T^(n+1)+T_B^(n+1)))/(1+2*k_1+2*k_2)`

Note: With every iteration, the `n` term does not change while the`(n+1)` terms change.

The Implicit form of the Jacobi method for the Transient Heat conduction equation is represented as,

`T_P^(n+1) = T_P^n + (k_1*(T_L^(n)+T_R^(n)) + k_2*(T_T^(n)+T_B^(n)))/(1+2*k_1+2*k_2)`

The Implicit form of the Gauss-seidel method is represented as,

`T_P^(n+1) = T_P^n + (k_1*(T_L^(n+1)+T_R^(n)) + k_2*(T_T^(n)+T_B^(n+1)))/(1+2*k_1+2*k_2)`

The Implicit form of the Successive Over-Relaxation method is represented as,

`T_P^(n+1) = T_P^n*(1-lamda) + lambda*[T_P^n + (k_1*(T_L^(n+1)+T_R^(n)) + k_2*(T_T^(n)+T_B^(n+1)))/(1+2*k_1+2*k_2)]`

Explicit Form

An explicit function is one which is given in terms of the independent variable.

In the case of explicit functions, we arrange the `n` terms and simplify to obtain

`T_P^(n+1) = T_P^n*(1-(2*k_1)-(2*k_2)) + (k_1*(T_L^n+T_R^n) + (k_2(T_B^n+T_T^n)))`

The above equation is the Explicit form of the Jacobi method for the Transient Heat conduction equation.

Similarly, The Explicit form of the Gauss-seidel method is represented as,

`T_P^(n+1) = T_P^n*(1-(2*k_1)-(2*k_2)) + (k_1*(T_L^(n+1)+T_R^n) + (k_2(T_B^(n+1)+T_T^n)))`

The Explicit form of the Successive Over-Relaxation method is represented as,

`T_P^(n+1) = T_P^n*(1-lambda) + lambda*[T_P^n*(1-(2*k_1)-(2*k_2)) + (k_1*(T_L^(n+1)+T_R^n) + (k_2(T_B^(n+1)+T_T^n)))]`

The above derived equations will be implemented in the matlab code.

Code:

The main code is the given below:

% Solving Steady and Unsteady states for 2-D conduction equation using
% iterative methods

clear all
close all
clc

% Initial Conditions
L = 1; % Unit square length domain
nx = 10; % grid points along X axis
ny = nx; % grid points along Y axis

% Grid spacing

x = linspace(0, L, nx);
dx = L/(nx-1);
y = linspace(0, L, ny);
dy = L/(ny-1);

error = 9e9;
tol = 1e-5;
t = 5;
dt = 1e-2;
alpha = 0.3;

k = 2*(1/(dx^2)+1/(dy^2)); % Assumption for Steady state constants
% Assumption for Unsteady state constants
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);

% Boundary Conditions
T = ones(nx, ny);
T(:, 1) = 400; % Left
T(:, end) = 800; % Right
T(1, :) = 600; % Top
T(end, :) = 900; % Bottom

% Edge Refining
T(1,1) = (600 + 400)/2;
T(1, end) = (600 + 800)/2;
T(end, 1) = (900 + 400)/2;
T(end, end) = (900 + 800)/2;

% Steady State Analysis
for iter_method = [1 2 3]
    if iter_method == 1
       time = jacobi( nx, ny, x, y, dx, dy, T, tol, k );
    else if iter_method == 2
            time = gauss_seidel( nx, ny, x, y, dx, dy, T, tol, k );
        else if iter_method == 3
                time = sor( nx, ny, x, y, dx, dy, T, tol, k );
            end
        end
    end
end

%Unsteady or Transient State Analysis
%Implicit Approach
for iter_method = [1 2 3]
    if iter_method == 1
       time = implicit_jacobi( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 );
    else if iter_method == 2
            time = implicit_gauss_seidel( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 );
        else if iter_method == 3
                time = implicit_sor( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 );
            end
        end
    end
end

%Explicit Approach
for iter_method = [1 2 3]
    if iter_method == 1
        time = explicit_jacobi( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 );
    else if iter_method == 2
            time = explicit_gauss_seidel( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 );
        else if iter_method == 3
                time = explicit_sor( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 );
            end
        end
    end
end

The code for the iterative methods are given below:

Steady State

function [ out ] = jacobi( nx, ny, x, y, dx, dy, T, tol, k )
% Jacobi method for steady state

meshgrid(x,y);

T_old = T;
jac_iter = 1;
error = 9e9;

tic
while (error > tol)
    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)));
    T_old = T;
    jac_iter = jac_iter + 1;
end
iter_time = toc;

figure(1)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Jacobi Method Iterations = %d', jac_iter);
txt2 = sprintf('Jacobi Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Steady State';txt1;txt2;txt3});

out = iter_time;
end


function [ out ] = gauss_seidel( nx, ny, x, y, dx, dy, T, tol, k )
% Gauss Seidel Method for Steady state

meshgrid(x,y);

T_old = T;
gs_iter = 1;
error = 9e9;

tic
while (error > tol)
    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)));
    T_old = T;
    gs_iter = gs_iter + 1;
end
iter_time = toc;

figure(2)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Gauss-Seidel Method Iterations = %d', gs_iter);
txt2 = sprintf('Gauss-Seidel Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Steady State';txt1;txt2;txt3});

out = iter_time;
end


function [ out ] = sor( nx, ny, x, y, dx, dy, T, tol, k )
% Successive Over-Relaxation Method for Steady state

lambda = 1.5; % Over Relaxation factor
meshgrid(x,y);

T_old = T;
sor_iter = 1;
error = 9e9;

tic
while (error > tol)
    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)));
    T_old = T;
    sor_iter = sor_iter + 1;
end
iter_time = toc;

figure(3)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Successive Over-Relaxation Method Iterations = %d', sor_iter);
txt2 = sprintf('Successive Over-Relaxation Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Steady State';txt1;txt2;txt3});

out = iter_time;
end

Transient State

Implicit Form

function [ out ] = implicit_jacobi( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 )
% Implicit Jacobi method for Transient state

meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
jac_iter = 1;


tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = (T_old(i-1, j) + T_old(i+1, j));
                term2 = (T_old(i, j-1) + T_old(i, j+1));
                T(i,j) = (To(i,j) + (k1*term1) + (k2*term2))/(1 + 2*k1 + 2*k2);
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        jac_iter = jac_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(4)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Implicit Jacobi Method Iterations = %d', jac_iter);
txt2 = sprintf('Implicit Jacobi Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Transient State';txt1;txt2;txt3});

out = iter_time;

end


function [ out ] = implicit_gauss_seidel( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 )
% Implicit Gauss Seidel method for Transient state

meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
gs_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = T(i-1, j) + T_old(i+1, j);
                term2 = T(i, j-1) + T_old(i, j+1);
                T(i,j) = (To(i,j) + (k1*term1) + (k2*term2))/(1 + 2*k1 + 2*k2);
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        gs_iter = gs_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(5)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Implicit Gauss-Seidel Method Iterations = %d', gs_iter);
txt2 = sprintf('Implicit Gauss-Seidel Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Transient State';txt1;txt2;txt3});

out = iter_time;

end


function [ out ] = implicit_sor( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 )
% Implicit Successive Over-Relaxation Method for Transient State

lambda = 1.09; % Over Relaxation factor
meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
sor_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = T(i-1, j) + T_old(i+1, j);
                term2 = T(i, j-1) + T_old(i, j+1);
                T(i,j) = (1-lambda)*T_old(i,j) + (lambda*((To(i,j) + (k1*term1) + (k2*term2))/(1 + 2*k1 + 2*k2)));
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        sor_iter = sor_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(6)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Implicit Successive Over-Relaxation Method Iterations = %d', sor_iter);
txt2 = sprintf('Implicit Successive Over-Relaxation Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Transient State';txt1;txt2;txt3});

out = iter_time;

end

Explicit Form

function [ out ] = explicit_jacobi( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 )
% Explicit Jacobi method for Transient state

meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
jac_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = (T_old(i-1, j) + T_old(i+1, j));
                term2 = (T_old(i, j-1) + T_old(i, j+1));
                T(i,j) = (To(i,j)*(1 - 2*k1 - 2*k2) + (k1*term1) + (k2*term2));
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        jac_iter = jac_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(7)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Explicit Jacobi Method Iterations = %d', jac_iter);
txt2 = sprintf('Explicit Jacobi Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Transient State';txt1;txt2;txt3});

out = iter_time;

end


function [ out ] = explicit_gauss_seidel( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 )
% Explicit Gauss Seidel method for Transient State

meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
gs_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = T(i-1, j) + T_old(i+1, j);
                term2 = T(i, j-1) + T_old(i, j+1);
                T(i,j) = ((To(i,j)*(1 - 2*k1 - 2*k2)) + (k1*term1) + (k2*term2));
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        gs_iter = gs_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(8)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Explicit Gauss-Seidel Method Iterations = %d', gs_iter);
txt2 = sprintf('Explicit Gauss-Seidel Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Transient State';txt1;txt2;txt3});

out = iter_time;

end


function [ out ] = explicit_sor( nx, ny, x, y, dx, dy, T, t, dt, tol, k, k1, k2 )
% Explicit Successive Over-Relaxation Method for Transient State

lambda = 1.3; % Over Relaxation factor
meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
sor_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = T(i-1, j) + T_old(i+1, j);
                term2 = T(i, j-1) + T_old(i, j+1);
                T(i,j) = (1-lambda)*T_old(i,j) + (lambda*((To(i,j)*(1 - 2*k1 - 2*k2)) + (k1*term1) + (k2*term2)));
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        sor_iter = sor_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(9)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X Length Domain');
ylabel('Y Length Domain');
txt1 = sprintf('Explicit Successive Over-Relaxation Method Iterations = %d', sor_iter);
txt2 = sprintf('Explicit Successive Over-Relaxation Method Simulation Time = %f s', iter_time);
txt3 = sprintf('Final error = %f', error);
colorbar
title({'Transient State';txt1;txt2;txt3});

out = iter_time;

end

Results and Observations:

Steady State Analysis

From the above plots, we can infer that the Successive Over Relaxation method is the most efficient as it has lower simulation time and number of iterations when compared to the Jacobi and Gauss-Seidel methods. It is also very important to understand that to obtain the most optimum or an efficient result in the SOR method, the over relaxation factor `(lambda)`, should be tweaked i.e according to the given problem parameters, the over relaxation factor should be experimented to obtain an efficient solution. For the above state analysis, the most efficient lambda value is 1.5.

Transient State Analysis

Implicit Approach

From the above plots, we observe a similar trend in the iterative method statistics i.e. In the above case of Transient state (Implicit), we infer that the SOR method is the most efficient as it has lower simulation time and iterations. As discussed in the previous plots, the over relaxation factor `(lamda)`is experimented to find the best solution possible. The over relaxation factor for the above state analysis is 1.09. We also observe that the final error is 0 which depicts that the solution method is accurate.

Explicit Approach

From the above plots, we observe a similar trend in the iterative method statistics i.e. In the above case of Transient state (Explicit), we infer that the SOR method is the most efficient as it has lower simulation time and iterations. As discussed in the previous plots, the over relaxation factor `(lamda)`is experimented to find the best solution possible. The over relaxation factor for the above state analysis is 1.3. We also observe that the final error is 0 which depicts that the solution method is accurate.

Inference

The above table shows the overall statistics of the methods implemented in the current analysis. We can successfully infer that of all the iterative methods used, the Successive Over-Relaxation method has the best results for the analysis as it has the lowest possible simulation time and iterations. Of all the approaches present, the Explicit Transient state analysis has the best possible solution statistics with respect to the simulation time. The steady state also has the lowest stats, but it is important to understand that the steady state is time independant and does not rely on Thermal Conductivity.

As discussed, it is very important to understand that the SOR method also depends on a lot of factors like the Over Relaxation Factor, Thermal Conductivity, grid size and the time parameters. Particular to the SOR method, the over relaxation factor gives an edge when compared to the other methods as it can be tweaked easily to find the best possible accurate solution.