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:

- the just-mentioned Input Library files;
- the Fortran file GXCCVG.F, in d_phoe21/d_earth/d_opt/d_numalg; and
- the following Encyclopaedia article.

NOTES ON NON- STAGGERED CALCULATIONS OF FLUID FLOW by S.V. Zhubrin

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

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.

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.

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 :

Continuity.

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.

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.

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.

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

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.

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

Pp-PE Ue= He + ----- * Ae,u ( 2.15 ) (aP)e where, (sum(ai*Ui))p+(Su)p (sum(ai*Ui))E+(Su)E He= 0.5*[ ------------------- + ------------------- ] (aP)p (aP)E and 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.

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

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

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

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

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

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,

- To allocate storage, declare the names, initialize and activate the solutions of CCV variables.
- To introduce the boundary conditions and special sources for co-located velocities.
- To activate the desired low dispersion convection scheme.
- To activate the calculation of the pressure gradients as additional source terms for UC1, VC1 and WC1.
- To specify RG and LG variables for transmission of special data to CCV program.
- To switch the CCV algorithm on, or off.
- 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.

SOLVE(UC1,VC1,WC1)

STORE(EPRP,NPRP,HPRP)

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.

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

COVAL(SOUTH,UC1,1.0/PRNDTL(UC1),UC1WALL)

COVAL(SOUTH,WC1,1.0/PRNDTL(WC1),WC1WALL) - Laminar zero- slip non- moving wall conditions (BFC=T)
WALL (SOUTH,SOUTH,1,NX,1,1,1,NZ,1,1)

COVAL(SOUTH,UC1,1.0/PRNDTL(UC1),GRND)

COVAL(SOUTH,VC1,1.0/PRNDTL(VC1),GRND)

COVAL(SOUTH,WC1,1.0/PRNDTL(WC1),GRND) - 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)

COVAL(SOUTH,UC1,GRND2,UC1WALL)

COVAL(SOUTH,WC1,GRND2,WC1WALL) - Turbulent zero- slip non- moving wall conditions (BFC=T)
WALL (SOUTH,SOUTH,1,NX,1,1,1,NZ,1,1)

COVAL(SOUTH,UC1,GRND,GRND)

COVAL(SOUTH,VC1,GRND,GRND)

COVAL(SOUTH,WC1,GRND,GRND) - 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

- Buoyancy forces
PATCH(BUOY,PHASEM,1,1,1,NY,1,NZ,1,1) COVAL(BUOY,WC1,FIXFLU,GRND3)

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:

- The Rhie and Chow interpolation is used to evaluate the cell- face velocities followed by the linear interpolation for pressure is used to calculate the cell- face cartesian velocities and pressure differences.
- One alternative to the linear pressure interpolation,Symmetrical Momentum Source Scheme, SMS, /2/ is provided for CARTES=T.
- The zero- slip conditions, for the case of rigid walls, should be set by means of PATCH and COVAL commands ( No CONPOR facilities for automatical activation of wall function for co- located velocities near the blockages are provided so far).
- For BFC problems the K-E turbulence model is the only one which supplied by correct implementation of the wall functions and kinetic energy generation term.
- Just part of existing GREX3 facilities are activated for application on non-staggered grid ( see Section 5 below ) automatically. But the user is allowed to fill the standard GROUND file by problem specific extensions in conventional manner guided by GXCCVG subroutine codings.
- Node porosities should be specified and initialized by user along with the specifications for cell-face ones.
- No XCYCLE option for co-located velocities is provided.
- Linear relaxation is provided for cell- face velocities.
- SBS- solver should be used for the co- located velocities solutions.

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.

wbs