Solving 1D Linear convection equation using MATLAB

Objective

> Write a Matlab code to simulate a One-dimensional Linear Convection equation.

Conditions

> 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

> length is L = 1m

> First-order forward differencing for the time derivative

> First-order backwards/rearward differencing for the space term.

> Constant convection coefficient, c = 1 m/s

> Time step, dt = 0.01

> Simulation run time, t = 0.4 seconds.

> To Perform the simulation for the Number of nodes, n=20,40,80 and 160.

Governing equation

> One-dimensional Linear Convection equation is 

    `(delu)/(delt)+c(delu)/(delx) = o`

    where,

    `u` = Velocity

    `c` = Convection coefficient

> As per the condition, Forward differencing is done for time derivative and Rearward differencing is used for space derivatives

    `(delu)/(delt) = (u_i^(n+1) - u_i^n)/(Deltat)`

    `(delu)/(delx) = (u_i^n - u_(i-1)^n)/(Deltax`

> Substituting the time and space derivatives in the 1D linear convection equation we get

    `(u_i^(n+1) - u_i^n)/(Deltat)+c(u_i^n - u_(i-1)^n)/(Deltax) = 0`

    `u_i^(n+1) = u_i^n - (c(Deltat)/(Deltax))*((u_i^n - u_(i-1)^n)/(Deltax))`

> Courant-Friedrichs-Lewy condition to check stability of the above equation

    `C = c*(dt/dx)<=1`

    where, 

    `C` = Courant number

    `c` = Constant convection coefficient

    `dt` = Time step

    `dx` = Size of mesh

Programs

Function program

% Function to get the node points of the step function
function[ns,nst] = convec_func(L,n,xs,xst)
dx = L/(n-1);
ns = (xs/dx)+1;
nst = (xst/dx)+1;
end

Program explanation

> When the mesh(x) changes with the number of nodes(n), the start and stop positions of the square wave also changes.

> This function program gives the start and stop positions of the square wave despite the change in the number of nodes.

> The program takes the length(L), nodes(n), start(xs) and stop(xst) positions of the square wave and calculates the nodal start(ns) and nodal stop(nst) values.

> The difference(dx) for the mesh is calculated. The difference value changes according to the number of nodes in the mesh.

> As the value of 'n' get higher the 'dx' reduces

> Then the 'ns' and 'nst' values are acquired using the formula `(x/d)+1`

Main program

% Simulation of One-dimensional Linear Convection equation with different number of nodes
clear all
close all
clc

% Inputs
L = 1;
n = [20 40 80 160];
xs = 0.1;
xst = 0.3;
t = 0.4;
dt = 1e-2;
t_m = t/dt;
c = 1;

for i = 1:length(n)
    % Creating domain length 'x' and dx
    x = linspace(0,L,n(i));
    dx = L/(n(i)-1);    
    % Creating velocity variable 'u'    
    u = ones(1,n(i));    
    % Using function to get the start stop values respective to nodes    
    [start,stop] = convec_func(L,n(i),xs,xst);
    u(start:stop) = 2;
    % Storing the velocity values
    u_old = u;
    u_in = u;
    % Calculating Courant number
    C(i) = c*(dt/dx);
    
    % Time loop
    for m = 1:t_m        
        % Space loop
        for k = 2:n(i)
            u(k) = u_old(k) - (c*(dt/dx)*(u_old(k)-u_old(k-1)));
        end
        
        % Updating velocity profile
        u_old = u;
        
        % Plotting the time marching solution      
        figure(i)
        plot(x,u_in,'b','linewidth',2)
        hold on
        plot(x,u,'k')
        xlabel('Domain Length [m]')
        ylabel('Velocity [m/s]')
        title(sprintf('Simulation of 1D Linear Convection equation with n = %d',n(i)))
        
        if m == t_m
            figure(i+4)
            plot(x,u_in,'b','linewidth',2)
            hold on
            plot(x,u,'-o')
            xlabel('Domain Length [m]')
            ylabel('Velocity [m/s]')
            title(sprintf('Original and final velocity profiles with n = %d',n(i)))
        end
    end
end

Program  explanation

> The length of the step function(L), number of nodes(n), the start(xs) and stop(xst) positions of the step function, the run time(t), the time step(dt), time marching period(t_m) and the constant convection coefficient(c) are given as inputs.

> A loop is created for each value of n

> Inside the loop, a mesh(x) for the length(L) is created according to the nodes(n).

> The distance between nodes(dx) is calculated.

> Then the velocity(u) is defined as 1m/s at all points.

> The function is called to get the start and stop positions of the square wave which is then used to designate the initial velocity condition.

> The initial velocity condition is stored in another variable (u_in) for plotting purposes.

> The Courant number C is calculated.

> A time loop is defined for the period of time marching

> Within the time loop the space loop is specified from 2 to n.

> The derived 1D linear convection equation is used to calculate the solution profile.

> The convected velocity profile acquired is stored for the next iteration till the end of the loop.

> The simulation for each n with its initial condition is plotted accordingly.

> The original and final time marched profile is plotted in separate plots for comparison. 

Results

At n = 20

At n = 40 

At n = 80

At n = 160

Courant number C

for - n=20, n=40, n=80, n=160 respectively are

Conclusion

> From the above plots, it is seen that as the number of nodes increases, the accuracy of the convected solution increases.

> As per the Courant-Friedrichs-Lewy condition, for an equation to be stable the courant number 'C' must be less than or equal to unity. If the Courant number is greater than 1, it means that the information propagates through more than one grid cell at each time step, making the solution inaccurate and potentially resulting in divergence of the solution

> It is clearly seen that when n=80, the value of C is closest to unity and when the value of n=160 the value of C>1.

> Therefore, at n=80 the solution is most accurate and it becomes unstable when n=160.

> A dissipative behaviour occurs in the solutions due to the truncation error caused by the even-order derivatives. The phenomenon is known as numerical dissipation interchangeably artificial viscosity.

> Due to the change in the grid size the peak velocity when n=20 and n=40 get numerically diffused or dissipated rapidly compared to when n=80.

> When the grid size is increased further, the courant number exceeds unity and the solution becomes unstable.