In all finite-volume CFD codes for which cell-centre values of variables are stored, the question arises: what value should be ascribed to the fluid which crosses the cell face?

In the diagram below, values of the variable phi are known for the cell centres W, P and E; but what is the value of phi in the fluid at face w, which travels from cell W to cell P, or from P to W?

            ----------------------------------------------
            |              |              |              |
            |      W      w|      P       |      E       |
            |         -------->           |              |
            |              |              |              |
            ----------------------------------------------

A safe answer, ie one which ensures fairly good solution, is:

phi at w = phiW when the flow is from W to P, but phiP when the flow is from P to W.

In other words, phiw equals the phi on the UPWIND side of the cell face.

This, or rather the so-called "hybrid" variant of it, is what is used in PHOENICS unless any other scheme is switched on by the user.

The hybrid variant uses the upwind scheme only when the Peclet number ( = normal-to-face velocity times inter-node distance divided by diffusivity) exceeds 2.0 . Otherwise, the arithmetic average of phiW and phiP is taken.

The two schemes are called respectively the "upwind discretization scheme" (UDS) and the "hybrid discretization scheme".

Numerical diffusion

However, there is an objection: when the flow direction is diagonal to the grid, cell P receives fluid from both the west and the south cells and so takes up an intermediate value. This intermediate value is then passed on to cell N; and so on.

            |              |      N       |              |
            |              |              |              |
            |              |      ^       |              |
            ----------------------|----------------------
            |              |      |       |              |
            |      W      w|      P       |      E       |
            |           ------->          |              |
            |              |      ^s      |              |
            ----------------------|-----------------------
            |              |      |       |              |
            |              |      S       |              |

The result is that physically-present discontinuities become "smeared" by the numerical procedure.

Remedies

Many means have been proposed for reducing the magnitude of this effect, which for obvious reasons is called "numerical diffusion" or "false diffusion".

Some, such as Raithby's "skew-upwind scheme" , address directly the influence of the diagonality of the flow.

3.1 The general form

All of the formulae will be expressed in terms of the face value, phiw, and the cell-centre values phiWW, phiW and phiP, thus:

phiw = phiW + 0.5 * B * (phiW - phiWW) where B is a function of r, defined as (phiP - phiW)/(phiW - phiWW).

It is assumed that the flow is from left to right. Thus W is on the upwind side of P, and WW is on the upwind side of W. If the flow direction were reversed, phiE would be involved instead of phiWW.

       --------------------------------------------------------
       |              |              |              |
       |      WW      |      W      w|      P       |     E
       |              |        ---------->          |
       |              |              |              |
       --------------------------------------------------------

3.2 Particular formulae

(a) Origin of the compilation

The first part of present article is based on the literature review and general study published by Waterson and Deconinck [1995].

These authors distinguish the formulae for phiw into linear and non-linear. This classification is retained here.

The numbers on the left are the values of ISCHEM, the index which is used in PHOENICS to distinguish the schemes.

1. LUS Linear-upwind scheme
2. FROMM FROMM's scheme
3. CUS Cubic-upwind scheme
4. QUICK Quadratic upwind scheme
5. CDS Central-differencing scheme
6. SMART Bounded QUICK
7. KOREN Bounded Cubic, Total-Variation Diminishing
8. VANL1 van Leer MUSCL, bounded FROMM, TVD
9. HQUICK Harmonic QUICK, smooth
10. OSPRE Based on FROMM, smooth/continuous
11. VANL2 van Leer Harmonic, smooth, TVD
12. VANALB Based on FROMM, smooth/continuous
13. MINMOD = Zhu-Rodi SOUCUP scheme, TVD
14. SUPBEE Compressive limiter for sharp gradients, TVD
15. UMIST bounded QUICK, TVD
16. HCUS Harmonic CUS, smooth
17. CHARM Based on QUICK, polynomial ratio

(c) The linear schemes

The five linear schemes, or rather six when the upwind scheme is added, can be compactly expressed as follows, B and r being defined as in 3.1 above.

[Strictly speaking, PHOENICS uses Waterson's (1994) generalization of the definition of r, which is valid for non-uniform and BFC grids.]

UDS upwind difference B(r) = 0
CDS central difference B(r) = r

The twelve non-linear schemes can be expressed as follows.

SMART bounded QUICK B(r) = max(0, min(2*r, 0.75*r+0.25, 4))
KOREN B(r) = max(0, min(2*r, 2*r/3 + 1/3, 2))
VANL1 (or MUSCL) B(r) = max(0, min(2*r, 0.5 + 0.5*r, 2))
HQUICK harmonic QUICK B(r) = 2*(r+|r|)/(r+3)
OSPRE B(r) = 3*(r**2 + r)/{2.*(r**2 + r + 1)}
VANL2 (or VANLH) B(r) = (r + |r|)/(r + 1)
VANALB B(r) = (r**2 + r)/(r**2 + 1)
MINMOD B(r) = max(0, min(r, 1))
SUPBEE Superbee B(r) = max(0, min(2*r, 1), min(r, 2))
UMIST bounded QUICK B(r) = max(0, min( 2*r, 0.25 + 0.75*r, 0.75 + 0.25*r, 2))
HCUS B(r) = 1.5*(r + |r|)/(r + 2)
CHARM bounded QUICK B(r) = r*(3r + 1)/(r + 1)**2 for r > 0 B(r) = 0. for r <= 0

4. Implementation in PHOENICS

(a) What the user must do

When the user has decided which discretization scheme he wishes to use for each solved-for variable, he should either:-

• make his selection via the General Menu, or
• insert the appropriate SCHEME command in the Q1 file as follows:

SCHEME(scheme_name, variable_name 1, variable_name 2, ...etc.)

That is all, because everything else is performed automatically, including, to mention those things which become visible, the creation of the patch command:

PATCH(HOCS,CELL,1,NX,1,NY,1,NZ,1,LSTEP)

and corresponding COVALs, for example:

The coding which is activated thereby, when EARTH execution starts, is in the open-source GXHOCS.F file, which forms part of the advanced-numerical-algorithms option.

The file contains subroutine GXHOCS and its ancillary subroutines GXHOEN, GXHOHL, GXFLPS and BLKSLD. (The acronym HOCS stands for "higher-order convection schemes".)

This coding ensures that extra "source" terms are supplied, at each relevant cell and for each relevant variable, to make up the differences between the effects of the in-EARTH default convection fluxes, and those which are appropriate to the selected schemes.

One detail to note is that, when ANY use of schemes is made, the in-EARTH default becomes the upwind rather than the hybrid scheme. This is effected by the automatically-made setting: DIFCUT = 0.0 .

In Group 1 Section 1, GXHOCS forms the array INLCS(NPHI) and allocates storage for several geometrical quantities. The array INLCS contains an integer which defines for each variable the scheme selected by the user in Q1 via the SCHEME PIL command.

In Group 1 Section 2, GXHOCS computes and stores the position vectors for all cells in the solution domain, as these are required in Group 13 for BFC calculations.

In Group 13, the deferred-correction source terms are calculated, in Section 13 for linear schemes, and in Section 14 for non-linear schemes. Subroutine GXHOEN computes the east/west and north/south cell-face contributions to the source term, and Subroutine GXHOHL computes the low and high cell-face contributions.

Subroutine GXFLPS is called directly from Subroutines GXHOEN and GXHOHL so as to calculate the flux-limiter function for the non- linear schemes.

Function BLKSLD is used by subroutines GXHOEN and GXHOHL so as to detect the presence of blocked faces and solid cells within the stencil used to calculate the higher-order correction for a single-cell face.

(a) An example

The flow which has just been discussed, although very simple, is an important test of the schemes; for sources of one kind or another abound in fluid-flow phenomena.

The schemes which, on the basis of the test, should NOT be used when sources are present, are, of the linear schemes:-

LUS, FROMM, CUS, CDS and QUICK, ie all of them; and of the non-linear ones:-

OSPPRE and VANALB.

These disqualify themselves primarily because they produce out-of- bounds values, which often lack physical meaning.

When the schemes are applied to problems WITHOUT sources, all of them perform better than the UDS, if the still-present out-of-bounds values are tolerable.

PHOENICS Library Case N101, and others, show this.

Some of the non-linear schemes have been reported as sometimes leading to convergence difficulties [Waterson and Deconinck, 1995]. MINMOD has been particularly mentioned in this regard.

1. Since sources are so prevalent in CFD, it appears wise NEVER to use any of the linear schemes.

   They may be satisfactory in some no-source circumstances; but there is nothing to suggest that they are better than those non-linear schemes which also handle sources well.

2. Use SUPERBEE only for variables for which interfaces must remain sharp.
3. Unless other problem-specific advice points in another direction, use SMART.
4. Remember that, precisely because they ARE non-linear, what works well in some cases may work badly in others.
5. Take note of the fact that CLDA and its successor X-Cell, are likely to furnish more-satisfactory long-term solutions.

7. References

1. R.K.Agrawal, 'A third-order accurate upwind scheme for Navier- Stokes at high Reynolds numbers', AIAA 81-0112, (1981).
2. R.Courant, E.Isaacson & M.Rees, 'On the solution of non-linear hyperbolic differential equations by finite differences', Comm. Pure Appl. Maths, 5, p243, (1952).
3. P.H.Gaskell and A.K.C.Lau, 'Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm, Int.J.Num.Meth.Fluids, Vol.8, p617, (1988).
4. A.Harten, 'On a class of High Resolution Total-Variation -Stable Finite-Difference Schemes', SIAM J.Num.Analysis, 21, p1, (1984).
5. J.E.Fromm, 'A method for reducing dispersion in convective difference schemes', J.Comp.Phys., Vol.3, p176, (1968).
6. C.Hirsch, 'Numerical computation of internal and external flows', Computational Methods for Inviscid and Viscous Flows, Vol.2,Wiley Interscience, (1990).
7. B.Koren, 'A robust upwind discretisation method for advection, diffusion and source terms', In: Numerical Methods for Advection-Diffusion Problems, Ed. C.B.Vreugdenhil & B.Koren, Vieweg, Braunschweig, p117, (1993).
8. P.K.Khosla & S.G.Rubin, 'A diagonally dominant second-order accurate implicit scheme', Computers & Fluids, Vol.2, p207, (1974).
9. B.P.Leonard, M.A.Leschziner & J.McGuirk, 'The QUICK algorithm: a uniformly 3rd-order finite-difference method for highly- convective flows', Proc. 1st Conf. on Numerical Methods in Laminar & Turbulent Flow, Swansea, p807, (1978).
10. B.P.Leonard, 'A stable and accurate convective modelling procedure based on quadratic upstream interpolation', Comp. Methods Appl. Mech. Eng., Vo.19, p59, (1979).
11. B.P.Leonard, 'Locally-modifed QUICK scheme for highly convective 2D and 3D flows', Proc. 5th Conf. Num. Meth. in Laminar and Turbulent Flow, Montreal, p35, (1987).
12. B.P.Leonard, 'Simple high-accuracy resolution program for convective modelling of discontinuities', Int.J.Num.Meth. Fluids, Vol.8, p1291, (1988).
13. F.S.Lien and M.A.Leschziner, 'Upstream monotonic interpolation for scalar transport with application to complex turbulent flows', Int.J.Num.Meth. Fluis, Vol.19, p527, (1994).
14. B.Noll, 'Evaluation of a bounded high-resolution scheme for combustor flow computations', AIAA J., Vol.30, No.1, p64, (1992).
15. H.S.Price, R.S.Varga and J.E.Warren, 'Applications of oscillation matrices to diffusion-convection equations', J. Math. & Phys., Vol.45, p301, (1966).
16. P.L.Roe, 'Finite-volume methods for the compressible Navier- Stokes equations', Proc. 5th Int. Conf. Num. Methods in Laminar and Turbulent Flow, Montreal, Vol. V, p2088, (1987).
17. P.L.Roe, 'Characteristic-based schemes for the Euler equations', Ann.Rev.Fluid Mech., 18, p337, (1986).
18. P.L.Roe, 'A survey of upwind differencing techniques', CFD Lecture Series 1989-04, VKI, Rhode-Saint-Genese, Belgium, (1989).
19. S.G.Rubin and P.K.Khosla, 'Polynomial interpolation method for viscous flow calculations', J.Comp.Phys., Vol.27, p153, (1982).
20. D.B.Spalding, 'A novel finite-difference formulation for differential expressions involving both first and second derivatives', Int.J.Num.Meth.Eng., Vol.4, p551, (1972).
21. S.P.Spekreijse, 'Multigrid solution of the steady Euler equations', PhD Thesis, TUD, Delft, The Netherlands, (1986).
22. P.K.Sweby, 'High resolution schemes using flux-limiters for hyperbolic conservation laws', SIAM J.Num.Anal., 21, p995, (1984).
23. G.D.Van Albada, B.Van Leer & W.W.Roberts, 'A comparative study of computational methods in cosmic gas dynamics', Astron. Astrophysics, Vol.108, p76, (1982).
24. B.Van Leer, 'Towards the ultimate conservative difference scheme II', J.Comp.Phys., Vol.14, p361, (1974).
25. B.Van Leer, 'Towards the ultimate conservative difference scheme V', J.Comp.Phys., Vol.32, p101, (1979).
26. B.Van Leer, 'Upwind difference methods for aerodynamic problems governed by the Euler equations', Lectures in Applied Maths, 22, p327, Am. Math. Soc., (1985).
27. N.P.Waterson & H.Deconinck, 'A unified approach to the design and application of bounded higher-order convection schemes', VKI Preprint 1995-21, (1995).
28. N.P.Waterson, 'Development of bounded higher-order convection scheme for general industrial applications', VKI Project Report 1994-33, (1994).
29. G.Zhou, 'Numerical simulations of physical discontinuities in single and multi-fluid flows for arbitrary Mach numbers', PhD Thesis, Chalmers Univ. of Tech., Goteborg, Sweden (1995).
30. J.Zhu & W.Rodi, 'A low dispersion and bounded convection scheme', Comp.Meth. Appl. Mech. & Engng, Vol.98, p345, (1991).
31. J.Zhu, 'On the higher-order bounded discretisation schemes for finite-volume computations of incompressible flows', Comp. Meth. Appl. Mech. & Engng., Vol.98, p345, (1992).