1D - Linear convection varying over time step (dt)

Objective:

Numerical solution to the one-dimensional Linear Convection equation.

1. Assume that the domain length is L = 1m

2. The initial velocity profile is a step function. It is equal to 2m/s between x= 0.1 and 0.3 and 1m/s everywhere else

3. Use first-order forward differencing for the time derivative

4. Use first-order backward/rearward differencing for the space term.

5. Time step should be variables in your calculation

6. Constant Velocity, C = 1

7. Set grid points n = 80

Expected output

1. For C = 1 and for time = 0.4 seconds, compare the original and final velocity profiles.

2. Make the comparison for time step = 1e-4, 1e-3, 1e-2, and 1e-1

Formula of 1D linear convection:

`(delu)/(delt) = -C*(delu)/(delx)`

Firstly, i have created a function to access the formula of 1D linear convection

Matlab code with explaination:

function [x,u_old]=oned_linear_convection(L,n,n_start,n_end,Change_in_velocity,tm,c,dt,z)

x = linspace(0,L,n);    
dx = x(2)-x(1);            %grid spacing
n1 = ((n_start*L)/dx)+1    %starting node where the velocity is changed
n2 = ((n_end*L)/dx)+1      %ending node where the velocity is changed
u = ones(1,n);             %assigining all the values of u to 1
u(n1:n2) = Change_in_velocity;  %assigining the values to u where there is a change in velocity

u_old = u;                 
u_initial =u;              %u_initial is just to show the initial profile

for j=1:tm                  %time loop
for i = 2:n                 %space loop
        u(i) = u_old(i) - (c*(dt/dx)*(u_old(i) - u_old(i-1)));  
end
    u_old = u;              %to update the velocity to u_old
    figure(1)
end
    plot(x,u) 
    hold on
    axis([0 L+0.3 0 Change_in_velocity+0.3])
    xlabel('Domain length')
    ylabel('velocity')
    legend('1e-1','1e-2','1e-3','1e-4')
    
end

now we need to enter the parameters and the function is called

close all
clear all
clc

%The initial velocity profile is a step function. 
%It is equal to 2m/s between x= 0.1 and 0.3 and 1m/s everywhere else
%1D linear convection = (du/dt) + C* (du/dx)

L = 1;                 %length of domain
c = 1;                 %velocity constant
t=0.4;                 %time
dt = [0.1 0.01 0.001 0.0001];             %time step
n = 80;                %number of nodes
n_start=0.1;
n_end=0.3;      
Change_in_velocity=2

for z=1:length(dt)
   tm=t/dt(z);            %time marching steps

   [x1,u1]=oned_linear_convection(L,n,n_start,n_end,Change_in_velocity,tm,c,dt(z),z);

end

Output:

It is found that as a timestep increases diffusion errors in the solution decreases. But when timestep is increasing too much solutions becomes unstable. It can be concluded that making the time step too small makes the simulation run for too long and as a result diffusion error increases with each iteration and the final solution has a larger diffusion error.

Reasons behind such results:

For the simple problem of linear convection of a quantity u, the first order explicit upwind scheme becomes:

where a is the velocity magnitude, Δt is the timestep and Δx is the length between mesh elements.

By performing a Taylor series expansion on this scheme, this scheme introduces what is called a numerical diffusion, or numerical viscosity, equal to:

As expected from a diffusion coefficient, the numerical diffusion coefficient must be positive, otherwise, the solution will grow indefinitely with respect to time, making the numerical scheme unstable.

The right term inside the parenthesis of the above expression is commonly referred to as the Courant number, which is a dimensionless quantity. Therefore, the Courant number can be stated as follows:

It follows from the numerical diffusion coefficient discussion, that for any explicit simple linear convection problem, the Courant number must be equal or smaller than 1, otherwise, the numerical viscosity would be negative:

Therefore, the CFL condition expresses that the distance that any information travels during the timestep length within the mesh must be lower than the distance between mesh elements. In other words, information from a given cell or mesh element must propagate only to its immediate neighbours.

This discussion may only hold true for explicit schemes. Where as, implicit linear convection is unconditionally stable, and also lacks a CFL condition. This happens because implicit schemes use the entire domain to calculate each timestep. However, implicitly calculating (i.e., solving an implicit linear system) each timestep can be very costly. In addition, time accuracy decreases as the timestep increases. Therefore, it is logical is to use an explicit scheme with a lower timestep.

From the above case consider:- a=1; dx=0.0127; dt=1e-1

Courant number = 7.87 (i.e. greater than 1)

As the CFL number is greater than 1, the solution must go unstable and same can be seen from the above plot. The oscilation pattern seen the plot is due to unstablity.

For dt=1e-2, 1e-3, 1e-4

CFL number is less than 1 and the solution is stable as shown in figure.

Though the solution is stable, there is a possibillty of existing of some errors.

If we observe the above plot, though the profile with time step 1e-2,1e-3,1e-4 are stable but they are not of same shape which is due to the artificial diffusion error.