1. Introduction

2. Finite-Volume Equations

3. Basic Discretisation Schemes

3.1 Central-Differencing Scheme

3.2 Upwind-Differencing Scheme

3.3 Hybrid-Differencing Scheme

4. Higher Order Discretisation Schemes

4.1 Classification of Schemes

4.2 Flux-Limiter Formulation

4.3 Linear Higher-Order Schemes

4.4 Non-linear higher-order Schemes

5. Recommendations

6. Scheme Activation

7. Schemes in BFCs

8. Applications

9. Concluding Remarks

1. INTRODUCTION

An important consideration in CFD is the discretisation of the

convection terms in the finite-volume equations.

The accuracy, stability and boundedness of the solution depends

on the numerical scheme used for these terms.

The default scheme used in PHOENICS for all variables is the

hybrid-differencing scheme (HDS), which employs:

- the 1st-order upwind-differencing scheme (UDS) in high-convection

regions; and - the 2nd-order central-differencing scheme (CDS) in low-convection

regions.

The UDS is bounded and highly stable, but highly diffusive when

the flow direction is skewed relative to the grid lines.

The HDS is only marginally more accurate than the UDS, as the 2nd-

order CDS will be restricted to regions of low Peclet number.

The 2 approaches commonly used to remedy the problem of numerical

diffusion are:

- mesh refinement; and

- the adoption of schemes with an higher order of accuracy than UDS.

For engineering problems, the necessary degree of grid refinement is

generally impractical as the UDS & HDS are sluggish to grid-

refinement tests.

Thus, schemes with higher-order truncation errors than UDS have

been proposed in an attempt to improve resolution.

PHOENICS V3.4 now has an extensive set of higher-order convection

schemes which are valid in all coordinate systems for the staggered-

grid option.

Two classes of higher-order scheme are now available:

LINEAR SCHEMES:

- these include a number of classical schemes, eg QUICK & CDS;

- they offer good resolution but do not guarantee boundedness;

- the boundedness problem means that they may generate unwanted

oscillations around steep gradients, or unacceptable negative

values.

NON-LINEAR SCHEMES:

- these schemes, which include SMART and van Leer's MUSCL, also

offer improved resolution, but employ a non-linear flux limiter

to secure boundedness.

- the flux-limiter may reduce the numerical accuracy of the

solution to some extent.

These schemes may be used for both single- and two-phase flows,

although the facility has yet to be extended to include:

- - the volume fraction equations R1, R2 and RS; and

- - the energy variables TEM1 and TEM2.

The other limitation is that there is no special treatment of the

domain and internal boundaries, at which all higher-order schemes

are reduced to the UDS.

The purposes of this lecture are:

- to outline the existing convection schemes in PHOENICS;
- to describe the new schemes briefly;
- to present some library results validating their use;
- to point out how they were activated;
- to indicate the location of the relevant GROUND coding.

2. FINITE-VOLUME EQUATIONS

The discretised equation for a general specific variable F is:

J_{h} - J_{l} + J_{n} - J_{s} + J_{e} - J_{w}

+ D_{h} - D_{l }+ D_{n} - D_{s} + D_{e} - Dw =
S_{p} (1)

where S_{p} is the source term for the control volume p, and J_{f} and

D_{f} represent, respectively, the convective and diffusive fluxes

of F across the control-volume face f (f=h,l,n,s,e or w).

The convection fluxes through the cell faces are calculated as:

J_{f }= C_{f}*F_{f} (2)

where C_{f} is the mass flow rate across the cell face f.

The convected variable F_{f} is stored at the cell
centres, and so

its value must be determined by interpolation, i.e. from a scheme.

In PHOENICS, the various higher-order schemes are implemented by

maintaining the UDS discretisation and adding extra higher-order

terms in the form of a source term.

F_{f} can be explicitly formulated in terms of its
neighbouring nodal

values by a functional relationship of the form:

F_{f} = P(F_{nb})
(3)

where F_{nb} denotes the neighbouring-node F values.

From equations (1) through (3), the discretised equation becomes:

{ D_{h} + C_{h}*[P(F_{nb})]_{h}
} - { D_{l} + C_{l}*[P(F_{nb})]_{l}
} +

{ D_{n} + C_{n}*[P(F_{nb})]_{n}
} - { D_{s} + C_{s}*[P(F_{nb})]_{s}
} +

{ D_{e} + C_{e}*[P(F_{nb})]_{e}
} - { D_{w} + C_{w}*[P(F_{nb})]_{w}
} = S_{p} (4)

The higher-order (HO) schemes are introduced into (4) by using

the deferred-correction procedure, whereby:

F_{f} = F_{f}(U)
+ F_{f'} (5)

with

F_{f'} = F_{f}(H) - F_{f}(U) (6)

Here F_{f'} is a HO correction which represents the
difference between

the UDS value F_{f}(U) and the HO value F_{f}(H).

If eqn (5) is substituted into eqn (4), the resulting discretised

equation is:

{ D_{h} + C_{h}*F_{h}(U) } - { D_{l}
+ C_{l}*F_{l}(U) } +

{ D_{n} + C_{n}*F_{n}(U) } - { D_{s}
+ C_{s}*F_{s}(U) } +

{ D_{e} + C_{e}*F_{e}(U) } - { D_{w}
+ C_{w}*F_{w}(U) } = S_{p} + B_{p}
(7)

The deferred-correction source term, Bp, is given by:

B_{p} = C_{l}*F_{l'} - C_{h}*F_{h'} + C_{s}*F_{s'}
- C_{n}*F_{n'} + C_{w}*F_{w'} - C_{e}*F_{e'}
(8)

This treatment leads to a diagonally-dominant coefficient matrix

since it is formed using the UDS.

If eqn (8) is expanded in terms of nodal values, the final form of

the discretised equation is:

a_{p}*F_{p} = sum (a_{nb}*F_{nb}) + S_{p} + B_{p} (9)

where: ap and anb are the convection-diffusion coefficients obtained

from the UDS: F_{p} is the cell-average value of F stored at the cell

centre; and the summation is over the immediate neighbouring nodes

nb (=L,H,S,N,W and E).

All of the schemes provided in PHOENICS calculate the cell-face

values F_{f} using at most two upwind cell-centre
values and one

downwind cell-centre value, as shown in this Figure .

This stencil involves the upstream, central and downstream grid

points, designated by u, c and d respectively.

The most natural interpolation assumption for the cell-face

value of the convected variable Ff would appear to be the CDS,

which calculates the cell-face value from:

F_{f} = 0.5*(F_{c} + F_{d}) (10)

This scheme is 2nd-order accurate, but is unbounded so that

unphysical oscillations appear in regions of strong convection

and also in the presence of discontinuities, such as shocks.

The CDS may be used directly in very low Reynolds-number flows

where diffusive effects dominate over convection.

The UDS (see Courant et al [1952]) assumes that the convected

variable at the cell face f is the same as the upwind cell-centre

value:

F_{f} = F_{c}
(11)

The UDS is unconditionally bounded and highly stable, but it is

only 1st-order accurate in terms of truncation error.

The scheme is therefore highly diffusive when the flow direction

is skewed relative to the grid lines.

3.3 Hybrid-Differencing Scheme

The HDS (Spalding [1972]) switches between CDS and UDS according

to the local cell Peclet number, as follows:

F_{f} = 0.5*(F_{c}
+ F_{d}) for Pe < 2 (12)

F_{f} = F_{c} for
Pe > 2

The cell Peclet number is defined as:

Pe = r*abs(U_{f})*A_{f}/D_{f }(13)

in which A_{f}=cell-face area and D_{f}=physical diffusion coefficient.

When Pe > 2, CDS calculations tend to become unstable so that the

HDS reverts to the UDS. Physical diffusion is ignored when Pe > 2.

The HDS scheme is marginally more accurate than the UDS, because

the 2nd-order CDS will be used in regions of low Peclet number.

4.1 Classification of Higher-Order Schemes

LINEAR SCHEMES are those whose coefficients are not direct

functions of the convected variable when applied to a linear

convection equation.

Linear convection schemes of 2nd-order accuracy or higher may

suffer from unboundedness and are not unconditionally stable.

NON-LINEAR SCHEMES analyse the solution within the stencil and

adapt the discretisation to avoid any unwanted behaviour, such

as unboundedness.

These two types of scheme may be presented in a unified way by

use of the FLUX-LIMITER formulation.

4.2 Flux-Limiter Formulation

The FLUX-LIMITER formulation calculates the cell-face value of

the convected variable from:

F_{f} = F_{c} +
0.5*B(r)*(F_{c}-F_{u})
(14)

where B(r) is termed a limiter function, and the gradient ratio r

is defined as:

r = (F_{d}-F_{c})/(F_{c}-F_{u}) (15)

The generalisation of this approach to handle non-uniform meshes

has been given by Waterson [1994].

From equation (14) it can be seen that B(r)=0 gives the UDS, and

B(r)=r gives the CDS.

4.3 Linear Higher-Order Schemes

PHOENICS provides the following linear higher-order schemes:

- CDS - Linear-upwind scheme (LUS)
- Quadratic-upwind scheme (QUICK) - Fromm's upwind scheme
- Cubic upwind scheme (CUS)

These are unified for implementation purposes as members of the

Kappa class of schemes:

B(r) = 0.5*{(1+K)*r+(1-K)} (16)

where: K = 1 gives CDS, K = 0.5 gives QUICK, K = -1 gives LUS,

K = 0 gives Fromm and K = 1/3 gives CUS.

These schemes are plotted in the Flux-Limiter Diagram (FLD) in the

next panel, which takes the form of a plot of B(r) against r.

The two main regions of the flux-limiter diagram are given by r<0,

indicating an extremum, and r>0 indicating monotonic variation.

Linear flux limiter diagram

4.4 Non-Linear Higher-Order Schemes

The following QUICK-based
non-linear schemes have been provided:

- SMART (bounded QUICK, piecewise linear):

B(r) = max(0,min(2*r,0.75*r+0.25,4)) (17)

- H-QUICK (harmonic based on QUICK, smooth):

B(r) = 2*(r+|r|)/(r+3) (18)

- UMIST (bounded QUICK, piecewise linear):

B(r) = max(0,min(2*r,0.25+0.75*r,0.75+0.25*r,2)) (19)

- CHARM (bounded QUICK, smooth):

B(r) = r*(3r+1)/(r+1)**2 for r > 0

B(r) = 0. for r <= 0 (20)

PHOENICS provides the following Fromm-based
non-linear schemes:

- MUSCL:

B(r) = max(0,min(2*r,0.5+0.5*r,2)) (21)

- van-Leer harmonic :

B(r) = (r+|r|)/(r+1) (22)

- OSPRE :

B(r) = 3*(r^{2}+r)/{2.*(r^{2}+r+1)} (23)

- van Albada :

B(r) = (r^{2}+r)/(r^{2}+1) (24)

The remaining non-linear non-linear schemes in PHOENICS are:

- Superbee : B(r) = max(0,min(2*r,1),min(r,2)) (25)
- Minmod : B(r) = max(0,min(r,1)) (26)

and the following CUS-based flux limiters:

- H-CUS:

B(r) = 1.5*(r+|r|)/(r+2) (27)

- Koren:

B(r) = max(0,min(2*r,2*r/3+1/3,2)) (28)

Classical and CUS-based non-linear
schemes

Most limiters fall into the following 2 categories:

- Polynomial ratio (PR) limiters, which offer the possibility of

smooth, continuous limiter functions without discontinuous

switching, thereby aiding convergence; and

- Piecewise-linear (PL) limiters which switch between linear

schemes so as to produce bounded versions of existing linear

schemes. The disadvantage is that their discontinuous nature may

induce convergence problems.

Limiter functions are designed to fulfil particular boundedness

criteria, usually either the Total Variation Diminishing (TVD)

or Positivity conditions.

All flux-limited schemes provided in PHOENICS are positive, and

the following schemes are TVD: Koren, MUSCL, van Leer harmonic,

Minmod, Superbee and UMIST.

The LINEAR HO schemes offer good resolution, but are unbounded

and may produce unphysical oscillations, which can lead to severe

convergence problems. Therefore, it is recommended that NON-LINEAR

(bounded) schemes always be applied to:

- the momentum equations whenever physical discontinuites are

present, as for example in shock waves; - those turbulence transport equations for which negative values

are unacceptable; and - the species concentrations and enthalpy equations for which

bounded solutions are essential.

For incompressible flows, it is recommended that a LINEAR HO scheme

is applied to the momentum equations, unless severe convergence

difficulties are encountered.

LINEAR HO schemes:

- the CUS is formally the most accurate although QUICK gives

similar results. - the LUS is somewhat less accurate than these schemes, but gives

much better numerical stability.

NON-LINEAR schemes:

- the piecewise-linear kappa-based schemes: SMART, Koren and MUSCL

are likely to give the highest levels of accuracy. - the smooth limiters, e.g OSPRE and van Leer Harmonic, are likely

to give much better convergence at somewhat reduced accuracy. - the classical limiters Minmod and Superbee are not recommended

for general use: Minmod is diffusive and slow to converge;

Superbee is over-compressive, which is excellent for free-

surface scalar markers.

5. RECOMMENDATIONS

It should always be possible to obtain convergence when using

the HO schemes from the very start of the calculation.

However, typically, the false time-step requirements are between

0.01 and 0.1 of the value required by the UDS or HDS.

If convergence proves particularly problematic, then it is

suggested that the user try restarting the calculation from a UDS

or HDS solution.

6. SCHEME ACTIVATION

The default scheme for all variables is the HDS, which is

activated by the setting DIFCUT=0.5.

The UDS is activated for all variables by setting DIFCUT=0.

The HO schemes can be activated from the MENU or Q1; the MENU

has a default HO setting of van-Leer harmonic for all variables.

The PIL command SCHEME is used to select the required HO scheme

for particular variables.

The syntax of the SCHEME command is:

**SCHEME(NAME,variable name 1,variable name 2,...etc.)**

The 1st argument NAME identifies the required scheme, and the

2nd argument permits the user to specify those SOLVEd variables

which will use the selected scheme.

The 1st argument NAME identifies the required scheme, as follows:

(a) Linear schemes

NAME = LUS, FROMM, CUS, QUICK or CDS

(b) Non-linear schemes

NAME = SMART, HQUICK, UMIST, KOREN, SUPBEE, MINMOD, OSPRE,

VANALB, VANL1 (or MUSCL), VANL2 (or VANLH), CHARM,

or HCUS.

For example, the PIL commands:

SCHEME(QUICK,U1,V1)

SCHEME(SMART,H1,C1,C2)

select QUICK for U1 and V1, and SMART for H1,C1 and C2, and UDS

for any SOLVEd variables which do not appear in a SCHEME command.

If ALL is entered as the 2nd argument, then the selected scheme is

applied to all SOLVEd-for variables.

The schemes described thus far are also available for use with the **default
staggered-mesh BFC option.**

The SCHEME command and the associated FORTRAN coding in gxhocs.for are

implemented only for the staggered-grid option in PHOENICS.

The **CCM BFC** **option**
uses its own set of higher-order convection schemes for

single and multi-block meshes. These are activated by using the READ Q1

facility, i.e. by inserting the following commands in the Q1 file:

LSG7=T

SCHMBEGIN

VARNAM UC1 SCHEME MINMOD

VARNAM VC1 SCHEME MINMOD

VARNAM KE SCHEME MINMOD

VARNAM EP SCHEME MINMOD

SCHMEND

**CCM **supports the following schemes:

QUICK - quadratic upwind scheme

SUPERB - superbee scheme

MINMOD

SMART

Please note that the conventional SCHEME COMMAND is not
compatible with

the CCM option.

The **GCV option** has
its own set of higher-order convection schemes for both

single-block and multi-block meshes. They are activated by inserting the

following type of command in the Q1 file:

SPEDAT(SET,GCVSCH,UC1,C,MINMOD)

SPEDAT(SET,GCVSCH,VC1,C,MINMOD)

SPEDAT(SET,GCVSCH,KE,C,MINMOD)

SPEDAT(SET,GCVSCH,EP,C,MINMOD)

etc. The 3rd argument is the variable name and the 5th argument can be

one of:

CDS - central differencing scheme

SOUP - second-order upwind scheme? i.e linear upwind

QUICK - quadratic upwind scheme

MINMOD

MUSCL - van Leer scheme, VANL1 in staggered-mesh schemes

OSHER

SUPERB - superbee scheme

SMART

Please note that the conventional SCHEME COMMAND is not compatible with

the GCV option.

8. APPLICATIONS

The 'Numerical-Algorithms' library contains a number of Q1 files

exemplifying the use of the HO schemes, including:

2d Diagonal Scalar Convection N101

2d Laminar Wall-Driven Cavity N102

2d Laminar Backward Facing Step N103

2d Scalar Convection with Recirculation N104

2d Turbulent Backward Facing Step N105

2d Laminar Flow Over a Thin Fence N108

2d Turbulent Flow through an Orifice Plate N110

1d Shocked transonic flow in a laval nozzle N131

2d Transonic underexpanded free jet

Any of these cases may be loaded from the LIBRARY MENU, or from

the SATELLITE by using,for example, LOAD(N110).

The results of some of these cases may be found in the POLIS

Applications Album, and some will now be presented here.

- Diagonal convection of a scalar
step profile

- Laminar backward facing step
flow

- Laminar wall-driven cavity flow

- Convection of a scalar step in a
recirculating flow

9. CONCLUDING REMARKS

Some recent PHOENICS developments have been presented for HO

convection-discretisation schemes on staggered meshes.

The schemes comprised 5 linear schemes and 12 non-linear schemes.

Some successful applications of these schemes were reported.

More development work needs to be done, including:

- extensions to include R1, R2, RS, TEM1 & TEM2;
- extensions to handle cells adjacent to external and internal

boundaries; - the unification of these schemes with those currently provided

for the co-located multi-block (CCM & GCV) options.

Sources of further information are:

- FORTRAN coding: The GX-file GXHOCS.FOR
- POLIS entry: SCHEMES FOR CONVECTION DISCRETISATION
- Report: N.P.Waterson, VKI Dip. Rep. 1994-33, Belgium, (1994).