Improving the robustness of the control volume finite element method with application to multiphase porous media flow

Control volume finite element methods (CVFEMs) have been proposed to simulate flow in heterogeneous porous media because they are better able to capture complex geometries using unstructured meshes. However, producing good quality meshes in such models is nontrivial and may sometimes be impossible, especially when all or parts of the domains have very large aspect ratio. A novel CVFEM is proposed here that uses a control volume representation for pressure and yields significant improvements in the quality of the pressure matrix. The method is initially evaluated and then applied to a series of test cases using unstructured (triangular/tetrahedral) meshes, and numerical results are in good agreement with semianalytically obtained solutions. The convergence of the pressure matrix is then studied using complex, heterogeneous example problems. The results demonstrate that the new formulation yields a pressure matrix than can be solved efficiently even on highly distorted, tetrahedral meshes in models of heterogeneous porous media with large permeability contrasts. The new approach allows effective application of CVFEM in such models.


INTRODUCTION
The classical approach for simulating multiphase flow in porous media discretises space using fixed grids and solves the governing equations using finite-volume methods and the 2-point flux approximation. 1,2 Although widely used, complex geometries are often poorly represented in this approach because of the limitations of the 2-point flux approximation. Alternative methods use unstructured meshes to discretise space and control volume finite element methods (CVFEMs) or variants thereof to solve the governing equations (see the references herein [3][4][5][6][7][8][9][10][11][12][13][14] ). An unstructured mesh better captures the often complex, multiscale, and time-dependent solution fields such as pressure, saturation, or composition, especially when used in conjunction with dynamic adaptive mesh optimisation. 11,13,14 In most CVFEM schemes, pressure and velocity are discretised finite element (FE)-wise, whereas saturation and composition are discretised control volume (CV)-wise to enforce mass conservation. However, it is well known, although rarely reported in the literature, that solving the pressure equation in such schemes can be very challenging if there are elements with large angles. The pressure solution may fail to converge, or the results may be physically unrealistic. [15][16][17] A modified CVFEM was proposed by Larreteguy 18 to deal with obtuse elements, but author states that the method should only be used in local regions of the mesh; it should not be used as a global fix for a poor mesh. Since mesh quality is of such importance, a great deal of effort has been devoted to creating good quality meshes (eg, the previous studies [19][20][21][22][23][24] ). However, unless an extremely fine mesh resolution is used, large-angle elements may be inevitable when discretising the very high aspect ratio domains frequently found in geological models (eg, the previous studies [25][26][27][28]. Large-angle elements may also be created at displacement fronts when using dynamic adaptive mesh optimisation. 11,13 Here, a variation of the classical CVFEM is presented, termed the double control volume finite element method (DCVFEM). In this new method, pressure is discretised CV-wise rather than FE-wise. All other variables are discretised following the classical CVFEM. The approach conserves mass and, in the implementation reported here, time is discretised using an adaptive implicit method that provides efficient convergence using highly anisotropic meshes without reducing the time-step size. [29][30][31] However, it is likely that the new DCVFEM could be implemented with other time-stepping approaches.
The new method is first tested to prove that it can reproduce multiphase flow in a number of simple test cases. Its robustness is then demonstrated using example models containing complex geometries with very large aspect ratio domains and large permeability contrasts, discretised using highly distorted triangular and tetrahedral meshes. The example models represent various complex geological reservoirs and were constructed using a surface-based approach in which heterogeneous domains are represented as volumes bounded by surfaces. 11,[32][33][34] However, the new formulation is not restricted to geological problems. It is suitable for any multiphase porous media flow applications in which heterogeneities, or approaches such as dynamic mesh optimisation, give rise to highly distorted meshes.
The presented results demonstrate that the new DCVFEM creates a good quality pressure matrix that can be solved using widely available solvers even when the mesh contains a significant proportion of elements with very large angles. The new method is novel and significant because it provides a solution to a commonly encountered problem in the use of CVFEM with unstructured meshes to model the large group of problems that include multiphase flow in highly heterogeneous domains.

A DOUBLE CONTROL VOLUME FINITE ELEMENT METHOD
This paper builds on the formulation presented in Gomes et al 35 and focuses on the element pair P n DGP n+1 (FE), in which the velocity is discontinuous across elements and is represented with 1 • less than the continuously represented pressure. However, in the classical CVFEM, pressure is represented element-wise and varies continuously across an element ( Figure 1A), whereas in this paper, the pressure is represented using CVs and is piecewise constant within a CV ( Figure 1B). The notation used for the new method reported here is therefore P n DGP n+1 (CV). As this is the sole difference between the classical CVFEM of Gomes et al 35 and the new DCVFEM reported here, only a summary of the equations and classical formulation will be presented prior to explaining the new DCVFEM with CV pressure representation. We report here how the CV-based pressure system is formulated; other details of the method can be found in Gomes et al. 35 Darcy's law for multiphase flow can be written as where v is the force density, u is the phase saturation-weighted Darcy velocity of phase . p is the pressure and s u α is a source term, which may include gravity and/or capillary pressure. is defined as where K is the permeability tensor, K r α , , and S are the relative permeability, viscosity, and saturation of phase , respectively. In classical CVFEM, Equation 1 is discretised using an FE-based representation for velocity and pressure. Here, only the velocity is discretised using FE; the pressure is discretised using CV shape functions: where Q j and P j are FE Lagrangian basis functions and CV Lagrangian basis functions, respectively. N u and N p are the number of degrees of freedom for the velocity and pressure, respectively. Each component of the weak form of the force balance Equation 1 is tested with the FE Lagrangian basis function to obtain: ∑ where n is the outward-pointing normal vector, Ω E and Γ E are the volume and boundary of element E, respectively, and Γ Ω is the boundary of the computational domain. The last term in Equation 4 is used to weakly enforce the pressure boundary condition, p bc , on the domain boundary.
In matrix form, Equation 4 is where v and p are the solution vectors of velocity and pressure, respectively, defined as N α denotes the number of phases; v x , v y , and v z are the components of the velocity associated with each Cartesian dimension (x, y, and z). The mass matrix M, gradient matrix C, and source vector s u are defined as [ in which k and l are the degrees of freedom of element E and i and j are dimensions. n (E,l) is the normal to the boundary of element E. Assuming the incompressible flow, the saturation equation can be written as where is the porosity and s cty is a source term. This equation is discretised as described in Gomes et al 35 using CVs for the space discretisation and the implicit − method for the time discretisation. Equation 6 is discretised and summed over all the phases, obtaining Equation 7 is also bounded by the constraint: Solving for u n+1 in matrix form, Equation 7 becomes The pressure matrix is obtained by manipulating the matrix systems formed by Equations 5 and 9. First, Equation 5 is multiplied by The mass matrix (M) is block diagonal. Therefore, its inverse is not very expensive to calculate. Moreover, if a lumped mass matrix is used, the calculation of the inverse is straightforward as the matrix becomes diagonal. Next, Equation 10 is multiplied by B T : Combining Equation 11 with Equation 9, to remove the dependency from the velocity, the pressure system is obtained: The pressure matrix is, therefore, formed by B T M −1 C. This pressure matrix needs to be solved using a linear solver. Here, the package PETSc is used as the solver tool. 36 More specifically, to solve the pressure matrix, a GMRES solver with 30 restart iterations and a multigrid preconditioner 37 is used. The preconditioner is either hypre 38 or GAMG 39 both with their default PETSc configurations.

NUMERICAL EXPERIMENTS
The presented formulation is first tested against the semianalytical Buckley-Leverett solution (Model 1). In Model 2, the efficiency to solve the pressure matrix of the DCVFEM is tested against the classical CVFEM in a 3D homogeneous domain with single phase flow. Next, the efficiency to solve the pressure matrix is tested using two 3D test cases with very distorted meshes (see Table 1), one with isotropic (Model 3) and one with anisotropic permeability (Model 4). Models 1, 3, and 4 simulate 2-phase flow and the Brooks-Corey model for relative permeability 40 is used: where S w and S nw are the wetting and nonwetting phase saturations, respectively. The irreducible wetting phase saturation (S wirr ), Brooks-Corey exponents (n w and n nw ), viscosity ratio (M 0 ) defined with the viscosity of the wetting phase in the numerator, end-point relative permeability (k w max and k o max ) and residual nonwetting phase saturation (S nwr ) are defined in Table 2. Units are SI except for permeability that is reported in miliDarcy (1mD = 9.869233 × 10 −16 m 2 ). For all multiphase test cases, the domain is initially saturated with the nonwetting phase at (1 − S wirr ). The wetting phase is injected through the inlet boundary, displacing the nonwetting phase towards the outlet boundary. All other boundaries have no  flow. For Models 1, 2, and 4, the wetting phase is injected at a constant rate u in ; for Model 3 injection is at constant pressure yielding a fixed pressure drop Δp. Porosity and horizontal (K h ) and vertical (K v ) permeabilities are specified in Table 3 for the different test cases and domains R1 to R3 of the heterogeneous test cases.

Model 1: 1D immiscible displacement
In the first numerical experiment, the DCVFEM is tested against the 1D Buckley-Leverett (BL) solution, 41 in which a nonwetting phase is displaced by a wetting phase. First, the element pair P 0 DGP 1 (CV) is tested. The saturation profile using different meshes is shown in Figure 2, and error metrics are shown in Figure 3. Likewise, the element pair P 1 DGP 2 (CV) is tested. Figure 4 shows the saturation profile using different meshes. Both element pairs are able to capture the shock front (sharp saturation change), the rarefaction (smooth saturation change behind the shock front), and conserve mass. Figure 3 shows the corresponding errors compared with the classical CVFEM. Figure 3A shows that for a given resolution and using the element pair P 0 DGP 1 , the classical formulation CVFEM yields a more accurate result, although the difference in accuracy between the methods is small for a given mesh resolution. Figure 3B, on the other hand, shows that for a given resolution and using the element pair P 1 DGP 2 ,    the new DCVFEM formulation yields a more accurate result although the difference is again small. Note that flux limiters (see Gomes et al 35 ) revert to a low accuracy (first-order method) advection scheme around the shock front, so in this example problem, the DCVFEM yields slightly higher accuracy than the classical CVFEM for quadratic elements. The most significant finding here is that the new DCVFEM method has a similar quality of match to a known solution as the conventional CVFEM. These experiments show that the new DCVFEM formulation can correctly model multiphase flow. As in Gomes et al, 35 the results show close-to-linear convergence (Figure 3), which is the best that can be expected due to the requirement of flux-limiting to correctly capture the shock-front. 42 This result is important because the close-to-linear convergence obtained in Gomes et al 35 is very uncommon in the literature. Usually the observed convergence is 0.5 (see, for example, Schmid et al, 10 Abushaikha et al, 12 and Hoteit and Firoozabadi 43 ). Discretising the pressure using CVs does not affect the order of convergence of the method.

Model 2: 3D homogeneous case
In this case, the new DCVFEM is tested against the classical CVFEM in a test case in which both numerical methods can solve the pressure system using the iterative solver GMRES with 30 restart iterations with hypre as preconditioner. The domain has homogeneous material properties and only one fluid phase is considered. The aspect ratio of the problem is very high (1:1000), and the mesh has been distorted intentionally to obtain a poor quality mesh with a 33.33% of the elements with angles above 175 • (Table 1).
Three different meshes are tested to evaluate the efficiency of the proposed method. Mesh A is defined by 3000 elements. Meshes B and C are created by refining-by-splitting mesh A. Mesh B is defined by 24 000 elements and mesh C by 192 000 elements. Figure 5 shows the residual reduction versus the number of iterations for the different meshes and element pairs using the classical CVFEM and the novel DCVFEM. The DCVFEM can always achieve convergence and the number of iterations to  reduce the residual; 10 orders of magnitude is much smaller than when using the classical CVFEM, requiring 1 or 2 orders of magnitude fewer iterations. Moreover, Figure 6 displays the asymptotic convergence rates for both methods. The DCVFEM shows a significantly better asymptotic convergence rate than the classical CVFEM.

Model 3: 3D test case representing a common subsurface geological reservoir containing channelized sandbodies
This numerical experiment uses a heterogeneous 3D model that captures channelized features typical of many subsurface geological reservoirs, including hydrocarbon reservoirs, aquifers, and targets for geological CO 2 storage. The model contains 3 types of channelized sandbody, each with a different value of (isotropic) permeability, in a zero permeability mudstone background (not shown) (see Figure 7A and Jacquemyn et al 44 ). The aspect ratio of the problem is very large ( Figure 7A). Moreover, the subdomains created because of the intersections between different channels have even higher aspect ratios ( Figure 7B). Consequently, high aspect ratio ("sliver") elements with large angles are introduced in the mesh unless a very large number of small elements are used. Moreover, there are large contrasts in permeability between domains (Table 3). Such features are typical of subsurface reservoirs it is important that models can capture their geometries and correctly simulate multiphase flow.
Three different meshes are tested to evaluate the efficiency of the proposed method. Mesh A is defined by 12 854 elements (Figure 7). Meshes B and C are created by refining-by-splitting mesh A. Mesh B is defined by 102 832 elements and mesh C by 824 464 elements. The mesh quality in this numerical experiment is characterized by the angles defined in Table 1. It can be seen that the elements are very distorted, with the majority having an angle >120 • and a significant number >170 • . Figure 8 shows the residual reduction versus the number of iterations for the different meshes and element pairs using GMRES with 30 restart iterations together with hypre ( Figure 8A) or GAMG ( Figure 8B) preconditioners. It can be seen that the number of iterations is very low in all cases, especially using hypre as the preconditioner. Moreover, Figure 9 shows the asymptotic convergence rate for the different methods/resolutions. Using hypre, the convergence rate is almost independent of the mesh  used and there is only a small variation using GAMG. The convergence rate is especially good for the P 0 DGP 1 (CV) element pair. However, even the high-order element pair P 1 DGP 2 (CV) shows mesh independent performance.
It is important to note that using the classical CVFEM approach, no convergence could be achieved using this set of iterative solvers. With the new method presented here, the asymptotic convergence rates obtained are very good considering the quality of the mesh and the significant permeability contrasts between domains.

Model 4: 3D test case representing a common subsurface geological reservoir deposited in a shallow marine environment
We finish testing the new method by using a second heterogeneous 3D model of complex but common geological features to demonstrate that the method is robust. In this case, the model contains 3 different domains with different material properties, representing different types of sandstone deposited in a shallow marine environment. 26,27,32 The geometry of the domains is partially controlled by curved, inclined surfaces termed clinoforms. These clinoform surfaces intersect and create wedge-shaped, high aspect-ratio features. Meshing algorithms typically create large-angle elements around these features. The contrasts in material properties across domains are again very large (see Table 3). Moreover, in this example, the permeability in the different domains is anisotropic. Figure 10 shows the numerical domain with the corresponding dimensions, regions and mesh (36 160 elements).
Three different meshes are tested to evaluate the efficiency of the proposed method. Mesh A is defined by 36 160 elements ( Figure 10). Meshes B and C are created by refining-by-splitting mesh A. Mesh B is defined by 289 280 elements and mesh C by 2 314 240 elements. The mesh quality is again characterized by the angles defined in Table 1. Similar to the previous test case, the elements are very distorted. In this case, there are more elements with angles >150 • and also more elements with angles >170 • , making this numerical experiment more challenging. Figure 11 shows the residual reduction versus the number of iterations for the different meshes and element pairs using GMRES with 30 restart iterations together with hypre ( Figure 11A) and GAMG ( Figure 11B). In all cases, the number of  iterations is below 200. Moreover, when using hypre as a preconditioner, the number of iterations is always below 60. However, contrary to the previous numerical experiment, the convergence here does depend on the mesh. Nonetheless, considering the difficulty of the numerical experiment, the results obtained are very satisfactory, as the number of iterations is very low even when the mesh has more than 2 million elements. Figure 12 shows the asymptotic convergence rate for the different FIGURE 13 A, View from below of the pressure field using the classical control volume finite element method approach. A nonphysical pressure drop can be seen close to the right boundary. B, View from below of the pressure field obtained using the double control volume finite element method approach. C, 20× vertically exaggerated view of the model from below displaying the high aspect ratio wedge-shaped domains defined by the intersecting clinoforms [Colour figure can be viewed at wileyonlinelibrary.com] methods/resolutions. The asymptotic convergence is not as good as in Model 2 and, in some cases, is above 0.9. However, as seen in Figure 11, the convergence rate to reduce the residual by 10 orders of magnitude, which is usually sufficient to obtain the correct numerical solution, is reasonable. Moreover, it is important to note that using the classical CVFEM approach, no convergence could be achieved using this set of iterative solvers.
As an example of the robustness of the new DCVFEM approach, the CVFEM was used with mesh A and the low-order element pair P 0 DGP 1 (CV) to obtain the pressure field using a direct solver. Figure 13A shows that the result is clearly nonphysical and therefore incorrect (notice the nonphysical pressure drop near the right boundary). Using the new DCVFEM, the pressure field is as expected ( Figure 13B). Figure 13C shows the connectivity between the different layers of the numerical domain. The new DCVFEM presented here is able to properly calculate the pressure field using an iterative solver with no special modifications. Classical CVFEM formulations do not converge using an iterative solver and yield incorrect solutions using a direct solver.

CONCLUSIONS
A novel, high-order, DCVFEM with n − 1th-order representation for velocity and nth-order for pressure is presented. The new formulation is based on the well-known CVFEM and retains desirable features such as mass conservation. Its implementation in a conventional CVFEM code requires the use of CV shape functions to discretise the pressure. The novelty of the formulation is that it provides significantly better pressure matrices in complex heterogeneous problems. This is important because many problems of interest, particularly those focused on subsurface geological reservoirs, contain complex geometries with very high aspect ratio domains. These produce mesh elements with large angles unless very high resolution is used which is computationally expensive. The new method enables solutions to be obtained using iterative numerical solvers for meshes where classical CVFEM fails. It also allows convergence to be achieved in fewer iterations in models where both methods work.
The new formulation yields accurate solutions with first-order convergence when tested against the classical Buckley-Leverett benchmark. Simulations of heterogeneous test cases demonstrate that new method shows better performance than classical CVFEM, as it is always able to solve the pressure system using a black box iterative solver and good convergence rates are achieved even in models with highly distorted meshes. Moreover, the new method correctly calculates the pressure field in cases where the classical CVFEM fails even using a direct solver.