Numerical solution to the 1D-Linear Convection equation for varying grid points using MATLAB

OBJECTIVE:

To solve the one-dimensional Linear Convection equation and compare the original and final velocity profiles for various grid points using MATLAB.

Project Conditions:

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 and number of grid points should be variables in your calculation

6. Constant Velocity, C = 1

7. Set time step = 0.01

THEORY:

Linear Convection equation is the most accessible equation in CFD, from the Navier Stokes equation we kept only the accumulation and convection terms for the component of the velocity - as we already know, in CFD the variables to be computed are velocities. To make things even simpler, the coefficient of the first derivative of the velocity is constant, making the equation linear.

where,

`(∂u)/(∂t)` = time derivative

`(∂u)/(∂x)` = convective terms

This equation represents the propagation of a wave without change of shape; the speed of the wave is and with the initial condition given by:

The exact solution is

u(x,t)=u0(x−c⋅t)">u(x,t)=u0(x−c⋅t)u(x,t)=u0(x−c⋅t)

This means that on an axis with the positive direction to the right, the wave will move from left to right with the speed c.

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.

The partial derivative can be approximated using the derivative definition, removing the limit operator:

`(∂u)/(∂x)≈(u(x+Δx)−u(x))/(Δx)`

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:

where n+1 and n are two consecutive steps in time and i and i−1 are two neighboring 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:

Here we follow the method of Time Marching, which is well illustated in the image below

PROGRAM:

clear all
close all
clc

%Numerical solution to the one-dimensional Linear Convection equation.
%INPUT Parameter

L = 1;                  %Length of Domain
C = 1;                  %Constant Velocity
n = [20, 40, 80, 160];  %No. of grid points
dt = 0.01;              %Time step
t = 0.4;                %Simulation Time

%No. of time steps
nt = t/dt;

%Define the start and end positions when velocity is 2m/s
x_start = 0.1;
x_end = 0.3;


%Calculation
%Grid Loop
for i = 1:length(n)
    
    %Mesh Generation
    x = linspace(0, L, n(i));
    dx = L/(n(i) - 1);
    
    %Function of start and end values
    n_start(i) = (x_start + dx)/dx;
    n_end(i) = (x_end + dx)/dx;
    
    %Boundary Condition when velocity is 1m/s
    u = ones(1,n(i));
    %Boundary condition when velocity is 2m/s
    u(n_start(i):n_end(i)) = 2;
    
    %Assigning u to u_old
    u_old = u;
    u_initial = u;
    
    %TIME MARCHING
    %TIME LOOP
    for k = 1:nt 
        
        %SPACE LOOP
        for j = 2:n(i)
            %Difference Equation
            u(j) = u_old(j)-(((C*dt)/dx)*(u_old(j) - u_old(j - 1)));
        end
        
        %Updating Old values
        u_old = u;

        
        %Plotting Movement and change in velocity profile 
        figure(i)
        plot(x,u,'b')
        hold on
        plot(x,u_initial,'r')
        grid on
        
        xlabel('Length of Domain')
        ylabel('Velocity')
        title('Movement and change in Velocity profile for various Grid point')
        legend('Final Velocity Profile','Initial Velocity Profile') 
        axis([0 1 0.8 2.2])
        pause(dt)
    end
    
end

Explanation:

• Program is started with commands clear all, close all and clc to establish a fresh start without the manipulation of previous values.
• Inputs parameters are assigned for length, Constant Velocity, Grid points, Time step and Simulation time period.
• Next step we have defined the start and end positions denoting 0.1 & 0.3 at 2m/s and rest of others at 1m/s.
• for loop has been generated to define the various grid values (i.e) 20,40,80,160.
• The values of x, dx, n_start, n_end, boundary condition for u at 2 & 1m/s are assigned here.
• Since the array u contains new values it can be reused for the time & space loop, so we assign u again with old values (i.e) u_old, u_initial.
• Space Loop is now inputed between Time loop in a for loop command with the difference equation within the grid loop.
• Old value of u is been inputed again for plotting our results.
• Plotting is done for both initial and final velocity profile for four of the grid values and given appropriate labels.

OUTPUTS:

Grid Point = 20

Grid Point = 40

Grid Point = 80

Grid Point = 160

CONCLUSION:

From the above plots we can conclude that

At grid point 20

• We can observe that the profile is more like Trapezoidal and as it moves along the right to final velocity it losses its sharpness, gradient and peak velocity decreases looking like a sine wave.
• The reason maybe because of the truncation error occuring due to the solving of difference equation causing diffusion in the plot.

At grid point 40

• As we have increased the grid point the plot looks much better than the previous plot, yet not a perfect square.
• Here the truncation error is less and hence we could clearly see the solutions getting stable.

At grid point 80

• At this grid point we could say that solution is nearly a square as the shape of the wave and amplitude remains constant.
• Diffusion is less causing very less truncation error while the wave moves rightside with constant velocity.

At grid point 160

• Here the initial velocity is a perfect square but as the solution moves right it gets blown up by oscillating above the line of velocities.
• This could have happened because of the wrong choice of scheme leading to a ignored diffusion terms in numerical calculation which may be termed as Artificial Diffusion.