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:
NOTES ON NON- STAGGERED CALCULATIONS OF FLUID FLOW by S.V. Zhubrin
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 :
PGCp = 1./ ( sum (ai) + SP ) ( 2.17 )
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.
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,
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.
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)
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)
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)
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)
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
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:
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