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 shall comprehend the idea of the 2D heat conduction equation under both steady-state and transient flow conditions. Before moving on to the encoding phase, we must first comprehend the derivation phase in which the concepts are incrementally discretized. Before continuing, there are a few points to keep in mind. For the steady state, iterative solvers are required, and for transients, we must solve both of the implicit and explicit schemes from the previous section.
We will now comprehend how the steady-state condition was derived.

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 start the coding portion in MATLAB, where these terms are all taken into account by using for loops and while loops.

% 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.