A multipoint stress-flux mixed finite element method for the Stokes-Biot model

The interaction of a free fluid with a deformable porous medium, referred to as fluid-poroelastic structure interaction (FPSI), is a challenging multiphysics problem. It has applications to predicting and controlling processes arising in gas and oil extraction from naturally or hydraulically fractured reservoirs, modeling arterial flows, and designing industrial filters, to name a few. For this physical phenomenon, 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 Biot system of poroelasticity. In the latter, the volumetric deformation of the elastic porous matrix is complemented with the Darcy equation that describes the average velocity of the fluid in the pores. 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. The model exhibits features of both coupled Stokes-Darcy flows and fluid-structure interaction (FSI).

To the authors’ knowledge, one of the first works in analyzing the Stokes-Biot coupled problem is [53], where well-posedness for the fully dynamic problem is established by developing an appropriate variational formulation and using semigroup methods. One of the first numerical studies is presented in [12], where monolithic and iterative partitioned methods are developed for the solution of the coupled system. A non-iterative operator splitting scheme with a non-mixed Darcy formulation is developed in [21]. Finite element methods for mixed Darcy formulations, where the continuity of normal flux condition becomes essential, are considered in [20] using Nitsche’s coupling and in [8] using a pressure Lagrange multiplier. More recently, a nonlinear quasi-static Stokes–Biot model for non-Newtonian fluids is studied in [3]. The authors establish well-posedness of the weak formulation in Banach space setting, along with stability and convergence of the finite element approximation. In [25], the fully dynamic coupled Navier-Stokes/Biot system with a pressure-based Darcy formulation is analyzed. Additional works include optimization-based decoupling method [24], a second order in time split scheme [44], various discretization methods [56, 13, 23], dimensionally reduced model for flow through fractures [22], and coupling with transport [4]. All of the above mentioned works are based on displacement formulations for the elasticity equation. In a recent work [47], the first mathematical and numerical study of a stress-displacement mixed elasticity formulation for the Stokes-Biot model is presented.

The goal of the present paper is to develop a new fully mixed formulation of the quasi-static Stokes-Biot model, which is based on dual mixed formulations for all three components - Darcy, elasticity, and Stokes. In particular, we use a velocity-pressure Darcy formulation, a weakly symmetric stress-displacement-rotation elasticity formulation, and a weakly symmetric stress-velocity-vorticity Stokes formulation. This formulation exhibits multiple advantages, including local conservation of mass for the Darcy fluid, local poroelastic and Stokes momentum conservation, and accurate approximations with continuous normal components across element edges or faces for the Darcy velocity, the poroelastic stress, and the free fluid stress. In addition, dual mixed formulations are known for their locking-free properties and robustness with respect to the physical parameters, including the regimes of almost incompressible materials, low poroelastic storativity, and low permeability [45, 60].

Our five-field dual mixed Biot formulation is based on the model developed in [45] and studied further in [7]. It is also considered in [47] for the Stokes-Biot problem. Our analysis also extends to the strongly symmetric mixed four-field Biot formulation developed in [59]. Our three-field dual mixed Stokes formulation is based on the models developed in [35, 34]. In particular, we introduce the stress tensor and subsequently eliminate the pressure unknown, by utilizing the deviatoric stress. In order to impose the symmetry of the Stokes stress and poroelastic stress tensors, the vorticity and structure rotation, respectively, are introduced as additional unknowns. The transmission conditions consisting of mass conservation, conservation of momentum, and the Beavers–Joseph–Saffman slip with friction condition are imposed weakly via the incorporation of additional Lagrange multipliers: the traces of the fluid velocity, structure velocity and the poroelastic media pressure on the interface. The resulting variational system of equations is then ordered so that it shows a twofold saddle point structure. The well-posedness and uniqueness of both the continuous and semidiscrete continuous-in-time formulations are proved by employing some classical results for parabolic problems [52, 54] and monotone operators, and an abstract theory for twofold saddle point problems [33, 1]. In the discrete problem, for the three components of the model we consider suitable stable mixed finite element spaces on non-matching grids across the interface, coupled through either conforming or non-conforming Lagrange multiplier discretizations. We develop stability and error analysis, establishing rates of convergence to the true solution. The estimates we establish are uniform in the limit of the storativity coefficient going to zero.

Another main contribution of this paper is the development of a new mixed finite element method for the Stokes-Biot model that can be reduced to a positive definite cell-centered pressure-velocities-traces system. We recall the multipoint flux mixed finite element (MFMFE) method for Darcy flow developed in [40, 57, 58, 19], where the lowest order Brezzi-Douglas-Marini $\mathbb{BDM}_{1}$ velocity spaces [17, 48, 18] and piecewise constant pressure are utilized. An alternative formulation based on a broken Raviart-Thomas velocity space is developed in [43]. The use of the vertex quadrature rule for the velocity bilinear form localizes the interaction between velocity degrees of freedom around mesh vertices and leads to a block-diagonal mass matrix. Consequently, the velocity can be locally eliminated, resulting in a cell-centered pressure system. In turn, the multipoint stress mixed finite element (MSMFE) method for elasticity is developed in [5, 6]. It utilizes stable weakly symmetric elasticity finite element triples with $\mathbb{BDM}_{1}$ stress spaces [16, 30, 46, 6, 10, 11]. Similarly to the MFMFE method, an application of the vertex quadrature rule for the stress and rotation bilinear forms allows for local stress and rotation elimination, resulting in a cell-centered displacement system. We also refer the reader to the related finite volume multipoint stress approximation (MPSA) method for elasticity [49, 50, 41]. Recently, combining the MSMFE and MFMFE methods, a multipoint stress-flux mixed finite element (MSFMFE) method for the Biot poroelasticity model is developed in [7]. There, the dual mixed finite element system is reduced to a cell-centered displacement-pressure system. The reduced system is comparable in cost to the finite volume method developed in [51].

In this paper we note for the first time that the MSMFE method for elasticity can be applied to the weakly symmetric stress-velocity-vorticity Stokes formulation from [35, 34] when $\mathbb{BDM}_{1}$-based stable finite element triples are utilized. With the application of the vertex quadrature rule, the fluid stress and vorticity can be locally eliminated, resulting in a positive definite cell-centered velocity system. To the best of our knowledge, this is the first such scheme for Stokes in the literature.

Finally, we combine the MFMFE method for Darcy with the MSMFE methods for elasticity and Stokes to develop a multipoint stress-flux mixed finite element for the Stokes-Biot system. We analyze the stability and convergence of the semidiscrete formulation. We further consider the fully discrete system with backward Euler time discretization and show that the algebraic system on each time step can be reduced to a positive definite cell-centered pressure-velocities-traces system.

The rest of this work is organized as follows. The remainder of this section describes standard notation and functional spaces to be employed throughout the paper. In Section 2 we introduce the model problem and in Section 3 we derive a fully-mixed variational formulation, which is written as a degenerate evolution problem with a twofold saddle point structure. Next, existence, uniqueness and stability of the solution of the weak formulation are obtained in Section 4. The corresponding semidiscrete continuous-in-time approximation is introduced and analyzed in Section 5, where the discrete analogue of the theory used in the continuous case is employed to prove its well-posedness. Error estimates and rates of convergence are also derived there. In Section 6, the multipoint stress-flux mixed finite element method is presented and the corresponding rates of convergence are provided, along with the analysis of the reduced cell-centered system. Finally, numerical experiments illustrating the accuracy of our mixed finite element method and its applications to coupling surface and subsurface flows and flow through poroelastic medium with a cavity are reported in Section 7.

We end this section by introducing some definitions and fixing some notations. Let $\mathcal{O}\subset\mathrm{R}^{n}$, $n\in\{2,3\}$, denote a domain with Lipschitz boundary. For $\mathrm{s}\geq 0$ and $\mathrm{p}\in[1,+\infty]$, we denote by $\mathrm{L}^{\mathrm{p}}(\mathcal{O})$ and $\mathrm{W}^{\mathrm{s},\mathrm{p}}(\mathcal{O})$ the usual Lebesgue and Sobolev spaces endowed with the norms $\|\cdot\|_{\mathrm{L}^{\mathrm{p}}(\mathcal{O})}$ and $\|\cdot\|_{\mathrm{W}^{\mathrm{s},\mathrm{p}}(\mathcal{O})}$, respectively. Note that $\mathrm{W}^{0,\mathrm{p}}(\mathcal{O})=\mathrm{L}^{\mathrm{p}}(\mathcal{O})$. If $\mathrm{p}=2$ we write $\mathrm{H}^{\mathrm{s}}(\mathcal{O})$ in place of $\mathrm{W}^{\mathrm{s},2}(\mathcal{O})$, and denote the corresponding norm by $\|\cdot\|_{\mathrm{H}^{\mathrm{s}}(\mathcal{O})}$. Similar notation is used for a section $\Gamma$ of the boundary of $\mathcal{O}$. By $\mathbf{M}$ and $\mathbb{M}$ we will denote the corresponding vectorial and tensorial counterparts of a generic scalar functional space $\mathrm{M}$. The $\mathrm{L}^{2}(\mathcal{O})$ inner product for scalar, vector, or tensor valued functions is denoted by $(\cdot,\cdot)_{\mathcal{O}}$. The $\mathrm{L}^{2}(\Gamma)$ inner product or duality pairing is denoted by $\left<\cdot,\cdot\right>_{\Gamma}$. For any vector field ${\mathbf{v}}=(v_{i})_{i=1,n}$, we set the gradient and divergence operators, as

For any tensor fields ${\boldsymbol{\tau}}:=(\tau_{ij})_{i,j=1,n}$ and ${\boldsymbol{\zeta}}:=(\zeta_{ij})_{i,j=1,n}$, we let $\mathbf{div}({\boldsymbol{\tau}})$ be the divergence operator $\mathrm{div}$ acting along the rows of ${\boldsymbol{\tau}}$, and define the transpose, the trace, the tensor inner product, and the deviatoric tensor, respectively, as

where $\mathbf{I}$ is the identity matrix in $\mathrm{R}^{n\times n}$. In addition, we recall the Hilbert space

equipped with the norm $\|{\mathbf{v}}\|^{2}_{\mathbf{H}(\mathrm{div};\mathcal{O})}:=\|{\mathbf{v}}\|^{2}_{\mathbf{L}^{2}(\mathcal{O})}+\|\mathrm{div}({\mathbf{v}})\|^{2}_{\mathrm{L}^{2}(\mathcal{O})}$. The space of matrix valued functions whose rows belong to $\mathbf{H}(\mathrm{div};\mathcal{O})$ will be denoted by $\mathbb{H}(\mathbf{div};\mathcal{O})$ and endowed with the norm $\|{\boldsymbol{\tau}}\|^{2}_{\mathbb{H}(\mathbf{div};\mathcal{O})}:=\|{\boldsymbol{\tau}}\|^{2}_{\mathbb{L}^{2}(\mathcal{O})}+\|\mathbf{div}({\boldsymbol{\tau}})\|^{2}_{\mathbf{L}^{2}(\mathcal{O})}$. Finally, given a separable Banach space $\mathrm{V}$ endowed with the norm $\|\cdot\|_{\mathrm{V}}$, we let $\mathrm{L}^{\mathrm{p}}(0,T;\mathrm{V})$ be the space of classes of functions $f:(0,T)\to\mathrm{V}$ that are Bochner measurable and such that $\|f\|_{\mathrm{L}^{\mathrm{p}}(0,T;\mathrm{V})}<\infty$, with

Let $\Omega\subset\mathrm{R}^{n}$, $n\in\{2,3\}$, be a Lipschitz domain, which is subdivided into two non-overlapping and possibly non-connected regions: fluid region $\Omega_{f}$ and poroelastic region $\Omega_{p}$. Let $\Gamma_{fp}=\partial\Omega_{f}\cap\partial\Omega_{p}$ denote the (nonempty) interface between these regions and let $\Gamma_{f}=\partial\Omega_{f}\setminus\Gamma_{fp}$ and $\Gamma_{p}=\partial\Omega_{p}\setminus\Gamma_{fp}$ denote the external parts on the boundary $\partial\Omega$. We denote by ${\mathbf{n}}_{f}$ and ${\mathbf{n}}_{p}$ the unit normal vectors that point outward from $\partial\Omega_{f}$ and $\partial\Omega_{p}$, respectively, noting that ${\mathbf{n}}_{f}=-{\mathbf{n}}_{p}$ on $\Gamma_{fp}$. Let $(\mathbf{u}_{\star},p_{\star})$ be the velocity-pressure pair in $\Omega_{\star}$ with $\star\in\{f,p\}$, and let ${\boldsymbol{\eta}}_{p}$ be the displacement in $\Omega_{p}$. Let $\mu>0$ be the fluid viscosity, let $\mathbf{f}_{\star}$ be the body force terms, and let $q_{\star}$ be external source or sink terms.

We assume that the flow in $\Omega_{f}$ is governed by the Stokes equations, which are written in the following stress-velocity-pressure formulation:

where ${\boldsymbol{\sigma}}_{f}$ is the stress tensor, ${\mathbf{e}}(\mathbf{u}_{f}):=\dfrac{1}{2}\,\left(\nabla\mathbf{u}_{f}+(\nabla\mathbf{u}_{f})^{\mathrm{t}}\right)$ stands for the deformation rate tensor, $\Gamma_{f}=\Gamma^{\mathrm{N}}_{f}\cup\Gamma^{\mathrm{D}}_{f}$, and $T>0$ is the final time. Next, we adopt the approach from [34, 1], and include as a new variable the vorticity tensor ${\boldsymbol{\gamma}}_{f}$,

In this way, owing to the fact that $\mathrm{tr}({\mathbf{e}}(\mathbf{u}_{f}))=\mathrm{div}(\mathbf{u}_{f})=q_{f}$, we find that (2.1) can be rewritten, equivalently, as the set of equations with unknowns ${\boldsymbol{\sigma}}_{f},{\boldsymbol{\gamma}}_{f}$ and $\mathbf{u}_{f}$, given by

Notice that the fourth equation in (2.2) has allowed us to eliminate the pressure $p_{f}$ from the system and provides a formula for its approximation through a post-processing procedure. For simplicity we assume that $|\Gamma^{\mathrm{N}}_{f}|>0$, which will allow us to control ${\boldsymbol{\sigma}}_{f}$ by ${\boldsymbol{\sigma}}_{f}^{\mathrm{d}}$. The case $|\Gamma^{\mathrm{N}}_{f}|=0$ can be handled as in [34, 35, 36] by introducing an additional variable corresponding to the mean value of $\mathrm{tr}({\boldsymbol{\sigma}}_{f})$.

In turn, let ${\boldsymbol{\sigma}}_{e}$ and ${\boldsymbol{\sigma}}_{p}$ be the elastic and poroelastic stress tensors, respectively, satisfying

where $0<\alpha_{p}\leq 1$ is the Biot–Willis constant, and $A$ is the symmetric and positive definite compliance tensor, which in the isotropic case has the form, for all tensors ${\boldsymbol{\tau}}$,

satisfying

In this case, ${\boldsymbol{\sigma}}_{e}\,:=\,\lambda_{p}\,\mathrm{div}({\boldsymbol{\eta}}_{p})\,\mathbf{I}+2\,\mu_{p}\,{\mathbf{e}}({\boldsymbol{\eta}}_{p})$, and $0<\lambda_{\min}\leq\lambda_{p}(\mathbf{x})\leq\lambda_{\max}$ and $0<\mu_{\min}\leq\mu_{p}(\mathbf{x})\leq\mu_{\max}$ are the Lamé parameters. The poroelasticity region $\Omega_{p}$ is governed by the quasi-static Biot system [14]:

where $\Gamma_{p}=\Gamma^{\mathrm{N}}_{p}\cup\Gamma^{\mathrm{D}}_{p}=\tilde{\Gamma}^{\mathrm{N}}_{p}\cup\tilde{\Gamma}^{\mathrm{D}}_{p}$, $s_{0}>0$ is a storativity coefficient and $\mathbf{K}(\mathbf{x})$ is the symmetric and uniformly positive definite rock permeability tensor, satisfying, for some constants $0<k_{\min}\leq k_{\max}$,

To avoid the issue with restricting the mean value of the pressure, we assume that $|\Gamma^{\mathrm{D}}_{p}|>0$. We also assume that $\Gamma^{\mathrm{D}}_{f}$, $\Gamma^{\mathrm{D}}_{p}$, and $\tilde{\Gamma}^{\mathrm{D}}_{p}$ are not adjacent to the interface $\Gamma_{fp}$, i.e., $\exists\ s>0$ such that $\mathrm{dist}\,(\Gamma^{\mathrm{D}}_{f},\Gamma_{fp})\geq s$, $\mathrm{dist}\,(\Gamma^{\mathrm{D}}_{p},\Gamma_{fp})\geq s$, and $\mathrm{dist}\,(\tilde{\Gamma}^{\mathrm{D}}_{p},\Gamma_{fp})\geq s$. This assumption is used to simplify the characterization of the normal trace spaces on $\Gamma_{fp}$.

Next, we introduce the following transmission conditions on the interface $\Gamma_{fp}$ [53, 12, 20, 8]:

where ${\mathbf{t}}_{f,j}$, $1\leq j\leq n-1$, is an orthogonal system of unit tangent vectors on $\Gamma_{fp}$, $\mathbf{K}_{j}=(\mathbf{K}\,{\mathbf{t}}_{f,j})\cdot{\mathbf{t}}_{f,j}$, and $\alpha_{\mathtt{BJS}}\geq 0$ is an experimentally determined friction coefficient. The first and second equations in (2.8) correspond to mass conservation and conservation of momentum on $\Gamma_{fp}$, respectively, whereas the third one can be decomposed into its normal and tangential components, as follows:

representing balance of normal stress and the Beaver–Joseph–Saffman (BJS) slip with friction condition, respectively.

Finally, the above system of equations is complemented by the initial condition $p_{p}(\mathbf{x},0)=p_{p,0}(\mathbf{x})$ in $\Omega_{p}$. We stress that, similarly to [47], compatible initial data for the rest of the variables can be constructed from $p_{p,0}$ in a way that all equations in the system (2.2)–(2.8), except for the unsteady conservation of mass equation in the first row of (2.6), hold at $t=0$. This will be established in Lemma 4.9 below. We will consider a weak formulation with a time-differentiated elasticity equation and compatible initial data $({\boldsymbol{\sigma}}_{p,0},p_{p,0})$.

In this section we proceed analogously to [3, Section 3] (see also [34]) and derive a weak formulation of the coupled problem given by (2.2), (2.3)–(2.6), and (2.8).