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

Aim:Solving the steady and unsteady 2D heat conduction problem

objective:

1. Steady state analysis & Transient State Analysis

Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

You will implement the following methods for solving implicit equations.

• Jacobi
• Gauss-seidel
• Successive over-relaxation

2. You will write code for both implicit and explicit schemes (This is for the transient analysis only)

3. Compare the convergence rate [Number of iterations] of the steady state and transient [explicit and implicit] simulations and justify the numerics.

Solution:

Here we will understand the concept of 2d equation of heat conduction in both steady state and transient state of flow. We need to understand the derivation part where the terms are discretized step by step and then we move on to the encoding part. There are a few things to consider before going through this, that for the steady state we need to use iterative solvers and for the transients we need to solve both the schemes we saw in the previous section as IMPLICIT and EXPLICIT.

Now we are going to understand the derivation part of the steady-state condition.

 Solving this by the methods 

1. Jacobi

2. Gauss-seidel

3. Successive Over Relaxation SOR

since d2t/dx2 + d2T/dy2 = 0, its has elliptic nature properties which we an discretize them using central differencing scheme

now what happens when represent them as jacobi

for Gauss seidel we represent

SOR is enhanced gauss-seidel method where nth time step and n+1 time step is all same.

 where omega is the relaxation factor as 1.4.

now we begin the coding part which is done in MATLAB where these all terms are considered by applying for loop and while loop.

% Solving Steady and unsteady 2D heat conduction problem
% Steady-state conduction 
% Iterative solver - Jacobi, Gauss -seidel, Relaxation methods

clear all
close all
clc

%initial conditions
nx = 10;
ny = 10;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);


% given data of boundary condition
% initial temp condition

t = 300*ones(nx,ny);  
t(:,1) = 400;
t(1,:) = 600;
t(:,ny) = 800;
t(nx,:) = 900;

% defining the temperature at corner points

t(1,1) = 500;
t(nx,ny) = 850;
t(1,ny) = 700;
t(nx,1) = 650;

dt = 1e-2; %time step
told=t;

error = 100;
tol = 1e-4;
k = 2*(dx^2+dy^2)/(dx^2*dy^2);

method = input('Iterative method = '); %iterative method selection

%JACOBI ITERATIVE METHOD

if (method==1)
    tic;
    jacobi = 1;
    while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
            t(i,j) = (((told(i-1,j)+told(i+1,j))/(dx^2))+((told(i,j+1)+told(i,j-1))/(dy^2)))/k;
            end
        end
    error = max(max(abs(t-told)));
    told = t;
    jacobi = jacobi+1;
    time_taken = toc;
    end

%plotting    
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
set(gca,'YDIR','reverse')
text1 = sprintf('number of iteration = %d', jacobi);
text2 = sprintf('simulation time = %f', time_taken);
title({'Steady state heat condustion profile by jacobi'; text1; text2});
end


%GAUSS-SEIDEL ITERATIVE METHOD

if (method==2)
    tic;
    gauss = 1;
    while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
            t(i,j) = (((t(i-1,j)+told(i+1,j))/(dx^2))+((told(i,j+1)+t(i,j-1))/(dy^2)))/k;
            end
        end
    error = max(max(abs(t-told)));
    told = t;
    gauss = gauss+1;
    time_taken = toc;
    end

%plotting
figure(2)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
set(gca,'YDIR','reverse')
text1 = sprintf('number of iteration = %d', gauss);
text2 = sprintf('simulation time = %f', time_taken);
title({'Steady state heat condustion profile by Gauss-Seidel'; text1; text2});
end


%SUCCESSIVE OVER RELAXATION ITERATIVE METHOD

if (method==3)
    tic;
    sor = 1;
    omega = 1.4; %over relaxation factor
    while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
            term1 =(t(i-1,j)+told(i+1,j))/(dx^2);
            term2 =(told(i,j+1)+t(i,j-1))/(dy^2);
            t(i,j) = ((told(i,j)*(1-omega)))+(omega*(term1+term2))/k;
            end
        end
    error = max(max(abs(t-told)));
    told = t;
    sor = sor+1;
    time_taken = toc;
    end

%plotting
figure(3)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
set(gca,'YDIR','reverse')
text1 = sprintf('number of iteration = %d', sor);
text2 = sprintf('simulation time = %f', time_taken);
text3 = sprintf('over relaxation factor = %f', omega);
title({'Steady state heat condustion profile by Successive Over Relaxation'; text1; text2; text3});
end

Output

By Jacobi Method

By Gauss-seidel Method

By Successive Over Relaxation method

From above output results we can see that the SOR method gives a better in terms of iterations carried out by them and no. of simulations time compare to other methods 

Simulation of a transient (unstable) state:

We will move on to the transient state simulation, where we will work on implicit and explicit schemes.

1. Implicit schemas

2. Explicit schemes

What are these schematics?

Explicit and implicit methods are approaches used in numerical analysis to obtain numerical approximations to solving time-dependent ordinary and PDR as required in computer simulations of physical processes. Explicit methods calculate the state of the system at a later time from the state of the system at the current time, while implicit methods find the solution by solving an equation involving both the current state of the system and the later state. Mathematically, if

 is the current system state and  is the state at the later time( is a small time step),

then, for an explicit method

while for an implicit method one solves an equation

By derivation we have :

Applying the first heat conduction equation in (1) to node  at the time moment of , the equation can be rewritten as

The partial differential in the two sides of (2) can be approximated by difference quotient. The temperature item in the right of the equal sign can be approximated by first-order time forward difference quotient.

The second-order partial differential in the left of the equal sign can be approximated by the central difference quotient.

Substituting the value, we can get the difference equation of 

Equation is the difference equation of heat conduction equation, and the truncation error [34] is O(DT + x2. y2).

Assuming that deltax = deltay = delta  and substituting it into the latest eq., we can obtain

While we discretize the whole equation we get here.

By substituting the discrtize values to Jacobi, Gauss-seidel and SOR. we get

 Jacobi

  Gauss-seidel 

 Over Relaxation

Below is Transient State Implicitly,  MATLAB code:

% Solving Steady and Transient 2D heat conduction problem
% Transient state solution with IMPLICIT method
% Iterative solver are - Jacobi, Gauss -seidel, Relaxation methods


clear all
close all
clc

%initial conditions
%L = 1
nx = 10;
ny = 10;
nt = 300; % No.of steps
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);


% given data of boundary condition
% initial temp condition

t = 300*ones(nx,ny);  
t(:,1) = 400;
t(1,:) = 600;
t(:,ny) = 800;
t(nx,:) = 900;

% defining the temperature at corner points

t(1,1) = 500;
t(nx,ny) = 850;
t(1,ny) = 700;
t(nx,1) = 650;

dt = 0.01; %time step
told=t;
alpha = 0.3; % Thermal diffusivity
tprev_dt = told;

tol = 1e-4;
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);

method = input('Iterative method = '); %iterative method selection

%JACOBI ITERATIVE METHOD

if (method==1)
    tic;
    jacobi = 1;
    for p = 1:nt
        error =100;
         while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                term1 = (1+(2*k1)+(2*k2))^-1;
                term2 = k1*term1;
                term3 = k2*term1;
                h = (told(i-1,j)+told(i+1,j));
                v = (told(i,j+1)+told(i,j-1));
                t(i,j) = ((tprev_dt(i,j)*term1)+(h*term2)+(v*term3));
            end
        end
        
    error = max(max(abs(t-told)));
    told = t;
    jacobi = jacobi+1;
         end
         tprev_dt = t ;
    end
    time_taken = toc;

%plotting    
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
set(gca,'YDIR','reverse')
text1 = sprintf('number of iteration = %d', jacobi);
text2 = sprintf('simulation time = %f', time_taken);
title({'Transient state heat conduction profile Implicit-Jacobi'; text1; text2});
end


%GAUSS-SEIDEL ITERATIVE METHOD

if (method==2)
    tic;
    gauss = 1;
    for p = 1:nt
        error =100;
         while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                term1 = (1+(2*k1)+(2*k2))^-1;
                term2 = k1*term1;
                term3 = k2*term1;
                h = (t(i-1,j)+told(i+1,j));
                v = (told(i,j+1)+t(i,j-1));
                t(i,j) = ((tprev_dt(i,j)*term1)+(h*term2)+(v*term3));
            end
        end
        
    error = max(max(abs(t-told)));
    told = t;
    gauss = gauss+1;
         end
         tprev_dt = t; 
    end
    time_taken = toc;

%plotting    
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
set(gca,'YDIR','reverse')
text1 = sprintf('number of iteration = %d', gauss);
text2 = sprintf('simulation time = %f', time_taken);
title({'Transient state heat conduction profile Implicit-Gauss-seidel'; text1; text2});
end

%SUCCESSIVE OVER RELAXATION ITERATIVE METHOD

if (method==3)
    tic;
    sor = 1;
    omega = 1.2; % Over relaxation factor
    for p = 1:nt
        error =100;
         while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                term1 = (1+(2*k1)+(2*k2))^-1;
                term2 = k1*term1;
                term3 = k2*term1;
                h = (t(i-1,j)+told(i+1,j));
                v = (told(i,j+1)+t(i,j-1));
                t(i,j) = (told(i,j)*(1-omega))+(omega*((tprev_dt(i,j)*term1)+(h*term2)+(v*term3)));
            end
        end
        
    error = max(max(abs(t-told)));
    told = t;
    sor = sor+1;
         end
         tprev_dt = t ;
    end
    time_taken = toc;

%plotting    
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
set(gca,'YDIR','reverse')
text1 = sprintf('number of iteration = %d', sor);
text2 = sprintf('simulation time = %f', time_taken);
title({'Transient state heat conduction profile Implicit-SOR'; text1; text2});
end

 Output:

The conclusion of the above results are just as silimar as what was in the steady state, and it was observed that SOR methods gives better results with no. of simulation and time elapsing on the run case.

Transient State with Explicit Scheme 

after discretizing the values for explicit schemes we have equation.

code for Explicit schemes of Transient state 

% Solving Steady and Transient 2D heat conduction problem
% Transient state solution with EXPLICIT


clear all
close all
clc

%initial conditions
%L = 1
nx = 10;
ny = 10;
nt = 300; % No.of steps
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);

% initial temp condition

t = 300*ones(nx,ny);  
t(:,1) = 400;
t(1,:) = 600;
t(:,ny) = 800;
t(nx,:) = 900;

% defining the temperature at corner points

t(1,1) = 500;
t(nx,ny) = 850;
t(1,ny) = 700;
t(nx,1) = 650;

dt = 0.01; %time step
told=t;
alpha = 0.3; % Thermal diffusivity

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
CFL = alpha*dt/(dx^2);

%EXPLICIT METHOD
% running for loop

tic;
    explicit = 1;
    for p = 1:nt
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                term1 = (told(i-1,j)-2*told(i,j)+told(i+1,j));
                term2 = (told(i,j-1)-2*told(i,j)+told(i,j+1));
                t(i,j) = ((told(i,j))+(k1*term1)+(k2*term2));
            end
        end
        
    told = t;
    explicit = explicit+1;
    end
    
  time_taken = toc; 
  
 %Plotting   
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
set(gca,'YDIR','reverse')
text1 = sprintf('timen = %d', dt*nt);
text2 = sprintf('simulation time = %f', time_taken);
title({'Transient state heat conduction profile Explicit'; text1; text2})

Output:

• From the above simulation, we can see that the explicit method converges faster than the implicit methods, which take less simulation time.
• For the steady state, we can opt for the SOR method because it reaches convergence faster than any other methods.
• Explicit methods produce greater accuracy with less computational time required to run it.

References:

Solving of Two-Dimensional Unsteady-State Heat-Transfer Inverse Problem Using Finite Difference Method and Model Prediction Control Method (hinda