Linear Convection - Effect of time step size

Here we need to solve the 1-D linear convection Equation which is given by:

`frac(partialu)(partialt) + cfrac(partialu)(partialx) =0`

Now to solve it numerically, we need to implement some differencing schemes of the numerical techniques to discretize the equation.

The discretization for the first-order time derivative is given by:

`frac(partialu)(partialt) = frac(u^(n+1)(i) - u^n(i))(deltat)`

where, `frac(partialu)(partialt)` = First order time derivative of the velocity u,

           `u^(n+1)(i)` = Velocity profile at (n+1)th time step

           `u^n(i)` = Velocity profile at nth time step`

The discretization of first order spacial derivative in x is given by:

`frac(partialu)(partialx) = frac(u^n(i) - u^n(i-1))(deltax)`

where, `frac(partialu)(partialx)` = First order space derivative of the velocity u`

           `u^n(i)` = velocity at the grid point ( i ) and time step n`

           `u^n(i-1)` = velocity at the grid point ( i -1 ) and time step n`

Discretizing the equation we get,

`frac(u^(n+1)(i) - u^n(i))(deltat) + cfrac(u^n(i) - u^n(i-1))(deltax) = 0`

`frac(u^(n+1)(i) - u^n(i))(deltat) = -cfrac(u^n(i) - u^n(i-1))(deltax)`

`u^(n+1)(i) = u^n(i) -cxxdeltatfrac(u^n(i) - u^n(i-1))(deltax)`

Now the domain Initialisations include :

L = length of the domain = 1m

c = diffusion velocity = 1m/s

Also assume, `c frac(deltat)(deltax)` = `lambda`

This Lambda also termed as the CFL constant plays a very important in governing the stability and convergence of the numerical scheme employed.

The stability of the numerical scheme is governed by CFL Stability Criteria. 

For a 1-D case, the value of the constant lambda, `lambda <= c` to ensure convergence,

                                                                         where C = velocity of the propagation.

This value of the CFL constant also determines the rate of convergence and the accuracy of the numerical scheme. 

The time step and the number of grid points of the domain are kept as user input in the Matlab code :

In this project, the number of the grid points is kept constant and is set to n=80. The Effect of different time step sizes is seen in the stability of the solution.

The code to solve the equation numerically in MATLAB is given below:

clear all;
close all;
clc;

% Defining the input parameters for the linear Convection Equation
L = 1;
c = 1;


n  = input(' Input the number of grid points you want : ');
dt = input(' Input the time step size ( dt ) : ');

% Calculations

x = linspace(0,L,n);
dx = x(2)-x(1);

% x_start = 0.2m and x_end = 0.3m
a = x;
b = x;

% I will do the subtraction of the elements of x with 0.1 and the minimum
% will correspond to the start of the n_index for 1m/s velocity profile
p = abs(x - 0.1);
[ min_value_start , n_start_index ] = min(p);

% I will do the subtraction of the elements of x with 0.3 and the minimum
% will correspond to the end of the n_index for 2m/s velocity profile
q = abs(x-0.3);
[ min_value_end , n_end_index ] = min(q);

% initialising the velocity profile
u = ones(1,n);
u(n_start_index:n_end_index ) = 2;

% creating a copy of u for reference
uold1 = u;
uold  = u;
% Inputing the time till I want the time-steps
t = input('Input the time till you want the velocity profile evaluation : ');

% number of time steps
nt = t/dt;

% first order forward differencing scheme
% creating a video file
v = VideoWriter('Linear challenge 1 with time step dt and n grid points.avi');

% changing the framerate of the video so as to have a deeper insight
v.FrameRate = 10 ;
open(v);

tic
% time loop
for k=1:nt
    
    % Space loop
    for i = 2:n
        
        % Discretization done to solve the equation
        u(i) = uold(i) - c*(dt/dx)*(uold(i) - uold(i-1));
        
    end
    
    % Updating the old velocity profile for the next time step
    uold = u;
    
    % plotting the velocity profile
    plot(x,u,'k','linew',2);
    
    % Specifying the length of the axis
    axis([0 1 -1 3]);
    
    % making a video file and then joining the images together
    frame = getframe(gcf);
    writeVideo(v,frame);

    % Pausing the loop to 0.02 seconds for smooth animation
    pause(0.02)
    
end

close(v);
  toc

 % plotting the final figures
 figure()
 plot(x,uold1,'b','linew',2);
 hold on;
 plot(x,u,'k','linew',2);
 legend('Initial Profile',[ 'Profile after 0.4 seconds for' num2str(n) ' grid points' ])
 title('Linear Convection Plots');
 xlabel('X');
 ylabel(' Velocity ');
 

Also as per the assumption of the perfect square profile, the velocity profile should just travel in the x-direction in the analytical solution which is not happening here because of the initial profile, not being exactly a step function. The reason for this is the way of initialization of the profile. At the jump, the profile should have 2 values that are not done in the MATLAB code. This causes the calculation at the grid points to be inaccurate because of the neighboring points and hence it deforms the initial square profile. The rate of the deformation is directly related to the CFL constant. As the grid points are increased, it causes more accuracy in the calculations as more information is available. Consequently, this results in an increase of the CFL constant. This is because as the number of the grid points is increased the grid spacing between them is decreased and it results in the increase of the constant. The same is true with the increase in the time step size is increased it causes `deltat` to increase and hence increase the CFL constant, `lambda`.

This result is verified in the numerical scheme employed. The extent of deformation of the initial profile is seen with the decrease in the time steps values.

To calculate the elapsed time tic - toc command is used in MATLAB which calculates the time elapsed in between the start of the tic and the end of toc.

1) For the first Input n = number of grid points = 80 and time step size dt = 1e-1

On running the MATLAB code for n = 80 grid points and time step size dt = 1e-1 the simulation is obtained as shown below:

The graph comparing the velocity profile after 0.4 seconds and the initial profile is shown below:

For this scheme `deltat = 0.1` and `deltax = 0.0127`

In this case the CFL constant, `lambda` = `c frac(deltat)(deltax)` = `1 xx 0.1 / 0.0127` = `7.87`

Clearly in this case the CFL constant is much greater than c(which is 1), and hence the numerical scheme is not stable which is seen clearly in the plots also. The solution has blown up as is seen from the values on the y-axis. It has reached an order of 10000. This type of selection of the CFL constant should be avoided as it would result in an unstable solution.

2) For the Second Input n = number of grid points = 80 and time step size dt = 1e-2 

On running the MATLAB code for n = 80 grid points and time step size dt = 1e-2 the simulation is obtained as shown below:

The graph comparing the velocity profile after 0.4 seconds and the initial profile is shown below:

For this scheme `deltat = 0.01` and `deltax = 0.0127`

In this case the CFL constant `lambda` = `c frac(deltat)(deltax)` = `1 xx 0.01 / 0.0127` = `0.787`

Clearly in this case the CFL constant is less than c(which is 1), and hence the numerical scheme is stable which is seen clearly in the plots also. The solution has not blown up as is seen from the values on the y-axis. Also, it is seen that the profile after 0.4 seconds has still retained its shape to some extent as the CFL constant is more in this case.

3) For the first Input n = number of grid points = 80 and time step size dt = 1e-3

 On running the MATLAB code for n = 80 grid points and time step size dt = 1e-3 the simulation is obtained as shown below:

The graph comparing the velocity profile after 0.4 seconds and the initial profile is shown below:

For this scheme `deltat = 0.001` and `deltax = 0.0127`

In this case the CFL constant `lambda` = `c frac(deltat)(deltax)` = `1 xx 0.001 / 0.0127` = `0.0787`

Clearly in this case the CFL constant is less than c(which is 1), and hence the numerical scheme is stable which is seen clearly in the plots also. The solution has not blown up as is seen from the values on the y-axis. Also, it is seen that the profile after 0.4 seconds has lost its shape to some extent as CFL constant is less in this case as compared to the above case.

4) For the first Input n = number of grid points = 80 and time step size dt = 1e-4

On running the MATLAB code for n = 80 grid points and time step size dt = 1e-4 the simulation is obtained as shown below:

The graph comparing the velocity profile after 0.4 seconds and the initial profile is shown below:

For this scheme `deltat = 0.0001` and `deltax = 0.0127`

In this case the CFL constant `lambda` = `c frac(deltat)(deltax)` = `1 xx 0.0001 / 0.0127` = 0.00787

Clearly in this case the CFL constant is less than c(which is 1), and hence the numerical scheme is stable which is seen clearly in the plots also. The solution has not blown up as is seen from the values on the y-axis. Also, it is seen that the profile after 0.4 seconds has still lost its shape to quite some extent as CFL constant is less in this case as compared to the other cases.

THE VELOCITY PROFILE AFTER 0.4 SEC IN A SINGLE PLOT :

1) For dt = 0.01 , 0.001 , 0.0001

2.) For all the dt in a single plot

Some Key Points:

1.) One interesting phenomenon to note here is the deformation of the square wave while propagating forwards. This effect can be attributed to the fact of Numerical Diffusion, a major error source in CFD studies. Numerical diffusion is the tendency for transported variables to diffuse more than they should. Theoretically, the diffusion shouldn't occur as there is no diffusion term. The square wave which was initialized in MATLAB isn't perfectly square, there is some slope in the jump of the wave. Because of this, the solution algorithm diffuses more than it should and we get a deformed profile. If we solve the problem numerically with pen and paper and account for the perfect square wave we will be able to solve with high physical consistency.

2.) Another important fact to note here is that as the time step size is decreased we see more deformation in the square wave. This is because as there are more time grid points numerical diffusion will occur in the provided time.

dt vs simulation time

It is seen that the simulation time increases as the timestep decreases and it makes sense also as more computations are done in order to reach a particular value with smaller increments.

Conclusion: The effect of the time step sizes is studied on the numerical scheme and the results are in agreement with the mathematical formulas and the conditions derived for governing them. The CFL constant also determines the rate of convergence and the accuracy of the numerical scheme. Hence the proper selection of the parameters is a must so that the CFL stability criterion is satisfied in order to solve the problem accurately and efficiently.