Stability analysis in the unsteady problem

For both implicit and explicit methods look at the effect of time step on the stability of the solution.

1. How does an unstable solution look like?

2. At what CFL number does your solution become unstable ?

Answer :

clear all 
close all 
clc

% Inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
alpha = 1.4;

x = linspace(0,L,nx);
y = linspace(0,L,ny);

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

% BC
T_l = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);

T(2:ny-1,1) = T_l;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;

T(1,1) = (T_T + T_l)/2;
T(nx,ny) = (T_B + T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_l + T_B)/2;
dt = 20;


k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 1000;
CFL = alpha*((1/dx^2) + (1/dy^2))*(dt);

[X,Y] = meshgrid(x,y);
told = T;
t_prev = told;
tol = 1e-4;
error = 9e9;


for iterative_solver = 1
    % Jacobi method
  if iterative_solver ==1
     tic;
     jacobi_iter = 1;
     % Time Loop
       for g = 1:nt
           while (error > tol)
               % Nodal Loop
             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) = ((t_prev(i,j)*term1) + (H*term2) + (V*term3));
                end
             end
                error = max(max(abs(told - T)));
                told = T;
                jacobi_iter = jacobi_iter + 1;      
           end
           t_prev = T;
       end
         p = toc;
         figure(1)
         [M,N] = contourf(X,Y,T);
         clabel(M,N);
         colormap(jet);
         colorbar
         title_text = sprintf('No. of iterations(Jacobi) = %d  simulation time = %fs',jacobi_iter,p);
         xlabel('x axis');
         ylabel('y axis');
         title(title_text);
         set(gca,'ydir','reverse')
  end
end

output of Jacobi method

• In implict method higher dt values give stable solution. As some implicit methods are unconditionally stable (i.e. no matter how high is dt it will give stable solution).
• Here we can see that the dt values are higher, whereas in explicit method dt values needs to be smaller, which states that in implicit method the dt value don't matter much. That means stability can be maintained at higher values of dt using fewer time steps and calculated in vey less computational time, as shown in figure below.
• In implicit method, as the dt increases so does the truncation error, which may show less accuracy. Also as the time step increasesit brings in errors and reduces the accuracy.
• For the output shown below, where dt = 40, this output still shows stability with the increment in time taken for the calculation than the one with dt = 20 which is shown above using same iterative method (i.e. Jacobi ).

Explicit method :-

clear all 
close all 
clc

% Inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
alpha = 1.4;
dt = 1/3000;

x = linspace(0,L,nx);
y = linspace(0,L,ny);

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

% BC
T_l = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);

T(2:ny-1,1) = T_l;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;

T(1,1) = (T_T + T_l)/2;
T(nx,ny) = (T_B + T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_l + T_B)/2;

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 800;
CFL = alpha*((1/dx^2) + (1/dy^2))*(dt)

[X,Y] = meshgrid(x,y);
told = T;
tol = 1e-4;

tic;
% Time Loop
for g = 1:nt
    % Nodal Loop
    for i = 2:nx-1
        for j = 2:ny-1
            term1 = told(i,j);
            term2 = k1*(told(i+1,j) - 2*told(i,j) + told(i-1,j));
            term3 = k2*(told(i,j+1) - 2*told(i,j) + told(i,j-1));
            T(i,j) = (term1 + term2 + term3);
        end
    end
    told = T;
end
p = toc;

figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('No. of iterations = %d \t time taken = %fs',nt,p);
xlabel('x axis');
ylabel('y axis');
title(title_text);
set(gca,'ydir','reverse')

Output of Explicit Method :

• Here as we using explicit method, so we need to make sure that CFL (i.e. alpha*dt/dx^2) should be < or = 0.5 to make the solution stable. Here CFL = k1 + k2 because we are dealing with 2d heat conduction equation, as k1 = alpha*dt/dx^2 & k2 = alpha*dt/dy^2.
• As nt = 800 and dt = 1/3000 we get the output stable as CFL is < 0.5.
• In explicit method, time step should be very less to maintain stability which results in higher computing time to make calculations than usual.
• Also the stability depends upon numerical errors i.e. if it amplifies from one martching step to next then the calculations will become unstable and if the numerical error decreases after every next marching step then the calculations will be stable.
• CFL plays a vital role in explicit method for stability.

The following outputs shows the unstable solutions for several CFL values above 0.5 :

• From the above 4 outputs the first one shows the start of unstability as CFL rises above 0.5 (i.e. 0.510).
• Same way it gradually starts getting more unstable with increment in CFL in above pictures that are (CFL = 0.511, 0.513, 0.7).
• So based on above outputs we say that the solution starts to get unstable when CFL > 0.5
• Also the CFL number relies on time step, so if dt increases so does the CFL which makes solution unstable for Explicit Method.