Encyclopaedia Index

NUMERICAL CONVECTION SCHEMES IN PHOENICS

Contents

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.

3. BASIC DISCRETISATION SCHEMES

3.1 Central-Differencing Scheme

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.

3.2 Upwind-Differencing Scheme

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. HIGHER-ORDER DISCRETISATION SCHEMES

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²+r)/{2.*(r²+r+1)} (23)

• van Albada :
  

B(r) = (r²+r)/(r²+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.

7. SCHEMES IN BFCs

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