Effect Of Time Step On Stability Of The Solution For 2D Heat Conduction Equation

Aim - Effect of time step on stability of the solution for both explicit and implicit methods.

In discretizing the governing equation forward time central spacing has been used.

Explicit method is a conditionally stable i.e for a given value of `Delta x , Delta t` is chosen such that it satisfies the courant-Friedrichs-Lewy condition. On performing von-neumann stability analysis on the discretized difference equation with equal grid spacing in both x and y directions, it has been found that `C <= 0.25`.

Code for explicit method is given below. Change the time step to know the effect.

% Effect of time step for Explicit method to 2D transient heat conduction equation
clear; close all; clc
nx = 101; % No of Grid points in x-dir
dx = 1/(nx-1);
ny = 101; % No of Grid points in y-dir
dy = 1/(ny-1);
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,-y);
tm = 50; % Simulating for 50 sec
a = 1.e-3; % Thermal diffusivity
dt = 0.02; % Time step
c = a*dt/dx^2; % Courant number
nt = tm/dt; % Number of time steps

u = ones(nx,ny)*298; % Initialising the solution 
% Boundary conditions
u(:,1) = 400; % Left wall
u(:,nx) = 800; % Right wall
u(1,:) = 600; % Top wall
u(ny,:) = 900; % Bottom wall

% Copying the u
ut = u; % For time increment

% Constants in the discritized equation
k1 = a*dt/dx^2;
k2 = a*dt/dy^2;
k = 1-2*k1-2*k2;

% Time loop 
for t = 1:nt
    % Nodal loop
    for i = 2:nx-1
        for j = 2:ny-1
            u(i,j) = k*ut(i,j) + k1*(ut(i-1,j)+ut(i+1,j)) + k2*(ut(i,j-1)+ut(i,j+1));
        end
    end
    ut = u; % Updating for next time step
end
%% plotting
contourf(xx,yy,u,'ShowText','on')
colorbar
yticklabels({'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'})
title(['Explicit, ','Time step : ',num2str(dt),', CFL no : ',num2str(c)])

Plots at different time steps -

If the `C <= 0.25` errors will not grow up as we march forward in time and the solution will be stable. This can be seen in the below two plots.

If time step is further increased such that it exceeds the stability criteria then the solution becomes unstable.

With a small increment in `Delta t` i.e from 0.025 to 0.0251, solution has started becoming unstable.

With further increase in `Delta t` i.e from 0.0251 to 0.0252, solution has blown up.

Implicit method is unconditionally stable i.e for even larger value of `Delta t` implicit method results in a stable solution.

Code for implicit method is given below.

% Effect of time step for Implicit method to 2D transient heat conduction equation
clear; close all; clc
nx = 101; % No of Grid points in x-dir
dx = 1/(nx-1);
ny = 101; % No of Grid points in y-dir
dy = 1/(ny-1);
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,-y);
tm = 50; % Simulating for 50 sec
dt = 0.02; % Time step size
nt = tm/dt; % Number of time steps
u = ones(nx,ny)*298; % Initialising the solution 
% Boundary conditions
u(:,1) = 400; % Left wall
u(:,nx) = 800; % Right wall
u(1,:) = 600; % Top wall
u(ny,:) = 900; % Bottom wall

% Copying the u
uges = u; % Guess values for iteration
ut = u; % For time increment

% Constants in the discritized equation
a = 1.e-3; % Thermal diffusivity
k1 = a*dt/dx^2;
k2 = a*dt/dy^2;
k = 1+2*k1+2*k2;

tol = 1.e-4; % Tolerance

%% Jacobi 
% Time loop
for t = 1:nt
    err = 1; % initializing error at each time step
    % convergence loop
    while err > tol
        % Nodal loop
        for i = 2:nx-1
            for j = 2:ny-1
                u(i,j) = (1/k)*( ut(i,j) + k1*(uges(i-1,j)+uges(i+1,j)) + k2*(uges(i,j-1)+uges(i,j+1)) );
            end
        end
        err = max(max(abs(u-uges)));
        uges = u; % updating for next iteration
    end
    ut = u; % updating for next time step
end

%% plotting
contourf(xx,yy,u,'ShowText','on')
colorbar
yticklabels({'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'})
title(['Implicit, ','Time step : ',num2str(dt)])

Plots at different time steps -

Implicit method with time step of 0.02 sec.

Implicit method with time step of 0.025 sec.

Implicit method with time step of 1 sec.

Implicit method with time step of 2 sec.

It can be seen that implicit method is stable even for larger time step value. Disadvantage with implicit method is it requires higher computation time compared with explicit method.