Fluid flow in porous media has been extensively studied, both theoretically and computationally, because of its broad applications in different branches of science and engineering. The most popular model of flow of an incompressible fluid in rigid porous media is the Darcy model, which is based on the assumption that the domain contains only one pore-network. Due to the restricting assumptions in the classical Darcy model [Rajagopal, 2007; Nakshatrala and Rajagopal, 2011; Chang et al., 2017], its application has been limited and several modifications and alternative models have been proposed that predict more realistic flow behaviors. In particular, due to the complexity of the pore-structure in many geo-materials such as shale, many studies have focused on developing mathematical models and computational frameworks that consider the presence of two (or more) dominant pore-networks exhibiting different hydro-mechanical properties. Some of the recent studies on multiple pore-networks include [Borja and Koliji, 2009; Choo et al., 2015].
The mathematical models pertaining to the flow in porous media with multiple pore-networks are complex and involve numerous field variables. It is not always possible to derive analytical solutions to these mathematical models, and one has to resort to numerical solutions for realistic problems. Different approaches are available for developing formulations for multi-field mathematical models. Mixed finite element formulations, which offer the flexibility of using different approximations for different field variables, are particularly attractive for multi-field problems. Accurate numerical solutions have been obtained using mixed finite element for various porous media models; for example, see [Masud and Hughes, 2002; Badia and Codina, 2010; Nakshatrala et al., 2006; Nakshatrala and Rajagopal, 2011; Choo and Borja, 2015]. Moreover, many of the mathematical models pertaining to the multiple pore-networks, and in particular, the mathematical model considered in this paper, cannot be written in terms of a single-field variable. Although mixed methods are considered a powerful tool, especially for modeling flow problems in porous media, they suffer from some restrictions. To obtain stable and convergent solutions, a mixed formulation should satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB) stability condition [Babuška, 1973; Brezzi and Fortin, 1991]
. Numerical instability of the solution and probable spurious oscillations in the profile of unknown variables are the main consequences of the violation of this condition. Such drawbacks are observed in many of the existing formulations and highlight the need for developing more robust computational frameworks. In order to resolve numerical instabilities resulting from violation of the LBB condition, computational approaches are divided broadly into two classes[Franca and Hughes, 1988]: those that satisfy the LBB condition and those that circumvent it.
In the former approach, elements are developed by placing restrictions on the interpolation spaces so as to satisfy the LBB condition under the classical mixed (Galerkin) formulation. Such elements are collectively referred to as the H(div) elements [Brezzi and Fortin, 1991; Brezzi et al., 2008]. Two popular works of this type are Raviart-Thomas (RT) spaces [Raviart and Thomas, 1977], and Brezzi-Douglas-Marini (BDM) spaces [Brezzi et al., 1985, 1987]. The class of stabilized methods, which falls under the latter approach, is an attractive way of circumventing the LBB condition. In a stabilized formulation, stabilization terms are augmented to the classical mixed formulation to avoid a saddle-point problem as well as mathematical instabilities [Hughes et al., 2004]. Various stabilized formulations have been published for flow problems (e.g., see [Badia and Codina, 2009; Brooks and Hughes, 1982; Hughes et al., 2000; Turner et al., 2009]) and for flow problems in porous media, in particular, (e.g., see [Masud and Hughes, 2002; Badia and Codina, 2010; Nakshatrala et al., 2006; Choo and Borja, 2015]).
Herein, we develop a stabilized mixed formulation of the double porosity/permeability model proposed by [Nakshatrala et al., 2016]. The stabilization terms and the stabilization parameter are derived in a mathematically consistent manner by appealing to the variational multiscale formalism [Hughes, 1995]. It is noteworthy that the nodal-based equal-order interpolation for all the field variables is stable under the proposed stabilized mixed formulation. Such a feature for interpolations is particularly desirable for studies in porous media for two reasons. The obvious reason is that the equal-order interpolation is computationally the most convenient. The second reason is that, in many porous media applications, the flow and transport equations are coupled (Section 9 of this paper deals with such a coupled problem). But many existing formulations (including the stabilized formulations) produce non-physical negative solutions for the transport equations (i.e., a negative value for concentration fields), especially when the diffusion/dispersion is anisotropic [Nagarajan and Nakshatrala, 2011]. The known robust non-negative finite element based formulations for the transport equations are nodal-based (e.g., refer to [Nagarajan and Nakshatrala, 2011; Mudunuru and Nakshatrala, 2016]). By choosing nodal-based unknowns even for the flow problem, one can avoid projections from nodal to non-nodal interpolation spaces and vice-versa.
To determine whether a computational framework is robust, systematic convergence and error analyses are required. To this end, we first perform a mathematically rigorous stability analysis of the proposed stabilized mixed formulation. Since the proposed formulation is residual-based, consistency is shown quite easily. We also present patch tests and representative numerical results to show that the obtained numerical results are stable. After establishing the stability of the proposed formulation, we perform a thorough accuracy assessment of the approximations by estimating the error associated with the numerical solutions. Specifically, we perform both a priori and a posteriori error estimations, which individually serve different purposes [Babuška et al., 2010]. A posteriori error estimations monitor different forms of the error in the numerical solution [Becker and Rannacher, 2001; Babuška and Strouboulis, 2001] and using the computed approximate solution, they provide an estimate of the form , where is the solution, is the finite element solution for a mesh with mesh size , denotes an appropriate norm, and is a constant (real) number. On the other hand, a priori error estimations provide us with the order of convergence of a given finite element method [Ainsworth and Oden, 1997].
[Shabouei and Nakshatrala, 2016] have shown that porous media models such as those defined by the Darcy and Darcy-Brinkman equations satisfy certain mechanics-based properties, and they have utilized these properties to construct solution verification procedures. Recently, [Nakshatrala et al., 2016] have shown that the double porosity/permeability model also enjoys properties with strong mechanics underpinning. These include the minimum dissipation theorem and a reciprocal relation. Herein, we utilize these mechanics-based properties to construct a posteriori solution verification procedures to assess the accuracy of numerical solutions obtained under the proposed formulation for the double porosity/permeability model.
Another type of numerical instability, known as Gibbs phenomenon, can also be observed in the numerical solutions of problems associated with flow through porous media with disparate properties. In layered porous domains, conventional continuous finite element methods are not capable of capturing abrupt changes in material properties and result in overshoots and undershoots in the profiles of numerical solutions along the interface of layers where there are jump discontinuities. In order to eliminate such erroneous oscillations, one possible approach is discontinuous Galerkin (DG) methods. DG methods have been successfully employed by[Hughes et al., 2006] for the case of Darcy equations. An extension of the proposed framework using discontinuous Galerkin method for double porosity/permeability model can be obtained using a method similar to the one proposed by [Hughes et al., 2006]. However, obtaining such an extension and comparison between the performance of continuous and discontinuous formulations for capturing abrupt changes in material properties are beyond the scope of this paper and will be addressed in a subsequent one.
A common assumption in models of flow in porous media is that of steady-state conditions. However, many flows occurring in porous media such as aquifers and oil-bearing strata are transient or unsteady in nature. In this paper, we extend the proposed stabilized mixed formulation for the double porosity/permeability mathematical model to the transient case, and we illustrate this extension can accurately capture the transient flow characteristics.
Recently, it has been shown that some stabilized methods (which are primarily designed to suppress numerical instabilities) when applied to solve problems with physical instabilities, suppress both types of instabilities [Shabouei and Nakshatrala, ; Shabouei, 2016]. Therefore, a good test of the proposed stabilized mixed formulation for a coupled flow and transport problem involves a problem that exhibits a physical instability similar to the classical Saffman-Taylor instability [Saffman and Taylor, 1958]. Using numerical simulations we show that the proposed formulation suppresses only the spurious numerical instabilities while capturing the underlying physical instability.
The rest of this paper is organized as follows. After an outline of the governing equations of the double porosity/permeability model in Section 2, the corresponding stabilized mixed formulation is presented in Section 3 with a derivation provided in Appendix A. The theoretical convergence analysis for the proposed stabilized mixed formulation is presented in Section 4, followed by the numerical convergence behavior of the elements presented in Section 5 where patch tests in one- and three-dimensional spaces are described. The representative numerical results are used to showcase the performance of the proposed mixed formulation in Section 6. Section 7 provides the mechanics-based assessment of the numerical accuracy. The transient analysis and the capability of the computational framework for modeling coupled problems and capturing well-known physical instabilities in fluid mechanics are discussed in Sections 8 and 9. Finally, conclusions are drawn in Section 10.
Throughout this paper, repeated indices do not imply summation. The terms classical mixed formulation and Galerkin formulation are used interchangeably.
2. Governing Equations for Double Porosity/permeability
For convenience to the reader and for future referencing, we document the equations that govern the double porosity/permeability mathematical model considered in [Nakshatrala et al., 2016]. Let be a bounded domain, where “” denotes the number of spatial dimensions. The boundary of the domain is assumed to be piecewise smooth. Mathematically, , where the superposed bar denotes the set closure [Evans, 1998]. A spatial point is denoted by . The gradient and divergence operators with respect to are denoted by and , respectively. The unit outward normal to the boundary is denoted by .
The porous domain is assumed to consist of two dominant pore-networks, which will be referred to as the macro-pore and micro-pore networks and are, respectively, denoted by subscripts and . These pore-networks are connected with the possibility of mass exchange between them. The pressure field and the discharge (or Darcy) velocity in the macro-pore network are, respectively, denoted by and , and the corresponding ones in the micro-pore network are denoted by and . The governing equations under the double porosity/permeability model take the following form:
where is the specific body force. The true density and coefficient of viscosity of the fluid are, respectively, denoted by and .
denotes the permeability tensor for macro-pore () and micro-pore () networks. denotes the part of the boundary on which the normal component of the velocity is prescribed in the macro-pore () and micro-pore () networks. Similarly, is that part of the boundary on which the pressure is prescribed in the macro-pore () and micro-pore () networks. and denote the prescribed pressures on and , respectively. and denote the prescribed normal components of the velocities on and , respectively. is the rate of volume exchange of the fluid between the two pore-networks per unit volume of the porous medium, and we model it as follows [Barenblatt et al., 1960]:
where is a dimensionless characteristic of the porous medium. In the rest of the paper, as is commonly done in the literature, will be simply referred to as the mass transfer. For mathematical well-posedness, we assume that
3. A Stabilized Mixed Weak Formulation
In this section, we present the proposed stabilized mixed formulation for the double porosity/permeability model. A derivation of the proposed formulation is provided in Appendix A. The proposed formulation is built upon the stabilization ideas put forth in a pioneering paper by [Masud and Hughes, 2002]. The proposed formulation for double porosity/permeability model can be obtained by adding a stabilization term, similar to the one proposed by [Masud and Hughes, 2002] for the case of single-pore network Darcy equations, to each pore-network. The stabilization terms are based on the residual of the balance of linear momentum in each pore-network. The stability can be achieved without adding residual-based stabilization terms related to the mass balance equations for any of the pore-networks. We also present an extension of the proposed formulation for enforcing the velocity boundary conditions weakly, which will be convenient for problems involving curved boundaries. This extension is achieved by employing a procedure similar to the one proposed by [Nitsche, 1971].
We define the relevant function spaces, which will be used in the rest of this paper. We denote the set of all square-integrable functions on by . For mathematical well-posedness, we assume that
where is a non-integer Sobolev space [Adams and Fournier, 2003]. The function spaces for the velocity and pressures fields are defined as follows:
where is a standard Sobolev space, and is the dual space corresponding to [Adams and Fournier, 2003]. The standard inner-product over a set is denoted as
For convenience, the subscript will be dropped if
. Moreover, the action of a linear functional on a vector from its associated vector space is denoted by.
A few remarks are needed regarding the following condition on the pressures in the function spaces and :
This condition of vanishing mean pressure in one of pore-networks is a mathematically elegant way of fixing the datum for the pressure. Without fixing the datum for the pressure (which will be the case when only the velocity boundary conditions are prescribed on the entire boundary), one can find the pressures only up to an arbitrary constant, which will be the case even under Darcy equations [Nakshatrala et al., 2006]. Herein, we introduced the vanishing mean pressure condition into the function spaces to ensure uniqueness of the solutions, which will be established later in this paper. However, it should be emphasized that vanishing mean pressure in one of the pore-networks is not necessary for all the problems under the double porosity/permeability model. One can fix the datum for the pressure under the double porosity/permeability model by prescribing the pressure in at least one of the pore-networks on a portion of the boundary, which is a set of non-zero measure. To put it differently, for problems with pressure boundary conditions, the datum for the pressure is automatically fixed through the prescribed boundary condition, and hence, for those problems, one does not include the zero mean pressure condition in the function spaces and . For example, see the problem in subsection 5.1, which deals with prescribed pressure boundary conditions.
The classical mixed formulation, which is based on the Galerkin formalism, reads as follows: Find and such that we have
where the bilinear form and the linear functional are, respectively, defined as follows:
In a subsequent section, we will show that the equal-order interpolation for all the variables, which is computationally the most convenient, is not stable under the classical mixed formulation. Of course, one could use divergence-free elements (e.g., Raviart-Thomas spaces [Raviart and Thomas, 1977]) but they need special data structures and computer implementations. We, therefore, present a stabilized mixed formulation, which is stable under the equal-order interpolation for all the field variables.
[breakable] The proposed stabilized mixed formulation reads as follows: Find and such that we have
where the bilinear form and the linear functional are, respectively, defined as follows:
In subsequent sections, we show that the proposed stabilized mixed formulation is consistent, stable and accurate.
3.1. Weak enforcement of velocity boundary conditions
In the previous derivations made earlier in this section, the pressure boundary conditions (i.e., equations (2.1g) and (2.1h)) are enforced weakly under the proposed stabilized mixed formulation and the classical mixed formulation. However, the velocity boundary conditions in which the normal components of the velocities are prescribed (i.e., equations (2.1e) and (2.1f)) are enforced strongly. For domains with curved boundaries, which are commonly encountered in subsurface modeling, it is desirable to even prescribe the velocity boundary conditions weakly. We, therefore, provide a possible extension of the proposed stabilized mixed formulation for weak enforcement of the velocity boundary conditions. To this end, we follow the approach proposed by [Nitsche, 1971]. The Nitsche’s method is a powerful tool for weakly enforcing Dirichlet boundary conditions without the use of Lagrange multipliers, and has been utilized by several works such as [Bazilevs and Hughes, 2007; Embar et al., 2010; Annavarapu et al., 2014; Schillinger et al., 2016]. The Nitsche’s method is sometimes referred to as a variationally consistent penalty method to enforce Dirichlet boundary conditions [Hansbo, 2005]. We extend the Nitsche’s method to the proposed four-field stabilized formulation to enforce the prescribed normal components of the velocities in the macro- and micro-pore networks.
The stabilized mixed formulation that enforces the velocity boundary conditions weakly can be obtained as follows: Find and such that we have
where the bilinear form and the linear functional are, respectively, defined as follows:
where is the mesh size and is the penalty parameter. In this paper, we have taken to be the maximum edge length in the mesh, and have taken the penalty parameter to be 10. In the above statement of the weak formulation, since the velocity boundary conditions are enforced weakly, the appropriate function space for the velocities and the associated weighting functions will be , which can be mathematically defined as follows:
The function spaces for the pressures and their weighting functions, however, remain same as before (i.e., the space).
4. A Theoretical Analysis of the Proposed Mixed Formulation
In this section, we present a systematic mathematical analysis (i.e., existence, uniqueness and well-posedness) and error analysis (i.e., consistency, stability, order of convergence) of the proposed stabilized mixed formulation. For convenience, we define the following product spaces:
We group the field variables as follows:
Then, the proposed mixed formulation in equation (3.7) can be compactly written as: Find such that we have
We shall establish the stability of the formulation under the following norm:
where denotes the norm corresponding to the standard inner-product. We need to first show that is in fact a norm on and . To this end, the following lemma will be used.
(A property of semi-norms) If and are semi-norms, then is also a semi-norm.
The homogeneity of directly stems from the homogeneity of the semi-norms and . To wit,
The non-negativity of is straightforward; that is, . The triangle inequality for the semi-norms and implies that
These inequalities imply that
We have employed the AM-GM inequality in obtaining equation (4.7), which further implies that
This establishes the triangle inequality for . The homogeneity, non-negativity and triangle inequality imply that is a semi-norm. ∎
(Stability norm) is a norm on and .
We first note that and are symmetric and positive definite tensors. The square root of a symmetric and positive definite tensor exists, and is itself a symmetric and positive definite tensor [Gurtin, 1981]. This implies that the following individual terms form semi-norms on and :
Then, Lemma 4.1 implies that is a semi-norm. It is easy to show that implies that
where is a constant. Noting that and utilizing the following condition in the definition of :
we conclude that . With this, we have established that = 0 implies that . Hence, is a norm. ∎
(Uniqueness of weak solutions) The weak solution under the proposed mixed formulation is unique.
On the contrary, assume that and are both (weak) solutions of the weak formulation. This implies that
By subtracting the above two equations and noting the linearity in the second slot, we obtain
Since , we can choose . This particular choice implies that
Using Proposition 4.1 (which establishes that is a norm on ) we conclude that . ∎
(Boundedness) The bilinear form is bounded. That is,
where is a constant.
A direct application of the triangle inequality of the absolute value on real numbers implies that
Cauchy-Schwartz inequality on inner-product implies that
By applying Cauchy-Schwartz inequality on -tuple real numbers (i.e., on Euclidean spaces) we obtained the following:
That is, we have established that
which completes the proof. ∎
(Coercivity) The bilinear form is coercive. That is, the bilinear form is bounded below.
The coercivity of the bilinear form can be established from the definition of and Proposition 4.1 (i.e., is a norm on ) as
Given the coercivity and boundedness of the bilinear form and the continuity of the linear functional, one can conclude that the proposed mixed weak formulation is well-posed by invoking the Lax-Milgram theorem [Brenner and Scott, 1994].
4.1. Convergence and error analysis of the finite element formulation
We decompose the computational domain into “” subdomains (which will be the elements in the context of the finite element method) such that
where a superposed bar indicates the set closure. We denote the finite element solution by . That is,
If we denote the set of all polynomials up to and including -th order over a set by , and the set of all continuous functions defined on (which is the set closure of ) by , then the following finite-dimensional spaces can be defined:
We define the corresponding product spaces as follows:
It is important to note that and are closed linear subspaces of and , respectively. The finite element formulation corresponding to the proposed stabilized mixed formulation reads: Find such that we have
In a given coordinate system, we denote . For a given multi-index (i.e., tuple) of non-negative integers, , with order , the corresponding partial derivative of a scalar field, , can be written as follows:
Using the above notation, the Sobolev semi-norm, , for scalar and vector fields can be compactly written as follows:
where denotes the summation over all the possible tuples of non-negative integers with order , and denotes the characteristic length of the domain.