Linear Convection

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

To discretize this equation in space and time we will use the Forward Difference scheme for time derivative and Backward Difference scheme for space derivative.

This approach will be applied both for space and time; for time steps we will use the superscript n and for space steps (grid cells) we will use the subscript i

`(underset(i)u^(n+1)-underset(i)u^n)/(Deltat)-c.(underset(i)u^n-underset(i-1)u^n)/(Deltax)=0`

where n+1 and n are two consecutive steps in time and 'i ' and 'i−1' are two neighbouring points of the discretized x coordinate (grid points). The only unknown in the equation written in this way is un+1i. We must solve for this unknown to get an equation that allows us to advance in time, as follows:

INPUTS:

1) 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) Constant Velocity, C = 1

4) Set time step = 0.01

5) Simulation end time = 0.4seconds

Using forward differencing for time derivative and backwards derivative for space terms.

• We will compare the original and final velocity profiles for various grid scales and points.

1) For n=20.

Final Velocity profile:

INITIAL PROFILE:

Due to the less number of grid points, the velocity profile seems to be decreasing from 2m/s to 1.6m/s. Also the shape of the square wave changes to a sine wave.

2) For n=40.

Final Velocity Profile:

INITIAL PROFILE:

As the number of grid points is increased the final velocity profile improves, but still is decreasing to 1.8m/s from 2 m/s. In this case, too the square wave changes to a sine wave.

3) For n=80

Final Velocity Profile:

INITIAL VELOCITY PROFILE:

As seen from the pictures above, it seems that the number of grid points n=80 has been suitable to maintain the velocity profile throughout the simulation. Still, we can see that the wave profile does not remain the same.

4) For n=160:

Final Velocity profile:

INITIAL VELOCITY PROFILE:

At n=160 grid points, the simulation blows up as it progresses. It shows that there is a threshold to be maintained where, below or above it, the results me not support the outcome of the simulation.

Also, the waves changing their shapes may be due to the elimination of the diffusion terms from the equation and the first-order backward differencing scheme in space creates false diffusion. Therefore, it may be said that to maintain the velocity behaviour over the simulation diffusion terms needed to be introduced.

Following is the code for the Calculation of the equation at n=20. For n=40,80,160 only the 'n' values are changed.

clear all
close all
clc

%Solving 1D linear convection
%Input for the equation
L=1;
c=1;
n=20;

%GRID specification
x=linspace(0,L,n);
dx= L/(n-1);

%Time specification
dt=0.01;
nt=(0.4/dt);

%defining the start and end of wave profile
x_start=0.1;
x_end=0.3;

%grid start and grid end of wave profile
n_start=(x_start/dx)+1;
n_end=(x_end/dx)+1;

%Initialize the wave profile
u= ones(1,n);
u(n_start:n_end)=2;
%Applying BC's
u(1)=1;
uold=u;

%Time loop
for k=1:nt
    
    %space loop
    for i=2:n
        
    u(i)= uold(i)-((c*dt/dx)*(uold(i)-uold(i-1)));
    end
    
    %Updating the velcoties
    uold=u;
    
    plot(x,u,'r');
    grid on
    axis([0 1 1 2]);
    xlabel('Grid space');
    ylabel('VELOCITY');
    title('Linear Convection Problem');
    pause(0.001)
     saveas(gca,'Linear convection_20.jpg');
end