Encyclopaedia Index



This article explains the motivation and nature of the so-called higher-order schemes which are available in PHOENICS for the simulation of convective transport.

Results are also shown of their application to the simple problem of flow down-stream of a point source in a uniform stream flowing obliquely to the grid.

There are seventeen schemes, in addition to the default "hybrid" one; but consideration of their performance suggests that seven of them should not be used.


  1. The question considered
  2. The upwind answer
  3. Higher-order schemes
    3.1The general form
    3.2Particular formulae

    (a) Origin of the study
    (b) the list of schemes
    (c) The linear schemes
    (d) Non-linear schemes

  4. Implementation in PHOENICS

    (a) What the user must do
    (b) What happens thereafter
    (c) Some details of the GXHOCS coding

  5. Performance of the schemes

    (a) An example
    (b) Discussion
    (c) Further information

  6. Recommendations
  7. References

See also Applications Album entries

1. The question considered

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

The answer is given to this question influences the balance equations for both cell W and cell P.

2. The upwind answer

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.


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.

The present topic

Other authors have sought to find formulae for cell-face values in simpler ways, involving only: the phi values on either side of the face, and one still farther upstream.

It is these formulae which are the subject of the present article.

Following the jargon of the specialist literature, they are called the "higher-order schemes".

3. Higher-order schemes in PHOENICS

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.

There are five linear schemes and twelve non-linear ones.

(b) the list of schemes

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

B(r) = 0.5*{ (1+K)*r + (1-K) } where
LUS linear-upwind K = -1
QUICK quadratic upwind K = 0.5
FROMM Fromm's upwind scheme K = 0
CUS cubic upwind scheme K = 0.3333

(d) Non-linear schemes

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:-

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:


and corresponding COVALs, for example:


(b) What happens thereafter

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 .

The treatment is valid for non-uniform meshes and for body-fitting coordinates [Waterson, 1994].

(c) Some details of the GXHOCS coding

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.

If BLKSLD is .TRUE., then no higher-order correction is calculated and the scheme reduces to the UDS.

5. Performance of the schemes

(a) An example

Case N121 from the PHOENICS Input-File Library enables the performance of all the schemes to be compared.

It concerns the computation of the value of a scalar variable downstream of a point source in a stream of uniform velocity which has a direction diagonal to the grid, from south- west to north-east.

The concentrations of the scalar will be displayed.

The exact solution would consist of a diagonal band, with the value of the scalar equal to 10.5. with zero-value regions on either side,

The two points to examine are:

  1. does the band spread unduly? and
  2. do the values lie within the range from 10.5 to 0?

UDS: The expected numerical diffusion is evident, but the range of values, from 0 to 10.5 , is correct UDS

LUS: Numerical diffusion is low; but a negative value of -0.2 appears, and the maximum is only 7.8 LUS

FROMM: Numerical diffusion is low; but a negative value of - 1.5 appears. The maximum value is nearly correct. FROMM

CUS: Numerical diffusion is low; but a negative value of -2.7 appears, and maximum of 12.7 . CUS

QUICK: Numerical diffusion is low; but a negative value of -3.8 appears, and a maximum of 14.6 . QUICK

CDS: As expected of the central-difference scheme, these results are unsatisfactory from all points of view CDS

SMART: Numerical diffusion is low, and the upper and lower limits are correct. SMART

KOREN: Numerical diffusion is low, and the upper and lower limits are correct. KOREN

VLEER1: Numerical diffusion is low, and the upper and lower limits are correct. VLEER1

HQUICK: Numerical diffusion is low, and the upper and lower limits are correct. hquick

OSPRE: Numerical diffusion is low, but the upper limit is 14.0 . OSPRE

VLEER2: Numerical diffusion is low, and the upper and lower limits are correct. VLEER2

VANALB: Numerical diffusion is low, but the upper limit is 11.7 . VANALB

MINMOD: Numerical diffusion is low, and the upper and lower limits are correct. MINMOD

SUPBEE: Numerical diffusion is very low, and the upper and lower limits are correct. SUPBEE

UMIST: Numerical diffusion is low, and the upper and lower limits are correct. UMIST

HCUS: Numerical diffusion is low, and the upper and lower limits are correct. HCUS

CHARM: Numerical diffusion is low, and the upper and lower limits are correct. charm

(b) Discussion

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:-


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

Therefore, the only acceptable linear scheme is the UDS, despite its numerical-diffusion tendency; whereas most of the non-linear schemes are acceptable.

(c) Further information

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.

SUPERBEE has been commented upon in the same reference as being "over-compressive", in that it tends to undo the effects of physical diffusion as well as that of numerical diffusion.

6. Recommendation

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