Domain decomposition and partitioning methods for mixed finite element discretizations of the Biot system of poroelasticity

In this paper we study several non-overlapping domain decomposition methods for the Biot system of poroelasticity [14], which models the flow of a viscous fluid through a poroelastic medium along with the deformation of the medium. Such flow occurs in many geophysics phenomena like earthquakes, landslides, and flow of oil inside mineral rocks and plays a key role in engineering applications such as hydrocarbon extraction through hydraulic or thermal fracturing. We use the classical Biot system of poroelasticity with the quasi-static assumption, which is particularly relevant in geoscience applications. The model consists of an equilibrium equation for the solid medium and a mass balance equation for the flow of the fluid through the medium. The system is fully coupled, with the fluid pressure contributing to the solid stress and the divergence of the solid displacement affecting the fluid content.

The numerical solution of the Biot system has been extensively studied in the literature. Various formulations have been considered, including two-field displacement–pressure formulations [27, 45, 51], three-field displacement–pressure–Darcy velocity formulations [48, 52, 62, 32, 41, 60, 49, 59], and three-field displacement–pressure–total pressure formulations [42, 47]. More recently, fully-mixed formulations of the Biot system have been studied. In [61], a four-field stress–displacement–pressure–Darcy velocity mixed formulation is developed. A posteriori error estimates for this formulation are obtained in [2]. In [39], a weakly symmetric stress–displacement–rotation elasticity formulation is considered, which is coupled with a mixed pressure–Darcy velocity flow formulation. Fully-mixed finite element approximations carry the advantages of local mass and momentum conservation, direct computation of the fluid velocity and the solid stress, as well as robustness and locking-free properties with respect to the physical parameters. Mixed finite element (MFE) methods can also handle discontinuous full tensor permeabilities and Lamé coefficients that are typical in subsurface flows. In our work we focus on the five-field weak-stress-symmetry formulation from [39], since weakly symmetric MFE methods for elasticity allow for reduced number of degrees of freedom. Moreover, a multipoint stress–multipoint flux mixed finite element approximation for this formulation has been recently developed in [7], which can be reduced to a positive definite cell-centered scheme for pressure and displacement only, see also a related finite volume method in [46]. While our domain decomposition methods are developed for the weakly symmetric formulation from [39], the analysis carries over in a straightforward way to the strongly symmetric formulation from [61].

Discretizations of the Biot system of poroelasticity for practical applications typically result in large algebraic systems of equations. The efficient solution of these systems is critical for the ability of the numerical method to provide the desired resolution. In this work we focus on non-overlapping domain decomposition methods [57, 50]. These methods split the computational domain into multiple non-overlapping subdomains with algebraic systems of lower complexity that are easier to solve. A global problem enforcing appropriate interface conditions is solved iteratively to recover the global solution. This approach naturally leads to scalable parallel algorithms. Despite the abundance of works on discretizations of the Biot system, there have been very few results on domain decomposition methods for this problem. In [28], a domain decomposition method using mortar elements for coupling the poroelastic model with an elastic model in an adjacent region is presented. In that work, the Biot region is not decomposed into subdomains. In [26, 25], an iterative coupling method is employed for a two-field displacement–pressure formulation, and classical domain decomposition techniques are applied separately for the elasticity and flow equations. A monolithic domain decomposition method for the two-field formulation of poroelasticity combining primal and dual variables is developed in [31]. To the best of our knowledge, domain decomposition methods for mixed formulations of poroelasticity have not been studied in the literature.

In this paper we study both monolithic and split non-overlapping domain decomposition methods for the five-field fully mixed formulation of poroelasticity with weak stress symmetry from [39, 7]. Monolithic methods require solving the coupled Biot system, while split methods only require solving elasticity and flow problems separately. Both methods have advantages and disadvantages. Monolithic methods involve solving larger and possibly ill-conditioned algebraic systems, but may be better suitable for problems with strong coupling between flow and mechanics, in which case split or iterative coupling methods may suffer from stability or convergence issues and require sufficiently small time steps. Our methods are motivated by the non-overlapping domain decomposition methods for MFE discretizations of Darcy flow developed in [29, 22, 8] and the non-overlapping domain decomposition methods for MFE discretizations of elasticity developed recently in [35].

In the first part of the paper we develop a monolithic domain decomposition method. We employ a physically heterogeneous Lagrange multiplier vector consisting of displacement and pressure variables to impose weakly the continuity of the normal components of stress and velocity, respectively. The algorithm involves solving at each time step an interface problem for this Lagrange multiplier vector. We show that the interface operator is positive definite, although it is not symmetric in general. As a result, a Krylov space solver such as GMRES can be employed for the solution of the interface problem. Each iteration requires solving monolithic Biot subdomain problems with specified Dirichlet data on the interfaces, which can be done in parallel. We establish lower and upper bounds on the spectrum of the interface operator, which allows us to perform analysis of the convergence of the GMRES iteration using field-of-values estimates.

In the second part of the paper we study split domain decomposition methods for the Biot system. Split or iterative coupling methods for poroelasticity have been extensively studied due to their computational efficiency. Four widely used sequential methods are drained split, undrained split, fixed stress split, and fixed strain split. Decoupling methods are prone to stability issues and a detailed stability analysis of the aforementioned schemes using finite volume methods can be found in [36, 37], see also [20] for stability analysis of several split methods using displacement–pressure finite element discretizations. Iterative coupling methods are based on similar splittings and involve iterating between the two sub-systems until convergence. Convergence for non-mixed finite element methods is analyzed in [44], while convergence for a four-field mixed finite element discretization is studied in [63]. An accelerated fixed stress splitting scheme for a generalized non-linear consolidation of unsaturated porous medium is studied in [18]. Studies of the optimization and acceleration of the fixed stress decoupling method for the Biot consolidation model, including techniques such as multirate or adaptive time stepping and parallel-in-time splittings have been presented in [1, 13, 17, 56, 3].

In our work we consider drained split (DS) and fixed stress (FS) decoupling methods in conjunction with non-overlapping domain decomposition. In particular, at each time step we solve sequentially an elasticity and a flow problem in the case of DS or a flow and an elasticity problem in the case of FS splitting. We perform stability analysis for the two splittings using energy estimates and show that they are both unconditionally stable with respect to the time step and the physical parameters. We then employ separate non-overlapping domain decomposition methods for each of the decoupled problems, using the methods developed in [29, 8] for flow and [35] for mechanics.

In the numerical section we present several computational experiments designed to verify and compare the accuracy, stability, and computational efficiency of the three domain decomposition methods for the Biot system of poroelasticity. In particular, we study the discretization error and the number of interface iterations, as well as the effect of the number of subdomains. We also illustrate the performance of the methods for a physically realistic heterogeneous problem with data taken from the Society of Petroleum Engineers 10th Comparative Solution Project.

The rest of the paper is organized as follows. In Section 2 we introduce the mathematical model and its MFE discretization. The monolithic domain decomposition method is developed and analyzed in Section 3. In Section 4 we perform stability analysis of the DS and FS decoupling methods and present the DS and FS domain decomposition methods. The numerical experiments are presented in Section 5, followed by conclusions in Section 6.

Let $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ be a simply connected domain occupied by a linearly elastic porous body. We use the notation $\mathbb{M}$, $\mathbb{S}$, and $\mathbb{N}$ for the spaces of $d\times d$ matrices, symmetric matrices, and skew-symmetric matrices respectively, all over the field of real numbers, respectively. Throughout the paper, the divergence operator is the usual divergence for vector fields, which produces vector field when applied to matrix field by taking the divergence of each row. For a domain $G\subset\mathbb{R}^{d}$, the $L^{2}(G)$ inner product and norm for scalar, vector, or tensor valued functions are denoted $\left(\cdot,\cdot\right)_{G}$ and $\|\cdot\|_{G}$, respectively. The norms and seminorms of the Hilbert spaces $H^{k}(G)$ are denoted by $\|\cdot\|_{k,G}$ and $|\cdot|_{k,G}$, respectively. We omit $G$ in the subscript if $G=\Omega$. For a section of the domain or element boundary $S\subset\mathbb{R}^{d-1}$ we write $\langle\cdot,\cdot\rangle_{S}$ and $\|\cdot\|_{S}$ for the $L^{2}(S)$ inner product (or duality pairing) and norm, respectively. We will also use the spaces

	$\displaystyle H(\operatorname{div};\Omega)=\{q\in L^{2}(\Omega,\mathbb{R}^{d}):\operatorname{div}q\in L^{2}(\Omega)\},$
	$\displaystyle H(\operatorname{div};\Omega,\mathbb{M})=\{\tau\in L^{2}(\Omega,\mathbb{M}):\operatorname{div}\tau\in L^{2}(\Omega,\mathbb{R}^{d})\},$

equipped with the norm

$$\|\tau\|_{\operatorname{div}}=\left(\|\tau\|^{2}+\|\operatorname{div}\tau\|^{2}\right)^{1/2}.$$

The partial derivative operator with respect to time, $\frac{\partial}{\partial t}$, is often abbreviated to $\partial_{t}$. Finally, $C$ denotes a generic positive constant that is independent of the discretization parameters $h$ and $\Delta t$.

Given a vector field $f$ representing body forces and a source term $g$, the quasi-static Biot system of poroelasticity is [14]:

	$\displaystyle-\operatorname{div}\sigma(u)$	$\displaystyle=f,$	$\displaystyle\text{in $\Omega\times(0,T]$},$		(2.1)
	$\displaystyle{K^{-1}}z+\nabla p$	$\displaystyle=0,$	$\displaystyle\text{in $\Omega\times(0,T]$},$		(2.2)
	$\displaystyle\frac{\partial}{\partial t}(c_{0}p+\alpha\operatorname{div}u)+\operatorname{div}z$	$\displaystyle=g,$	$\displaystyle\text{in $\Omega\times(0,T]$},$		(2.3)

where $u$ is the displacement, $p$ is the fluid pressure, $z$ is the Darcy velocity, and $\sigma$ is the poroelastic stress, defined as

$$\sigma=\sigma_{e}-\alpha pI.$$

(2.4)

Here $I$ is the $d\times d$ identity matrix, $0<\alpha\leq 1$ is the Biot-Willis constant, and $\sigma_{e}$ is the elastic stress satisfying the stress-strain relationship

$$A\sigma_{e}=\epsilon(u),\quad\epsilon(u):=\frac{1}{2}\left(\nabla u+(\nabla u)^{T}\right),$$

(2.5)

where $A$ is the compliance tensor, which is a symmetric, bounded and uniformly positive definite linear operator acting from $\mathbb{S}\to\mathbb{S}$, extendible to $\mathbb{M}\to\mathbb{M}$. In the special case of homogeneous and isotropic body, $A$ is given by,

$$A\sigma=\frac{1}{2\mu}\left(\sigma-\frac{\lambda}{2\mu+d\lambda}\operatorname{tr}(\sigma)I\right),$$

(2.6)

where $\mu>0$ and $\lambda\geq 0$ are the Lamé coefficients. In this case, $\sigma_{e}(u)=2\mu\epsilon(u)+\lambda\operatorname{div}u\,I$. Finally, $K$ stands for the permeability tensor, which is symmetric, bounded, and uniformly positive definite, and $c_{0}\geq 0$ is the mass storativity. To close the system, we impose the boundary conditions

	$\displaystyle u$	$\displaystyle=g_{u}$	on $\Gamma_{D}^{u}\times(0,T]$	,	$\displaystyle\sigma n$	$\displaystyle=0$	$\displaystyle\text{on $\Gamma_{N}^{\sigma}\times(0,T]$},$		(2.7)
	$\displaystyle p$	$\displaystyle=g_{p}$	on $\Gamma_{D}^{p}\times(0,T]$	,	$\displaystyle z\cdot n$	$\displaystyle=0$	$\displaystyle\text{on $\Gamma_{N}^{z}\times(0,T]$},$		(2.8)

where $\Gamma_{D}^{u}\cup\Gamma_{N}^{\sigma}=\Gamma_{D}^{p}\cup\Gamma_{N}^{z}=\partial\Omega$ and $n$ is the outward unit normal vector field on $\partial\Omega$, along with the initial condition $p(x,0)=p_{0}(x)$ in $\Omega$. Compatible initial data for the rest of the variables can be obtained from (2.1) and (2.2) at $t=0$. Well posedness analysis for this system can be found in [53].

We consider a mixed variational formulation for (2.1)–(2.8) with weak stress symmetry, following the approach in [39]. The motivation is that MFE elasticity spaces with weakly symmetric stress tend to have fewer degrees of freedom than strongly symmetric MFE spaces. Moreover, in a recent work, a multipoint stress–multipoint flux mixed finite element method has been developed for this formulation that reduces to a positive definite cell-centered scheme for pressure and displacement only [7]. Nevertheless, the domain decomposition methods in this paper can be employed for strongly symmetric stress formulations, with the analysis carrying over in a straightforward way. We introduce a rotation Lagrange multiplier $\gamma:=\frac{1}{2}\left(\nabla u-\nabla u^{T}\right)\in\mathbb{N}$, which is used to impose weakly symmetry of the stress tensor $\sigma$. We rewrite (2.5) as

$$A\left(\sigma+\alpha pI\right)=\nabla u-\gamma.$$

(2.9)

Combining (2.5) and (2.4) gives $\operatorname{div}u=\operatorname{tr}(\epsilon(u))=\operatorname{tr}(A\sigma_{e})=\operatorname{tr}A(\sigma+\alpha pI)$, which can be used to rewrite (2.3) as

$$\partial_{t}(c_{0}p+\alpha\operatorname{tr}A\left(\sigma+\alpha pI\right))+\operatorname{div}z=g.$$

(2.10)

The combination of (2.9), (2.1), (2.2), and (2.10), along with the boundary conditions (2.7)–(2.8), leads to the variational formulation: find $(\sigma,u,\gamma,z,p):[0,T]\to\mathbb{X}\times V\times\mathbb{Q}\times Z\times W$ such that $p(0)=p_{0}$ and for a.e. $t\in(0,T)$,

	$\displaystyle\left(A(\sigma+\alpha pI),\tau\right)+\left(u,\operatorname{div}{\tau}\right)+\left(\gamma,\tau\right)=\langle g_{u},\tau\,n\rangle_{\Gamma_{D}^{u}},$	$\displaystyle\forall\tau\in\mathbb{X},$		(2.11)
	$\displaystyle\left(\operatorname{div}{\sigma},v\right)=-\left(f,v\right),$	$\displaystyle\forall v\in V,$		(2.12)
	$\displaystyle\left(\sigma,\xi\right)=0,$	$\displaystyle\forall\xi\in\mathbb{Q},$		(2.13)
	$\displaystyle\left({K^{-1}}z,q\right)-\left(p,\operatorname{div}{q}\right)=-\langle g_{p},q\cdot n\rangle_{\Gamma_{D}^{p}},$	$\displaystyle\forall q\in Z,$		(2.14)
	$\displaystyle c_{0}\left(\partial_{t}{p},w\right)+\alpha\left(\partial_{t}A(\sigma+\alpha pI),wI\right)+\left(\operatorname{div}{z},w\right)=\left(g,w\right),$	$\displaystyle\forall w\in W,$		(2.15)

where

	$\displaystyle\mathbb{X}=\big{\{}\tau\in H(\operatorname{div};\Omega,\mathbb{M}):\tau\,n=0\text{ on }\Gamma_{N}^{\sigma}\big{\}},\quad V=L^{2}(\Omega,\mathbb{R}^{d}),\quad\mathbb{Q}=L^{2}(\Omega,\mathbb{N}),$
	$\displaystyle Z=\big{\{}q\in H(\operatorname{div};\Omega):q\cdot n=0\text{ on }\Gamma_{N}^{z}\big{\}},\quad W=L^{2}(\Omega).$

It was shown in [7] that the system (2.11)–(2.15) is well posed.

Next, we present the MFE discretization of (2.11)–(2.15). For simplicity we assume that $\Omega$ is a Lipshicz polygonal domain. Let $\mathcal{T}_{h}$ be a shape-regular quasi-uniform finite element partition of $\Omega$, consisting of simplices or quadrilaterals, with $h=\text{max}_{E\in\mathcal{T}_{h}}\text{diam}(E)$. The MFE method for solving (2.11)–(2.15) is: find $(\sigma_{h},u_{h},\gamma_{h},z_{h},p_{h}):[0,T]\to\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}\times Z_{h}\times W_{h}$ such that, for a.e. $t\in(0,T)$,

	$\displaystyle\left(A(\sigma_{h}+\alpha p_{h}I),\tau\right)+\left(u_{h},\operatorname{div}{\tau}\right)+\left(\gamma_{h},\tau\right)=\langle g_{u},\tau\,n\rangle_{\Gamma_{D}^{u}},$	$\displaystyle\forall\tau\in\mathbb{X}_{h},$		(2.16)
	$\displaystyle\left(\operatorname{div}{\sigma_{h}},v\right)=-\left(f,v\right),$	$\displaystyle\forall v\in V_{h},$		(2.17)
	$\displaystyle\left(\sigma_{h},\xi\right)=0,$	$\displaystyle\forall\xi\in\mathbb{Q}_{h},$		(2.18)
	$\displaystyle\left({K^{-1}}z_{h},q\right)-\left(p_{h},\operatorname{div}{q}\right)=-\langle g_{p},q\cdot n\rangle_{\Gamma_{D}^{p}},$	$\displaystyle\forall q\in Z_{h},$		(2.19)
	$\displaystyle c_{0}\left(\partial_{t}{p_{h}},w\right)+\alpha\left(\partial_{t}A(\sigma_{h}+\alpha p_{h}I),wI\right)+\left(\operatorname{div}{z_{h}},w\right)=\left(g,w\right),$	$\displaystyle\forall w\in W_{h},$		(2.20)

with discrete initial data obtained as the elliptic projection of the continuous initial data. Here $\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}\times Z_{h}\times W_{h}\subset\mathbb{X}\times V\times\mathbb{Q}\times Z\times W$ is a collection of suitable finite element spaces. In particular, $\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}$ could be chosen from any of the known stable triplets for linear elasticity with weak stress symmetry, e.g. [55, 5, 9, 10, 11, 21, 30, 15, 16, 24, 6, 40], satisfying the inf-sup condition

$\displaystyle\forall v\in V_{h},\xi\in\mathbb{Q}_{h},\quad\|v\|+\|\xi\|$

$\displaystyle\leq C\sup_{0\neq\tau\in\mathbb{X}_{h}}\frac{\left(v,\operatorname{div}{\tau}\right)+\left(\xi,\tau\right)}{\|\tau\|_{\text{div}}}.$

(2.21)

For the flow part, $Z_{h}\times W_{h}$ could be chosen from any of the known stable velocity-pressure pairs of MFE spaces such as the Raviart-Thomas ($\mathcal{RT}$) or Brezzi-Douglas-Marini ($\mathcal{BDM}$) spaces, see [19], satisfying the inf-sup condition

$\displaystyle\forall w\in W_{h},\quad\|w\|$

$\displaystyle\leq C\sup_{0\neq q\in Z_{h}}\frac{(\operatorname{div}q,w)}{\|q\|_{\text{div}}}.$

(2.22)

Let $\Omega=\cup_{i=1}^{m}\Omega_{i}$ be a union of non-overlapping shape-regular polygonal subdomains, where each subdomain is a union of elements of $\mathcal{T}_{h}$. Let $\Gamma_{i,j}={\partial\Omega}_{i}\cap{\partial\Omega}_{j},\,\Gamma=\cup_{i,j=1}^{m}\Gamma_{i,j},$ and $\Gamma_{i}={\partial\Omega}_{i}\cap\Gamma={\partial\Omega}_{i}\setminus{\partial\Omega}$ denote the interior subdomain interfaces. Denote the restriction of the spaces $\mathbb{X}_{h}$, $V_{h}$, $\mathbb{Q}_{h}$, $Z_{h}$, and $W_{h}$ to $\Omega_{i}$ by $\mathbb{X}_{h,i}$, $V_{h,i}$, $\mathbb{Q}_{h,i}$, $Z_{h,i}$, and $W_{h,i}$, respectively. Let $\mathcal{T}_{h,i,j}$ be a finite element partition of $\Gamma_{i,j}$ obtained from the trace of $\mathcal{T}_{h}$, and let $n_{i,j}$ be a unit normal vector on $\Gamma_{i,j}$ with an arbitrarily fixed direction. In the domain decomposition formulation we utilize a vector Lagrange multiplier $\lambda_{h}=(\lambda_{h}^{u},\lambda_{h}^{p})^{T}$ approximating the displacement and the pressure on the interface and used to impose weakly the continuity of the normal components of the poroelastic stress tensor $\sigma$ and the velocity vector $z$, respectively. We define the Lagrange multiplier space on $\mathcal{T}_{i,j}$ and $\cup_{i<j}\mathcal{T}_{i,j}$ as follows:

$$\Lambda_{h,i,j}:=\begin{pmatrix}\Lambda_{h,i,j}^{u}\\ \Lambda_{h,i,j}^{p}\end{pmatrix}:=\begin{pmatrix}\mathbb{X}_{h}\,n_{i,j}\\ Z_{h}\cdot n_{i,j}\end{pmatrix},\quad\Lambda_{h}^{u}:=\bigoplus_{1\leq i<j\leq m}\Lambda_{h,i,j}^{u},\quad\Lambda_{h}^{p}:=\bigoplus_{1\leq i<j\leq m}\Lambda_{h,i,j}^{p},\quad\Lambda_{h}:=\begin{pmatrix}\Lambda_{h}^{u}\\ \Lambda_{h}^{p}\end{pmatrix}.$$

The domain decomposition formulation for the mixed Biot problem in a semi-discrete form reads as follows: for $1\leq i\leq m$, find $(\sigma_{h,i},u_{h,i},\gamma_{h,i},z_{h,i},p_{h,i},\lambda_{h}):[0,T]\to\mathbb{X}_{h,i}\times V_{h,i}\times\mathbb{Q}_{h,i}\times Z_{h,i}\times W_{h,i}\times\Lambda_{h}$ such that, for a.e. $t\in(0,T)$,

	$\displaystyle\left(A(\sigma_{h,i}+\alpha p_{h,i}I),\tau\right)_{\Omega_{i}}+\left(u_{h,i},\operatorname{div}{\tau}\right)_{\Omega_{i}}+\left(\gamma_{h,i},\tau\right)_{\Omega_{i}}$
	$\displaystyle\qquad\quad=\langle g_{u},\tau\,n_{i}\rangle_{{\partial\Omega}_{i}\cap\Gamma_{D}^{u}}+\langle\lambda_{h}^{u},\tau\,n_{i}\rangle_{\Gamma_{i}},$	$\displaystyle\forall\tau\in\mathbb{X}_{h,i},$		(3.1)
	$\displaystyle\left(\operatorname{div}{\sigma_{h,i}},v\right)_{\Omega_{i}}=-\left(f,v\right)_{\Omega_{i}},$	$\displaystyle\forall v\in V_{h,i},$		(3.2)
	$\displaystyle\left(\sigma_{h,i},\xi\right)_{\Omega_{i}}=0,$	$\displaystyle\forall\xi\in\mathbb{Q}_{h,i},$		(3.3)
	$\displaystyle\left({K^{-1}}z_{h,i},q\right)_{\Omega_{i}}-\left(p_{h,i},\operatorname{div}{q}\right)_{\Omega_{i}}=-\langle g_{p},q\cdot n_{i}\rangle_{{\partial\Omega}_{i}\cap\Gamma_{D}^{p}}-\langle\lambda_{h}^{p},q\cdot n_{i}\rangle_{\Gamma_{i}},$	$\displaystyle\forall q\in Z_{h,i},$		(3.4)
	$\displaystyle c_{0}\left(\partial_{t}{p_{h,i}},w\right)_{\Omega_{i}}+\alpha\left(\partial_{t}A(\sigma_{h,i}+\alpha p_{h,i}I),wI\right)_{\Omega_{i}}+\left(\operatorname{div}{z_{h,i}},w\right)_{\Omega_{i}}=\left(g,w\right)_{\Omega_{i}},$	$\displaystyle\forall w\in W_{h,i},$		(3.5)
	$\displaystyle\sum_{i=1}^{m}\langle\sigma_{h,i}\,n_{i},\mu^{u}\rangle_{\Gamma_{i}}=0,$	$\displaystyle\forall\mu^{u}\in\Lambda_{h}^{u},$		(3.6)
	$\displaystyle\sum_{i=1}^{m}\langle z_{h,i}\cdot n_{i},\mu^{p}\rangle_{\Gamma_{i}}=0,$	$\displaystyle\forall\mu^{p}\in\Lambda_{h}^{p},$		(3.7)

where $n_{i}$ is the outward unit normal vector field on $\Omega_{i}$. We note that both the elasticity and flow subdomain problems in the above method are of Dirichlet type. It is easy to check that (3.1)–(3.7) is equivalent to the global formulation (2.16)–(2.20) with $(\sigma_{h},u_{h},\gamma_{h},z_{h},p_{h})|_{\Omega_{i}}=(\sigma_{h,i},u_{h,i},\gamma_{h,i},z_{h,i},p_{h,i})$.