Use of Proper Closure Equations in Finite Volume Discretization Schemes

In the case of thermo-fluid problems, decision regarding the choice of unknown variables needs to be taken at the continuous level of formulation. A popular choice, among many available options, is the pressure-based primitive variable formulation (Acharya, 2007), in which the flow continuity equation is considered as a constraint on fluid pressure. Pressure-based finite volume methods, in turn, can be categorized based on the computational grid arrangement used for the discretization of solution domain. In a staggered grid arrangement, velocity components and other unknown scalars are defined at different nodal points (Patankar, 1980). On the contrary, all unknown variables are defined at the same nodal locations in a co-located grid arrangement. The discussion of finite volume discretization schemes in this article is limited to pressure-based formulations on co-located grids.


Introduction
Finite Volume Method (FVM) is a popular field discretization approach for the numerical simulation of physical processes described by conservation laws. This article presents some advances in finite volume discretization of heat and fluid flow problems.
In the case of thermo-fluid problems, decision regarding the choice of unknown variables needs to be taken at the continuous level of formulation. A popular choice, among many available options, is the pressure-based primitive variable formulation (Acharya, 2007), in which the flow continuity equation is considered as a constraint on fluid pressure. Pressure-based finite volume methods, in turn, can be categorized based on the computational grid arrangement used for the discretization of solution domain. In a staggered grid arrangement, velocity components and other unknown scalars are defined at different nodal points (Patankar, 1980). On the contrary, all unknown variables are defined at the same nodal locations in a co-located grid arrangement. The discussion of finite volume discretization schemes in this article is limited to pressure-based formulations on co-located grids.
At the continuous level of formulation in a finite volume method, an intensive variable  is constrained by an integral balance equation for a control volume CV bounded by the control surface CS as follows: In Eq. (1) J   represents the flux of  across the control surface, S   stands for the rate of generation of  within the control volume and n  is the outward unit normal vector to the control surface. In most of the engineering problems of interest, there are two mechanisms for the transport of  , i.e. advection and diffusion, and the flux function is mathematically described as follows: In Eq.
(2) V   is the advective mass flux and   is the diffusion coefficient associated with variable  . By properly defining  , J  mathematically describe the transport of mass, momentum and energy in fluid flow problems (Patankar, 1980).
A typical two dimensional discrete domain covered by contiguous control volumes (cells) is shown in Fig. 1. Each cell represents a nodal point and is bounded by a number of panels which comprise the control surface. Typically, each panel includes one integration point. There isn't any limitation on the shape and the number of panels of a cell. Here we use a background finite element mesh to define the control volumes. Figure 1 shows structured and unstructured two-dimensional grids in an element-based finite volume method (EBFVM) context (Ashrafizadeh et al., 2011).
The objective of any finite volume discretization scheme is to use the CV balance equations, similar to Eq. (1), to obtain a set of algebraic equations which constrains the unknown nodal values throughout the domain. This objective is achieved by carrying out the discretization in two stages.
At the first stage of discretization, here called level-1 approximation, Eq. (1) is converted to a semi-discrete form comprising of both nodal and integration point variables. Level-1 approximate form of Eq. (1) for an internal finite volume in Fig. 1 can be written as follows: In Eq. (3), subscripts np and ip stand for nodal and integration point variables respectively. Closure equations are now needed to relate ip variables to np variables.
At the second stage of discretization, here called level-2 approximation, Eq. (3) is converted to a fully-discrete form, i.e. an algebraic equation which constrains the variable  at nodal point P . The computational molecule obtained at this stage of discretization is written in the following convenient form: Here, subscript nb stands for the neighbor nodes and superscript nb N refers to the number of influential neighbor nodes. Influence coefficients P a and nb a and the source term S contain terms and parameters which model all physical effects relevant to the value of  at the nodal point P .
If Eq. (1) is a nonlinear equation, linearization also becomes necessary and should be carried out during the discretization. In principle, the equation can be linearized before or after each of the approximation levels and this might affect the properties of the final discrete model.
In most engineering problems, including heat and fluid flow ones, a coupled set of nonlinear partial differential equations needs to be discretized and solved. In such cases, the inter-equation couplings should also be numerically modeled. Two-dimensional, incompressible Navier-Stokes equations, for example, can be written in a form similar to Eq.
(1) in which  is equal to 1 in the mass balance equation, u in the x-momentum balance equation, v in the y-momentum balance equation and e (specific energy) in the energy www.intechopen.com balance equation. In general, the computational molecule associated with each equation of a coupled set includes more than one variable to reflect the fact that the equations and the variables are inter-related. Therefore, in the case of a set of equations with two coupled variables  and  , the following modified form of Eq. (4) In this article, we focus on level-2 approximation of Navier-Stokes equations in the context of co-located, pressure-based finite volume method. This level of approximation is particularly important because it deals with the modeling of physical influences such as upwind and downwind effects and the couplings between variables. A number of famous level-2 approximations are reviewed and the method of proper closure equations (MPCE), proposed by the first author, is introduced. To avoid the complexities associated with multi-dimensional problems and to focus on the nature of the approximations in a simple setting, one dimensional equations and test cases are discussed in details. Two-dimensional results are only briefly presented to show the applicability of the MPCE in multi-dimensional problems.

The semi-discrete form of equations in a 1D context
A one-dimensional variable-area duct is shown in Fig. 2. The complete set of governing equations for unsteady 1D-inviscid compressible flow in this duct, using the momentum variable formulation, is as follows (Chterental, 1999): Where  is time, fu   is the momentum variable, u is velocity, p is pressure and a is the cross sectional area. Note that no artificial energy source term is used in the energy equation.
In order to close the above set of equations, an auxiliary equation is needed. Here, the ideal gas equation of state is used: Where t is temperature and R is the gas constant. Neglecting the potential energy term, the total energy per unit mass ( e ) of an ideal gas can be written as follows: Where v c is the specific heat at constant volume.
By integrating Eq. (6) over the P control volume in Fig. 2 and after using divergence theorem, the semi-discrete forms of the governing equations are obtained as follows: Note that capital letters (like F and P ) are used for the nodal values and small letters (like f and p ) are used for integration point variables. Here, P V stands for the volume of cell P. Closure equations for the integration point variables are discussed next.

Closure equations for integration point variables
After obtaining the semi-discrete form of governing equations and linearizing them, closure equations are needed to approximate integration point (ip) quantities in terms of the nodal values. By employing these closure equations, the semi-discrete equations are converted to fully-discrete equations. These ip-equations play a critical role in the robustness and accuracy of a collocated scheme. Therefore, different solution strategies for defining ipequation are briefly described in this section. Here, we only discuss the closure equations for "east" integration point of a typical cell around node P. The closure equations at the "west" ip are obtained similarly.

Methods based on geometrical interpolation
Simple symmetric interpolation is the most trivial choice for an ip-value in a uniform grid: On non-uniform grids, weighted interpolation is used as follows: (1 ) Where e  is the weight factor for "east" ip, calculated by the following formula: This scheme is 2nd-order accurate, but is unbounded so that non-physical oscillations may appear in regions with high gradients. To show the problem associated with this scheme, consider the 1-D inviscid incompressible flow in a constant area duct shown in Fig. 3. The semi-discrete form of continuity and momentum equations in this case are as follows: Where A is the cross sectional area of the duct and m  is the mass flow rate based on the most recent available value of velocity. Using simple symmetric interpolation, Eqs. (13) and (14) are converted to the following fully-discrete forms: Combining the above equations, one obtains the following constraint on the pressure at nodal point P: Note that pressure at node P is absent in this equation and this constraint cannot distinguish between a uniform pressure field and the wiggly pressure field shown in Fig. 4.

Methods based on taylor-series expansion
It is logical to assume that convected quantities in a thermo-fluid problem are influenced by the upstream condition. Therefore, assuming the flow direction from P to E, the following truncated Taylor series expansions are both valid approximations for e  : Equation (18) (13) and (14), would not be able to detect the checkerboard pressure fields. Such numerical schemes are prone to unphysical pressurevelocity decoupling symptom.

Methods based on momentum interpolation
Rhie and Chow (Rhie, 1983) employed a co-located grid to solve the flow field around an airfoil. To avoid nonphysical, wiggly numerical solutions, they came to the conclusion that different closure equations had to be used for the convecting or mass-carrying (û ) and convected or transported ( u ) velocity components at the integration points. The convecting velocity, used in the semi-discrete continuity equation, is linked to the pressure through a momentum interpolation procedure. Hence, this discretization scheme is called the Momentum Interpolation Method (MIM). As compared to the classical schemes on staggered grids, one may assume that nodal momentum balances are conceptually staggered to obtain closure equations for convecting velocity at integration points. The discrete momentum balance for node P can be written as follows: The convecting velocity at "e" is obtained by staggering the stencil of the "P" computational molecule to the "e" integration point: The closure for the convected velocity, and any other transported variable, comes from an upwind scheme. The simplest option at the integration point e would be the FOU: The interpolation formula, Eq. (22), is not unique as explained in (Acharya, 2007) and one might successfully implement the idea using modified forms of this approach. The key point, however, is to use different closure equations for the convected and convecting velocity components. The convecting velocity should be properly related to the pressure field.
The closure equation for e p has not been a matter of concern, at least for incompressible flows. Taking into account the elliptic behavior of the pressure, the following closure equation is used in the MIM: This technique guarantees the required physical pressure-velocity coupling in the numerical model and prevents any wiggly solutions (Ashrafizadeh et al., 2009).

Methods based on physical influences
The pleasure of working with all variables on a single nodal position, as opposed to the pain of working with staggered grids limited to simple geometries, has made the MIM very popular. Therefore, attempts have been made to employ and improve the MIM for the solution of flows at all speeds. Schneider and Raw (Schneider & Raw, 1987a, 1987b proposed the physical influence scheme to unify the integration point velocity interpolation formulas in co-located grids. They retained a unified definition for the integration-point velocity components, but argued that the closure for e u should not be obtained solely by purely mathematical upwind schemes. They suggested the non-conservative form of the momentum equation to provide a closure for e u : Where the over bar refers to the most recent available value in the nonlinear iterations. Equation (25) is used as the closure for e p in this scheme.
While successfully implemented in a number of multidimensional viscous problems, the physical influence scheme fails to prevent pressure checkerboard problem in inviscid incompressible flows (Bakhtiari, 2008). This failure reinforces the belief that two definitions for convected and convecting velocities at integration points are necessary when the calculations are carried out on a co-located grid.

Methods based on modified physical influences
To fix the problem associated with the physical influence scheme, Karimian and Schneider (Karimian, 1994a) accept the notion of dual velocities at the integration points and employ a combination of momentum and continuity equations to provide a closure for the convecting velocity ( They also use Eq. (27) for convected velocity ( e u ) and Eq. (25) as the closure for e p . Numerical experiments with the 1-D test case show that this method works well and successfully suppresses any oscillatory numerical solution even in time-depended all speed flows (Karimian, 1994b). This method has also been extended to 2D flows on structured (Alisadeghi et al., 2011a) and unstructured (Alisadeghi et al., 2011b) grids.
Darbandi (Darbandi et al., 1997) proposed a pressure-based all speed method similar to the Karimian's method, in which the momentum components ( fu   ) are used as the flow dependent variables instead of the velocity components ( u ). They used the following form of the momentum equation as the physical constraint on the integration point convected velocity: In this method, like all other methods discussed so far, Eq. (25) is used as the closure for e p . Numerical results have shown that this set of closure equations removes the possibility of a checkerboard solution and provides strong coupling between pressure and velocity fields. Another analysis has shown that the use of momentum-variable formulation, instead of primitive variable formulation, improves the stability and accuracy of the solution especially in the solution of compressible flows with shock waves (Darbandi, 2004).

One-dimensional incompressible flow
Ashrafizadeh et al. (Ashrafizadeh et al., 2009) have shown that it is possible to develop a colocated numerical scheme without resorting to the convecting and convected velocity concepts, originally proposed by Rhie and Chow. The proposed method, called the method of proper closure equations (MPCE), employs a proper set of physically relevant equations to constrain the velocity and pressure at integration points. It has already been successfully implemented and used in the solution of steady 1-D inviscid incompressible and compressible flow problems (Ashrafizadeh et al., 2008(Ashrafizadeh et al., , 2009) and steady 2-D viscous incompressible flow problems on both structured and unstructured grids (Alinia, 2011;Bakhtiari, 2008).
In this method, following the physical-based approach proposed by Schneider and Raw (Schneider et al., 1987) and modifications proposed by Karimian and Darbandi (Darbandi, 1997;Karimian, 1994a), a combination of non-conservative form of momentum and continuity equations is used to obtain an interpolation formula for the ip-velocity: which results in: Here u S and m S are artificial momentum and mass source terms, and a and u  are the ipcross sectional area and a scheme control parameter respectively. There are options available to the user at this stage. Using UDS or CDS Schemes for terms 1 and 2 and 0, 1, 2 u   , as discussed in (Rezvani, 2008), different formulas for e u can be obtained. The most common formula is as follows: However, only the pressure gradient appears in the momentum equation and CDS-based discretization, physically proper for the pressure, results in wiggly solutions. As discussed in (Rezvani, 2010), by taking the divergence of the momentum equation, a pressure poisson equation for the pressure is obtained. The pressure Poisson equation is an appropriate constraint equation for the pressure discretization. Therefore, as compared to Eq. (34), the following closure equation is proposed for the integration point pressure: Where 0 P   is a scheme control parameter. Numerical test results have shown that the method works well with any nonzero positive value for P  (Rezvani-2010).

One-dimensional compressible flow
Extension of the MPCE to compressible flow problems requires two additional tasks. First, the energy equation should also be solved in order to find the temperature filed. Second, depending on the chosen unknown variables, proper linearization is crucially important. Following approaches presented in (Darbandi et al., 2007;Karimian, 1994a), the Newton-Raphson linearization strategy is employed to linearize the nonlinear terms in the balance equations. For example, term pu in the energy equation is linearized as follows: The non-conservative form of the continuity equation in Eq. (6), in the absence of its source term, is proposed as a proper closure equation for the density (Darbandi et al., 1997;Karimian, 1994a). Here, in order to control the stability and accuracy of the code in both subsonic and supersonic regimes, a smart blending factor (  ) is used which makes the equation hyperbolic in supersonic flow regions [Karimian, 1994a]: where e M is local Mach number at the "east" integration point and m=2 is suggested for the superscript m.
The final required closure equation is an interpolation formula for the integration point temperature. Here, the non-conservative form of the energy equation in the absence of the source term is a natural candidate (Darbandi et al., 1997)  are right hand side constant terms. Subscript nb stands for the immediate neighbors of node P (i.e. W and E ). The first superscript refers to the relevant equation and the second superscript points to the relevant physical variable. The assembled linear equations can be solved using any direct or iterative solver. Here, the Gauss elimination method is used to solve the set of simultaneous equations. The coupled set of the algebraic equations might also be solved using a semi-implicit (segregated) solution strategy. www.intechopen.com

Two-dimensional incompressible flow
The MPCE has also been used to solve 2D, incompressible laminar flows on both structured and unstructured grids (Bakhtiari, 2008;Alinia, 2011). The implementation follows the basic philosophy already described in a 1D context. However, more complex interpolation formulas are required to take into account the effect of cell shape and orientation on balance equations. Implementation details are not discussed here for the sake of brevity.

Incompressible flow
Performance of the MPCE for solving incompressible and compressible flows is examined here through two different test cases. Steady incompressible flow in a constant area duct is the first test case which is used to examine the pressure-velocity coupling in MPCE. Quasi 1D flow in a convergent-divergent nozzle is another test case used to check the accuracy of the proposed method in inviscid incompressible flow problems.

Steady incompressible flow in a constant area duct
The 1D incompressible flow in a constant area duct is used to see if oscillatory solutions appear when the MPCE equations govern the discrete flow field. For this purpose, artificial element and volume sources/sinks are used in continuity and momentum equations. In a domain with 30 nodes, the artificial element source terms are located at the 10th and 20th cell-faces and control volume source terms are activated at the 10th and 20th nodal positions.
Numerical results show that the effects of sharp variations in velocity and pressure due to the artificial source terms are quickly damped out and MPCE works well even with 0 p   . Figures 5 and 6 show the effects of control-volume artificial sources on the numerical solution.
These results show that numerical oscillations are developed in both velocity and pressure fields for 0 p   . However, oscillations are strongly damped for any non-zero value of p  .
Note that the notation "N-D" in these figures stands for "Non-Dimensional" and velocity is normalized with the inlet velocity and pressure is normalized using the exit pressure.

The quasi-1-D test case
Ideal flow in a converging-diverging nozzle is another test case to examine the performance of MPCE. Figure 7 compares the non-dimensional velocity and pressure distributions along the nozzle with the analytical solutions. Excellent agreement is observed.

Compressible flow
For compressible flow, performance of MPCE is examined in three different test problems. Steady subsonic compressible flow in convergent-divergent nozzle, steady compressible flow with normal shock in convergent-divergent nozzle and unsteady compressible flow with normal shock and expansion waves in a shock tube.

Steady subsonic compressible flow in a convergent-divergent nozzle
Results of solving subsonic compressible flow in a symmetric convergent-divergent nozzle with aspect ratio 2 using 51 nodes is shown in this section. For the boundary condition implementation, mass flow rate and the static temperature are specified at the inlet and the www.intechopen.com static pressure is provided at the outlet. In order to assess the accuracy of the proposed method, numerical solutions corresponding to inlet Mach numbers of 0.05, 0.1, 0.15, 0.2, and 0.25 are presented in Fig. 8. Flow with the inlet Mach number of 0.3059, for which the nozzle is choked, is also considered. In this test case, 1 p   is used and pressure, temperature, and density are nondimensionalized with the outlet pressure, inlet temperature, and inlet density, respectively.

Steady flow with shock waves
In the divergent section of the previous test case, normal shock waves appear for certain inlet boundary conditions while the position and the strength of the shock can be controlled by regulating the back pressure. Figure 9 presents the Mach number, non-dimensional pressure, temperature and density distributions along the nozzle for the pressure ratios ( outlet inlet PP ) of 0.65, 0.75 and 0.85 using 81 nodes. It is seen that the normal shock is well captured using a few nodes in all cases and no numerical oscillation is developed along the nozzle. Note that this results are computed using 2 m  in Eq.(43) and p  takes values lower than 1 in the supersonic parts of the solution domain and values higher than 1 in the subsonic regions.

Unsteady compressible flow in a shock tube
As the final unsteady 1D test case with both compression and expansion waves, shock tube problem is discussed. The geometry shown in Fig.10 is a constant area duct divided into two regions separated by a diaphragm at the middle of the duct. At 0   the gas temperature is 298 K, the pressure at the right side of the membrane is 100 KPa and the pressure at the left side of the membrane is 1000 KPa. The time step for the marching numerical solution is equal to 6 41 0    seconds. Figure 11 shows the numerically calculated distributions of Mach number, non-dimensional pressure, temperature, and density after 3 0.42 10   seconds. Here, the Courant number is 0.229 which is computed using the maximum magnitude of velocity in the domain. The pressure, temperature and density are non-dimensionalized using the low pressure side values. Good agreement with the exact solution is observed.

Two-dimensional incompressible flow test cases
The MPCE has also been used to solve many 2D benchmark test problems. Here, only the results for two standard test cases are presented.

Flow in a 2D lid-driven cavity
Flow in a lid-driven cavity, shown in Fig. 12, is a widely used test case for 2D steady incompressible flows governed by Navier-Stokes equations. Numerical solution via MPCE on an unstructured grid with 12656 triangular elements, shown in Fig. 13, is compared to the numerical results of Erturk (Erturk, 2005). Flow patterns at different Reynolds numbers are shown in Fig. 14. Excellent agreement between the calculated velocities at cavity centerlines and the benchmark results in Fig. 14 is also observed.

Flow over a backward facing-step
The backward step flow is another commonly used test case for numerical solution algorithms. The geometry of this test case is shown in Fig. 15. Forty nodes are used across the channel, 300 nodes are used from the step to 30h and 50 nodes are used from 30h to 50h along the channel. Nodes are distributed uniformly. Figure 16 shows the calculated streamlines at Re=800, in which a recirculation zone appears at the upper wall. Normalized locations of the re-attachment point at the lower wall (X1), the separation point at the upper www.intechopen.com wall (X2), and the re-attachment point at the upper wall (X3) are compared with the results obtained by Gartling (Gartling, 1990) for flow at Re = 800. Comparison results are shown in Table 1 and excellent agreement is obvious (Erturk, 2008). A comparison has also been made between numerical results via MPCE and the benchmark results in (Gartling, 1990) and shown in Fig. 18. Again, excellent agreement is observed. Similar calculations and comparisons have also been made for a backward-step flow model with an inlet channel. Excellent results, not shown here, are also obtained in this test case.

Conclusion
Different pressure-based finite volume discretization schemes for the solution of fluid flow problems are explained and compared in the context of co-located grid arrangements and a new scheme, called the Method of Proper Closure Equations (MPCE), is proposed. The proposed method is the only discretization scheme on co-located grids which employs only one definition for the integration point velocity and successfully solves incompressible as well as compressible flow problems. Simple one-dimensional test cases are presented to discuss different available schemes and the building blocks of MPCE. Implementation of the MPCE on 2D structured and unstructured grids has also been carried out but the details are not reported here for the sake of brevity. Numerical solutions via MPCE, however, are presented for both 1D and 2D problems which clearly show the satisfactory performance of the MPCE in the numerical solution of fluid flow problems.