Effect of Time Step on Stability of Numerical Solution for 2D Heat Conduction

Objective

The objective of the problem is to analyze the effect of time step on the stability of the solution for both implicit and explicit schemes of the 2D Transient Heat Conduction Equation.

The stability of a numerical method for analysing an explicit time integration problem depends on the CFL Number or Cournat Number associated with the method in use.

Courant Number (or) CFL Number

CFL Number or Courant Number is a dimensionless number that is derived from the Cournat-Friedrichs-Levy (CFL) condition which states that, the distance that any information travels during the timestep length through the mesh should be less than the distance between mesh elements. In other words it is a measure of how much information travels through a grid cell in a given time step.

Effect of CFL Number on Explicit Scheme

For an explicit linear convection problem, we know that,

`u_(i)^(n+1) = u_(i)^(n) - ((c*dt)/(dx))(u_i-u_i-1)^(n)`

where,

u = Velocity of a particular element (m/s)

c = Convecting Velocity (m/s)

dt = Time step value

dx = Grid size

From the above equation, the CFL Number for a linear explicit convection problem can be written as,

`C = (c*dt)/(dx)`

C = Courant (or) CFL number

The condition for stability for this type of problem is given as,

`C <= 1`

For an explicit diffusion problem (as in the case of the 2D Heat Conduction Problem at hand), the expression changes to,

`T_(i,j)^(n+1) = T_(i,j)^n +alpha*dt((T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(dx^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(dy^2))`

From the above expression,

`CFL_x = (alpha*dt)/ (dx^2)`

and

`CFL_y = (alpha*dt)/(dy^2)`

Thus, for an explicit diffusion problem, the total CFL number should adhere to the given condition as follows,

`CFL = CFL_x + CFL_y <=0.5`

So, `CFL_x <=0.25` and `CFL_y <=0.25`

Thus, the above conditions must be satisfied to obtain a stable solution for the 2D Explicit Transient Heat Conduction Problem.

MATLAB Code to analyze effect of timestep on the stability of solution for the 2D Explicit Transient Heat Conduction Problem

The code below gives a stable solution for the Explicit scheme. Here, the dt values are changed to analyze the effect of variation in time step.

clear all
close all
clc 
% 2D Transient Heat Conduction using Explicit Methods 
%dT/dt = alpha*(((d^2T)/(dx^2))+((d^2T)/(dy^2)))
%Length of Domain in x and y directions (unit square)
Lx=1;
Ly=1;

%No. of grid points
nx=50;
ny=nx;
t=10;

%Creating the mesh
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

%Creating timestep dt
dt=0.01;
nt=t/dt;

%Initial Conditions
T=298*ones(nx,ny);
%Boundary Conditions
T(1,:) = 900; %First row or Bottom row of nodes in the matrix
T(end,:) = 600; %Last row or Top row of nodes in the matrix
T(:,1) = 400; %First column or Left side column of nodes in the matrix
T(:,end) = 800; %Last column or Right side column of nodes in the matrix
%Corner Node Temperatures
T(1,1) = 650; %Bottom left node
T(1,end) = 850 ; %Bottom right node
T(end,1) = 500; %Top left node
T(end,end) = 700; %Top right node
T_old = T;
T_initial = T;
alpha = 1e-2; 
tol = 1e-4;
error = 9e9;
K1 = (alpha*dt)/(dx^2);
K2 = (alpha*dt)/(dy^2);
count=1;
%Code for solving explicitly

for m=1:nt 
while(error>tol)
    for i=2:(nx-1)
        for j=2:(ny-1)
            h = K1*(T_old(i+1,j) - (2*T_old(i,j)) + T_old(i-1,j));
            v = K2*(T_old(i,j+1) - (2*T_old(i,j)) + T_old(i,j-1));
            T(i,j) = T_old(i,j)+h+v;
        end
    end
    error = max(max(abs(T_old-T)));
    T_old = T;
    count = count +1;
    
end
end
%Plotting the Results
[a,b] = contourf(x,y,T,'showText','on');
clabel(a,b);
colorbar;
colormap(jet);
xlabel('X Axis');
ylabel('Y Axis');
text = sprintf('2D Transient Heat Conduction (Explicit) \n Time Step: %d \n No. of iterations: %d \n CFL Number: %d',dt,count,(K1+K2));
title(text);

Plotting the Results

Case 1: dt = 0.001 secs

CFL Number = 0.048

Case 2: dt = 0.01 secs

CFL Number = 0.4802

Case 3: dt = 0.010425

CFL Number = 0.5

Case 4: dt = 0.01045

CFL Number = 0.501

Inference

Thus, the solution becomes unstable when the CFL Number becomes greater than 0.5, which occurs at time step dt = 0.01045 seconds, for a grid space of 0.0204

From the plots and the formula for the CFL Number it is clear that, for a constant value of Thermal Diffusivity (alpha) and Grid Spacing/Size (dx and dy), the time step size is directly proportional to the CFL Number. Thus, increasing the time step size, will evidently lead to an increase in the CFL Number, which leads to an unstable solution.

Effect of CFL Number on Implicit Scheme

MATLAB Code to analyze effect of timestep on the stability of solution for the 2D Implicit Transient Heat Conduction Problem

The code below gives a stable solution for the Explicit scheme. Here, the dt values are changed to analyze the effect of variation in time step.

clear all
close all 
clc

% 2D Transient Heat Conduction using Implicit Methods 
%dT/dt = alpha*(((d^2T)/(dx^2))+((d^2T)/(dy^2)))
%Length of Domain in x and y directions (unit square)
Lx=1;
Ly=1;

%No. of grid points
nx=50;
ny=nx;
t=10;

%Creating the mesh
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

%Creating timestep dt
dt=0.01;
nt=t/dt;

%Initial Conditions
T=298*ones(nx,ny);
%Boundary Conditions
T(1,:) = 900; %First row or Bottom row of nodes in the matrix
T(end,:) = 600; %Last row or Top row of nodes in the matrix
T(:,1) = 400; %First column or Left side column of nodes in the matrix
T(:,end) = 800; %Last column or Right side column of nodes in the matrix
%Corner Node Temperatures
T(1,1) = 650; %Bottom left node
T(1,end) = 850 ; %Bottom right node
T(end,1) = 500; %Top left node
T(end,end) = 700; %Top right node
T_old = T;
T_initial = T;
alpha = 0.001; 
tol = 1e-4;
error = 9e9;
K1 = (alpha*dt)/(dx^2);
K2 = (alpha*dt)/(dy^2);
a = (1+(2*K1)+(2*K2))^(-1);
t1 = K1*a;
t2 = K2*a;

%Code for solving Implicitly
method = input ('Solver Type : (1 - Jacobi 2 - Gauss Seidel 3 - SOR Jacobi 4 - SOR Gauss Seidel)');
if(method == 1) %Solving using Jacobi Method
    jcount = 1;
    while (error>tol)
        for m=1:nt
            for i=2:(nx-1)
                for j=2:(ny-1)
                    T(i,j) = a*T_old(i,j)+t1*(T_old(i-1,j)+T_old(i+1,j))+t2*(T_old(i,j+1)+T_old(i,j-1));
                end
            end
       
            error = max(max(abs(T_old-T)));
            T_old = T;
            jcount = jcount+1;
        end
        T_initial = T;
    end
%Plotting the results - Implicit Jacobi Method
[a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
xlabel('X Axis');
ylabel('Y Axis');
text = sprintf('2D Transient Heat Conduction using Implicit Jacobi Method \n CFL Number: %d \n Time Step: %d \n No. of iterations: %d',(K1+K2),dt,jcount);
title(text);

else if(method == 2) %Solving using Gauss Seidel Method
        gscount =1;
        while(error>tol)
            for m=1:nt
                for i=2:(nx-1)
                    for j=2:(ny-1)
                         T(i,j) = a*T_old(i,j)+t1*(T(i-1,j)+T_old(i+1,j))+t2*(T_old(i,j+1)+T(i,j-1));
                    end
                end
                 error = max(max(abs(T_old-T)));
                 T_old = T;
                 gscount = gscount+1;
            end
            T_initial = T;
        end
%Plotting the Results - Implicit Gauss Seidel Method
[a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
xlabel('X Axis');
ylabel('Y Axis');
text = sprintf('2D Transient Heat Conduction using Implicit Gauss Seidel Method \n No. of iterations: %d',gscount);
title(text);

else if (method == 3) % Solving using Implicit SOR Jacobi Method
        w = input ('Relaxation Parameter:');
  sor_jcount=1;
  while(error>tol)
      for m=1:nt
          for i=2:(nx-1)
              for j=2:(ny-1)
                  h =  a*T_old(i,j)+t1*(T_old(i-1,j)+T_old(i+1,j))+t2*(T_old(i,j+1)+T_old(i,j-1));
                  T(i,j) = (1-w)*T_old(i,j) + w*h;
              end
          end
          error = max(max(abs(T_old - T)));
          T_old = T;
          sor_jcount=sor_jcount+1;
      end
      T_initial = T;
  end
  %Plotting the Results - Implicit SOR Jacobi Method
  [a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
xlabel('X Axis');
ylabel('Y Axis');
text = sprintf('2D Transient Heat Conduction using Implicit SOR Jacobi Method \n Relaxation Parameter: %d \n No. of iterations: %d',w,sor_jcount);
title(text);

    else %Solving using SOR Gauss Seidel Method
        w = input('Relaxation Parameter :');
        sor_gscount = 1;
        while(error>tol)
            for m=1:nt
                for i=2:(nx-1)
                    for j=2:(ny-1)
                        h = a*T_old(i,j)+t1*(T(i-1,j)+T_old(i+1,j))+t2*(T_old(i,j+1)+T(i,j-1));
                        T(i,j) = T_old(i,j)*(1-w) + w*h;
                    end
                end
                error = max(max(abs(T_old-T)));
                T_old = T;
                sor_gscount = sor_gscount+1;
            end
            T_initial = T;
        end
         %Plotting the Results - Implicit SOR Gauss Seidel Method
  [a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
xlabel('X Axis');
ylabel('Y Axis');
text = sprintf('2D Transient Heat Conduction using Implicit SOR Gauss Seidel Method \n Relaxation Parameter: %d \n No. of iterations: %d',w,sor_gscount);
title(text);
    end
    end
end

   

Plotting the Results

Case 1: dt = 0.01

CFL Number = 0.048

Case 2: dt = 0.1

CFL Number = 0.48

Case 3: dt = 1

CFL Number = 4.8

Inference

From the above plots, we can see that increasing the time step does not have any effect on the stability of the solution. This implies that Implicit Methods are inherently stable irrespective of variation in time steps or CFL number. Thus, larger time steps can be used in combination with implicit schemes to reduce computational efforts without affecting the stability of the solution.