Encyclopaedia Index


This alternative to the staggered-grid / velocity-resolute formulation, which is still the default hydrodynamic algorithm of PHOENICS, is activated by setting CCV=T in the Q1 file.

It forms is part of the Advanced Numerical Algorithms option.

It is employed in several cases supplied with the Input File Library of that option.

Its nature and mode of operation can be learned by study of:



  1. Introduction
  2. Mathematical basis of the solution of Navier-Stokes equations on non- staggered grid.
  3. Practical details

1. Introduction

These notes describe the non - staggered calculations methods of fluid flow phenomena which have been used and developed for PHOENICS. It also demonstrates a main points for implementing non- staggered calculations in computer code. The specific aims of these notes are: to describe the existing non-staggered technique, to explain how it can be used in general purpose code and to present the main points of implementation procedure.

2. Mathematical basis of the solution of Navier Stokes equations on non- staggered grid.

2.1 Introduction

In this section, the mathematical basis for the non- staggered calculations is described. The terms in partial differential equations for co- located cartesian velocities for each velocity component are laid out and the corresponding difference expressions are represented following the PHOENICS conventions.

The Rhie and Chow velocity/ pressure interpolation treatment is used to close the set of governing equations.

2.2 Model equations.

For non- staggered calculations the momentum conservation equtions are the same as for staggered grid arrangments. The main differences lie in the way by which the pressure gradients and cell- face velocities are calculated.

To clarify the problem and solutions let us write down the governing equations distinguishing the nodal velocities from cell- face velocities. For two- dimensional steady flow in cartesian coordinate system the set of momentum conservation equations may be written as follows :


      d            d
     --(rho*U) + --(rho*V) = Sm                          ( 2.1 )
     dX           dY

X- direction momentum conservation.

d d --(rho*U*Uc) + --(rho*V*Uc)= dX dY

dPc d dUc d dUc = - --- + -- ( emu* ---) + -- ( emu* --- ) ( 2.2 ) dX dX dX dY dY

Y- direction momentum conservation.

d d --(rho*U*Vc) + --(rho*V*Vc)= dX dY

dPc d dVc d dVc = - --- + -- ( emu* ---) + -- ( emu* --- ) ( 2.3 ) dY dX dX dY dY

Similarly, a 3D simulation using the non- staggered grid will solve finite- volume approximations to partial- differential equations for the following variables:

     : VARIABLE                 : PHOENICS name  :   VARIABLE symbol :
     : Pressure                 : P1             :   Pc              :
     : Node velocity components : UC1, VC1, WC1  :   Uc, Vc, Wc      :
     : Cell face velocity       : U1, V1, W1     :   U, V, W         :
     : components               :                :                   :

This equation set has 4 equations for seven uknowns. So, the problem of closure arises. For two-dimensional cases, only three equations are usually solved but there are still five uknowns ( Pc, U, V, Uc, Vc ).

These equations can be converted into a set of algebraic equations by control- volume based finite- differencing. The set of resulting equations can be written in a generalised form as :

aP*PHIp = aE*PHIE+aW*PHIW+aS*PHIS+aN*PHIN+Sphi ( 2.4 )

where to preserve clarity coefficient pertaining to upwind form on can be given as :

aE = max( de, de - me ) ( 2.5 )

aW = max( dw, dw + mw ) ( 2.6 )

aN = max( dn, dn - mn ) ( 2.7 )

aS = max( ds, ds + ms ) ( 2.8 )

The diffusion ( d's ) and convection ( m's ) fluxes are given as in TR/99. The source terms for each PHI are given below
Equation PHI Sphi
Continuity 1 Sm
X-momentum Uc (Pw-Pe)*Ac,u + Su
Y-momentum Vc (Pn-Ps)*Ac,v + Sv
Z-momentum Wc (Ph-Pl)*Ac,w + Sw

In general case, another three unknown cell- face pressures, Pw, Pn and Ph, are added to the problem followed by the node free surface area ( Ac's) calculations. Su, Sv and Sw are the extra momentum sources, if any.

2.3 Analysis of the model.

To close the mathematical problem, relationships of the following form are required:

U = Lu ( Uc, Vc, Wc, Pc ) ;
V = Lv ( Uc, Vc, Wc, Pc ) ;
W = Lw ( Uc, Vc, Wc, Pc ) ;
P = Lp ( Pc ) . ( 2.9 )

The node-centred free-surface areas should also be determined either in terms of known cell face ones or as independent set of values. The latter practice is , obviously, preferable and will be used throughout below.

The functions Lu, Lv, Lw and Lp represent the required relationships which can be treated as interpolation formulae defining U, V, W and P in terms of nodal velocities and pressures.

2.4 Linear pressure gradient interpolation.

In this case, the cell- face pressures are written as linearly interpolated values from neighbouring nodal values and for uniform grid spacing they are:

Pw=(PW+Pp)/2.; Pe=(PE+Pp)/2.;
Pn=(PN+Pp)/2.; Ps=(PS+Pp)/2. ( 2.10 )

The above practice removes appearance of Pp in balance equations for Up and Vp. This results in the velocity- pressure decoupling because the difference equations for nodal velocities do not contain the pressure in node in question.

2.5 Cell face velocity interpolation.

The coupling can, however, be restored by calculation of the cell- face velocities.

2.5.1 Linear velocity interpolation.

In this case, the linear interpolation is applied to calculate cell face velocities as, for uniform grid spacing, arithmetic mean between neighbouring nodal values, as follows

Uw=(UW+Up)/2.; Ue=(UE+Up)/2.;
Un=(UN+Up)/2.; Us=(US+Up)/2. ( 2.11 )

Experience, however, has shown that this coupling is very weak and may results in zig- zag variations of pressure, while the predictable velocities are often nearly accurate.

2.5.2 The Rhie and Chow "momentum interpolation".

This interpolation technique is currently the most popular for non-staggered calculations since it was first proposed by Rhie and Chow (1983) followed by several researches over the last years.

It can be interpreted as follows: to evaluate the velocity component at the cell face, use the discretised momentum equations for the two neighbouringnodes as the interpolation formulae, but replace the terms representing pressure gradient across the cell face, by one which is centered about it.

An example below will clarify this procedure. On the non- staggered grid nodal velocity, say U, are represented as :

         (sum(ai*Ui))p    Pw-Pe            (Su)p
     Up= ------------- + ------- * Ap,u + -------         ( 2.12 )
            (aP)p         (aP)p            (aP)p

In the method considered the linear interpolation is still used to obtain the interface (cell- face) velocities. Substituting the last equations in arithmetic mean formulae we have:

               (sum(ai*Ui))p+(Su)p   (sum(ai*Ui))E+(Su)E
     Ue= 0.5*[ ------------------- + ------------------- ]
                   (aP)p                 (aP)E

Pw-Pe Pe-PeE + 0.5*( ----- * Ap,u + ------ * AE,u ) (aP)p (aP)E ( 2.13 )

The key idea is to replace the last two terms in RHS of this equation by one which contain the pressure at the node P in an explicit manner, namely

                   1        1
         + 0.5*( -----  + ------ )*( Pp-PE ) * Ae,u      ( 2.14 )
                 (aP)p     (aP)E

This modification provides the required pressure- velocity coupling

                 Ue= He +  ----- * Ae,u                   ( 2.15 )
               (sum(ai*Ui))p+(Su)p   (sum(ai*Ui))E+(Su)E
     He= 0.5*[ ------------------- + ------------------- ]
                    (aP)p                 (aP)E
          1              1        1
        ----- =  0.5*( -----  + ------ )
         (aP)e         (aP)p     (aP)E                   ( 2.16 )

Similar equations can be derived for other cell- face velocities. Then they can be used to calculate mass fluxes for finite- difference coefficients. Now, by adding linear pressure interpolation for calculation of pressure force differences the set of momentum equations is closed and can be solved to get nodal velocities.

2.5.3 Summary of the Rhie and Chow algorithm.

The technique just described divides the calculations of the cell- face velocities in the stages as follows :

  1. Calculate the pressure gradient coefficient distribution as a reciprocal of left hand side coefficient of (2.4),i.e

    PGCp = 1./ ( sum (ai) + SP ) ( 2.17 )

  2. Calculate the pressure- free velocities

    Hp = (sum(Ui*ai) + SU) * PGCp ( 2.18 )

    In (2.17, 2.18) SP and SU contain the contributions of all sources except pressure differences.

  3. Calculate the cell- face velocities by difference equations assembling linearly interpolated pressure gradient coefficients and pressure- free velocities together with actual pressure force differences on both sides of the face in question, i.e. for uniform grid spacing

    Ue = (Hp+HE)/2 + (PGCp+PGCE)/2*(Pp-PE)*Au,e ( 2.19 )

    Vn = (Hp+HN)/2 + (PGCp+PGCN)/2*(Pp-PN)*Av,n ( 2.20 )

    Wh = (Hp+HH)/2 + (PGCp+PGCH)/2*(Pp-PH)*Aw,h ( 2.21 )

    Then the mass fluxes needed for Up, Vp and Wp calculations are easily followed by linear interpolation of pressures for their differences as for uniform grid ( 2.10 ).

2.6 Solution procedure.

The co- located velocity component conservation equations are solved using the standard PHOENICS solver as normal scalar variables at the cell centres. As outlined above, additional sources are added. To solve each co- located velocity component, boundary conditions are specified. The source terms are added before solving the equation using the EARTH internal solver. For the turbulent calculations the turbulent viscosity is stored and calculated, and is used for building up the diffusion coefficients.

In this implementation of non- staggered algorithm U1, V1 and W1 will still be solved for; but the values of these quantities will be over- written, just after they have been solved by PHOENICS, with values which have been obtained from the cell- face velocity interpolation formulae.

Within each iteration, all sources except pressure differences are assembled first and then pressure gradient sources are added. After this, variables are solved one by one in the following order :

Co- located velocities UC1, VC1, WC1 Cell- face velocities U1, V1, W1 Mass balance P

3 Practical details

3.1 Activation

The input parameters required by the CCV program are, at present, obtained from standard Q1 file. The Q1 file is used conventionally to define the problem for PHOENICS built-in staggered fluid flow calculations. In addition, it must perform seven tasks particular to CCV.

These are in turn,

  1. To allocate storage, declare the names, initialize and activate the solutions of CCV variables.
  2. To introduce the boundary conditions and special sources for co-located velocities.
  3. To activate the desired low dispersion convection scheme.
  4. To activate the calculation of the pressure gradients as additional source terms for UC1, VC1 and WC1.
  5. To specify RG and LG variables for transmission of special data to CCV program.
  6. To switch the CCV algorithm on, or off.
  7. To call to the CCV Ground-station.

The first task is achieved by adding the following statements to related groups of the Q1 file followed by initialization commands.


The conventional Q1 commands can still be used to introduce the certain kind of boundary conditions and special sources for co-located cartesian velocities. Thus, one can retain practically unchanged the specification of INLET and OUTLET fluxes. To provide zero-slip wall conditions for laminar flows on Cartesian grids no extra efforts are needed to specify obvious near wall sources for co-located velocity components as can be seen, say, from U201 CCV library case.

In general case, however, it is, at present, the user or CCV developer responsibility to specify the boundary values and to set required source terms for co- located velocities. It can be done partly by rewriting the relevant GX subroutines of GREX3 adopting them to co-located velocities and partly by creation of the completely new codings particular to CCV.

The examples of both kinds can be found in subroutine GXCCVG (see Section 6.1 for more information). The latter works for CCV -EARTH in the same manner as GREX3 works for EARTH. It is the GXCCVG contents which is mainly responsible for present CCV limitations listed above.

A list of the sample PATCH and COVAL statements which activate the number of conditions extracted from the existing GXCCVG collection are shown below.

  1. Laminar zero- slip wall conditions (BFC=F)

    They are activated by following standard statements:

    WALL (SOUTH,SOUTH,1,NX,1,1,1,NZ,1,1)

  2. Laminar zero- slip non- moving wall conditions (BFC=T)

    WALL (SOUTH,SOUTH,1,NX,1,1,1,NZ,1,1)

  3. Turbulent zero- slip wall conditions (BFC=F)

    They are activated by following standard statements:

    WALL (SOUTH,SOUTH,1,NX,1,1,1,NZ,1,1)

  4. Turbulent zero- slip non- moving wall conditions (BFC=T)

    WALL (SOUTH,SOUTH,1,NX,1,1,1,NZ,1,1)

  5. Rate of the turbulent energy production (BFC=T)

    No special statements are needed. When CCV=T the subroutine GENBGR is called from GXCCVG to re-write the rate of turbulent energy production, see Section 6.2

  6. Buoyancy forces


3.2 Current capabilities and limitations:

The non- staggered CCV calculations, as they are formulated here, may only be applied to the one- phase fluid flow problems which are solved on Cartesian and on a general non-orthogonal BFC grids for which the mass flux velocities are set in the calculation procedure by way of re-writing the standard PHOENICS staggered velocities by the values obtained from the momentum interpolation and coordinate transformation. ( However, problems with polar coordinate system can usually be solved for if the corresponding grid is created by means of BFC grid generation techniques.) A further known limitations on the use of the procedure are at present, as follows:

These limitations are essential for the present stage of the development. They will be removed one by one in course of further improvment of the coding and implementation style.

The authors of the coding are Sergei Zhubrin and Nicolai Pavitskiy at CHAM MEI. The coding was created, tested and developed at various times during 1992 and 1993.