Linear Convection for varying number of nodes

Objective:

1. Find Numerical Solution to the One-Dimensional Linear Convection equation.
2. Compare the original and final velocity profiles in a single figure window for Number of nodes, n = 20, 40, 80 and 160.
3. Explain the following phenomena observed from the solutions- numerical diffusion and instability and how they relate to the courant number and the accuracy of the solution. How does a change in grid size affect these phenomena.

Derivation:

For the solution, we have to use First-Order Forward Differencing for the Time Derivative and First-Order Backward/Rearward Differencing for the space term.

One-Dimensional Linear Convection Equation:

∂u∂t+c(∂u∂x)=0

Using Forward Differencing for Time Derivative- ∂u∂t=un+1i−uniΔt

Using Backward Differencing for Space Derivative- ∂u∂x=uni−uni−1Δx

un+1i−uniΔt+c(uni−uni−1Δx)=0

and, un+1i=uni−c.ΔtΔx.(uni−uni−1)

Here n denotes the Time Step and i denotes the node.

Main Code:

% Week 4.1 - Linear Convection for varying number of nodes

clear all
close all
clc

% Inputs for Linear Convection
L = 1;          % Length in meters
c = 1;          % Constant convection coefficient in meters/second
dt = 0.01;      % Time step in seconds
t = 0.4;        % Simulation time in seconds

% Number of nodes- to be used 20, 40, 80 and 160
n = 20;

% Node array
x = linspace(0, L, n);
% Element size
dx = L/(n-1);

% Defining the Initial Wave for time t = 1
x_start = 0.1;      % Start value for u=2m/s
x_end = 0.3;        % End value for u=2m/s
n_start = round((x_start/dx)+1);    % Finding starting node
n_end = round((x_end/dx)+1);        % Finding end node

% Initialising arrays and boundary conditions
u = 0.1 * ones(1, n);
u(n_start:n_end) = 0.2;
u_old = u;          % Assinging u to new variable
u_initial = u;

% Total number of time steps
nt = t / dt;

% for loop for time loop
for  k=1:nt
    % for loop for space loop
    for i = 2:n
        u(i) = u_old(i) - (c*dt/dx)*(u_old(i)-u_old(i-1));
    end
    u_old = u;      % Assigning u_old with new u values
    
    plot(x, u, 'b')
    hold on
    plot(x, u_initial, 'r', 'linewidth',2)
    axis([0 1 0 0.2])
    pause(0.003)
    xlabel('Domain Length')
    ylabel('Velocity')
end

Output plots:

1. Plot for n=20

2. Plot for n=40

3. Plot for n=80

4. Plot for n=160

Result Explanation:

It can be seen for n=20, the highest velocity is lowered to about 0.17m/s, and it goes on increasing for n=40, n=80, and n=160.

But for n=160 we get an unexpected plot due to Numerical Dispersion.

Numerical Dispersion:

Dispersion results in a distortion of the propagation of different phases of wave(as seen in the plot of n=160), which shows up as a wiggle at the start and end of the wave.

This kind of error in solutions occurs due to the approximation in numerical solution(First-Order, Second-Order, etc.)

Numerical Dissipation or Artificial Velocity:

By performing the above operations Numerical Diffusion/Viscosity(Artificial Viscosity) is introduced, which is equal to -

μNUM=c.(ΔxΔ2).(1−c.ΔtΔx)

Here the term, Courant Number

C=c.ΔtΔx

---