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

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



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.

    Navier-Stokes Equations

    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".
  3. The terms used is the time derivative of velocity of the wave propagating in (+ve) x-axis and the convection term of velocitycomponent.
  4. 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.


clear all
close all

% 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)));
   u_old = u;  
   hold on
   pause(0.003);  % pause at every time interval of 0.003 sec 
                  % axis([0 10 0 0.2]);
   grid on 
   xlabel('Grid Space');
   M(ct) = getframe(gcf);
   ct = ct +1;


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);

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);

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.

Projects by Ankit Chakrabarty

Title:- Geometry cleanup and meshing of an automotive hood for structural analysis   Objective:  Perform geometry cleanup of automotive hood (both inner and outer panel). Perform mid-surface extraction of both the panels and Hem both the surfaces . Perf Read more

Abstract: This project presents the Modelling of the  behaviour of  fluid and its properties when it passes through a convergent-divergent nozzle using a Finite Difference Method. However, to simplify the understanding of the dynamics of the flow in the nozz Read more

 Turbocharger: Literature Review: A turbocharger is mechanical rotating device which forcebly induces air at high pressure into the internal combustion engine, in turn increases the power and effeciency (thermal and volumetric effeciency) of the engine. A turb Read more

  (A) Problem defination:- Solve the following Linear equation using Gauss-Jacobi, Gauss-Seidal and SOR iterative method and understand the following parameters:- Spectral radius and the impact of spectral radius on the convergence of the solution. Effect of d Read more

CASE STUDY:-  A Thermodynamic Analysis Of A Tractor Engine For Use In India With Given Engine Performance Targets. A Brief Literature Review: A tractor is an engineered vehicle specifically designed to deliver Power and traction/tractive effort to haul heavy equi Read more

Solving 2-D 2nd Order Linear Transient &Steady  Heat Conduction Equation Using Both Explicit And Implicit Numerical Schemes In MATLAB.   Explicit Scheme : We say that a scheme is explit when the information at time level "n+1" depends on previous time Read more

  CASE STUDY: Effect of Burn Rate on engine performance Of A 6-Cyl-Turbocharged-DI Engine Using A Non-predictive Combustion model (DI WIEBE) Abstract: Analytical functions approximating the burn rate in internal combustion engines are useful and cost-effective to Read more

Deriving the 4th order approximation for 2nd order derivative of a function using Taylor table for centered, skewed right-sided stencil and skewed left_sided stencil and compare the error results for the same using MATLAB.   Given function, f(x)    =& Read more

Objective:- (1) Determine the engine parameters running at 1800 rpm and record the same. (2) Further increase the power output obtained at 3600 rpm by 10%.   Model:-  The GT model Of the 1-cylinder engine is shown below:- Observations (@1800rpm):-  Read more


The End