Solving 1-D Linear Convection Using First-Order Backward Difference And Forward Difference method Using MATLAB

LITERATURE REVIEW:

WHAT IS CONVECTION?

Convection is the sum of bulk transport of the fluid and brownian/ osmotic dispersion of fluid constituents from high density region to lower density region.

Convection = Advection + Diffusion.

• Advection: Advection is the motion of the particle along the bulk flow.The equation is:-
• Diffusion: Diffusion is the process when a particle comes in contact with another particle and dissipiates its momentum and energy to another particle while moving along the flow.

• Navier-Stokes Equations

  Navier-Stokes Equations

  In fluid dynamics, the Navier-Stokes equations are equations, that describe the three-dimensional motion of viscous fluid substances. These equations are named after Claude-Louis Navier (1785-1836) and George Gabriel Stokes (1819-1903). In situations in which there are no strong temperature gradients in the fluid, these equations provide a very good approximation of reality.

  The Navier-Stokes equations consists of a time-dependent continuity equation for conservation of mass, three time-dependent conservation of momentum equations and a time-dependent conservation of energy equation. There are four independent variables in the problem, the x, y, and z spatial coordinates of some domain, and the time t.

  As can be seen, the Navier-Stokes equations are second-order nonlinear partial differential equations, their solutions have been found to a variety of interesting viscous flow problems. They may be used to model the weather, ocean currents, air flow around an airfoil and water flow in a pipe or in a reactor. The Navier–Stokes equations in their full and simplified forms help with the design of aircraft and cars, the study of blood flow, the design of nuclear reactors and many other things

1. The 1-D Convection equation which is mentioned below is accesible directly from the NAVIER-STROKE equation.
2. This equation represents "The Propagation of wave without change in shape with a speed of c in 1-D".
4. The terms used is the time derivative of velocity of the wave propagating in (+ve) x-axis and the convection term of velocitycomponent.
5. The initial boundary conditions are:-

• U(x,t =0) = 1m/s
• Speed of Wave, c = 1

Before writing codes, we have considered the following parameters:-

• Domain length = 1m
• Time step, dt = 0.01sec
• Grid size, dx = [0.05 , 0.025, 0.0125, 0.00625].
• Solved on x-axis.
• The initialization is made with a velocity of 1ms/s all over the grid points, except those that lie between x=0.1 and x=0.3; on these points the velocity is 2m/s.

Now, to discretize the equation in space and time, we use "Forward Difference scheme" for the time derivative and "Backward Difference Scheme" for the space derivative.

• The "Forward difference Scheme" is implemented for the time derivative in order to determine the velocity for the current node in the next time step.
• The "Backward Difference Scheme" is implemented for the space derivative because of the direction of the" flow of information" which is flowing from the left (boundary condition) to right.

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

This approach will be implemented on both the time and space derivative.

For time step, we use the subscript (n) and for the grid spacing we use the subscript as (i).

n+1 and n represents two consecutive time steps.

i and i-1 represents two consecutive grid points. The equation is therefore:-

Now, rearranging the above equation can give us velocity for a node in the next time.

The above equation represents an "EXPLICIT-SOLVER" equation as all the term of RHS is known.

MATLAB MAIN SCRIPT:

clear all
close all
clc


% Input for linear convection

 L = 1;      % Length of domain.
 
 n = 40;    % grid points.
 
 c = 1;      % cconstant value in convection term.
 
 dt = 1e-2;  % time step

             % calculation
 
 x = linspace(0,L,n); % Length of mesh
 
 dx = L/(n-1);        % Grid spacing of each mesh.
 
 
                      % Initializing the value of computation.
 
 u =1*ones(1,n);      % velocity row matrix
 
 n_start_1 = 0.1;
 
 x_start =n_start_func(n_start_1,dx);
 
 n_end_2 = 0.3;
 
 x_stop  = n_end_func(n_end_2,dx);
 
 u(x_start:x_stop) = 2;
 
 %u(1) = 1;           % Assigning thevalue of 1st velocity element in the matrix 
                      % (Boundary Condition)
  
 u_old = u;
 
 ct =1;
               
 u_initial = u;
 
 % Time loop( time step)
 % Total time  = dt*40 = 0.4 sec
 
 for k = 1:40
     
        for i = 2:n
     
          
            u(i) = u_old(i) - (c*(dt/dx)*(u_old(i) - u_old(i-1)));
    
       
        end
     
   u_old = u;  
   
   plot(x,u_initial,'r','linewidth',1);
   
   hold on
   
   plot(x,u,'b');
 
   pause(0.003);  % pause at every time interval of 0.003 sec 
   
                  % axis([0 10 0 0.2]);
   grid on 
   
   xlabel('Grid Space');
   
   ylabel('Velocity');
   
   M(ct) = getframe(gcf);
   
   ct = ct +1;
     
 end

RESULTS (PLOTS):

This is the plot of an approximate square- wave profile which will be incorporated with the above 1-D convection equation to analyse and under the behaviour at different grid points over a time step. With more grid points (finer grid size), the profile gets more symmetrical.

The FUNCTION SCRIPTS in MATLAB for above profils are:-

Function Script name is n_start_func(n_start_1,dx)

function output_1 = n_start_func(n_start_1,dx);


output_1 = round((n_start_1)/(dx)+1);
end

This function takes the value of n_start and for each value of dx from the main script and solves the eqn. in the function script and gives the value of the output_1.

Function script name is n_end_func(n_end_2,dx)

function output_2 = n_end_func(n_end_2,dx);


output_2 = round((n_end_2)/(dx)+1);
end

This function takes the value of n_end and for each value of dx from the main script and solves the eqn. in the function script and gives the value of the output_2.

For Grid points = 20

For Grid points 40

For grid Points 80

For Grid Points 160

# Discussion

# From the above plots, it is clear that the wave propagates or is convected towards right hand side.

# The wave is convected due to constant wave speed, c = 1>0 (+ve).

# The wave seems to diffuses down as the time step increases which is known as false diffusion /artificial diffusion. This is caused due to the error generated after using 1st order numerical difference schemes for both the time and space derivatives.

#False diffusion is a type of error which is observed while using upwind schemes to discretize the hyperbolic partial differential equation into a linear equation and approximate the convection equation.

# The error generated is due to neglecting higher order term of 0(h2) in Taylor infinite series which is very similar to the Diffusion terms in Navier-Strokes Equation.

# The finer the grid size / higher mesh density, the lesser diffusion errors are observed and hence the wave profile retains more of it peak velocity about x-axis.

# The false-diffusion error can also be reduced by using power law schemes, quick schemes , exponential schemes.

# While modelling the equation, viscosity term is not considered and thus no shape change is observed near the wall.

# Increasing the nodes as well as time step will make the scheme unstable and will generate ripples as seen in fig-4. Hence, space steps and time steps are critical.