Direct numerical simulation of 1D linear convection equation using MATLAB / Octave

Objective:-> To compare velocity profile of square wave between original and final velocity at different number of nodes.

Theory:- We have 1D linear convection as-

`(delu)/(delt)+ C(delu)/(delx)=0`

The above equation is restricted to 1D only that's why the velocity propagate only in single direction i.e. in X-direction. For discretization we are using first order forward differencing method for time derivative (Time marching) and first order backward differencing scheme for space derivative.

As we know formulae,

First order forward differencing, `(f(i+1)-f(i))/dt`

Hence, for time derivative, first order forward differencing can be, `(u^(t+1)-u^t)/dt`

where, t is superscript and denotes timestep

First order backward differencing, `(f(i)-f(i-1))/dx`

Hence, for space derivative, first order backward differencing can be, `(u_i - u_(i-1))/dx`

Where, i is subscript and denotes number of nodes.

`:.` Our first order 1D linear convection equation can be written as,

`:.` `(u_(i)^(t+1)-u_(i)^t)/dt+C(u_(i)^t-u_(i-1)^t)/dx=0`

`:.` `(u_(i)^(t+1)-u_(i)^t)/dt=-C(u_(i)^t-u_(i-1)^t)/dx`

`:.` `u_(i)^(t+1)-u_(i)^t=-(C.(dt)/(dx)) xx u_(i)^t-u_(i-1)^t`

`:.` `u_(i)^(t+1) = u_(i)^t - C(dt)/(dx) xx (u_(i)^t-u_(i-1)^t)`

Here, `C(dt)/(dx)` is Courant number.

Given data:- Domain length is L = 1m

Constant propagation Velocity, C = 1m/s

Set time step, dt = 0.01 second

Time = 0.4 seconds

Velocity, u = 2m/s for x= 0.1m to 0.3m

Velocity, u = 1m/s for x = 0m to 0.1m and for x = 0.3m to 1m

Number of nodes, n = 20,40,80 and 160

Matlab code:-

clear all
close all
clc

L=1; %length in meter
n=[20 40 80 160]; % array for number of grids (nodes)
dx = L./(n-1); % difference between two successive grids (mesh size)
c=1; % constant propogation velocity in m/s
dt=1e-2; %timestep value
cfl = (c*dt)./dx; %To find courant number
% for 0.4 seconds at 0.01 timestep 
% t=(0.4/dt)
t=40; % number of timesteps 

%grid loop
for o=1:length(n) 
  x= linspace(0,L,n(o));
  u=rand(1,n(o)); %initialization

  %loop for velocity values at diffrent grid positions
  for k=1:length(x)
    
    %loop to define square wave
    if 0.1 <= x(k) && x(k)<= 0.3
      u(k) = 2;
    else
      u(k) = 1;
    end
  endfor 

  uold=u;
  uinit=u;

  %time loop
  for j=1:t

    %space loop  
    for i=2:n(o)
      u(i)= uold(i) - (c*(dt/dx(o))*(uold(i)-uold(i-1)));
    end

    %updating old velocity
    uold = u;
    
  end
  
    %plotting figure
    figure(o);
    
    plot(x,u,'color','b');
    hold on;
    plot(x,uinit,'color','r');
    
    xlabel('Length in meter');
    ylabel('Velocity in m/s^2');
    legend('Final velocity profile','Initial velocity profile')
    %loop for title 
    if o==1
      title(['Number of grids = 20 ; cfl = ',num2str(cfl(1))]);
      axis([0 1 0.7 2.1]);
    elseif o==2
      title(['Number of grids = 40 ; cfl = ',num2str(cfl(2))]);
      axis([0 1 0.7 2.1]);
    elseif o==3
      title(['Number of grids = 80 ; cfl = ',num2str(cfl(3))]);
      axis([0 1 0.7 2.1]);
    else
      title(['Number of grids = 160 ; cfl = ',num2str(cfl(4))]);
    end

end

The Courant–Friedrichs–Lewy i.e. CFL condition is essential condition for CFD problems. It is used for numerical methods that consist time marching phenomena. In short CFL condition expresses that information from a current cell or mesh element must delivered only to its immediate cells or mesh.

Above length v/s velocity plot is generated for 20 grids. For this condition 'dx' = 0.052632.

Above length v/s velocity plot is generated for 20 grids. For this condition 'dx' = 0.025641.

Above length v/s velocity plot is generated for 20 grids. For this condition 'dx' = 0.012658.

Courant number, `K=C(dt)/(dx)`

Courant number is a non dimensional number and can be manipulate by changing the time or space interval. When the space interval is reduced or time interval is increased the courant number increases. In general the explicit methods depend greatly on courant number for their stability. For 1D time marching problem Courant number must be less than one for stable output.

`K<=1`

Hence for above plot with number of nodes 160, the value of Courant number is `(1xx0.01/0.0062893)=1.59` i.e. greater than 1.

Hence we are getting unstable graph for above condition.

Conclusion:- From studying above graphs we can conclude that-

1. The amplitude of wave is continuously decreasing until solution is converged.
2. As the number of nodes increases, dx value decreases and the error between initial velocity curve and final velocity curve decreases.
3. Due to velocity 'C' the square wave is convected.
4. For Courant number < 1 solution is stable & for Courant number > 1 solution gets unstable.