An Alternative Finite Volume Discretization of Body Force Field on Collocated Grid

Collocated grids are more suitable for the implementation on general geometries than the staggered counterparts, but their use requires the enhancement of the pressure-velocity (p − v) field coupling. This is achieved by thoughtful interpolation of the velocities on finite volume faces. However, employing standard interpolation schemes, such as the well known Rhie-Chow scheme, can cause unphysical spikes in the velocity field when an abruptly changing body force field is present; an example is shown in Fig. 1. To understand the problem and find the remedy, proposed originally in (Mencinger & Žun, 2007), we should analyze the connection between p and v fields. This connection is highlighted in the following subsections.


Introduction
Collocated grids are more suitable for the implementation on general geometries than the staggered counterparts, but their use requires the enhancement of the pressure-velocity (p − v) field coupling. This is achieved by thoughtful interpolation of the velocities on finite volume faces. However, employing standard interpolation schemes, such as the well known Rhie-Chow scheme, can cause unphysical spikes in the velocity field when an abruptly changing body force field is present; an example is shown in Fig. 1. To understand the problem and find the remedy, proposed originally in (Mencinger & Žun, 2007), we should analyze the connection between p and v fields. This connection is highlighted in the following subsections.

The Navier-Stokes equation
The conservation of momentum, expressed with the Newton's second law, is in fluid dynamics represented through equation where ρ is the density of an infinitesimal fluid particle and D v/Dt is its acceleration due to the presence of body force field f and stress σ in the considered fluid. The stress tensor σ contains both the pressure p and the viscous stress which is (for Newtonian fluid) proportional to the rate of strain. For incompressible fluids (∇· v = 0), considered in this text, it is written as where η denotes the viscosity and I the identity tensor. If η does not change substantially in the fluid, then (1) becomes the Navier-Stokes equation (Landau & Lifshitz, 1987) The components of (3), i.e. the x-component in the Cartesian system (where v =(u, v, w) t ) can be written in the form of the general transport equation for a specific scalar quantity φ ρ Dφ Dt = ∇·(Γ∇φ) + S φ (5) by setting φ = u, Γ = η, S φ = −∂p/∂x + f x and can thus be straightforwardly discretized with the finite volume (FV) method. However, finding the solution of the obtained discretized momentum equations is not straightforward. Firstly, the unknown pressure field needs to be found. Secondly, the uninformed discretization of the terms in (4) leads to obtaining nonphysical numerical artifacts which can overwhelm the solution. One such typical artifact is the appearance of checkerboarding pressure field which originates from weak p − v coupling, explained in the next subsection.

Origin of weak p − v coupling
The finite volume (FV) discretization of (4), obtained by its integration 1 over control volume P sized V P (e.g., see Ferziger & Perić (2002)) and time interval Δt results in 1 Actually, the so-called conservative form of (4) is considered, where ρDu/Dt is replaced by ∂(ρu) ∂t + ∇·(ρ vu) on the left hand side.
where the values at the beginning and at the end of the time step are denoted as () 0 and written without the superscript, respectively. Fully implicit time integration is used for simplicity, although the following discussion is not limited to this choice of temporal scheme. The subscripts () P and () f in (6) denote the average values in FV and on FV face, respectively, and the summation ∑ f comprises all FV faces. Each face is represented by surface vector S f = S f n f , where S f is the face area and n the normal vector pointing out of the FV. The velocity has two roles in (6). First, it appears in the mass flux through f-th face following from the assumed incompressibility. Second, the x-component of the velocity, u f in (6), is also convected velocity i.e. the formal unknown of the transport equation (by setting φ = u). To obtain a solvable form of (6), it is written as  To concretize and simplify the discussion, we consider a 1-dimensional uniform grid, shown schematically in Fig. 2, so that S f = 1, F e = ρ e u e , F w = −ρ w u w , and Δx e = Δx w ≡ Δx. If linear interpolation and central differencing are used to approximate values and derivatives on FV faces, then (8) can be written as with the coefficients The solution of (9), written formally as is obtained throughout an iterative process which usually requires under-relaxation of the calculated solution with the one from previous iteration u * P , weighted with the under-relaxation factor α u . We can also notify that u P , as written above, depends only on p E and p W and not on p P .
To obtain the solution of (9), we need convecting velocities (appearing in the coefficients as factors in FV-mass fluxes) and the values of pressure. The first, u e for example, can be obtained simply as 1 2 (u P + u E ); using (11) for u P and a matching relation for u E then yields where the superscript () α is used as a α P ≡ a P /α u and A corresponding equation can be obtained for u w .
The pressure field, on the other hand, is calculated from the equation which is obtained from the discrete incompressibility condition (7). On the considered 1D grid it simplifies to or, using the above interpolation, to u E − u W = 0. This means that the conservativeness of mass in cell P, expressed through (7), does not depend on u P . Furthermore, the pressure field equation connects the values of p in cells P, EE, and WW. So, the discretized pressure field disintegrates in two mutually independent parts (the other part contains cells E and W), which in practice often results in an unphysical zigzagging pressure pattern. An analogous difficulty appears also on three-and two-dimensional grids. It is manifested in a checkerboarding pattern of the pressure field. One might expect that the problem would disappear on nonuniform grids; unfortunately, this does not happen in practical calculations because the connection between the values of p in neighboring cells remains too weak.

Staggered versus collocated grids
The origin of the problem described above lies in the mentioned fact that u P depends on (∂p/∂x) P which is calculated using the values at grid points at two alternate cells (E and W).
As the two cells are 2Δx apart, this means that the gradient is actually obtained from the grid that is twice the coarse than the one actually set. A well known solution to circumvent this unwanted situation is to move the grid belonging to u component for a distance 1 2 Δx, so the grids belonging to p and u are staggered to one another.
Staggering the grids belonging to the velocity field components in the direction of the corresponding component first appeared in the paper of Harlow & Welch (1965). It enables more accurate representation of the continuity equation and of the pressure gradient in the Navier-Stokes equation. More importantly, it insures strong pressure-velocity coupling, required to obtain realistic solution of the equation. On the other hand, collocated (i.e. nonstaggered) grids have some obvious advantages over the staggered ones; Perić et al. (1988) describe them as follows: "(i) all variables share the same location; hence there is only one set of control volumes, (ii) the convection contribution to the coefficients in the discretized equations is the same for all variables, (iii) for complex geometries Cartesian velocity components can be used in conjunction with nonorthogonal coordinates, yielding simpler equations than when coordinate oriented velocity components are employed, and (iv) there are fewer constraints on the numerical grid, since there is no need to evaluate the so-called curvature terms." In short, the collocated grids offer much simpler CFD code implementation than the staggered counterparts when the domain geometry is complex. This seems to be the main reason why the majority of the popular commercial codes use collocated grids. As the collocated grid arrangement does not inherently insure strong p − v coupling which prevents the appearance of a nonphysical checkerboard pressure field, the coupling has to be insured by other means than the grid staggering. The established method for the coupling enhancement is the employment of the momentum interpolation scheme of Rhie & Chow (1983). This scheme together with additional important corrections (Choi, 1999;Gu, 1991;Majumdar, 1988) is described below.

Corrections of convecting velocity
The interpolated convecting velocity u e in (12) contains weighted interpolation of the pressure gradient, which leads to the disintegration of pressure field. The idea of the interpolation scheme of Rhie and Chow is to replace the interpolated derivative with the one calculated directly. However, the latter is multiplied with the corresponding interpolated coefficient so that (12) is corrected as where V e /(a α P ) e ≡ 1 2 (V E /(a α P ) E + V P /a α P ) and a := b is read as "a becomes b." Inserting (16) and an equivalent relation for u w in (14) changes the pressure field equation to which defines more compact computational molecules than (15) and prevents the previously described breakup of the pressure field.
The equation (17) connects the value of p P with the values in adjacent cells p E and p W . Yet, p P is not directly related with f x P , but only with f x E and f x W . This situation can be resolved with an another important correction, proposed by Gu (1991) so that (17) becomes where f x e and f x w are the body forces on the corresponding FV faces.
Additional corrections and proposed by Majumdar (1988) and Choi (1999), respectively, are obtained in the same spirit. The first one prevents the dependence of u e on the under-relaxation factor α u in the converged solution, while the second one diminishes the dependence of u e on the time-step size. The dependence on Δt is not completely removed because Δt is still contained in a α P . The complete removal the dependence on Δt was proposed, for example, by Yu et al. (2002) and recently by Pascau (2011). The latter work shows that this topic is still actual.

The problem of nonphysical spikes in velocity field
It turns out that the Rhie-Chow interpolation scheme works well as long as the pressure field is sufficiently smooth, i.e. without any abrupt variations. As explained below, it produces nonphysical spikes in the velocity field such as shown in Fig. 1; the spikes appear near the abrupt variations of the pressure field. The latter are generally a consequence of the abrupt changes in the body force field as can also be understood from the Navier-Stokes equation. Namely, the abrupt variation of f in (3) is counter-balanced by such a variation in ∇p. This is more obvious when the fluid is quiescent ( v = 0) so that (3) becomes and the equation for the pressure field is obtained by calculating the divergence of (22)

Example of abrupt body force variation: multiphase flow
In most cases dealing with fluids, the body force field originates from gravity: f = ρ g where g is the gravitational acceleration. Thus, an abrupt variation in ρ results in such a variation in f and the former appears often when dealing with multiphase system such as gas-liquid. Multiphase systems present an added difficulty in the flow simulations and require additional modeling. The most widespread approach is the employment of an Eulerian model, which make (3) valid throughout the whole flow domain regardless of the phase; this is typically achieved by using a phase identifying scalar field C, defined as C( x, t)= 1, x occupied by fluid 1 (e.g. gas), 0, x occupied by fluid 2 (e.g. liquid).
The finite volume discretization transforms C( x, t) to volume averaged C P which represents the volume fraction of fluid 1 in volume P; thus, the cells where 0 < C P < 1 contain the interface.
A characteristic representative of Eulerian models is the Volume-of-Fluid (VOF) model (Hirt & Nichols, 1981) in which a two-phase system is treated as a single fluid with material properties defined as a linear combination of the phase specific properties ρ 1 , ρ 2 , η 1 , and η 2 where C is advected (passively) through the considered domain. This is described by The solution of (26) is far from trivial and requires special discretization methods which surpass the scope of this chapter; a comprehensive overview of such methods is written, for example, by Scardovelli & Zaleski (1999).
Returning back to the body force, using (25) f becomes so that the abrupt changes in f , proportional to the density difference, are present at the phase interface. An another example of the sudden change in the body force field, when employing VOF, is due to the surface tension. The latter is modeled by the continuum surface force (CSF) model (Brackbill et al., 1992) which 'converts' the surface tension to a body force acting in the vicinity of the interface f = σκ∇C (28) where σ and κ are the surface tension coefficient and the curvature of the interface, respectively.

Examination of the problem in 1D
Clearly, the obtained pressure field equation depends on the interpolation of the velocities on CV boundaries. To investigate the influence of different interpolations in the situations with the presence of the abrupt body force field, we setup a one-dimensional case with 0 < x < 0.1 m and discretize the defined domain with the uniform grid containing 40 elements. The body force field is defined as f x = (Cρ 1 +(1 − C)ρ 2 ) g where g = 10.0 m 2 /s. In the considered case, f x is determined with the phase discrimination function where x 1 and x 2 are set to 0.030 m and 0.062 m, respectively. Material properties of air (ρ 1 = 1.29 kg/m 3 , η 1 = 1.8 × 10 −5 Ns/m 2 ) and water (ρ 2 = 1.0 × 10 3 kg/m 3 , η 2 = 1.0 × 10 −3 Ns/m 2 ) are used.
By setting the velocity at the boundary to zero we expect uniform zero velocity field all over the domain. Also, linear (by parts) pressure field is expected since ∂p/∂x is counterbalancing the body force field. The resulting pressure field is shown in Fig. 3(a): the pressure field obtained without the Rhie-Chow interpolation scheme (16) exhibits a zigzagging pattern which in this case does not overwhelm the solution. Expectedly, this pattern disappears with the employment of (16). However, Fig. 3(b) shows that this scheme produces unwanted spikes near the discontinuity. Adding Gu's correction does not appear to have a notable effect on the pressure field ( Fig. 3(a)) whereas the spikes (Fig. 3(b)) appear even larger in the presented case.

The remedy in 1D
The velocity field in Fig. 3(b) is represented by the values of u inside CVs; the values on CV boundaries are, on the other hand, equal to zero (numerically) as required by (14). To obtain the zero velocity field also inside CVs must hold. This follows from (9) by setting the velocities to zero. It simply means that the discretized body force and the pressure gradient should be counterbalanced to obtain zero velocity inside FV. At the same time, when both the Rhie-Chow and Gu's scheme are employed, the pressure field is obtained from which follows from (19) for quiescent fluid. The solution of the above equation can be constructed in a rather simple manner: by setting p in a selected starting cell to an arbitrary value and then using the relations, following obviously from (31), to calculate the values in the neighboring cells. Inserting (32) into (30) results in which is the condition to obtain zero velocity field inside FV when the pressure field is calculated from (31). However, the above equation could also be interpreted as a rule how to discretize the body force field. It instructs to construct the body force field inside CVs as a linear combination of the corresponding average values on CV boundaries. The rule eliminates completely (i.e. to round-off error) the unwanted spikes as demonstrated in Fig. 3(b).
The proposed body force discretization rule (33) is obtained for a uniform 1D grid. For a nonuniform grid it becomes where I e =( x e − x P )/Δx e and I e =( x P − x w )/Δx e denote the interpolation factors (with reference to Fig. 2). The above equation is obtained by using which follows from the Gauss's theorem (39).

The generalization of the remedy
Using the same procedure, let us obtain the body force discretization rules for general (i.e. nonorthogonal and/or nonstructured) 2-and 3-dimensional grids. Again, we assume that the solution of the pressure field equation for the quiescent fluid on all CV faces satisfies when both Rhie-Chow and Gu's corrections are used. The dot product of (36) with vector Δ f results in where subscript () Nb(f) denotes the value at neighboring cell Nb(f), sharing face f with cell P, and vector Δ f points from point P to neighboring point Nb(f) as shown in Fig. 4. When dealing with a 1-dimensional grid, (37) can simply be used to construct the pressure field starting from an arbitrary grid point. On 2-and 3-dimensional grids, however, it has to be noted that condition (36) can be satisfied on all FV faces simultaneously only if the body force field is conservative on the discrete level (Mencinger & Žun, 2007). This is generally not true, and the pressure field can not be constructed simply by using (37). Furthermore, the quiescent solution can not be obtained, i.e. the so-called spurious currents appear. Nevertheless, relation (37) can be considered as a reasonably good approximation and the described methodology still valid.
As the rule is obtained from the requirement for the zero velocity field following from the discretized Navier-Stokes equation, its form depends on the discretization of ∇p used in the equation. For example, (∇p) P can be obtained by employing the Gauss's theorem If p f is written in terms of the neighboring cells using the interpolation coefficient I f as then it can be written using (36) as Inserting (41) and (39) in (38) results in where the first summation term on the right hand side equals zero as ∑ f S f = 0 must hold for any closed surface. Finally, is obtained.
The rule can be generalized even further; if (∇p) P is written in terms of p P and the values of p in the neighboring cells as where γ P and γ Nb are geometrical vector coefficients. The latter can be obtained, for example, with the least squares method. Inserting (44) and (37) in (38) now results in The term in square brackets in (45) should be zero from the definition (44): the effect of adding a constant field p 0 to p should vanish. Therefore, (45) simplifies to and presents the discretization rule when (44) is used to obtain (∇p) P in the discretized Navier-Stokes equation. As rule (43), it contains the surface averaged body forces and the geometrical factors. Fig. 5. The velocity field near the rising bubble indicated by the contour; the same calculation as in Fig.1 but with the alternative body force discretization.

Testing the alternative discretization
The proposed alternative body force discretization is obtained considering a quiescent fluid which is seldom of research interest. To compare its performance against the standard discretization for a moving fluid we consider two cases: (i) raise of a bubble in a rectangular cavity and (ii) natural convection in a square cavity. The first case demonstrates that the unwanted velocity field spikes are removed or at least largely diminished by using the alternative body force discretization. Whereas an abrupt change of body force field is dealt with in the first case, the body force field is smooth in the second case. Namely, we also want to check the effect of the new discretization in such cases.

Rising bubble
The case considers a two-dimensional cavity with no-slip boundary condition at all walls. The cavity is rectangular (width w, height h). Initially, it is filled with water and an air bubble of radius R, centered at (x 0 , y 0 ); both the water and the air are quiescent. To follow the rise of the bubble using the VOF model, (26) needs to be solved besides the Navier-Stokes equation. In the latter, both the buoyancy (27) and the surface force (28) are taken into account. Whereas standard central differencing is used to calculate advection and diffusion terms in the momentum equation, the discretization of (26) requires nonstandard differencing scheme in order to reduce numerical diffusion which results in undesirable interface smearing. The low-diffusive CICSAM scheme (Ubbink & Issa, 1999) is used in the presented case. The curvature κ, needed to obtain the forces due to the surface tension by using the CSF model are calculated as suggested by Williams et al. (1998). In the presented case we set w = 0.02 m, h = 0.04 m, x0 = w/2, y0 = h/4, and R = 0.002 m. Uniform grid with 40x80 CVs was used. Fig. 5 shows the calculated velocity field around bubble at t = 0.115 s when using the alternative discretization; this figure can be compared directly with Fig. 1 which presents the same calculation except for using the standard discretization. Obviously, more realistic results are obtained with the proposed discretization. Interestingly, despite large spikes in the first (i.e. standard) calculation, the bubble contours of the two calculations indicated in Fig. 6 does not differ as much as one would expect. This is perhaps due to the fact that CV face velocities, which appear in the advection terms, are corrected so that they satisfy the continuity equation. The second case deals with a buoyancy driven flow in a two-dimensional square cavity and presents a classical CFD code benchmark problem (De Vahl Davis, 1983) with well known solutions. The left and the right wall are set to fixed temperatures T h and T l (T h > T l ), respectively. Both horizontal walls are thermally insulated. The velocity components vanish at the walls. The velocity field is obtained by solving (3) with Boussinesq approximation

Natural convection in a square cavity
where β is the thermal expansion coefficient and ρ 0 is the density at reference temperature T 0 . The temperature field is determined with the enthalpy transport equation where c p and λ are the thermal conductivity and the specific heat at constant pressure, respectively. Actually, both (3) and (48) are solved in nondimensional form which read where, for brevity, the same notation is used for nondimensional and dimensional quantities. Obviously, the problem is determined with the two dimensionless parameters Pr and Ra denoting Prandtl and Rayleigh number, respectively. They are defined as where L is the size of the cavity. Following the work of de Vahl Davis, we set Pr = 0.71 and consider only the highest value of Ra used in their test: Ra = 10 6 .
The calculation is performed using relatively coarse uniform grids containing 10x10, 20x20, 40x40, and 80x80 cells. The central discretization scheme is implemented for the interpolation of the values on CV-faces for both (49) and (50). The obtained streamlines and isotherms are shown in Fig. 7. Fig. 8 shows the variation of the calculated velocity components and the temperature along the horizontal and the vertical centerline of the cavity. Obviously, the difference between the results obtained with the standard (dashed line) and the proposed (solid line) body force discretization vanishes with the increased grid density. The difference on the 80x80 grid is not noticeable on the presented scale and is therefore not drawn in Fig. 8. Both the new and the standard discretization converge to the same solution, thus we can assume that the proposed discretization is consistent.

Conclusions
Although it was obtained for a quiescent fluid, the proposed discretization works well also for moving fluid as demonstrated in the considered cases. One should keep in mind that the obtained rules are valid when both the Rhie-Chow and Gu's correction are used. Together with the corrections of Choi (1999) and Yu et al. (2002), they approach the calculations on collocated grids to those on staggered grids. It is shown that the proposed discretization of the body force field is more appropriate than the standard one when dealing with abruptly variable body force fields. In fact, it can be used generally as it does not significantly change the calculated solutions when the body force field is smooth.
The proposed discretization does not exactly follow the spirit of the FV method where simply the average body force within a FV is considered. Even so, it is consistent since it converges to the same solution as obtained with the standard discretization. Its form depends on the chosen discretization of ∇p in the Navier-Stokes equation; in the presented calculations, Gauss's divergence theorem was used. Other possibilities such as using the least squares method are also possible, so that the discretization rule changes its form. Nevertheless, the idea remains that FV average value of f is to be replaced with the linear combination of the average FV face values.

Acknowledgments
The author is grateful to the Slovenian Research Agency (ARRS) for its financial support throughout the P2-0162 project.  We hope that among these chapters you will find a topic which will raise your interest and engage you to further investigate a problem and build on the presented work. This book could serve either as a textbook or as a practical guide. It includes a wide variety of concepts in FVM, result of the efforts of scientists from all over the world. However, just to help you, all book chapters are systemized in three general groups: New techniques and algorithms in FVM; Solution of particular problems through FVM and Application of FVM in medicine and engineering. This book is for everyone who wants to grow, to improve and to investigate.