Abstract

In this paper, we study a finite element computational model for solving the interaction between a fluid and a poroelastic structure that couples the Stokes equations with the Biot system. Equilibrium and kinematic conditions are imposed on the interface. A mixed Darcy formulation is employed, resulting in continuity of flux condition of essential type. A Lagrange multiplier method is used to impose weakly this condition. With the obtained finite element solutions, the error estimators are performed for the fully discrete formulations.

1. Introduction

The paper presents a computation friendly finite element formulation of coupled fluid motions at a free-poroelastic interface, where Navier–Stokes equations are employed for free fluid and Darcy’s law is used for the permeable material. A reliable and efficient a posteriori error estimator is analysed. This is a challenging multiphysics problem that has a wide range of applications, including processes arising in gas and oil extraction from naturally or hydraulically fractured reservoirs, designing industrial filters, and blood-vessel interaction. In these applications, it is important to model properly the interaction between the free fluid with the fluid within the porous medium and to take into account the effect of the deformation of the medium. For example, geomechanical effects play an important role in hydraulic fracturing, as well as in modeling phenomena such as subsidence and compaction.

The free fluid region can be modeled by the Stokes or Navier–Stokes equations, while the flow through the deformable porous medium is modeled by the quasi-static Biot systems of poroelasticity [1]. The two regions are coupled via dynamic and kinematic interface conditions, including balance of forces, continuity of normal velocity, and a no slip or slip with friction tangential velocity condition [210]. These multiphysics models exhibit features of coupled Stokes–Darcy flows and fluid-structure interaction (FSI) [1115].

The well-posedness of the mathematical model based on the Stokes–Biot system for the coupling between a fluid and a poroelastic structure is studied in [16]. A numerical study of the problem, using a Navier–Stokes equations for the fluid, is presented in [11, 17], utilizing a variational multiscale approach to stabilize the finite element spaces. The problem is solved using both a monolithic and a partitioned approach, with the latter requiring subiterations between the two problems.

Finite element analysis of an arbitrary Lagrangian–Eulerian method for Stokes/parabolic moving interface problem with jump coefficients has been studied in [18]. The authors in [19] studied a numerical solution of the coupled system of the time-dependent Stokes and fully dynamic Biot equations. They established stability of the scheme and derived error estimates for the fully discrete coupled scheme. Numerical errors and convergence rates for smooth problems as well as tests on realistic material parameters have been presented. In [20], Jing Wen and Yinnian He considered a strongly conservative discretization for the rearranged Stokes–Biot model based on the interior penalty discontinuous Galerkin method and mixed finite element method. The existence and uniqueness of solution of the numerical scheme have been presented. Then, the analysis of stability and priori error estimates has been derived. The numerical examples under uniform meshes which well validate the analysis of convergence and the strong mass conservation are presented. A staggered finite element procedure for the coupled Stokes–Biot system with fluid entry resistance has been studied by Bergkamp et al. in [21] while Ambartsumyan et al. in [22] studied flow and transport in fractured poroelastic media using Stokes flow in the fractures and the Biot model in the porous media. In Section 6 in [23], fully discrete continuous approximation has been proposed for the weak coupled mixed formulation. For the discretization of the fluid velocity and pressure, the authors have used the finite elements which include the MINI-elements, the Taylor–Hood elements, and the conforming Crouzeix–Raviart elements. For the discretization of the porous medium problem, they choose the spaces that include Raviart–Thomas and Brezzi–Douglas–Marini elements. An a priori error analysis is performed with some numerical tests confirming the convergence rates.

A posteriori error estimators are computable quantities, expressed in terms of the discrete solution and of the data that measure the actual discrete errors without the knowledge of the exact solution. They are essential to design adaptive mesh refinement algorithms which aqui-distribute the computational effort and optimize the approximation efficiency. Since the pioneering work of Babuška and Rheinboldt [2427], adaptive finite element methods based on a posteriori error estimates have been extensively investigated.

In [28], we have proposed a family error indicator for semidiscrete approximation for the Stokes–Biot system. To the best of our knowledge, there is no a posteriori error estimation for the Stokes–Biot fluid-poroelastic structure interaction model for fully discrete finite element methods. Here, we develop such a posteriori error analysis for the fully discrete conforming finite element methods. We have got a new family of a local indicator error (see equation (59) in Definition 3) and global (equation (60)). The difference between this paper and our previous work [28] is that our discretization is fully discrete formulation. As an advantage, the error indicators in this work are more accessible to computation.

The schedule of the paper is as follows. Section 2 is devoted to notations and basic results that are used throughout the document. Our main results regarding a posteriori error analysis are stated in Section 3. We prove that our indicator errors are efficient, reliable, and then optimal. The global inf-sup condition is the main tool yielding the reliability. In turn, the local efficiency result is derived using the technique of bubble function introduced by R. Verfürth [29] and used in similar context by C. Carstensen [30]. Finally, this paper is summarized with further works in Section 4.

2. Preliminaries and Notations

2.1. Stokes–Biot Model Problem

We consider a multiphysics model problem for free fluid’s interaction with a flow in a deformable porous media, where the simulation domain , , is a union of nonoverlapping regions and . Here, is a free fluid region with flow governed by the Stokes equations and is a poroelastic material governed by the Biot system. For simplicity of notation, we assume that each region is connected. The extension to nonconnected regions is straightforward. Let (see Figure 1).

Let be the velocity-pressure pair in , , and let be the displacement in . Let be the fluid viscosity, let be the body force terms, and let be external source or sink terms. Let and denote, respectively, the deformation rate tensor and the stress tensor:

In the free fluid region , satisfies the Stokes equations:where is the final time. Let and be the elastic and poroelastic stress tensors, respectively:where and are the Lamé parameters and is the Biot–Willis constant. The poroelasticity region is governed by the quasi-static Biot system [23]:where is a storage coefficient and is the symmetric and uniformly positive definite rock permeability tensor, satisfying, for some constants .

Following [1], the interface conditions on the fluid-poroelasticity interface are mass conservation, balance of stresses, and the Beavers–Joseph–Saffman (BJS) condition [31] modeling slip with friction:where and are the outward unit normal vectors to and , respectively, , , is an orthogonal system of unit tangent vectors on , , and is an experimentally determined friction coefficient.

Here, equations (8) and (9) represent mass conservation, equation (10) is the balance of normal forces, and equation (11) is the Beavers–Joseph–Saffman conditions. We note that continuity of flux constraints the normal velocity of the solid skeleton, while the condition accounts for its tangential velocity.

The above system of equations needs to be complemented by a set of boundary and initial conditions. Let and . Let . We assume that for simplicity homogeneous boundary conditions,

To avoid the issue with restricting the mean value of the pressure, we assume that . We also assume that is not adjacent to the interface , i.e., . Nonhomogeneous displacement and velocity conditions can be handled in a standard way by adding suitable extensions of the boundary data. The pressure boundary condition is natural in the mixed Darcy formulation, so nonhomogeneous pressure data would lead to an additional boundary term. We further say that the initial conditions are as follows:

Equations (2)–(10) consist of the model of the coupled Stokes and Biot flows problem that we will study below.

2.2. Weak Formulation

In this part, we first introduce some Sobolev spaces [32] and norms. If is a bounded domain of and is a nonnegative integer, the Sobolev space is defined in the usual way with the usual norm and seminorm . In particular, , and we write for . Similarly, we denote by the or inner product. For shortness, if is equal to , we will drop the index , while for any , , , and , for . The space denotes the closure of in . Let be the space of vector valued functions with components in . The norm and the seminorm on are given by

For a connected open subset of the boundary , we write for the inner product (or duality pairing); that is, for scalar valued functions and , one defines

In the following, we derive a Lagrange multiplier type weak formulation of the system, which will be the basis for our finite element approximation. Letwhere is the space of -vectors with divergence in with a norm

We define the global velocity and pressure spaces aswith norms

The weak formulation is obtained by multiplying the equations in each region by suitable test functions, integrating by parts the second-order terms in space, and utilizing the interface and boundary conditions.

Let be the bilinear forms related to Stokes, Darcy, and the elasticity operator, respectively. Let

Integration by parts in (2) and the two equations in (5) leads to the interface term as follows:

Using the first condition for balance of normal stress in (9), we setwhich will be used as a Lagrange multiplier to impose the mass conservation interface condition (8). Utilizing the BJS condition (10) and the second condition for balance of stresses in (9), we obtainwhere

For the well-posedness of , we require that . According to the normal trace theorem, since , then . Furthermore, since on and , then , (see, e.g., [33]). Therefore, we take .

The Lagrange multiplier variational formulation is for , find , , , , , and , such that and , and for all , , , , , and ,where we used the notation .

The assumptions on the fluid viscosity and the material coefficients , , and imply that the bilinear forms , , and are coercive and continuous in the appropriate norms. In particular, there exist positive constants , , , , , and such thatwhere (29), (30), (33), and (34) hold true thanks to Poincaré inequality and (33) and (34) also relies on Korn’s inequality.

In summary, from Corollary 3.1 in [23] (p. 7), the following result holds.

Theorem 1. There exists a unique solution to problems (26)–(28).

2.3. Finite Element Discretization

Let and be shape-regular and quasi-uniform partition of and , respectively, both consisting of affine elements with maximal element diameter . The two partitions may be nonmatching at the interface . For the discretization of the fluid velocity and pressure, we choose finite element spaces and , which are assumed to be inf-sup stable. Examples of such spaces include the mini-elements, the Taylor–Hood elements, and the conforming Crouzeix–Raviart elements. For the discretization of the porous medium problem, we choose and to be any of well-known inf-sup stable mixed finite element spaces, such as the Raviart–Thomas or the Brezzi–Douglas–Marini spaces. The global spaces are as follows:

We employ a conforming Lagrangian finite element space to approximate the structure displacement. Note that the finite element spaces , , and satisfy the prescribed homogeneous boundary conditions on the external boundaries. For the discrete Lagrange multiplier space, we take

The semidiscrete continuous-in-time problem reads that given and , for , find , , , , , and such that for all , , , , , and ,

We will take and to be suitable projections of the initial data and .

We introduce the errors for all variables as follows:

The following results hold cf. [23].

Theorem 2 (a priori error estimation). There exists a unique solution in of the discrete problems (37)–(39), and if the solution of the continuous problems (26)–(28) is smooth enough, then we have

Remark 1. For , we can subtract (37)–(39) to (26)–(28) to obtain the Galerkin orthogonality relation for all :

2.4. Fully Discrete Formulation

For the time discretization, we employ the backward Euler method. Let be the time step, , and let , . Let be the first-order (backward) discrete time derivation, where . Then, the fully discrete model reads that given and , find , , , , , and , , such that for all , , , , , and ,

We introduce the discrete-in-time norms as follows:and the errors for all variables-in-time as follows:

Remark 2. For and , we can subtract (37)–(39) to (43) and (44) to obtain the Galerkin orthogonality relation:

3. Error Estimation

In order to solve the Stokes–Biot model problem by efficient adaptive finite element methods, reliable and efficient a posteriori error analysis is important to provide appropriated indicators. In this section, we first define the local and global indicators (Section 3.1) and then the lower and upper error bounds are derived (Sections 3.4 and 3.5).

3.1. Residual Error Estimators

The general philosophy of residual error estimators is to estimate an appropriate norm of the correct residual by terms that can be evaluated easier and that involve the data at hand. To this end, define the exact element residuals.

Definition 1. (exact element residuals). Let and be an arbitrary finite element function. The exact element residuals over a triangle or tetrahedra and over are defined byAs it is common, these exact residuals are replaced by some finite-dimensional approximation called approximate element residual , , , , . This approximation is here achieved by projecting and on the space of piecewise constant functions in and piecewise functions in ; more precisely for all , we takeWhile for all , we take and as the unique element of , respectively, such thatrespectively,Thereby, we define the approximate element residuals.

Definition 2. (approximate element residuals). Let and be an arbitrary finite element function. Then, the approximate element residuals are defined byNext, introduce the gradient jump in normal direction bywhere is the identity matrix of .

Definition 3. (residual error estimators). Let and be the finite element solution of the fully discrete problems (43)–(45) in .
We definewhereandThen, the residual error estimator is locally defined byThe global residual error estimator is given byFurthermore, denote the local approximation terms bywhereWe set

Remark 3. The residual character of each term on the right-hand sides of (55)–(58) is quite clear since if would be the exact solution of (2)–(10), then they would vanish.

3.2. Analytical Tools
3.2.1. Inverse Inequalities

In order to derive the lower error bounds, we proceed similarly as in [30, 34] (see also [35]), by applying inverse inequalities and the localization technique based on simplex-bubble and face-bubble functions. To this end, we recall some notation and introduce further preliminary results. Given and , we let and be the usual simplex-bubble and face-bubble functions (see (1.5) and (1.6) in [29]), respectively. In particular, satisfies , , , and . Similarly, , , and . We also recall from [36] that, given , there exists an extension operator that satisfies and . A corresponding vectorial version of , that is, the component-wise application of , is denoted by . Additional properties of , , and are collected in the following lemma (see [36]).

Lemma 1. Given , there exist positive constants depending only on and shape-regularity of the triangulations (minimum angle condition), such that for each simplex and , there hold

Lemma 2 (continuous trace inequality). There exists a positive constant depending only on such that

3.2.2. Clément Interpolation Operator

In order to derive the upper error bounds, we introduce the Clément interpolation operator that approximates optimally nonsmooth functions by continuous piecewise linear functions:

In addition, we will make use of a vector valued version of , that is, , which is defined component-wise by . The following lemma establishes the local approximation properties of (and hence of ), for a proof (see Section 3 in [37]).

Lemma 3. There exist constants , independent of , such that for all , there holdwhere and .

3.2.3. Helmholtz Decomposition

Lemma 4 (see [7]). There exists such that every can be decomposed as , where , , , and

3.3. Error Equation

We set and .

For , let and . We define now the operator by

Then, the continuous problems (26)–(28) are equivalent to the following: for , find such that we have

We define the discrete version by the same way: find such that

Since for all and , , then from (70), we obtain

Now, we set

Then, from (72), we can writewhere

3.4. Reliability of the A Posteriori Error Estimator

The first main result is given by the following theorem.

Theorem 3 (upper error bound). Let such that be the exact solution and with be the finite element solution. Then, there exist a positive constant such that the error is bounded globally from above by

Proof. Let . We consider the residual equation (73), and we set . We take with and . As , then by Lemma 4, admits the decomposition where and with and . We consider with and . Thus, . Now, we recallTherefore, integrating by parts element by element, we may writewhereThus, from equations (74) and (78), we obtainCoercivity of operator leads to inf-sup condition as follows:By consequently, identities (73) and (80), inf-sup condition of operator (81), Cauchy–Schwarz inequality, estimation of Lemma 4, and the approximation properties of Lemma 3 imply the required estimate and finish the proof. □

3.5. Efficiency of the A Posteriori Error Estimator

For , we setwhere

The family estimator is consider efficient if it satisfies the following theorem.

Theorem 4 (lower error bound). Let such that be the exact solution and with be the finite element solution. Then, there exist a positive constant such that the error is bounded locally from below for all bywhere is a finite union of neighbouring elements of and is the product norm.

Proof. The lower bound is proved using the standard elementwise integration by parts, namely, error equation of Section 3.3 (i.e., identity (73) and equation (80)) and some inverse estimates of Lemma 1 (cf. [38] for details).

4. Discussion

In this paper, we have discussed a posteriori error estimates for a finite element approximation of the Stokes–Biot system where homogeneous boundary conditions are employed. The approach utilizes a fully discrete conforming finite element method. A residual type a posteriori error estimator is provided, that is, both reliable and efficient.

In a future paper, we will study the influence of nonhomogeneous boundary conditions on the a posteriori error indicators presented in this work. Further, it is well known that an internal layer appears at the interface as the permeability tensor degenerates; in that case, anisotropic meshes have to be used in this layer (see, for instance, [10]). Hence, we intend to extend our results to such anisotropic meshes.

5. Nomenclatures

(i): bounded domain(ii): the poroelastic medium domain(iii)(iv)(v), (vi)(respectively, ): the unit outward normal vector along (respectively, )(vii): the fluid velocity in (viii): the fluid pressure in (ix): the fluid velocities in (x): the fluid pressure in (xi)In , the of a scalar function is given as usual by (xii)In , the of a vector function is given as usual by , namely, (xiii): the space of polynomials of total degree not larger than (xiv): triangulation of (xv): the corresponding induced triangulation of , (xvi)For any , is the diameter of and is the diameter of the largest ball inscribed into (xvii)and (xviii): the set of all the edges or faces of the triangulation(xix): the set of all the edges () or faces () of a element (xx)(xxi): the set of all the vertices of a element (xxii)(xxiii)For , (xxiv)For , we associate a unit vector such that is orthogonal to and equals to the unit exterior normal vector to , (xxv)For , is the jump across in the direction of (xxvi)In order to avoid excessive use of constants, the abbreviations and stand for and , respectively, with positive constants independent of , or (xxvii)(xxviii)(xxix)(xxx)(xxxi)(xxxii)(xxxiii): product norm on (xxxiv): local product norm on .

Data Availability

There are no data underlying the findings in this paper to be shared.

Disclosure

The results presented in this paper constitute a continuation of our work posted on arxiv at the following link: https://arxiv.org/abs/2004.10676.

Conflicts of Interest

The authors declare that they have no conflicts of interest.