Comparison of Two- and Three-field Formulation Discretizations for Flow and Solid Deformation in Heterogeneous Porous Media

We illustrate the advantages and disadvantages among twoand three-field formulations that are used to mimic the flow and solid deformation in heterogeneous porous media. Local mass conservation, flux approximation, average permeability alteration at each time step, and degrees of freedom (DOF) are utilized to evaluate each method. Our result presents that four out of six methods provide the local mass conservative while three out of those four methods produce approximately the same flux approximation and permeability alteration. Three-field formulation methods generally require a smaller time step to converge for solving the system of nonlinear equations. Besides, they have more DOF than that of the two-field formulation because they have one more primary variable, i.e. fluid velocity. The two-field formulation that is a combination of continuous and enriched Galerkin function space enjoys all the benefits while requires the least DOF among the methods that preserve local mass conservation property. 1 Statement of problem Let Ω ⊂ Rd be the domain of interest in d-dimensional space bounded by boundary, ∂Ω. ∂Ω can be decomposed to displacement and traction boundaries, ∂Ωu and ∂Ωt, respectively, for the solid deformation problem. For the fluid flow problem, ∂Ω is decomposed to pressure and flux boundaries, ∂Ωp and ∂Ωq, respectively. The time domain is denoted by T = (0, T ] with T > 0. The coupling between the fluid flow and solid deformation can be captured through the application of Biot’s equation of poroelasticity, which is composed of force and mass balance equations (Biot, 1941). The linear momentum balance equation is read as follows: ∇ · 2μlε−∇λl∇ · uI + α∇p = f in Ω× T, (1) u = uD on ∂Ωu × T and σ · n = σD on ∂Ωt × T, (2) where λl and μl are are Lamé constants, ε is strain, u is displacement vector, I is the secondorder identity tensor, α is Biot's coefficient defined as α = 1 − K KS , K is the bulk modulus of a rock matrix, Ks is bulk rock matrix material (e.g. solid grains), p is fluid pressure, n is a normal unit vector, f is body force vector, which is neglected in this study, σ is stress, uD and σD are prescribed displacement and traction at boundaries, respectively. The second equation is the mass balance equation, which is read as: ρ ( φcf + α− φ Ks ) ∂p ∂t + ρα ∂∇ · u ∂t −∇ · (ρv) = g in Ω× T, (3) v = ∇ · κ ρ (∇p− ρg) in Ω× T, (4) p = pD on ∂Ωp × T and −∇ · κ(∇p− ρg) · n = qD on ∂Ωq × T, (5) where ρ is fluid density, φ is initial porosity, cf is fluid compressibility, t is time at the specific point, v is Darcy velocity vector, g is a gravitational vector, g is sink/source, pD and qD are specified pressure and flux, respectively, and κ is hydraulic conductivity defined as κ = ρkm μ , km is matrix permeability tensor and μ is fluid viscosity. The rock displacement caused by changing in pore pressure can influence media conductivity by: km = km0 (1 + εv/φ) 3 / (1 + εv) , (6) where km0 and km represent initial and current rock matrix permeability, respectively, and εv is a volumetric strain (Du and Wong, 2007). There are six discretizations compared in this study. Three of discretizations are two-field formulation (u×p), and the rest are three-field formulation (u×v×p) as presented in Table 1. The details of two-field formulation discretization can be found in Lee et al. (2016), Choo and Lee (2018), and Kadeethum et al. (2019), while Haga et al. (2012) presents the details of threefield formulation discretization. The solving algorithm, solver settings, and coupling scheme can be found in Kadeethum et al. (2019). The Picard iteration is used to solve the nonlinear pressure-dependent km model, and the PETSc LU solver is used inside each iteration for solving the system of linear equations (Balay et al., 2018). Formulation Mixed space Displacement Velocity Pressure u×p CG2 × CG1 CG2 CG1 CG2 × EG1 CG2 EG1 CG2 ×DG1 CG2 DG1 u×v×p CG2 × RT1 ×DG0 CG2 RT1 DG0 CG2 × CG2 × CG1 CG2 CG2 CG1 CG2 × CG + 2 ×DG1 CG + 2 CG + 2 DG1 Table 1: Two(u×p) and three-field (u×v×p) formulations are presented their choices of function space. CGk, EGk, DGk, RTk, and CG + k refer to continuous Galerkin, enriched Galerkin, discontinuous Galerkin, Raviart-Thomas, and continuous Galerkin enriched by bubble function with polynomial degree k approximation function space. The classical form of CG2 × CG1 method is well-known to suffer from lack of the local mass conservation and producing unphysical (spurious) pressure oscillation (Haga et al., 2012; Choo and Lee, 2018; Kadeethum et al., 2019). This problem affects not only velocity calculation and coupling with the transport equation but also rock porosity and permeability alterations, which are determined by either the effective stress or the volumetric strain. This inaccuracy influences both the reservoir compaction and the conductivity alteration effects on the media's ability to deliver or withhold the desired fluid. Even though three-field formulation was developed to mitigate this problem, the unphysical pressure oscillation still occurs in CG2 × CG2 × CG1 space (Haga et al., 2012). In this study, we evaluate the performance of the six aforementioned methods by comparing: 1) the local mass conservation property, 2) the flux approximation, and 3) the computational cost, i.e. degrees of freedom (DOF) and maximum time step, ∆t. 2 Results and discussion The example presented in this section is adapted from Kadeethum et al. (2019), and its geometry and boundary conditions are presented in Figure 1. The first investigation focuses on the local mass conservation, rmass, which is calculated as: