3.4.1. Review of Control Volume Finite Element Methods

The earliest reference to control-volume finite-element methods is the 1980 work by Baliga and Patankar [133] for the convection-diffusion equation, a refinement of Baliga’s 1978 dissertation [134]. Baliga and Patankar [135] first apply their approach to the Navier-Stokes equations of fluid mechanics in 1983. At the same time, Schneider and Zedan [136] develop a control-volume finite-element for heat conduction. Schneider and Raw [137] then develop a control-volume finite-element method for fluid flow in 1986. The work of Baliga/Patankar and Raw/Schneider are the foundations for two of the main control-volume finite-element methods that are used today for fluid mechanics. A third control-volume finite-element method is adapted from Galerkin Least Squares (GLS) finite-element methods by Swaminathan and Voller [138] in 1992, but there is no evidence of widespread use.

There are three difficult issues that must be addressed in all numerical methods for the Navier-Stokes equations: 1) stability at high Reynolds number and Peclet numbers, where pure centered differencing for the convection terms, or the analogs in FEM and FVM, can lead to numerical oscillations, 2) coupling of the pressure and velocity field, where “checker-boarding” can occur when the variables are co-located and use similar interpolations, and 3) updating of the pressure field.

There are three main schools of thought in the CVFEM community for addressing the three issues above. With the Baliga/Patankar approach, upwinding is achieved with exponential shape functions on linear triangular and tetrahedral elements. Originally, pressure-velocity coupling was attained using mixed-order elements. Later, an equal-order scheme was developed that involved pressure terms in the interpolation functions. Convecting and convected velocities were maintained for pressure-velocity coupling. The pressure is solved using a projection method similar to the SIMPLER [139] algorithm. The method is practically limited to triangles and tetrahedra because of the form of the interpolation functions.

With the Raw/Schneider approach, upwinding is achieved using the skewed upwinding or positive influence coefficient approaches. The pressure and velocity are solved fully coupled using an approximation of the transport equations as an interpolation function. Two velocity fields are maintained, a convecting and a convected field. The method is applicable to all element forms and has been successfully implemented in a commercial computational fluid dynamics code, TASCflow [140].

With the Swaminathan/Voller approach, the methods of streamline upwinding and pressure stabilization are adapted from finite-element methods. There is only a small amount of literature on this particular CVFEM.

The following historical synopsis of CVFEM’s addresses research for solving the pressure-based incompressible Navier-Stokes equations and the more elementary convection-diffusion equations.

1966 Winslow [141] presents a control volume formulation for a Poisson equation based on linear triangular elements. This work is important because it is one of the first applications of the finite volume method on unstructured meshes.

1980 Ramadhyani and Patankar [142] compare the accuracy of the Galerkin finite element method with a control volume method for the Laplacian operator. They use bilinear shape functions and rectangular elements, where the control volume method uses the bilinear shape functions as interpolation functions. The numerical errors of the control volume method are half those of the finite element method.

1980 Baliga and Patankar [133] introduce a flow-oriented upwind interpolation for convection-diffusion problems on triangular elements, a refinement of 1978 dissertation work [134]. The upwinding is introduced through an interpolation function based on a locally analytic solution to the velocity-aligned transport equation, resulting in exponential shape functions. They solve both radial heat conduction in a rotating hollow cylinder for Peclet numbers up to $100$ , and the transport of a step scalar field, all with specified velocity fields. The directional upwinding provides better solutions than uniform first-order upwinding.

1983 Baliga and Patankar [135] develop a mixed-interpolation scheme for solving the Navier-Stokes equations with heat transfer on triangular elements. The mixed interpolation keeps the pressure from decoupling from velocity. The pressure is solved by applying the continuity equation over macro-triangles. The interpolation function for the convecting velocity contains the pressure gradient. Each macro-triangle is subdivided into four sub-triangles for the momentum and energy equations. The flow-oriented upwind scheme is used to interpolate velocity and temperature for their respective transport equations. The velocity is assumed to vary linearly over the element for computing mass flow rates. The equations are solved in a segregated manner using an approach similar to the SIMPLER method [139]. This work is the first application of the CVFEM for the Navier-Stokes equations.

1983 Baliga, Pham, and Patankar [143] apply the mixed-interpolation scheme [135] to fluid flow and heat transfer. They solve flow between rotating cylinders for Reynolds numbers up to $1000$ , fully developed flow in a square duct with a laterally imposed velocity for Reynolds numbers up to $100$ , natural convection in rectangular enclosures for Rayleigh numbers up to $10^5$ , and natural convection in a trapezoidal enclosure for Rayleigh numbers up to $10^6$ .

1985 Prakash and Patankar [144] solve the Navier-Stokes equations with an equal-order interpolation for velocity and pressure on triangular elements. The mass flow velocity, used for continuity, is different from that derived from momentum, thus avoiding staggering or mixed-interpolation. The flow-oriented upwind scheme is used to interpolate the convected velocity for the momentum equation while the pressure gradient is treated as an element-constant source term. The coefficient matrices for momentum are used to define the velocities for the continuity equation which include the now-unknown pressure gradient across control volume faces. They use a pressure correction approach similar to the SIMPLER algorithm to update velocity and accelerate convergence. The continuity and momentum equations are segregated in the solution process. They solve flow between rotating cylinders for Reynolds numbers up to $1000$ , and natural convection in a closed cavity with a Boussinesq-type buoyancy term for Grashof numbers up to $10^5$ . The solutions are more accurate than with the mixed interpolation scheme of Baliga [135]. They note problems with negative coefficients during the first iterations of a solution.

1985 Ramadhyani and Patankar [145] extend the flow-oriented upwind interpolation scheme from linear triangles to bilinear quadrilateral elements for convection-diffusion problems. Three-point quadratures (Simpsons Rule) are used to evaluate flux integrals, as in all the previously mentioned work. They argue that one-point quadratures are less accurate because of the nonlinear nature of the interpolation functions, but only at intermediate values of cell-Reynolds number. They present solutions for five different test cases, including the convection of scalar profiles and diffusion in rotating systems. After this article, there are no further publications for quadrilateral or hexahedral elements using methods developed by Baliga, Patankar, and Prakash.

1986 Schneider and Raw [137] develop a positive influence-coefficient extension to skewed upwind interpolation [146] for convection terms, based on 1985 dissertation work [147]. They apply the scheme to the convection-diffusion equation on quadrilateral elements. Diffusion terms are calculated by integrating the gradients of the isoparametric, bilinear interpolation functions. They solve several convected-scalar cases and claim smooth solutions where the methods of Baliga and Patankar exhibit oscillations. The skewed upwind method has less dependence on the element orientation than flow-oriented streamline upwinding.

1986 LeDain-Muir and Baliga [148] extend the flow-oriented upwind interpolation scheme to linear tetrahedral elements in three dimensions for the convection-diffusion problem. Each tetrahedron contains six control volume faces. A single unit normal is calculated for each control volume face. For integration, each face is subdivided into two triangles. A three-point quadrature is used on each triangular subface where the sample points are taken along the midpoints of the triangle edges. They solve radial heat conduction in a rotating hollow sphere, scalar transport of a step profile, and transport with radial convection between concentric spheres.

1986 Prakash [149] modifies the flow-oriented upwind interpolation to include source terms from the transport equations on triangular elements, with applications to the incompressible Navier-Stokes equations. The pressure gradient in the momentum equations is treated as a source term in the interpolation function for velocity, directly coupling the pressure to the velocity. The source term has a streamwise-linear influence on the interpolation function. The mass flux is calculated using the new interpolation function instead of assuming a linear variation. The pressure is then calculated through the continuity equation by directly applying the new velocity interpolation function, replacing the SIMPLER scheme but keeping the segregated approach. A pressure correction step is included to make sure the velocity interpolation function satisfies continuity. He solves flow between rotating cylinders up to a Reynolds number of $1000$ , the lid-driven cavity for Reynolds numbers up to $400$ , and natural convection in a square cavity for Grashof numbers up to $10^5$ . The solutions are more accurate than with the original collocated scheme of Praskash and Patankar [144].

1987 Schneider and Raw [150, 151] extend the positive-coefficient, skewed upwind interpolation [137] to the incompressible Navier-Stokes equations on quadrilateral elements. They use a element-local discretization of the transport equations to derive interpolation functions at control volume faces that couple the velocity and pressure. The convection terms are constructed with the positive-coefficient, skewed upwinding. The skewed upwinding couples all the control volume face values together within an element, so an internal matrix inversion must be applied to calculate individual face values. The momentum and continuity equations are solved all at once as a coupled system. They solve convection of a scalar field with a step profile, the lid-driven cavity for Reynolds numbers up to 1000, the inviscid forward-facing step to test the conservation of total pressure, and flow between rotating cylinders. Grid convergence studies suggest spatial accuracy near second order. They call their method Finite Element Difference Scheme (FIELDS).

1987 Schneider [152] extends their algorithm [150] to cylindrical, axisymmetric coordinates and presents solutions for the cylindrical driven cavity.

1987 Prakash [153] examines a donor-cell method for replacing flow-oriented upwind interpolation on triangular elements. The donor cell method provides positive coefficients, where the flow-oriented upwinding can yield negative coefficients, leading to oscillations. The donor-cell scheme is applied to several of the previous convected scalar problems and the thermally driven cavity. The scheme exhibits excessive diffusion and is not generally recommended.

1988 Hookey, Baliga, and Prakash [154] modify the treatment of the source term in the flow-oriented upwind interpolation for triangular elements relative to the previous source term modifications of Prakash [149]. A crossflow-quadratic multiplier is added for the source term in the the interpolation function. They apply the scheme to the convection-diffusion equation for radial heat conduction between rotating cylinders and radial heat conduction in radial flow between cylinders. The new source treatment proves better than the previous scheme of Prakash only when the flow has multidimensional features.

1988 Hookey and Baliga [155] apply the flow-oriented upwind interpolation with the modified source treatment [154] to the incompressible Navier-Stokes equations on triangles. Instead of calculating pressure by applying the interpolation functions directly to the continuity equation as was done by Prakash [149], a method similar to SIMPLEC [156] is used. The previous approach converged poorly at higher Reynolds numbers. A pressure correction approach is still employed to force the interpolation function for velocity to satisfy continuity, but with the penalty of an enlarged stencil for the pressure-correction equation. The continuity and momentum equations are solved simultaneously. They solve a polar lid-driven cavity for Reynolds numbers up to 380 and the natural convection for Rayleigh numbers up $10^6$ . Solutions are compared against results from the older methods of Prakash [149] and Baliga [135], and prove to be more accurate.

1988 Reviews of control volume finite element methods for fluid flow and heat transfer are given in the Handbook of Numerical Heat Transfer by both Baliga [157] and Schneider [158]. They provide implementation details for many of the methods published to date.

1992 Swaminathan and Voller [138, 159] extend of the ideas of the Streamline-Upwind Petrov-Galerkin (SUPG) method [160] to solving the convection-diffusion equation with quadrilateral elements. They solve several convected-scalar problems and compare the results to a FEM implementation of the SUPG scheme. The CVFEM analog of SUPG performs just as well, except for time accurate solutions where the phase error is larger. The SUPG method provides better solutions than the skew upwinding or flow-oriented upwinding for heat conduction between rotating cylinders, but worse for the scalar transport of a step profile.

1992 Baliga and Saabas [161] provide a critical review of control volume finite element methods. They criticize the schemes of Hookey [155] and Raw [150] for being too expensive, computationally. They introduce the mass advection weighted scheme of Saabas where he adapts the concept of positive influence-coefficients from Schneider and Raw to the formulation of Baliga and Prakash. They call the original flow-oriented upwind scheme of Baliga and Patankar FLO, the source-term modified scheme of Prakash FLOS, and the mass advection-weighted scheme of Saabas MAW. The FLO(S) schemes result in mixed-sign off-diagonal coefficients if the triangular elements are obtuse, potentially admitting oscillations. Additionally, many of the schemes developed to date, for CVFEM, over-specify the pressure boundary conditions, leading to poor convergence.

1992 Naterer and Schneider [162] extend the approach of Schneider and Raw [150] to compressible flow. An explicit predictor-corrector time integration is used for transient solutions. The influence-coefficient matrices are used to interpolate density, velocity, and internal energy at control volume faces at an intermediate time level. These values are then used to correct the state variables using a forward Euler integration. They solve a transient shock tube problem for an initial pressure ratio of 10, flow through a converging-diverging nozzle with an area ratio of 2, and Mach 3 supersonic flow over a forward-facing step.

1993 Swaminathan, Voller, and Patankar [163] extend the streamline-upwind Petrov-Galerkin method and the pressure-stabilized Petrov-Galerkin [160] method to a conservative form for the control volume finite element method. The streamline-upwind control-volume, pressure-stabilized control-volume (SUCV/PSCV) method is applied to the incompressible Navier-Stokes equations with quadrilateral elements. They evaluate the integrals using mid-point quadrature and solve the segregated equations using a SIMPLER approach. They solve the lid-driven cavity at a Reynolds number of $400$ , natural convection in a square enclosure for a Rayleigh number of $10^5$ , and natural convection in a cylindrical annulus at a Rayleigh number of $10^4$ .

1994 Saabas and Baliga [164, 165] adapt the positive influence-coefficient scheme of Schneider and Raw [150] to triangular and tetrahedral elements and call the method mass advection weighting (MAW). They introduce a new control volume construction for tetrahedral elements. Their tetrahedral element contains one four-point planar face and two three-point planar faces, whereas the control volume construction of LeDain-Muir [148] contained six four-point surfaces. The reduced number of control volume faces makes the MAW scheme less expensive to apply, but the element shape functions become dependent on the shape of each element. They solve for pressure using the original approach of Prakash [144] with a SIMPLER method. The solution technique is segregated. For solving practical problems, they recommend using the FLO scheme for the convection terms and switching to the MAW scheme only if there are problems with negative coefficients. They advise against using the FLOS schemes of Prakash [149] and Hookey [155] because they typically do not provide enough improvement in accuracy to justify their slower convergence properties. Additionally, they claim that carrying pressure gradient terms in the velocity interpolation function requires the boundary conditions for pressure to be over-specified for inflow/outflow problems. They solve the 2D lid-driven cavity for Reynolds numbers up to $1000$ , 2D turbulent flow over a backward-facing step using a $k-\epsilon$ turbulence model for a Reynolds number of about $10^6$ , 3D natural convection in a cavity for Rayleigh numbers up to $10^6$ , and a turbulent jet injection into crossflow for jet Reynolds numbers up to $53600$ . The MAW scheme is required for the jet problem because of negative coefficient problems with the FLO scheme.

1994 Masson, Saabas, and Baliga [166] extend the MAW scheme of Saabas [164] to axisymmetric flows with triangular elements. They solve developing pipe flow for a Reynolds number of $40$ , pipe flow with a step constriction up to a Reynolds number of $1000$ , natural convection in a cylindrical enclosure for a Grashof number of $2$ and Prandtl number of $2500$ , and flow through and arterial section for Reynolds numbers up to $350$ .

1994 Masson and Baliga [167] apply the MAW scheme of Saabas [164] to dilute-particle flows with triangular elements. They solve equations for the gas phase and the dispersed phase. They solve for flow through a constricted channel for a Reynolds number of $100$ and Stokes numbers between $10^{-3}$ and $10^{-1}$ , and for flow in a split inertial separator for a Reynolds number of $200$ .

1994 Karimian and Schneider [168] improve the velocity-pressure coupling of the original Schneider-Raw scheme [150]. The original scheme, referred to as FIELDS, has poor performance for inviscid flow. They improve the coupling by adding a discrete continuity relation to the interpolation functions for the convecting velocity. The additional terms help smooth oscillations that occur for a mass sink test problem. They verify the new interpolation function on the lid-driven cavity for Reynolds numbers up to $3200$ , and the backward-facing step for Reynolds numbers up to $230$ .

1994 Karimian and Schneider [169] extend the method of Schneider and Raw [150] to both compressible and incompressible flow for the quasi-one-dimensional Euler equations. They solve for flow through a converging-diverging nozzle with an area ratio of $2.035$ with and without a shock.

1994 Deng et al. [170] present a new flux reconstruction scheme to replace the FIELDS scheme of Schneider and Raw [151]. They note that the FIELDS scheme is similar to the original work of Rhie and Chow [171] who where some of the first researchers to solve incompressible flow on collocated grids. Deng takes features of both schemes to create a compact reconstruction that does not require matrix inversions to calculate the integration point values in terms of nodal values. Since they question the consistency of the FIELDS scheme, they call their new scheme consistent physical interpolation (CPI). They apply the scheme to two and three-dimensional Navier-Stokes calculations on structured Cartesian meshes. They solve the lid-driven cavity for Reynolds numbers up to $1000$ , a 3D lid-driven cavity for Reynolds numbers up to $1000$ , and turbulent vortex shedding over a square cylinder for a Reynolds number of $22000$ .

1995 Costa et al. [172] apply the MAW scheme of Saabas [164] to three-dimensional turbulent flows with tetrahedral elements. They solve a turbulent jet injected into a crossflow for jet Reynolds number up to $53600$ , and flow through a T-junction in ducts at Reynolds numbers near $90000$ .

1995 Karimian and Schneider [173] apply control-volume finite-element methods to a shock-tube problem.

1995 Karimian and Schneider [174] extend the FIELDS scheme and the convecting velocity corrections to compressible flow. They solve the lid-driven cavity for Reynolds numbers up to $3200$ , flow over a shallow bump in a channel with Mach numbers from $0.5$ to $1.65$ , and flow through a ramped inlet for a Mach number of $2.5$ .

1995 Padra and Larreteguy [175] develop an error estimator with mesh refinement for the convection-diffusion equation. They use the formulation of Baliga and Patankar [133] with triangles. Larreteguy [176] then extends the scheme to fluid flow with triangles.

1996 Harms et al. [177] introduce a simplified interpolation function for the control volume finite element method. They develop a method for applying analytic shape functions on nonorthogonal meshes. They apply the scheme to flow between rotating cylinders for Reynolds numbers up to $1000$ and the scalar transport of a step profile.

1996 Comini et al. [178] compare CVFEM and GFEM formulations for the convection-diffusion equations.

1996 Neises and Steinbach [179] develop a control volume finite element method based on a mixed interpolation approach to facilitate pressure-velocity coupling and artificial dissipation for convective stability. They begin with a Galerkin finite element method and then manipulate off-diagonal terms to force conservation. They solve a laminar, 3D obstructed channel flow for Reynolds number from $0.1$ to $50000$ .

1996 Volker, Burton, and Vanka [180] apply a multigrid solution technique to a control volume finite element method on triangular elements. Linear interpolation is used throughout the triangles with a three-point quadrature to integrate fluxes. The pressure is solved using the SIMPLE method. They solve natural convection problems in square, triangular, and semicircular cavities for Rayleigh numbers up to $10^6$ .

1996 Botta and Hempel [181] describe a finite-volume projection method for unstructured, triangular meshes with element-centered variables.

1997 Darbandi and Schneider [182] develop a scheme for both compressible and incompressible flow using a momentum variable formulation of the Schneider/Raw scheme [150, 151]. The interpolation formula for the convecting velocities is derived from an approximation of the momentum equation with an additional velocity-weighted continuity equation term. Solutions are demonstrated for velocities up to Mach 0.9.

1997 Baliga [183] gives an overview of the control volume finite element method as applied to fluid flow.

1998 O’Rourke and Sahota [184] develop an edge-based scheme in 3D for the convection operator. The convection operator is constructed from a multidimensional upwind scheme. Within each element, the quadrature points are associated with edge mid-points instead of sub-face mid-points, so the amount of work is reduced over the traditional CVFEM.

1998 Gresho and Sani [185] compare CVFEM methods to GFEM methods.

1998 Venditti and Baliga [186] describe an error estimation strategy for incompressible flow with CVFEM.

2001 Reyes, Rincon, and Damia [187] present a CVFEM approach for turbulent flow with wall functions.

2001 Campos Silva and de Moura [188] present a method for 9-noded quad elements with the mass advection weighted scheme.

2002 Zhao, Tai and Ahmed [189] implement a 2D CVFEM on triangles for micro flows. They use an upwind scheme where nodal gradient are used to reconstruct the high-order fluxes at the control volume faces.

With respect to fluid flow, the CVFEM methods have been developed primarily for triangular and tetrahedral elements [133, 135, 144, 148, 149, 154, 155, 164, 166]. Development focused on triangular and tetrahedral elements because the shape functions are linear and gradient terms become constant over the element. Constant first derivatives simplify the formulation of many of the schemes. Fewer articles have been published on the use of quadrilateral elements [137, 138, 145, 150, 163] in two dimensions and no articles have been published for CVFEM with hexahedral elements in three dimensions. In addition, there have been CVFEM formulations for the streamline-vorticity equations [38, 190, 191, 192, 193, 194], for the heat equation [136, 195, 196, 197, 198, 199, 200, 201, 202], for flow in porous media [203, 204, 205, 206, 207, 208, 209], for overland flows [210, 211], and for linear elasticity [212, 213].