2D HEAT CONDUCTION EQUATION EQUATION - STEADY STATE & UNSTEADY STATE

OBJECTIVE:

To simulate Steady state & Transient state analysis.

To Solve the 2D heat conduction equation by using the point iterative techniques by implementing the following methods;

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

The absolute error criteria is 1e-4

ABSTRACT:

The main objective of this challenge is to solve the 2-D heat conduction equation for steady and unsteady state using iterative techniques such as jacobi method of iteration, gauss seidal method of iteration, and successive over-relaxation method of iteration.

Suppose one has a function u that describes the temperature at a given location (x, y, z). This function will change over time as heat spreads throughout space. The heat equation is used to determine the change in the function u over time. The rate of change of u is proportional to the "curvature" of u.

THE MAIN CODE:

IMPLICIT

% steady state analysis
clear all
close all
clc

%inputs
nx = 10; % number of grid points along x axis

ny = nx ; % number of grid points along the y axis

x = linspace(0,1,nx); % 0 to 1 with nx values in between

dx = x(2)-x(1);

y = linspace(0,1,ny); % 0 to 1 with ny values in between

z = fliplr(y);

f = 1.5;

iter = 1;

t = ones([nx,ny])*298; % initializing temperature gradient

% Applying BC

t(1,1:nx) = 600; % top neighbouring point 
t(ny,1:nx) = 900; % bottom neighbouring point
t(2:ny-1,1) = 400; % left neighbouring point
t(2:ny-1,nx) = 800; % right neighbouring point

% Creating a copy
told = t;
r1 = t;

error = 9; % error value
tol = 1e-4; % tolerance
for iterative_solver = 1:3;
  while(error>tol); %executing the loop until error is lesser than tolerance keep executing the for loop
   for i = 2:nx-1;
    for j = 2:ny-1;
     if(iterative_solver == 1);
        t(i,j)= jacobi(t,told,i,j);
        elseif(iterative_solver == 2)
        t(i,j) = GS(t,told,i,j);
        elseif(iterative_solver == 3)
        t(i,j) = sor(t,told,i,j,f);
      end
     end
    end
 error = max(max(t-told));
        told = t; % to always keep updating the t values if not we will always have errors and while loop will never end
        [u1,v1] = meshgrid(x,z);
        %plot
        
        if(iterative_solver == 1);
        figure (1);
        contourf(u1,v1,t);
        [C,h] = contourf(u1,v1,t);
        colorbar;
        colormap(jet);
        clabel(C,h);
        xlabel('x-axis');
        ylabel('y-axis');
        title_text = sprintf('jacobi iterantion number=%d',iter); % including iter number in the plot title
        title(title_text);
        pause(0.003);
        elseif(iterative_solver == 2);
        
        figure (2);
        contourf(u1,v1,t);
        [C,h] = contourf(u1,v1,t);
        colorbar;
        colormap(jet);
        clabel(C,h);
        xlabel('x-axis');
        ylabel('y-axis');
        title_text=sprintf('GS iterantion number=%d',iter);
        title(title_text);
        pause(0.003);
        else
        figure (3);
        [C,h] = contourf(u1,v1,t);
        colorbar;
        colormap(jet);
        clabel(C,h);
        xlabel('x-axis');
        ylabel('y-axis');
        title_text=sprintf('sor iterantion number=%d',iter);
        title(title_text);
        pause(0.003);
     end
      iter = iter+1;
    end
error = 9;
iter = 1;
t = r1;
told = r1;
end

% Transient state problem(explicit)
clear all 
close all
clc

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

%convective coefficient
c = 1e-4

%NUMBRER OF TIME STPES
 nt = 200;

%TIME STEP
dt = 32;

%CFL NUMBER
k = (c*dt)/dx^2;

t = ones([nx,ny])*673;

%applying the boundary conditions
t(1,1:nx) = 600;
t(ny,1:nx) = 900;
t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
told = t;
[xx,yy] = meshgrid(x,y);
tic
iter = 1
  for p = 1:nt
      for i = 2:nx-1
          for j = 2:ny-1
            t(i,j) = k*(told(i-1,j) + told(i,j-1) + told(i+1,j) + told(i,j+1)) + (1-4*k)*told(i,j);
          end
       end
   
    error = max(abs(told-t));
    told = t;
    iter = iter+1;
  end
 
time = toc;
figure (1);
[C,h] = contourf(xx,yy,t);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unsteady state heat conduction (Explicit Method) \n Computation Time: %0.4f \n CFL Number: %0.4f \n iterations=%d',time,k,iter));
xlabel('xx');
ylabel('yy');

EXPLICIT

% steady state analysis
clear all
close all
clc

%inputs
nx = 10; % number of grid points along x axis

ny = nx ; % number of grid points along the y axis

x = linspace(0,1,nx); % 0 to 1 with nx values in between

dx = x(2)-x(1);

y = linspace(0,1,ny); % 0 to 1 with ny values in between

z = fliplr(y);

f = 1.5;

iter = 1;

t = ones([nx,ny])*298; % initializing temperature gradient

% Applying BC

t(1,1:nx) = 600; % top neighbouring point 
t(ny,1:nx) = 900; % bottom neighbouring point
t(2:ny-1,1) = 400; % left neighbouring point
t(2:ny-1,nx) = 800; % right neighbouring point

% Creating a copy
told = t;
r1 = t;

error = 9; % error value
tol = 1e-4; % tolerance
for iterative_solver = 1:3;
  while(error>tol); %executing the loop until error is lesser than tolerance keep executing the for loop
   for i = 2:nx-1;
    for j = 2:ny-1;
     if(iterative_solver == 1);
        t(i,j)= jacobi(t,told,i,j);
        elseif(iterative_solver == 2)
        t(i,j) = GS(t,told,i,j);
        elseif(iterative_solver == 3)
        t(i,j) = sor(t,told,i,j,f);
      end
     end
    end
 error = max(max(t-told));
        told = t; % to always keep updating the t values if not we will always have errors and while loop will never end
        [u1,v1] = meshgrid(x,z);
        %plot
        
        if(iterative_solver == 1);
        figure (1);
        contourf(u1,v1,t);
        [C,h] = contourf(u1,v1,t);
        colorbar;
        colormap(jet);
        clabel(C,h);
        xlabel('x-axis');
        ylabel('y-axis');
        title_text = sprintf('jacobi iterantion number=%d',iter); % including iter number in the plot title
        title(title_text);
        pause(0.003);
        elseif(iterative_solver == 2);
        
        figure (2);
        contourf(u1,v1,t);
        [C,h] = contourf(u1,v1,t);
        colorbar;
        colormap(jet);
        clabel(C,h);
        xlabel('x-axis');
        ylabel('y-axis');
        title_text=sprintf('GS iterantion number=%d',iter);
        title(title_text);
        pause(0.003);
        else
        figure (3);
        [C,h] = contourf(u1,v1,t);
        colorbar;
        colormap(jet);
        clabel(C,h);
        xlabel('x-axis');
        ylabel('y-axis');
        title_text=sprintf('sor iterantion number=%d',iter);
        title(title_text);
        pause(0.003);
     end
      iter = iter+1;
    end
error = 9;
iter = 1;
t = r1;
told = r1;
end

% Transient state problem(explicit)
clear all 
close all
clc

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

%convective coefficient
c = 1e-4

%NUMBRER OF TIME STPES
 nt = 200;

%TIME STEP
dt = 32;

%CFL NUMBER
k = (c*dt)/dx^2;

t = ones([nx,ny])*673;

%applying the boundary conditions
t(1,1:nx) = 600;
t(ny,1:nx) = 900;
t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
told = t;
[xx,yy] = meshgrid(x,y);
tic
iter = 1
  for p = 1:nt
      for i = 2:nx-1
          for j = 2:ny-1
            t(i,j) = k*(told(i-1,j) + told(i,j-1) + told(i+1,j) + told(i,j+1)) + (1-4*k)*told(i,j);
          end
       end
   
    error = max(abs(told-t));
    told = t;
    iter = iter+1;
  end
 
time = toc;
figure (1);
[C,h] = contourf(xx,yy,t);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unsteady state heat conduction (Explicit Method) \n Computation Time: %0.4f \n CFL Number: %0.4f \n iterations=%d',time,k,iter));
xlabel('xx');
ylabel('yy');

SOLVERS:

JACOBI SOLVER;

function out = jacobi(t,told,i,j)
  out = 0.25*((told(i,j-1)+told(i,j+1)+told(i-1,j)+told(i+1,j)));
end

GAUSS SIEDEL SOLVER;

function out = GS(t,told,i,j)
  out=0.25*(t(i,j-1)+told(i,j+1)+t(i-1,j)+told(i+1,j));
end

SUCESSIVE OVER RELAXATION METHOD;

function out = sor(t,told,i,j,f)
  out =((1-f)*told(i,j)+(f*0.25*(t(i,j-1)+told(i,j+1)+t(i-1,j)+told(i+1,j))));
end

OUTPUT FOR STEADY STATE:

OUTPUT FOR TRANSIENT STATE ANALYSIS:

Now, let us look at the Transient state outputs with explicit method.

The output of iteration with CFL number,

Now, the Transient state simulation with Implicit method,

Output,

RESULT & CONCLUSION:

STEADY STATE - The iteration of jacobi's method had most number of iterations while the SOR method had the least number of iterations.

The SOR method may not be very yielding in vast problems or industry standard problems. As far as industry standard problems go from the internet source it is mentioned that Jacobi and Gauss Siedel methods are fruitful. 

Steady state vs Unsteady state analysis.

In steady state simultion there is no need of the time criteria, as there is no time marching achieved. When there is no time marching achieved, the temperature gradient remains the same throughout the simulation. Also, the number of grid points are the same; nx = ny.

In the Transient state simulation, there is the time marching criteria which can be detected by the presence of a term that denotes time in the equation. Due to which the temperature gradient will keep changing according to the time. There are two ways to do the unsteady state analysis. One is the explicit method, the one with least number of iterations in comparision to the other method, that is the implicit method which takes longer and more number of iterations for the explicit method.

Explicit method is much faster than the Implicit Method.

However, Implicit Method is much accurate, as it uses various iterative techniques to converge solutions from the entire values of the domain.

We can also use the inverse matrix technique in the Implicit method.