Modelling Convection in OpenFOAM

As you may know by now, the Navier-Stokes equations are composed from three parts: a time component, a convection component and a diffusion component.
The time component accounts for the change in time of the given scalar or vector; the convection component accounts for the transport of the scalar due to the velocity of the fluid; and the diffusion components determines the transport due to scalar concentration and will make the scalar to spread from high concentration to low concentration even with zero velocity.
In the previous post we've seen that when computing only the convection we had some diffusion also, and said that it was a numerical error.

In this post we will see this in more detail; we'll use the same shockTube problem as before, changed to work for the scalar transport solver, and model the convection only problem for 3 seconds.
For the convection, which mathematically is a divergence we'll use different schemes and see their effects.

The divergence schemes available are documented in OpenFoam2.3.1 UserGuide:

Behavior of interpolation schemes used in divSchemes 

For this post we'll try the following schemes: upwind, linear, linearUpwind and QUICK.
To better understand how each scheme actually works you should read some CFD theory, but I'll try describe here in a few words.
Basically the OpenFOAM (as any other finite volume method code) computes the values of any property at grid cell centers; to determine the needed value actually it accounts on the property value at the cell faces based on the previous iteration; now, the way these face values are determined is the divergence scheme to chose.
Let's assume we have two cells and we want to determine the value for the face between them, as shown bellow (credit for the picture to Jozsef Nagy):

The cells are named by their centers P and N, the face will be f; let's assume that the computed velocity for the last iteration has the direction from P to N. Now we are ready to understand how is the face value determined for each scheme.
If we use the upwind scheme, the face property will have the same value as the property at cell P center.
If we use the linear scheme, the face value will be determined as a linear interpolation between P and N values, weighted by the distance between P and N.
The linearUpwind scheme will determine the face value based on the P value and the gradient between P and N.
And finnaly, the QUICK scheme which only works on hexahedral meshes, is applying a quadratic interpolation in which will be involved three cells.

This being said, lets see the results for the four cases; the order in which the plots are shown, from top to bottom, is the following: upwind, linear, linearUpwind and QUICK.

The first thing to notice is that the upwind scheme has the largest error; to investigate even more, we see that we have some negative values, that shouldn't be there; we may need to change the plotted range.
If we rescale the plots between -0.2 and 0.2 we obtain the following  plots:

This clearly shows that the linear scheme creates some numerical fluctuations that propagate in the domain and may affect the solution; special treatment must be applied to stop this propagation if this scheme is to be used.
As we can see in the first results plot all the schemes are affected by numerical diffusion; but in the last two have this error is smaller; therefore we can conclude that any of the last two can be used, of course knowing  that an error still exists.

This concludes for today, if you have any questions please leave a comment below.