MATLAB Scripting & Stability analysis of Linear Convection Equation

Objective

Write a MATLAB Code solving Linear Convection Equation using Finite Difference Method. Study the stability of numerical scheme and analyzing the numerical diffusion error.

Input

1D Domain Length - 1 m

Constant Velocity(C) - 1 m/s

Total Solution Time - 0.4 seconds

Time Step - 0.01 seconds

Number of Grid Points  - 20,40,80,101,120

Initial Velocity Profile - U = 2m/s between x= 0.1 and 0.3 and 1m/s everywhere else.

1D Linear Convection Equation

`(du)/dt + c(du)/dx = 0`

Here the quantity is transported with a constant velocity C.

Any arbitrary initial profile U(X,t=0) is transported with velocity C.

After time ‘t’, the quantity U(X) with the initial profile is displaced with distance C*t

This is a Pure Convection phenomenon whose physical solution is the propagation of the wave.

Discretized form of 1D linear convection Equation using forward differencing in time & explicit backward difference in space.

`((U^(n+1))_i - (U^(n))_i )/(deltat) = c*((U^n)_i - (U^n)_(i-1))/(deltax)`

By doing some rearrangement we can write it as:-

`(U^(n+1))_i = (U^(n))_i - c(deltat)/(deltax)*((U^n)_i - (U^n)_(i-1))`

Now we can solve for time n+1 

MATLAB CODE

We start with writing all the input variables such as the length of the domain, velocity, time step.

The number of grid points will be user input i.e. n =  20,40,80,120

Total Simulation Time is also user input i.e. 0.4 seconds.

Then we will generate a 1D grid using the linspace command.

We will initialize the complete domain with a velocity of 1 m/s.

Then we will use FOR loop and IF-ELSE statement. In FOR Loop we will march through each and every grid point. IF the location of that gridpoint is between 0.1 meters to 0.3 meters we will change the velocity to 2 m/s otherwise velocity will remain as 1 m/s. In this way, we can define our initial velocity profile.

Finally, we will use a nested FOR loop for time and space marching. Inside that, we will solve our discretize equation in the complete domain for the given number of time steps.

%Solving 1D Linear Convection Equation using Finite Backward Difference Scheme

clear al
close all
clc

%inputs
L = 1;        %1D Domain Length 1 m
c = 1;        %convecting velocity = 1 m/s
dt = 0.01;    %time step in seconds

n = input('Enter the number of grid points: ');    %grid/mesh points

total_simulation_time = input('Enter the total simulation time: '); %seconds

nt = total_simulation_time/dt;    %number of time step

%Generating a 1D Grid

x = linspace(0,L,n);
dx = L/(n-1);            %grid spacing

%initilizing velocity as 1 m/s in complete domain
u = ones(1,n);       

%Initial Condition
% u = 2 m/s between x = 0.1 to 0.3 m else u = 1 m/s

for i = 1:n
    if x(i) >= 0.1 && x(i) <= 0.3
        u(i) = 2;
    end
end


u_initial = u;
u_old = u;

CFL_No = c*(dt/dx);  %distance traveled by numerical solution by mesh size

fprintf("The CFL number is: %.2f",CFL_No)

%time marching
for k = 1:nt
    
    %space marching Backward Differencing scheme
    for i = 2:n   
        u(i) = u_old(i) - c*(dt/dx)*(u_old(i) - u_old(i-1));    
    end

    %update old velocities
    u_old = u;

end

%ploting

plot(x,u_initial,'k');

hold on

plot(x,u,'b');
axis([0 1 1 2.1])

xlabel('X-Axis [m]')
ylabel('Velocity [m/s]')

legend('Initial Solution','Final Solution','FontSize',13)

title("1D Linear Convection using FDM ["+ n + " Grid Points]")

text(0.8,1.91,sprintf("Solution Time: %.1fs",total_simulation_time),'FontSize',13)
text(0.8,1.80,sprintf("CFL No: %.2f",CFL_No),'FontSize',13)

if CFL_No > 1
    text(0.8,1.70,"Solution Diverged!",'FontSize',13)
end

Results & Observations

A significant amount of diffusion can be seen in the numerical solution.

With the increase in grid point, numerical diffusion is getting reduced.

Solution diverges for N =  120!

The propagation of the initial profile can be seen but it is also getting diffuse as it propagates further. This diffusion is known as numerical or artificial diffusion which is an error. This is an error because the physics of the Convection Equation should only propagate the wave without causing any diffusion but we can see the profile is getting diffused.

This error is coming because of the numerical approximation or truncating the Taylor series up to a certain order. Since in this case, it is a first-order scheme so the derivative associated with the highest order is a second-order derivative or even derivative whose physical phenomena is diffusive hence our numerical scheme is associated with a numerical diffusion type of error.

Another thing to notice, With an increase in the number of the grid points or by refining the mesh, Artificial diffusion is getting reduced. This is happening because the `deltax` is getting reduced with the increase in the number of grid points so the truncation error will come down. The effect of numerical diffusion will reduce as the product of `deltax` and second-order derivative (truncated term) will come down.

When n = 101, CFL no. become 1 and numerical diffusion goes zero, the solution depicts pure convection phenomena.

With the increase in the number of grid points, Solution accuracy is increasing.

For grid point n = 120, Numerical solutions become unstable and diverge.

Stability Analysis

To understand this we can define a non-dimensional number i.e. CFL no.

CFL No. `= (c*deltat)/(deltax)`

CFL No. describes how fast does the numerical solution is traveling in the mesh.

`c*deltat` is the distance traveled by numerical solution and `deltax` is the length of one cell.

When CFL No. is < or = 1, It means the Numerical solution only travels 1 cell per time step.

When CFL No. is > 1, It means the Numerical solution travels more than 1 cell per time step.

When n = 120, CFL no. is  1.19 hence this makes the numerical solution unstable and the solution blows up.

So with refining the mesh CFL No. is increasing and the numerical solution is becoming unstable. Time step has to be reduced to maintain stability.

By seeing the result, We can define Stability criteria for First Order Explicit Backward Difference Scheme.

When CFL No. < 1, the solution is stable.

When CFL No. > 1, the solution is unstable.