# Analysis of a multiphysics finite element method for a poroelasticity model

Analysis of a multiphysics finite element method for a poroelasticity model Abstract This article concerns finite element approximations of a quasi-static poroelasticity model in displacement–pressure formulation, which describes the dynamics of poro-elastic materials under an applied mechanical force on the boundary. To better describe the multiphysics process of deformation and diffusion for poroelastic materials, we first present a reformulation of the original model by introducing two pseudo-pressures. We then propose a time-stepping algorithm that decouples the reformulated partial differential equation (PDE) problem at each time step into two sub-problems: one of which is a generalized Stokes problem for the displacement vector field (of the solid network of the poroelastic material) along with one pseudo-pressure field, and the other is a diffusion problem for the other pseudo-pressure field (of the solvent of the material). To make this multiphysics approach feasible numerically, two critical issues must be resolved: the first one is the uniqueness of the generalized Stokes problem and the other is to find a good boundary condition for the diffusion equation, so that it also becomes uniquely solvable. To address the first issue, we discover certain conserved quantities for the PDE solution that provide ideal candidates for a needed boundary condition for the pseudo-pressure field. The solution to the second issue is to use the generalized Stokes problem to generate a boundary condition for the diffusion problem. A practical advantage of the time-stepping algorithm allows one to use any convergent Stokes solver (and its code) together with any convergent diffusion equation solver (and its code) to solve the poroelasticity model. In this article, the Taylor–Hood mixed finite element method combined with the $$P_1$$-conforming finite element method is used as an example to demonstrate the viability of the proposed multiphysics approach. It is proved that the solutions of the fully discrete finite element methods fulfill a discrete energy law, which mimics the differential energy law satisfied by the PDE solution and converges optimally in the energy norm. Moreover, it is showed that the proposed formulation also has a built-in mechanism to overcome so-called ‘locking phenomenon’ associated with the numerical approximations of the poroelasticity model. Numerical experiments are presented to show the performance of the proposed approach and methods, and to demonstrate the absence of ‘locking phenomenon’ in our numerical experiments. This article also presents a detailed PDE analysis for the poroelasticity model, especially it is proved that this model converges to the well-known Biot’s consolidation model from soil mechanics as the constrained specific storage coefficient tends to zero. As a result, the proposed approach and methods are robust under such a limit process. 1. Introduction A poroelastic material (or medium) is a fluid–solid interaction (FSI) system at pore scale, and poro-mechanics is a branch of continuum mechanics and acoustics that studies the behavior of fluid-saturated porous materials. If the solid is an elastic material, then the subject of the study is known as poroelasticity. Moreover, the elastic material may be governed by linear or nonlinear constitutive law, which then leads, respectively, to linear and nonlinear poroelasticity. Examples of poroelastic materials include soil, polymer gels and medicine pills. Poroelastic materials exhibit an important state of matter found in a wide variety of mechanical, biomedical and chemical systems (cf. Biot, 1955; Hamley, 2007; Terzaghi, 1943; Tanaka & Fillmore, 1979; Doi & Edwards, 1986; Coussy, 2004; Yamaue & Doi, 2004 and the references therein). They also possess some fascinating properties, in particular, they display thixotropy, which means that they become fluid when agitated, but resolidify when resting. In general, the behavior of a poroelastic material is described by a multiphysics FSI process at pore scale. However, for poroelastic materials to be considered in this article, one important physical phenomenon of the multiphysics process is not explicitly revealed in their mathematical model instead, it is hidden in the model. This is a different feature (and a source of extra difficulties) of poroelastic materials compared with standard (macroscopic) FSI systems. This article considers a general quasi-static model of linear poroelasticity, which is broad enough to contain the well-known Biot’s consolidation model from soil mechanics (cf. Murad & Loula, 1992; Schreyer-Bennethum, 2007) and the Doi’s model for polymer gels (cf. Yamaue & Doi, 2004; Feng & He, 2010). The quasi-static feature is due to the assumption that the acceleration of the solid (described by the second-order time derivative of the displacement vector field) is assumed to be negligible. For the Biot’s model, some locking-free formulations have recently been studied in Lee (2014, 2016), Oyarzúa & Ruiz-Baier (2016) and Lee et al. (2017). We refer the reader to Terzaghi (1943), Coussy (2004) and Phillips & Wheeler (2007a) for a derivation of the model and to Schowalter (2000) for its mathematical analysis. When the parameter $$c_0$$, called the constrained specific storage coefficient, vanishes in the model, it reduces to the Boit’s model and Doi’s model arising from two distinct applications. Their mathematical analysis can be found in Feng & He (2010), and their finite element numerical approximations based on two very different approaches were carried out in Murad & Loula (1992) and Feng & He (2010), respectively. In Phillips & Wheeler (2007a,b), the authors proposed and analysed a semidiscrete and a fully discrete mixed finite element method, which simultaneously approximate the pressure and its gradient along with the displacement vector field. Since the implicit Euler scheme is used for the time discretization, a combined linear system must be solved at each time step. It is observed in the numerical tests that the proposed fully discrete mixed finite method may exhibit a ‘locking phenomenon’ in the sense that the computed pressure oscillates and its accuracy deteriorates when a rapidly changed initial pressure is given. It is explained in Phillips & Wheeler (2009) that such a ‘locking phenomenon’ is caused by the difficulty of satisfying the nearly divergence-free condition of the displacement field for very small time $$t$$. As a comparison, we recall that the ‘volumetric locking’ is used in poromechanics to describe the phenomenon when a finite element solution deteriorates as the Poisson ratio approaches $$0.5$$, that is, when the material becomes incompressible or nearly incompressible. The goal of this article is to present a multiphysics approach for approximating the poroelasticity model. A key idea of this approach is to derive a multiphysics reformulation for the original model that clearly reveals the underlying multiple physics process (i.e., the deformation and diffusion) of the pore-scale FSI system. At end, two pseudo-pressures are introduced, one of them is shown to satisfy a diffusion equation, whereas the displacement vector field along with the other pseudo-pressure variable is shown to satisfy a generalized Stokes system. It should be noted that the original pressure is eliminated in the reformulation, and thus it is not approximated as a primary variable; instead, it is computed as a linear combination of the two pseudo-pressures. On the basis of this multiphysics reformulation, we propose a time-stepping algorithm that decouples (or couples) the reformulated PDE problem at each time step into two sub-problems, a generalized Stokes problem for the displacement vector field along with a pseudo-pressures and a diffusion problem for another pseudo-pressure field. To make this multiphysics approach feasible numerically, two critical issues must be resolved: the first one is the uniqueness of the generalized Stokes problem and the other is to find a good boundary condition for the diffusion equation so that it also becomes uniquely solvable. To overcome these difficulties, we discover certain conserved quantities for the PDE solution that can be imposed as needed boundary conditions for the subproblems. Moreover, we demonstrate that, regardless the choice of discretization methods, the proposed formulation has a built-in mechanism to overcome the ‘locking phenomenon’ associated with numerical approximations of the poroelasticity model. The remainder of this article is organized as follows. In Section 2, we present a complete PDE analysis of the poroelasticity model that emphasizes the energy law of the underlying model. Several conserved quantities are derived for the PDE solution. Moreover, it is proved that the poroelasticity model converges to the Biot’s consolidation model as the constrained specific storage coefficient $$c_0\to 0$$. In Section 3, we propose and analyse some fully discrete finite element methods based on the multiphysics reformulation. Both coupled and decoupled time stepping are considered and compared. The Taylor–Hood mixed finite element method combined with the $$P_1$$-conforming finite element method is chosen as an example for spatial discretization. It is proved that the solution of the fully discrete finite element methods fulfills a discrete energy law that mimics the differential energy law satisfied by the PDE solution. Optimal order error estimates in the energy norm are also established. Finally, in Section 4, several benchmark numerical experiments are provided to show the performance of the proposed approach and methods, and to demonstrate the absence of ‘locking phenomenon’ in our numerical experiments. 2. PDE model and its analysis 2.1 Preliminaries $${\it{\Omega}} \subset \mathbb{R}^d \,(d=1,2,3)$$ denotes a bounded polygonal domain with the boundary $${\partial}{\it{\Omega}}$$. The standard function space notation is adopted in this article, and their precise definitions can be found in Temam (1977), Ciarlet (1978) and Brenner & Scott (2008). In particular, $$(\cdot,\cdot)$$ and $$\langle \cdot,\cdot\rangle$$ denote, respectively, the standard $$L^2({\it{\Omega}})$$ and $$L^2({\partial}{\it{\Omega}})$$ inner products. For any Banach space $$B$$, we let $$\mathbf{B}=[B]^d$$ and use $$\mathbf{B}^\prime$$ to denote its dual space. In particular, we use $$(\cdot,\cdot)_{\small\rm dual}$$ to denote the dual product on $${\mathbf{H}}^1({\it{\Omega}})' \times {\mathbf{H}}^1({\it{\Omega}})$$, and $$\Vert\,{\cdot}\,\Vert_{L^p(B)}$$ is a shorthand notation for $$\Vert\,{\cdot}\,\Vert_{L^p((0,T);B)}$$. We also introduce the function spaces   L02(Ω):={q∈L2(Ω);(q,1)=0},X:=H1(Ω). It is well known (Temam, 1977) that the following so-called inf–sup condition holds in the space $${\mathbf{X}}\times L^2_0({\it{\Omega}})$$:   supv∈X(divv,φ)‖v‖H1(Ω)≥α0‖φ‖L2(Ω)∀φ∈L02(Ω),α0>0. (2.1) Let   RM:={r:=a+b×x;a,b,x∈Rd} denote the space of infinitesimal rigid motions. It is well known (Temam, 1977; Girault & Raviart, 1981; Brenner, 1993) that $${\mathbf{RM}}$$ is the kernel of the strain operator $${\varepsilon}$$, that is, $${\mathbf{r}}\in {\mathbf{RM}}$$ if and only if $${\varepsilon}({\mathbf{r}})=0$$. Hence, we have   ε(r)=0divr=0,∀r∈RM. (2.2) Let $${\mathbf{L}}^2_\bot({\partial}{\it{\Omega}})$$ and $${\mathbf{H}}^1_\bot({\it{\Omega}})$$ denote, respectively, the subspaces of $${\mathbf{L}}^2({\partial}{\it{\Omega}})$$ and $${\mathbf{H}}^1({\it{\Omega}}),$$ which are orthogonal to $${\mathbf{RM}}$$, that is,   H⊥1(Ω):={v∈H1(Ω);(v,r)=0∀r∈RM},L⊥2(∂Ω):={g∈L2(∂Ω);⟨g,r⟩=0∀r∈RM}. It is well known (Dautray & Lions, 1990) that there exists a constant $$c_1>0$$ such that   infr∈RM‖v+r‖L2(Ω)≤c1‖ε(v)‖L2(Ω)∀v∈H1(Ω). Hence, for each $${\mathbf{v}}\in {\mathbf{H}}^1_\bot({\it{\Omega}})$$ there holds   ‖v‖L2(Ω)=infr∈RM‖v+r‖L2(Ω)2−‖r‖L2(Ω)2≤c1‖ε(v)‖L2(Ω), (2.3) which, and the well-known Korn’s inequality (Dautray & Lions, 1990), yield that for some $$c_2>0$$  ‖v‖H1(Ω) ≤c2[‖v‖L2(Ω)+‖ε(v)‖L2(Ω)] ≤c2(1+c1)‖ε(v)‖L2(Ω)∀v∈H⊥1(Ω). (2.4) By Lemma 2.1 of Brenner (1993), we know that for any $$q\in L^2({\it{\Omega}})$$, there exists $${\mathbf{v}}\in {\mathbf{H}}^1_\bot({\it{\Omega}})$$ such that $$\text{div}\,\, {\mathbf{v}}=q$$ and $$\|{\mathbf{v}}\|_{H^1({\it{\Omega}})} \leq C\|q\|_{L^2({\it{\Omega}})}$$. An immediate consequence of this lemma is that there holds the following alternative version of the inf–sup condition:   supv∈H⊥1(Ω)(divv,φ)‖v‖H1(Ω)≥α1‖φ‖L2(Ω)∀φ∈L02(Ω),α1>0. (2.5) Throughout the article, we assume $${\it{\Omega}} \subset \mathbb{R}^d$$ is a bounded polygonal domain such that $${\it{\Delta}}: H^1_0({\it{\Omega}}) \cap H^2({\it{\Omega}}) \rightarrow L^2({\it{\Omega}})$$ is an isomorphism (cf. Girault & Raviart, 1981; Dauge, 1988). In addition, $$C$$ is used to denote a generic positive (pure) constant which may be different in different places. 2.2 PDE model and its multiphysics reformulation The quasi-static poroelasticity model to be studied in this article is given by Phillips & Wheeler (2007a)   −divσ(u)+α∇p =fin ΩT:=Ω×(0,T)⊂Rd×(0,T), (2.6)  (c0p+αdivu)t+divvf =ϕin ΩT, (2.7) where   σ(u) :=με(u)+λdivuI,ε(u):=12(∇u+∇uT), (2.8)  vf :=−Kμf(∇p−ρfg). (2.9) Here, $${\mathbf{u}}$$ denotes the displacement vector of the solid and $$p$$ denotes the pressure of the solvent. $${\mathbf{f}}$$ is the body force and $$\mathbf{g}$$ is the gravitational acceleration, which is assumed to be a constant vector. $$I$$ denotes the $$d\times d$$ identity matrix, and $${\varepsilon}({\mathbf{u}})$$ is known as the strain tensor. The parameters in the model are Lamé constants $${\lambda}$$ and $$\mu$$ (warning: for the sake of notation brevity later, we use $$\mu$$ instead of $$2\mu$$ in (2.8) as often seen in the literature, in other words, precisely, here $$\mu/2$$ is a Lamé constant); the permeability tensor $$K=K(x)$$, which is assumed to be symmetric and uniformly positive definite in the sense that there exist positive constants $$K_1$$ and $$K_2$$ such that $$K_1|\zeta|^2\leq K(x)\zeta\cdot \zeta \leq K_2 |\zeta|^2$$ for a.e. $$x\in{\it{\Omega}}$$ and any $$\zeta\in \mathbf{\mathbb{R}}^d$$; the solvent viscosity $$\mu_f$$, Biot–Willis constant $$\alpha$$, and the constrained specific storage coefficient $$c_0$$. In addition, $$\sigma({\mathbf{u}})$$ is called the (effective) stress tensor. $$\widehat{\sigma}({\mathbf{u}},p):=\sigma({\mathbf{u}})-\alpha p I$$ is the total stress tensor. $${\mathbf{v}}_f$$ is the volumetric solvent flux and (2.9) is the well-known Darcy’s law. We assume that $$\rho_f\not\equiv 0$$, which is a realistic assumption. To close the above system, suitable boundary and initial conditions must be prescribed. The following set of boundary and initial conditions will be considered in this article:   σ^(u,p)n=σ(u)n−αpn =f1on ∂ΩT:=∂Ω×(0,T), (2.10)  vf⋅n=−Kμf(∇p−ρfg)⋅n =ϕ1on ∂ΩT, (2.11)  u=u0,p =p0in Ω×{t=0}. (2.12) We note that in some engineering literature, the second Lamé constant $$\mu$$ is also called the shear modulus and denoted by $$G$$, and $$B:={\lambda} +\frac23 G$$ is called the bulk modulus. $${\lambda}, \mu$$ and $$B$$ are computed from the Young’s modulus$$E$$ and the Poisson ratio$$\nu$$ by the following formulas:   λ=Eν(1+ν)(1−2ν),μ=G=E2(1+ν),B=E3(1−2ν). Unlike the existing approaches in the literature (Murad & Loula, 1992; Phillips & Wheeler, 2007a), in this article we will not directly approximate the above original model; instead, we first derive a (multiphysics) reformulation for the model, we then approximate the reformulated model. This is a key idea of this article, and it will be seen in the later sections that this new approach is advantageous. To this end, we introduce new variables   q:=divu,η:=c0p+αq,ξ:=αp−λq. It is easy to check that   p=κ1ξ+κ2η,q=κ1η−κ3ξ, (2.13) where   κ1:=αα2+λc0,κ2:=λα2+λc0,κ3:=c0α2+λc0. (2.14) Then (2.6)–(2.9) can be written as   −μdivε(u)+∇ξ =fin ΩT, (2.15)  κ3ξ+divu =κ1ηin ΩT, (2.16)  ηt−1μfdiv[K(∇(κ1ξ+κ2η)−ρfg)] =ϕin ΩT, (2.17) where $$p$$ and $$q$$ are related to $$\xi$$ and $$\eta$$ through the algebraic equations in (2.13). The boundary and initial conditions (2.10)–(2.12) can be rewritten as   σ(u)n−α(κ1ξ+κ2η)n =f1on ∂ΩT:=∂Ω×(0,T), (2.18)  −Kμf(∇(κ1ξ+κ2η)−ρfg)⋅n =ϕ1on ∂ΩT, (2.19)  u=u0,p =p0in Ω×{t=0}. (2.20) It is now clear that $$({\mathbf{u}}, \xi)$$ satisfies a generalized Stokes problem for a given $$\eta$$, where $$\kappa_3\xi$$ in (2.16) acts as a penalty term and $$\eta$$ satisfies a diffusion problem for a given $$\xi$$. This new formulation reveals the underlying deformation and diffusion multiphysics process, which occurs in the poroelastic material. In particular, the diffusion part of the process is hidden in the original formulation, but is apparent in the new formulation. To make use of the above reformulation for computation, we need to address a crucial issue of the uniqueness for the generalized Stokes problem and the diffusion problem after they are decoupled. The difficulty will be overcome by discovering some invariant quantities for the solution of the PDE model and using them to impose some appropriate boundary conditions for both subproblems (cf. Lemma 2.9). 2.3 Analysis of the PDE model We start this section with a definition of weak solutions to problems (2.6)–(2.12). For convenience, we assume that $${\mathbf{f}}, {\mathbf{f}}_1,\phi$$ and $$\phi_1$$ all are independent of $$t$$ in the remaining of the article. We note that all the results of this article can be easily extended to the case of time-dependent source functions. Definition 2.1 Let $${\mathbf{u}}_0\in{\mathbf{H}}^1({\it{\Omega}}), {\mathbf{f}}\in{\mathbf{L}}^2({\it{\Omega}}), {\mathbf{f}}_1\in {\mathbf{L}}^2({\partial}{\it{\Omega}}), p_0\in L^2({\it{\Omega}}), \phi\in L^2({\it{\Omega}})$$ and $$\phi_1\in L^2({\partial}{\it{\Omega}})$$. Assume $$c_0>0$$ and $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$ for any $${\mathbf{v}}\in \mathbf{RM}$$. Given $$T > 0$$, a tuple $$({\mathbf{u}},p)$$ with   u∈L∞(0,T;H⊥1(Ω)),p∈L∞(0,T;L2(Ω))∩L2(0,T;H1(Ω)),pt,(divu)t∈L2(0,T;H1(Ω)′) is called a weak solution to (2.6)–(2.12), if it holds for almost every $$t \in [0,T]$$  μ(ε(u),ε(v))+λ(divu,divv)−α(p,divv) =(f,v)+⟨f1,v⟩∀v∈H1(Ω), (2.21)   ((c0p+αdivu)t,φ)dual+1μf(K(∇p−ρfg),∇φ) =(ϕ,φ)+⟨ϕ1,φ⟩|,∀φ∈H1(Ω), (2.22)  u(0)=u0,p(0)=p0. (2.23) Similarly, we can define weak solutions to problems (2.15)–(2.17) and (2.18)–(2.20). Definition 2.2 Let $${\mathbf{u}}_0\in {\mathbf{H}}^1({\it{\Omega}}), {\mathbf{f}} \in {\mathbf{L}}^2({\it{\Omega}}), {\mathbf{f}}_1 \in {\mathbf{L}}^2({\partial}{\it{\Omega}}), p_0\in L^2({\it{\Omega}}), \phi\in L^2({\it{\Omega}})$$ and $$\phi_1\in L^2({\partial}{\it{\Omega}})$$. Assume $$c_0>0$$ and $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$ for any $${\mathbf{v}}\in \mathbf{RM}$$. Given $$T > 0$$, a $$5$$-tuple $$({\mathbf{u}},\xi,\eta,p,q)$$ with   u∈L∞(0,T;H⊥1(Ω)),ξ∈L∞(0,T;L2(Ω)),η∈L∞(0,T;L2(Ω))∩H1(0,T;H1(Ω)′),q∈L∞(0,T;L2(Ω)),p∈L∞(0,T;L2(Ω))∩L2(0,T;H1(Ω)) is called a weak solution to (2.15)–(2.17) and (2.18)–(2.20) if it holds for almost every $$t \in [0,T]$$  μ(ε(u),ε(v))−(ξ,divv) =(f,v)+⟨f1,v⟩∀v∈H1(Ω), (2.24)  κ3(ξ,φ)+(divu,φ) =κ1(η,φ)∀φ∈L2(Ω), (2.25)  (ηt,ψ)dual+1μf(K(∇p −ρfg),∇ψ) =(ϕ,ψ)+⟨ϕ1,ψ⟩∀ψ∈H1(Ω), (2.26)  p:=κ1ξ+κ2η,q:=κ1η−κ3ξ, (2.27)  η(0)=η0: =c0p0+αq0, (2.28) where $$q_0:=\text{div}\,\, {\mathbf{u}}_0$$, $$u_0$$ and $$p_0$$ are same as in Definition 2.1. Remark 2.3 (a) We note that $$p$$ and $$q$$ are derivative variables in the system, after $${\mathbf{u}}$$, $$\xi$$ and $$\eta$$ are computed, $$p$$ and $$q$$ are simply updated by their algebraic expressions in (2.27). (b) Equation (2.26) implicitly imposes the following boundary condition for $$\eta$$:   κ2K∇η⋅n=Kρfg⋅n−κ1K∇ξ⋅n. (2.29) (c) By an Aubin–Lions Lemma (Temam, 1977), we conclude that $${\mathbf{u}}\in C^0([0,T];{\mathbf{L}}^2({\it{\Omega}}))$$, $$p\in C^0([0,T]; L^2({\it{\Omega}}))$$ and $$\eta\in C^0([0,T]; L^2({\it{\Omega}}))$$. Hence, the initial conditions in (2.23) and in (2.28) all make sense. (d) It should be pointed out that the only reason for introducing the space $${\mathbf{H}}_\perp^1({\it{\Omega}})$$ in the above two definitions is that the boundary condition (2.10) is a pure ‘Neumann condition’. If it is replaced by a pure Dirichlet condition or by a mixed Dirichlet–Neumann condition, there is no need to introduce this space. Neither do we need to impose the compatibility condition $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$$$\forall {\mathbf{v}}\in \mathbf{RM}$$. This fact will be used in our numerical experiments in Section 4. We also note that from the analysis point of view, the pure Neumann condition case is the most difficult case. Lemma 2.4 Every weak solution $$({\mathbf{u}},p)$$ of problems (2.21)–(2.23) satisfies the following energy law:   E(t)+1μf∫0t(K(∇p−ρfg),∇p)ds−∫0t(ϕ,p)ds −∫0t⟨ϕ1,p⟩,ds=E(0) (2.30) for all $$t\in [0,T]$$, where   E(t): =12[μ‖ε(u(t))‖L2(Ω)2+λ‖divu(t)‖L2(Ω)2+c0‖p(t)‖L2(Ω)2−2(f,u(t))−2⟨f1,u(t)⟩]. (2.31) Moreover,   ‖(c0p+αdivu)t‖L2(0.T;H1(Ω)′) ≤1μf‖K∇p−ρfg‖L2(ΩT)+‖ϕ‖L2(ΩT)+‖ϕ1‖L2(∂ΩT)<∞. (2.32) Likewise, weak solutions of (2.24)–(2.28) satisfy a similar energy law which is a rewritten version of (2.30) in the new variables. Lemma 2.5 Every weak solution $$({\mathbf{u}},\xi,\eta,p,q)$$ of problems (2.24)–(2.28) satisfies the following energy law:   J(t)+1μf∫0t(K(∇p−ρfg),∇p)ds−∫0t(ϕ,p)ds −∫0t⟨ϕ1,p⟩ds=J(0) (2.33) for all $$t\in [0,T]$$, where   J(t): =12[μ‖ε(u(t))‖L2(Ω)2+κ2‖η(t)‖L2(Ω)2+κ3‖ξ(t)‖L2(Ω)2−2(f,u(t))−2⟨f1,u(t)⟩]. (2.34) Moreover,   ‖ηt‖L2(0.T;H1(Ω)′) ≤1μf‖K∇p−ρfg‖L2(ΩT)+‖ϕ‖L2(ΩT)+‖ϕ1‖L2(∂ΩT)<∞. (2.35) We skip the proofs of the above two lemmas to save space and refer the interested reader to Feng et al. (xxxx) for their proofs. The above energy law immediately implies the following solution estimates. Lemma 2.6 There exists a positive constant $$C_1=C_1\bigl(\|{\mathbf{u}}_0\|_{H^1({\it{\Omega}})}, \|p_0\|_{L^2({\it{\Omega}})},$$$$\|{\mathbf{f}}\|_{L^2({\it{\Omega}})},\|{\mathbf{f}}_1\|_{L^2({\partial} {\it{\Omega}})},\|\phi\|_{L^2({\it{\Omega}})}, \|\phi_1\|_{L^2({\partial}{\it{\Omega}})} \bigr)$$ such that   μ‖ε(u)‖L∞(0,T;L2(Ω))+κ2‖η‖L∞(0,T;L2(Ω)) +κ3‖ξ‖L∞(0,T;L2(Ω))+K1μf‖∇p‖L2(0,T;L2(Ω))≤C1. (2.36)   ‖u‖L∞(0,T;L2(Ω))≤C1,‖p‖L∞(0,T;L2(Ω))≤C1(κ212+κ1κ3−12). (2.37)   ‖p‖L2(0,T;L2(Ω))≤C1,‖ξ‖L2(0,T;L2(Ω))≤C1κ1−1(1+κ212). (2.38) We note that (2.38) follows from (2.36), (2.3), the Poincaré inequality and (2.50) below, and the relation $$p=\kappa_1\xi +\kappa_2\eta$$. Furthermore, by exploiting the linearity of the PDE system, we have the following a priori estimates for the weak solution. Theorem 2.7 Suppose that $${\mathbf{u}}_0$$ and $$p_0$$ are sufficiently smooth, then there exist positive constants $$C_2=C_2\bigl(C_1,\|{\nabla} p_0\|_{L^2({\it{\Omega}})} \bigr)$$ and $$C_3=C_3\bigl(C_1,C_2, \|{\mathbf{u}}_0\|_{H^2({\it{\Omega}})},\|p_0\|_{H^2({\it{\Omega}})} \bigr)$$ such that there hold the following estimates for the solution to problems (2.15)–(2.17) and (2.10)–(2.12):   μ‖ε(ut)‖L2(0,T;L2(Ω))+κ2‖ηt‖L2(0,T;L2(Ω)) +κ3‖ξt‖L2(0,T;L2(Ω))+K1μf‖∇p‖L∞(0,T;L2(Ω))≤C2. (2.39)  μ‖ε(ut)‖L∞(0,T;L2(Ω))+κ2‖ηt‖L∞(0,T;L2(Ω)) +κ3‖ξt‖L∞(0,T;L2(Ω))+K1μf‖∇pt‖L2(0,T;L2(Ω))≤C3. (2.40)   ‖ηtt‖L2(H1(Ω)′)≤K2μfC3. (2.41) Proof. On noting that $${\mathbf{f}}, {\mathbf{f}}_1,\phi$$ and $$\phi_1$$ all are assumed to be independent of $$t$$, differentiating (2.24) and (2.25) with respect to $$t$$, taking $${\mathbf{v}}={\mathbf{u}}_t$$ and $$\varphi=\xi_t$$ in (2.24) and (2.25), respectively, and adding the resulting equations yield   μ‖ε(ut)‖L2(Ω)2=(qt,ξt)=κ1(ηt,ξt)−κ3‖ξt‖L2(Ω)2. (2.42) Setting $$\psi=p_t=\kappa_1\xi_t + \kappa_2\eta_t$$ in (2.26) gives   κ1(ηt,ξt)+κ2‖ηt‖L2(Ω)2+K2μfddt‖∇p−ρfg‖L2(Ω)2=ddt[(ϕ,p)+⟨ϕ1,p⟩]. (2.43) Adding (2.42) and (2.43) and integrating in $$t$$, we get for $$t\in[0,T]$$  K2μf‖∇p(t)−ρfg‖L2(Ω)2+∫0t[μ‖ε(ut)‖L2(Ω)2+κ2‖ηt‖L2(Ω)2+κ3‖ξt‖L2(Ω)2]ds =K2μf‖∇p0−ρfg‖L2(Ω)2+(ϕ,p(t)−p0)+⟨ϕ1,p(t)−p0⟩, which readily infers (2.39). To show (2.40), first differentiating (2.24) one time with respect to $$t$$ and setting $${\mathbf{v}}={\mathbf{u}}_{tt}$$, differentiating (2.25) twice with respect to $$t$$ and setting $$\varphi=\xi_t$$ and adding the resulting equations, we get   μ2ddt‖ε(ut)‖L2(Ω)2=(qtt,ξt)=κ1(ηtt,ξt)−κ32ddt‖ξt‖L2(Ω)2. (2.44) Secondly, differentiating (2.26) with respect $$t$$ one time and taking $$\psi=p_t=\kappa_1\xi_t + \kappa_2\eta_t$$ yield   κ1(ηtt,ξt)+κ22ddt‖ηt‖L2(Ω)2+Kμf‖∇pt‖L2(Ω)2=0. (2.45) Finally, adding the above two inequalities and integrating in $$t$$ give for $$t\in [0,T]$$  μ‖ε(ut(t))‖L2(Ω)2+κ2‖ηt(t)‖L2(Ω)2+κ3‖ξt(t)‖L2(Ω)2+2Kμf∫0t‖∇pt‖L2(Ω)2ds =μ‖ε(ut(0))‖L2(Ω)2+κ2‖ηt(0)‖L2(Ω)2+κ3‖ξt(0)‖L2(Ω)2. (2.46) Hence, (2.40) holds. Equation (2.41) follows immediately from the following inequality   (ηtt,ψ)=−1μf(K∇p,∇ψ)≤Kμf‖∇p‖L2(Ω)‖∇ψ‖L2(Ω)∀ψ∈H01(Ω), (2.40) and the definition of the $$H^{1}({\it{\Omega}})'$$-norm. The proof is complete. □ Remark 2.8 As expected, the above high order norm solution estimates require $$p_0\in H^1({\it{\Omega}}), {\mathbf{u}}_t(0)\in {\mathbf{L}}^2({\it{\Omega}}), \eta_t(0)\in L^2({\it{\Omega}})$$ and $$\xi_t(0)\in L^2({\it{\Omega}})$$. The values of $${\mathbf{u}}_t(0), \eta_t(0)$$ and $$\xi_t(0)$$ can be computed using the PDEs as follows. It follows from (2.17) that $$\eta_t(0)$$ satisfies   ηt(0)=ϕ+1μfdiv[K(∇p0−ρfg)]. Hence, $$\eta_t(0)\in L^2({\it{\Omega}})$$ provided that $$p_0\in H^2({\it{\Omega}})$$. To find $${\mathbf{u}}_t(0)$$ and $$\xi_t(0)$$, differentiating (2.15) and (2.16) with respect to $$t$$ and setting $$t=0$$, we get   −μdivε(ut(0))+∇ξt(0) =0in Ω,κ3ξt(0)+divut(0) =κ1ηt(0)in Ω. Hence, $${\mathbf{u}}_t(0)$$ and $$\xi_t(0)$$ can be determined by solving the above generalized Stokes problem. The next lemma shows that weak solutions of problems (2.24)–(2.28) preserve some ‘invariant’ quantities, it turns out that these ‘invariant’ quantities play a vital role in the construction of our time-splitting scheme to be introduced in the next section. Lemma 2.9 Every weak solution $$({\mathbf{u}},\xi,\eta,p, q)$$ to problems (2.24)–(2.28) satisfies the following relations:   Cη(t):=(η(⋅,t),1)=(η0,1)+[(ϕ,1)+⟨ϕ1,1⟩]tt≥0. (2.47)  Cξ(t):=(ξ(⋅,t),1)=1d+μκ3[μκ1Cη(t)−(f,x)−⟨f1,x⟩]. (2.48)  Cq(t):=(q(⋅,t),1)=κ1Cη(t)−κ3Cξ(t). (2.49)  Cp(t):=(p(⋅,t),1)=κ1Cξ(t)+κ2Cη(t). (2.50)  Cu(t):=⟨u(⋅,t)⋅n,1⟩=Cq(t). (2.51) We skip the proof of the above lemma again to save space and refer the interested reader to Feng et al. (xxxx) for its proof. Remark 2.10 We note that $$C_\eta, C_\xi, C_q$$ and $$C_p$$ all are (known) linear functions of $$t$$, and they become (known) constants when $$\phi\equiv 0$$ and $$\phi_1\equiv 0$$. With the help of the above lemmas, we can show the solvability of problems (2.6)–(2.12). Theorem 2.11 Let $${\mathbf{u}}_0\in{\mathbf{H}}^1({\it{\Omega}}), {\mathbf{f}}\in{\mathbf{L}}^2({\it{\Omega}}), {\mathbf{f}}_1\in {\mathbf{L}}^2({\partial}{\it{\Omega}}), p_0\in L^2({\it{\Omega}}), \phi\in L^2({\it{\Omega}})$$ and $$\phi_1\in L^2({\partial}{\it{\Omega}})$$. Suppose $$c_0>0$$ and $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$ for any $${\mathbf{v}}\in \mathbf{RM}$$. Then there exists a unique solution to problems (2.6)–(2.12) in the sense of Definition 2.1, likewise, there exists a unique solution to problems (2.15)–(2.17) and (2.18)–(2.20) in the sense of Definition 2.2. Proof. We only outline the main steps of the proof and leave the details to the interested reader. First, Because the PDE system is linear, the existence of weak solution can be proved by the standard Galerkin method and compactness argument (cf. Temam, 1977). We note that the energy laws established in Lemmas 2.4 and 2.5 guarantee the required uniform estimates for the Galerkin approximate solutions. Secondly, to show the uniqueness, suppose there are two sets of weak solutions, again by the linearity of the PDE system it is trivial to show that the difference of the solutions satisfy the same PDE system with zero initial and boundary data. The energy law immediately implies that the difference must be zero, and hence the uniqueness is verified. □ We conclude this section by establishing a convergence result for the solution of problems (2.15)–(2.17) and (2.18)–(2.20) when the constrained specific storage coefficient $$c_0$$ tends to $$0$$. Such a convergence result is useful and significant for the following two reasons. First, as mentioned earlier, the poroelasticity model studied in this article reduces to the well-known Biot’s consolidation model from soil mechanics (cf. Murad & Loula, 1992; Phillips & Wheeler, 2007a) and Doi’s model for polymer gels (cf. Yamaue & Doi, 2004; Feng & He, 2010). Secondly, it proves that the proposed approach and methods of this article are robust under such a limit process. Theorem 2.12 Let $${\mathbf{u}}_0\in{\mathbf{H}}^1({\it{\Omega}}), {\mathbf{f}}\in{\mathbf{L}}^2({\it{\Omega}}), {\mathbf{f}}_1\in {\mathbf{L}}^2({\partial}{\it{\Omega}}), p_0\in L^2({\it{\Omega}}), \phi\in L^2({\it{\Omega}})$$ and $$\phi_1\in L^2({\partial}{\it{\Omega}})$$. Suppose $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$ for any $${\mathbf{v}}\in \mathbf{RM}$$. Let $$({\mathbf{u}}_{c_0},\eta_{c_0},\xi_{c_0},p_{c_0},q_{c_0})$$ denote the unique weak solution to problems (2.15)–(2.17) and (2.18)–(2.20). Then there exists $$({\mathbf{u}}_*, \eta_*,\xi_*, p_*, q_*)\in {\mathbf{L}}^\infty(0,T;{\mathbf{H}}^1_\perp({\it{\Omega}}))\times L^\infty(0,T;L^2({\it{\Omega}}))\times L^2(0,T; L^2({\it{\Omega}}))\times L^2(0,T;H^1({\it{\Omega}})\times L^\infty(0,T;L^2({\it{\Omega}})))$$ such that $$({\mathbf{u}}_{c_0},\eta_{c_0},\xi_{c_0},p_{c_0},q_{c_0})$$ converges weakly to $$({\mathbf{u}}_*, \eta_*,\xi_*, p_*, q_*)$$ in the above product space as $$c_0\to 0$$. Proof. It follows immediately from (2.35)–(2.38) and Korn’s inequality that $${\mathbf{u}}_{c_0}$$ is uniformly bounded (in $$c_0$$) in $${\mathbf{L}}^\infty(0,T; {\mathbf{H}}^1_\perp({\it{\Omega}}))$$. $$\sqrt{\kappa_2} \eta_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^\infty(0,T; L^2({\it{\Omega}}))\cap H^1(0,T; H^{1}({\it{\Omega}})')$$. $$\sqrt{\kappa_3}\xi_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^\infty(0,T; L^2({\it{\Omega}}))$$. $$\xi_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^2({\it{\Omega}}_T)$$. $$p_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^2(0,T; H^1({\it{\Omega}}))$$. $$q_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^\infty(0,T; L^2({\it{\Omega}}))$$. On noting that $$\lim_{c_0\to 0}\kappa_1=\frac{1}{\alpha}$$, $$\lim_{c_0\to 0}\kappa_2=\frac{{\lambda}}{\alpha^2}$$ and $$\lim_{c_0\to 0}\kappa_3=0$$, by the weak compactness of reflexive Banach spaces and Aubin–Lions Lemma (Dautray & Lions, 1990), we have that there exist $$({\mathbf{u}}_*, \eta_*,\xi_*, p_*, q_*)\in {\mathbf{L}}^\infty(0,T;{\mathbf{H}}^1_\perp({\it{\Omega}})) \times L^\infty(0,T;L^2({\it{\Omega}}))\times L^2(0,T; L^2({\it{\Omega}}))\times L^2(0,T;H^1({\it{\Omega}}))\times L^\infty(0,T;L^2({\it{\Omega}}))$$ and a subsequence of $$({\mathbf{u}}_{c_0},\eta_{c_0},\xi_{c_0},$$$$p_{c_0},q_{c_0})$$ (still denoted by the same notation) such that as $$c_0\to 0$$ (a subsequence of $$c_0$$, to be exact) $${\mathbf{u}}_{c_0}$$ converges to $${\mathbf{u}}_*$$ weak $$*$$ in $${\mathbf{L}}^\infty(0,T; {\mathbf{H}}^1_\perp({\it{\Omega}}))$$ and weakly in $${\mathbf{L}}^2(0,T; {\mathbf{H}}^1_\perp({\it{\Omega}}))$$. $$\sqrt{\kappa_2} \eta_{c_0}$$ converges to $$\frac{\sqrt{{\lambda}}}{\alpha} \eta_*$$ weak $$*$$ in $$L^\infty(0,T; L^2({\it{\Omega}}))$$ and strongly in $$L^2({\it{\Omega}}_T)$$. $$\kappa_3\xi_{c_0}$$ converges to $$0$$ strongly in $$L^2({\it{\Omega}}_T)$$. $$\xi_{c_0}$$ converges to $$\xi_*$$ weakly in $$L^2({\it{\Omega}}_T)$$. $$p_{c_0}$$ converges to $$p_*$$ weakly in $$L^2(0,T; H^1({\it{\Omega}}))$$. $$q_{c_0}$$ converges to $$p_*$$ weak $$*$$ in $$L^\infty(0,T; L^2({\it{\Omega}}))$$ and weakly in $$L^2({\it{\Omega}}_T)$$. Then setting $$c_0\to 0$$ in (2.24)–(2.28) yields (note that the dependence of the solution on $$c_0$$ is suppressed there)   2μ(ε(u∗),ε(v))−(ξ∗,divv) =(f,v)+⟨f1,v⟩∀v∈H1(Ω),(divu∗,φ) =1α(η∗,φ)∀φ∈L2(Ω),((η∗)t,ψ)dual+1μf(K(∇p∗ −ρfg),∇ψ) =(ϕ,ψ)+⟨ϕ1,ψ⟩∀ψ∈H1(Ω),p∗:=1αξ∗+λα2η∗,q∗:=1αη∗in L2(ΩT),u∗(0)=u0,η∗(0)=η0: =αq0. Equivalently,   2μ(ε(u∗),ε(v))−(ξ∗,divv) =(f,v)+⟨f1,v⟩∀v∈H1(Ω),(divu∗,φ) =(q∗,φ)∀φ∈L2(Ω),α((q∗)t,ψ)dual+1μf(K(∇(λα−1q∗+α−1ξ∗) −ρfg),∇ψ) =(ϕ,ψ)+⟨ϕ1,ψ⟩∀ψ∈H1(Ω),p∗ :=1α(ξ∗+λq∗)in L2(ΩT),q∗(0)=q0 :=divu0. Hence, $$({\mathbf{u}}_*, \eta_*,\xi_*,p_*,q_*)$$ is a weak solution of Biot’s consolidation model (cf. Yamaue & Doi, 2004; Feng & He, 2010). By the uniqueness of its solutions, we conclude that the whole sequence $$({\mathbf{u}}_{c_0},\eta_{c_0},\xi_{c_0},$$$$p_{c_0},q_{c_0})$$ converges to $$({\mathbf{u}}_*, \eta_*,\xi_*,p_*,q_*)$$ as $$c_0\to 0$$ in the above sense. The proof is complete. □ 3. Fully discrete finite element methods The goal of this section is to design and analyse some fully discrete finite element methods for the poroelasticity model based on the above new formulation. 3.1 Formulation of fully discrete finite element methods In this section, we consider the space–time discretization that combines the above splitting algorithm with appropriately chosen spatial discretization methods. To this end, we introduce some notation. Assume $${\it{\Omega}}\in\mathbb{R}^d (d=2, 3)$$ is a polygonal domain. Let $$\mathcal{T}_h$$ be a quasi-uniform triangulation or rectangular partition of $${\it{\Omega}}$$ with mesh size $$h$$ and $$\bar{{\it{\Omega}}}=\bigcup_{K\in\mathcal{T}_h}\bar{K}$$. Also, let $$({\mathbf{X}}_h, M_h)$$ be a stable mixed finite element pair, that is, $${\mathbf{X}}_h\subset {\mathbf{H}}^1({\it{\Omega}})$$ and $$M_h\subset L^2({\it{\Omega}})$$ satisfy the inf–sup condition   supvh∈Xh(divvh,φh)‖vh‖H1(Ω)≥β0‖φh‖L2(Ω)∀φh∈M0h:=Mh∩L02(Ω), β0>0. (3.1) A number of stable mixed finite element spaces $$({\mathbf{X}}_h, M_h)$$ have been known in the literature (Brezzi & Fortin, 1992). A well-known example is the following so-called Taylor–Hood element (cf. Bercovier & Pironneau, 1979; Brezzi & Fortin, 1992; Roberts & Thomas, 1991):   Xh ={vh∈C0(Ω¯);vh|K∈P2(K)∀K∈Th},Mh ={φh∈C0(Ω¯);φh|K∈P1(K)∀K∈Th}. In the next subsection, we shall only present the analysis for the Taylor–Hood element, we remark that the analysis can be extended to other stable mixed elements. However, piecewise constant space $$M_h$$ is not recommended because that would result in no rate of convergence for the approximation of the pressure $$p$$ (see Section 3.3). Finite element approximation space $$W_h$$ for $$\eta$$ variable can be chosen independently, any piecewise polynomial space is acceptable provided that $$W_h \supset M_h$$. Especially, $$W_h\subset L^2({\it{\Omega}})$$ can be chosen as a fully discontinuous piecewise polynomial space, although it is more convenient to choose $$W_h$$ to be a continuous (respectively, discontinuous) space if $$M_h$$ is a continuous (respectively, discontinuous) space. The most convenient choice is $$W_h =M_h$$, which will be adopted in the remainder of this article. Recall that $${\mathbf{RM}}$$ denotes the space of the infinitesimal rigid motions (see Section 2), evidently, $${\mathbf{RM}}\subset {\mathbf{X}}_h$$. We now introduce the $$L^2$$-projection $$\mathcal{P}_R$$ from $${\mathbf{L}}^2({\it{\Omega}})$$ to $${\mathbf{RM}}$$. For each $${\mathbf{v}}\in {\mathbf{L}}^2({\it{\Omega}})$$, $$\mathcal{P}_R{\mathbf{v}}_h\in {\mathbf{RM}}$$ is defined by   (PRvh,r)=(vh,r)∀r∈RM. Moreover, we define   Vh:=(I−PR)Xh={vh∈Xh;(vh,r)=0∀r∈RM}. (3.2) It is easy to check that $${\mathbf{X}}_h={\mathbf{V}}_h\bigoplus {\mathbf{RM}}$$. It was proved in Feng & He (2010) that there holds the following alternative version of the above inf–sup condition:   supvh∈Vh(divvh,φh)‖vh‖H1(Ω)≥β1‖φh‖L2(Ω)∀φh∈M0h,β1>0. (3.3) Finite element algorithm (FEA): (i) Compute $${\mathbf{u}}^0_h\in {\mathbf{V}}_h$$ and $$q^0_h\in W_h$$ by   uh0 =Rhu0,ph0=Qhp0,qh0=Qhq0 (q0=divu0),ηh0 =c0ph0+αqh0,ξh0=αph0−λqh0, where $$\mathcal{Q}_h$$ denotes the $$L^2$$-projection operator (see (3.23)). (ii) For $$n=0,1,2, \cdots$$, do the following two steps. Step 1: Solve for $$({\mathbf{u}}^{n+1}_h,\xi^{n+1}_h, \eta^{n+1}_h)\in {\mathbf{V}}_h\times M_h \times W_h$$ such that   μ(ε(uhn+1),ε(vh))−(ξhn+1,divvh)=(f,vh)+⟨f1,vh⟩∀vh∈Vh, (3.4)  κ3(ξhn+1,φh)+(divuhn+1,φh)=κ1(ηhn+θ,φh)∀φh∈Mh; (3.5)   (dtηhn+1,ψh)+1μf(K(∇(κ1ξhn+1+κ2ηhn+1) −ρfg,∇ψh)=(ϕ,ψh)+⟨ϕ1,ψh⟩∀ψh∈Wh, (3.6) where $$\theta=0$$ or $$1$$. Step 2: Update $$p^{n+1}_h$$ and $$q^{n+1}_h$$ by   phn+1 =κ1ξhn+1+κ2ηhn+θ,qhn+1=κ1ηhn+1−κ3ξhn+1. (3.7) Remark 3.1 (a) At each time step, problems (3.4)–(3.5) is a generalized Stokes problem with a mixed boundary condition for $$({\mathbf{u}}_h^{n+1},\xi_h^{n+1})$$. When $$\theta=0$$, the well posedness of the generalized Stokes problem follows from the fact that (3.4)–(3.5) is not a saddle point problem; instead, it results in a symmetric and positive definite linear system because $$\kappa_3>0$$. When $$\theta=1$$, the well posedness of (3.4)–(3.6) is an immediate consequence of the discrete energy law given in Lemma 3.3 below. (b) When $$\theta=0$$, Step 1 consists of two decoupled sub-problems that can be solved independently. On the other hand, when $$\theta=1$$, it is a coupled system, all three unknown functions must be solved together. 3.2 Stability analysis of fully discrete finite element methods The primary goal of this subsection is to derive a discrete energy law, which mimics the PDE energy law (2.30). It turns out that such a discrete energy law only holds if $$h$$ and $${\it{\Delta}} t$$ satisfy the mesh constraint $${\it{\Delta}} t=O(h^2)$$ when $$\theta=0$$, but for all $$h,{\it{\Delta}} t>0$$ when $$\theta=1$$. Before discussing the stability of (FEA), we first show that the numerical solution satisfies all side constraints which are fulfilled by the PDE solution. Lemma 3.2 Let $$\{({\mathbf{u}}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0}$$ be defined by the (FEA), then there hold   (ηhn,1) =Cη(tn)for n=0,1,2,⋯, (3.8)  (ξhn,1) =Cξ(tn−1+θ)for n=1−θ,1,2,⋯, (3.9)  ⟨uhn⋅n,1⟩ =Cu(tn−1+θ)for n=1−θ,1,2,⋯. (3.10) We skip the proof of the above lemma again to save space and refer the interested reader to Feng et al. (xxxx) for its proof. The next lemma establishes an identity that mimics the continuous energy law for the solution of (FEA). Lemma 3.3 Let $$\{({\mathbf{u}}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0}$$ be defined by (FEA), then there holds the following identity:   Jh,θℓ+Sh,θℓ=Jh,θ0for ℓ≥1,θ=0,1, (3.11) where   Jh,θℓ:=12[μ‖ε(uhℓ+1)‖L2(Ω)2+κ2‖ηhℓ+θ‖L2(Ω)2+κ3‖ξhℓ+1‖L2(Ω)2−2(f,uhℓ+1)−2⟨f1,uhℓ+1⟩],Sh,θℓ:=Δt∑n=1ℓ[μΔt2‖dtε(uhn+1)‖L2(Ω)2+1μf(K∇phn+1−Kρfg,∇phn+1) +κ2Δt2‖dtηhn+θ‖L2(Ω)2+κ3Δt2‖dtξhn+1‖L2(Ω)2−(ϕ,phn+1)−⟨ϕ1,phn+1⟩ −(1−θ)κ1Δtμf(Kdt∇ξhn+1,∇phn+1)].phn+1:=κ1ξhn+1+κ2ηhn+θ. Proof. Since the proof for the case $$\theta=1$$ is exactly the same as that of the PDE energy law, so we omit it and leave it to the interested reader to explore. Here, we only consider the case $$\theta=0$$. Based on (3.5), we can define $$\eta_h^{-1}$$ by   κ1(ηh−1,φh)=κ3(ξh0,φh)+(divuh0,φh). (3.12) Setting $${\mathbf{v}}_h=d_t{\mathbf{u}}_h^{n+1}$$ in (3.4), $$\varphi_h=\xi_h^{n+1}$$ in (3.5) and $$\psi_h= p_h^{n+1}$$ in (3.6) after lowing the super-index from $$n+1$$ to $$n$$ on the both sides of (3.6), we get   μ2dt‖ε(uhn+1)‖L2(Ω)2+μ2Δt‖dtε(uhn+1)‖L2(Ω)2 =dt(f,uhn+1)+dt⟨f1,uhn+1⟩+(ξhn+1,divdtuhn+1), (3.13)  κ3(dtξhn+1,ξhn+1)+(divdtuhn+1,ξhn+1)=κ1(dtηhn,ξhn+1), (3.14)   (dtηhn,phn+1)+1μf(K(∇(κ1ξhn+κ2ηhn)−ρfg),∇phn+1) =(ϕ,phn+1)+⟨ϕ1,phn+1⟩. (3.15) The first term on the left-hand side of (3.15) can be rewritten as   (dtηhn,phn+1) =(dtηhn,κ1ξhn+1+κ2ηhn) =κ1(dtηhn,ξhn+1)+κ2Δt2‖dtηhn‖L2(Ω)2+κ22dt‖ηhn‖L2(Ω)2. (3.16) Moreover,   1μf(K∇(κ1ξhn+κ2ηhn)−Kρfg,∇phn+1) =1μf(K∇phn+1−Kρfg,∇phn+1)−κ1Δtμf(Kdt∇ξhn+1,∇phn+1). (3.17)  κ3(dtξhn+1,ξhn+1)=κ32dt‖ξhn+1‖L2(Ω)2+κ3Δt2‖dtξhn+1‖L2(Ω)2. (3.18) Adding (3.13)–(3.15), using (3.16)–(3.18) and applying the summation operator $${\it{\Delta}} t\sum_{n=1}^{\ell}$$ to the both sides of the resulting equation yield the desired equality (3.11). The proof is complete. □ In the case $$\theta=1$$, (3.11) gives the desired solution estimates without any mesh constraint. On the other hand, when $$\theta=0$$, since the last term in the expression of $$S^\ell_{h,\theta}$$ does not have a fixed sign, and hence it needs to be controlled to ensure the positivity of $$S^\ell_{h,\theta}$$. Corollary 3.4 Let $$\{({\mathbf{u}}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0}$$ be defined by (FEA) with $$\theta=0$$, then there holds the following inequality:   Jh,0ℓ+S^h,0ℓ≤Jh,00for ℓ≥1, (3.19) provided that $${\it{\Delta}} t=O(h^2)$$. Here,   S^h,0ℓ:=Δt∑n=1ℓ[μΔt4‖dtε(uhn+1)‖L2(Ω)2+12μf‖∇phn+1‖L2(Ω)2−1μf(ρfg,∇phn+1)+κ2Δt2‖dtηhn‖L2(Ω)2+κ3Δt2‖dtξhn+1‖L2(Ω)2−(ϕ,phn+1)−⟨ϕ1,phn+1⟩]. Proof. By Schwarz inequality and inverse inequality (3.22), we get   κ1K1Δtμf(dt∇ξhn+1,∇phn+1) ≤κ12K22μf‖K12∇ξhn+1−K12∇ξhn‖L2(Ω)2 +12μfK2‖K12∇phn+1‖L2(Ω)2 ≤K22κ122μfh2‖ξhn+1−ξhn‖L2(Ω)2+12μf‖∇phn+1‖L2(Ω)2. (3.20) To bound the first term on the right-hand side of (3.20), we appeal to the inf–sup condition and get   ‖ξhn+1−ξhn‖L2 ≤1β1supvh∈Vh(divvh,ξhn+1−ξhn)‖∇vh‖L2(Ω) ≤μβ1supvh∈Vh(ε(un+1−un),ε(vh))‖∇vh‖L2(Ω) ≤μβ1Δt‖dtε(uhn+1)‖L2. (3.21) Substituting (3.21) into (3.20) and combining it with (3.11) imply (3.19) provided that $${\it{\Delta}} t \leq C(\mu_f\beta_1^2)(2\mu c_1^2\kappa_1^2)^{-1} h^2$$. The proof is complete. □ 3.3 Convergence analysis The goal of this section is to analyse the fully discrete finite element algorithm proposed in the previous subsection. Precisely, we shall derive optimal order error estimates for (FEA) in both $$L^\infty(0,T;L^2({\it{\Omega}}))$$ and $$L^2(0,T; H^1({\it{\Omega}}))$$-norm. To this end, we first list some facts, which are well known in the literature (Brezzi & Fortin, 1992; Brenner & Scott, 2008) about finite element functions. We first recall the following inverse inequality for polynomial functions Ciarlet (1978):   ‖∇φh‖L2(K)≤c1h−1‖φh‖L2(K)∀φh∈Pr(K),K∈Th. (3.22) For any $$\varphi\in L^2({\it{\Omega}})$$, we define its $$L^2$$-projection $$\mathcal{Q}_h: L^2\rightarrow W_h$$ as   (Qhφ,ψh)=(φ,ψh)ψh∈Wh. (3.23) It is well known that the projection operator $$\mathcal{Q}_h: L^2\rightarrow W_h$$ satisfies (cf. Brenner & Scott, 2008), for any $$\varphi\in H^s({\it{\Omega}}) (s\geq1)$$,   ‖Qhφ−φ‖L2(Ω)+h‖∇(Qhφ−φ)‖L2(Ω)≤Chℓ‖φ‖Hℓ(Ω)ℓ=min{2,s}. (3.24) We like to point out that when $$W_h\not\subset H^1({\it{\Omega}})$$, the second term on the left-hand side of (3.24) has to be replaced by the broken $$H^1$$-norm. Next, for any $$\varphi\in H^1({\it{\Omega}})$$, we define its elliptic projection $$\mathcal{S}_h\varphi$$ by   (K∇Shφ,∇φh) =(K∇φ,∇φh)∀φh∈Wh, (3.25)  (Shφ,1) =(φ,1). (3.26) It is well known that the projection operator $$\mathcal{S}_h: H^1({\it{\Omega}})\rightarrow W_h$$ satisfies (cf. Brenner & Scott, 2008), for any $$\varphi\in H^s({\it{\Omega}}) (s>1)$$,   ‖Shφ−φ‖L2(Ω)+h‖∇(Shφ−φ)‖L2(Ω)≤Chℓ‖φ‖Hℓ(Ω)ℓ=min{2,s}. (3.27) Finally, for any $${\mathbf{v}}\in {\mathbf{H}}^1_\perp({\it{\Omega}})$$, we define its elliptic projection $$\mathcal{R}_h{\mathbf{v}}$$ by   (ε(Rhv),ε(wh))=(ε(v),ε(wh))wh∈Vh. (3.28) It is easy to show that the projection $$\mathcal{R}_h{\mathbf{v}}$$ satisfies (cf. Brenner & Scott, 2008), for any $${\mathbf{v}}\in {\mathbf{H}}^1_\perp({\it{\Omega}})\cap {\mathbf{H}}^s({\it{\Omega}}) (s>1)$$,   ‖Rhv−v‖L2(Ω)+h‖∇(Rhv−v)‖L2(Ω)≤Chm‖v‖Hm(Ω)m=min{3,s}. (3.29) To derive error estimates, we introduce the following error notation   Eun :=u(tn)−uhn,Eξn:=ξ(tn)−ξhn,Eηn:=η(tn)−ηhn,Epn :=p(tn)−phn,Eqn:=q(tn)−qhn. It is easy to check that   Epn=κ1Eξn+κ2Eηn,Eqn=κ3Eξn+κ1Eηn. (3.30) Also, we denote   Eun =u(tn)−Rh(u(tn))+Rh(u(tn))−uhn:=Λun+Θun,Eξn =ξ(tn)−Sh(ξ(tn))+Sh(ξ(tn))−ξhn:=Λξn+Θξn,Eηn =η(tn)−Sh(η(tn))+Sh(η(tn))−ηhn:=Ληn+Θηn,Epn =p(tn)−Qh(p(tn))+Qh(p(tn))−phn:=Λpn+Θpn. Lemma 3.5 Let $$\{ ({\mathbf{u}}_h^n, \xi_h^n, \eta_h^n) \}_{n\geq0}$$ be generated by the (FEA) and $${\it{\Lambda}}_{{\mathbf{u}}}^n, {\it{\Theta}}_{{\mathbf{u}}}^n, {\it{\Lambda}}_{\xi }^n, {\it{\Theta}}_{\xi }^n, {\it{\Lambda}}_{\eta }^n$$ and $${\it{\Theta}}_{\eta }^n$$ be defined as above. Then there holds the following identity:   Ehℓ+Δt∑n=1ℓ[1μf(K∇Θ^pn+1−Kρfg,Θ^pn+1) +Δt2(μ‖dtε(Θun+1)‖L2(Ω)2+κ2‖dtΘηn+θ‖L2(Ω)2+κ3‖dtΘξn+1‖L2(Ω)2)] =Eh0+Δt∑n=1ℓ[(Λξn+1,divdtΘun+1)−(divdtΛun+1,Θξn+1)] +(Δt)2∑n=1ℓ(dt2ηh(tn+1),Θξn+1)+Δt∑n=1ℓ(Rhn+1,Θ^pn+1) +(1−θ)(Δt)2∑n=1ℓκ1μf(Kdt∇Θξn+1,∇Θ^pn+1), (3.31) where   Θ^pn+1 :=κ1∇Θξn+1+κ2∇Θηn+θ, (3.32)  Ehℓ :=12[μ‖ε(Θuℓ+1)‖L2(Ω)2+κ2‖Θηℓ+θ‖L2(Ω)2+κ3‖Θξℓ+1‖L2(Ω)2], (3.33)  Rhn+1 :=−1Δt∫tntn+1(s−tn)ηtt(s)ds. (3.34) Proof. Subtracting (3.4) from (2.24), (3.5) from (2.25), (3.6) from (2.26), respectively, we get the following error equations:   μ(ε(Eun+1),ε(vh))−(Eξn+1,divvh)=0∀vh∈Vh, (3.35)  κ3(Eξn+1,φh)+(divEun+1,φh)  =κ1(Eηn+θ,φh)+Δt(dtη(tn+1),φh)∀φh∈Mh, (3.36)   (dtEηn+1,ψh)+1μf(K∇(κ1Eξn+1+κ2Eηn+1)−Kρfg,∇ψh) =(Rhn+1,ψh)∀ψh∈Wh, (3.37)  Eu0=0,Eξ0=0,Eη−1=0. (3.38) Using the definition of the projection operators $$\mathcal{Q}_h, \mathcal{S}_h, \mathcal{R}_h$$, we have   μ(ε(Θun+1),ε(vh))−(Θξn+1,divvh)=(Λξn+1,divvh)∀vh∈Vh, (3.39)  κ3(Θξn+1,φh)+(divΘun+1,φh)=κ1(Θηn+θ,φh) −(divΛun+1,φh)+Δt(dtη(tn+1),φh)∀φh∈Mh, (3.40)   (dtΘηn+1,ψh)+1μf(K∇(κ1Θξn+1+κ2Θηn+1)−Kρfg,∇ψh) =(Rhn+1,ψh)∀ψh∈Wh, (3.41)  Eu0=0,Eξ0=0,Eη−1=0. (3.42) Equation (3.31) follows from setting $${\mathbf{v}}_h=d_t{\it{\Theta}}^{n+1}_{{\mathbf{u}}}$$ in (3.39), $$\varphi_h={\it{\Theta}}^{n+1}_\xi$$ (after applying the difference operator $$d_t$$ to the equation (3.40)), $$\psi_h=\hat{{\it{\Theta}}}_p^{n+1}=\kappa_1{\it{\Theta}}_\xi^{n+1} +\kappa_2{\it{\Theta}}_\eta^{n+\theta}$$ in (3.41), adding the resulting equations and applying the summation operator $${\it{\Delta}} t\sum^{\ell}_{n=1}$$ to both sides. □ Theorem 3.6 Let $$\{(u_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0}$$ be defined by (FEA), then there holds the error estimate for $$\ell \leq N$$   max0≤n≤ℓ[μ‖ε(Θun+1)‖L2(Ω)+κ2‖Θηn+θ‖L2(Ω)+κ3‖Θξn+1‖L2(Ω)] +[Δt∑n=0ℓ1μf‖Θ^pn+1‖L22]12≤C1(T)Δt+C2(T)h2, (3.43) provided that $${\it{\Delta}} t =O(h^2)$$ when $$\theta=0$$ and $${\it{\Delta}} t >0$$ when $$\theta=1$$. Here   C1(T) =C‖ηt‖L2((0,T);L2(Ω))2+C‖(η)tt‖L2((0,T);H1(Ω)′), (3.44)  C2(T) =C‖ξ‖L∞((0,T);H2(Ω))+C‖ξt‖L2((0,T);H2(Ω)) +C‖div(u)t‖L2((0,T);H2(Ω)). (3.45) Proof. To derive the above inequality, we need to bound each term on the right-hand side of (3.31). Using the fact $${\it{\Theta}}_{{\mathbf{u}}}^0=\mathbf{0}$$, $${\it{\Theta}}_\xi^0=0$$ and $${\it{\Theta}}_\eta^{-1}=0$$, we have   Ehℓ+Δt∑n=1ℓ[1μf(K∇Θ^pn+1−Kρfg,∇Θ^pn+1) +Δt2(μ‖dtε(Θun+1)‖L2(Ω)2+κ2‖dtΘηn+θ‖L2(Ω)2+κ3‖dtΘξn+1‖L2(Ω)2)] =(Δt)2∑n=1ℓ(dt2η(tn+1),Θξn+1)+Δt∑n=1ℓ(Rhn+1,Θ^pn+1)+μ2‖ε(Θu1)‖L2(Ω)2 +Δt∑n=1ℓ[(Λξn+1,divdtΘun+1)−(divdtΛun+1,Θξn+1)] +(1−θ)(Δt)2∑n=1ℓκ1μf(Kdt∇Θξn+1,∇Θ^pn+1). (3.46) We now estimate each term on the right-hand side of (3.46). To bound the first term on the right-hand side of (3.46), we first use the summation by parts formula and $$d_t\eta_h(t_0)=0$$ to get   ∑n=0ℓ(dt2η(tn+1),Θξn+1)=1Δt(dtη(tl+1),Θξl+1)−∑n=1ℓ(dtη(tn),dtΘξn+1). (3.47) Now, we bound each term on the right-hand side of (3.47) as follows:   1Δt(dtη(tℓ+1),Θξℓ+1)≤1Δt‖dtη(tℓ+1)‖L2(Ω)‖Θξℓ+1‖L2(Ω) ≤1Δt‖ηt‖L2((tℓ,tℓ+1);L2(Ω))⋅1β1supvh∈Vhμ(ε(Θuℓ+1),ε(vh))−(Λξℓ+1,divvh)‖∇vh‖L2(Ω) ≤Cμβ1Δt‖ηt‖L2((tℓ,tℓ+1);L2(Ω))[‖ε(Θuℓ+1)‖L2(Ω)+1μ‖Λξℓ+1‖L2(Ω)] ≤μ4(Δt)2‖ε(Θuℓ+1)‖L2(Ω)2+Cμβ12‖ηt‖L2((tℓ,tℓ+1);L2(Ω))2+Cβ12‖Λξℓ+1‖L2(Ω)2, (3.48)   ∑n=1ℓ(dtη(tn),dtΘξn+1)≤∑n=1ℓ‖dtη(tn)‖L2(Ω)‖dtΘξn+1‖L2(Ω) ≤∑n=1ℓ‖dtη(tn)‖L2(Ω)⋅1β1supvh∈Vhμ(dtε(Θun+1),ε(vh))−(dtΛξn+1,divvh)‖∇vh‖L2(Ω) ≤Cμβ1∑n=1ℓ‖dtη(tn)‖L2(Ω)[‖dtε(Θun+1)‖L2(Ω)+1μ‖dtΛξn+1‖L2(Ω)] ≤∑n=1ℓ[μ4‖dtε(Θun+1)‖L2(Ω)2+Cβ12‖dtΛξn+1‖L2(Ω)]+Cμβ12‖ηt‖L2((0,T);L2(Ω))2. (3.49) The second term on the right-hand side of (3.46) can be bounded as   |(Rhn+1,Θ^pn+1)| ≤‖Rhn+1‖H1(Ω)′‖∇Θ^pn+1‖L2(Ω) ≤K14μf‖∇Θ^pn+1‖L2(Ω)2+μfK1‖Rhn+1‖H1(Ω)′2 ≤K14μf‖∇Θ^pn+1‖L2(Ω)2+μfΔt3K1‖ηtt‖L2((tn,tn+1);H1(Ω)′)2, (3.50) where we have used the fact that   ‖Rhn+1‖H1(Ω)′2≤Δt3∫tntn+1‖ηtt‖H1(Ω)′2dt. The fourth term on the right-hand side of (3.46) can be bounded by   Δt∑n=1ℓ[(Λξn+1,divdtΘun+1)−(divdtΛun+1,Θξn+1)] =(Λξℓ+1,divΘuℓ+1)−Δt∑n=1ℓ[(dtΛξn+1,divΘun)+(divdtΛun+1,Θξn+1)]tag ≤12‖Λξℓ+1‖L2(Ω)2+12‖divΘuℓ+1‖L2(Ω)2+12Δt∑n=1ℓ[‖dtΛξn+1‖L2(Ω)2 +C‖ε(Θun+1)‖L2(Ω)2+‖divdtΛun+1‖L2(Ω)2+‖Θξn+1‖L2(Ω)2]. (3.51) When $$\theta=0$$ we also need to bound the last term on the right-hand side of (3.46), which is carried out below.    ∑n=1ℓ(dt∇Θξn+1,∇Θ^pn+1)≤∑n=1ℓ‖dtΘξn+1‖L2(Ω)‖∇Θ^pn+1‖L2(Ω) ≤∑n=1ℓsupvh∈Vhμ(dtε(Θun+1),ε(vh))−(dtΛξn+1,divvh)‖∇vh‖L2(Ω)‖∇Θ^pn+1‖L2(Ω) ≤∑n=1ℓ[μ2κ1Δth2β12‖dtε(Θun+1)‖L22+κ1Δth2β12‖dtΛξn+1‖L2+κ1−14Δt‖∇Θ^pn+1‖L2(Ω)2]. (3.52) Substituting (3.47)–(3.52) into (3.46) and rearranging terms, we get   μ‖ε(Θuℓ+1)‖L2(Ω)2+κ2‖Θηℓ+θ‖L2(Ω)2+κ3‖Θξℓ+1‖L2(Ω)2 +Δt∑n=1ℓ1μf‖K∇Θ^pn+1‖L2(Ω)2 ≤4μf(Δt)2K1‖ηtt‖L2((0,T);H1(Ω)′)2+4μ(Δt)2β12‖ηt‖L2((0,T);L2(Ω))2, +‖Λξℓ+1‖L2(Ω)2+Δt∑n=1ℓ‖dtΛξn+1‖L2(Ω)2+Δt∑n=1ℓ‖divdtΛun+1‖L2(Ω)2, (3.53) provided that $${\it{\Delta}} t\leq \mu_f\beta_1^2(4\mu\kappa_1^2K)^{-1} h^2$$ when $$\theta=0$$, but it holds for all $${\it{\Delta}} t>0$$ when $$\theta=1$$. Hence, (3.43) follows from using the approximation properties of the projection operators $$\mathcal{Q}_h, \mathcal{R}_h$$ and $$\mathcal{S}_h$$. The proof is complete. □ We conclude this section by stating the main theorem of the section. Theorem 3.7 The solution of the (FEA) satisfies the following error estimates:    max0≤n≤N[μ‖∇(u(tn)−uhn)‖L2(Ω)+κ2‖η(tn)−ηhn‖L2(Ω) +κ3‖ξ(tn)−ξhn‖L2(Ω)]≤C^1(T)Δt+C^2(T)h2. (3.54)   [Δt∑n=0N1μf‖∇(p(tn)−phn)‖L2(Ω)2]12≤C^1(T)Δt+C^2(T)h, (3.55) provided that $${\it{\Delta}} t=O(h^2)$$ when $$\theta=0$$ and $${\it{\Delta}} t>0$$ when $$\theta=1$$. Here,   C^1(T) :=C1(T),C^2(T) :=C2(T)+‖ξ‖L∞((0,T);H2(Ω))+‖η‖L∞((0,T);H2(Ω)) +‖∇u‖L∞((0,T);H2(Ω)). Proof. The above estimates follow immediately from an application of the triangle inequality on   u(tn)−uhn =Λun+Θun,ξ(tn)−ξhn=Λξn+Θξn,η(tn)−ηhn =Ληn+Θηn,p(tn)−phn=Λ^pn+Θ^pn and appealing to (3.24), (3.29) and Theorem 3.6. □ 4. Numerical experiments In this section, we shall present three two-dimensional numerical experiments to validate theoretical results for the proposed numerical methods, to numerically examine the performances of the approach and methods as well as to compare them with existing methods in the literature on two benchmark problems. One of these two problems was used to demonstrate the ‘locking phenomenon’ in Phillips & Wheeler (2009). Our numerical result shows that such a ‘locking phenomenon’ does not occur in our numerical methods, it confirms the fact that our approach and methods have a built-in mechanism to prevent the ‘locking phenomenon’. Test 1. Let $${\it{\Omega}}=[0,1]\times [0,1]$$, $${\it{\Gamma}}_1=\{(1,x_2); 0\leq x_2\leq 1\}$$, $${\it{\Gamma}}_2=\{(x_1,0); 0\leq x_1\leq 1\}$$, $${\it{\Gamma}}_3=\{(0,x_2); 0\leq x_2\leq 1\}$$, $${\it{\Gamma}}_4=\{(x_1,1); 0\leq x_1\leq 1\}$$ and $$T=0.001$$. We consider problems (2.26)–(2.28) with following source functions:   f =−(λ+μ)t(1,1)T+αcos⁡(x1+x2)et(1,1)T,ϕ =(c0+2κμf)sin⁡(x1+x2)et+α(x1+x2) and the following boundary and initial conditions:   p =sin⁡(x1+x2)eton ∂ΩT,u1 =12x12ton Γj×(0,T),j=1,3,u2 =12x22ton Γj×(0,T),j=,2,4,σn−αpn =f1on ∂ΩT,u(x,0)=0,p(x,0) =sin⁡(x1+x2)in Ω, where   f1(x,t)=μ(x1n1,x2n2)Tt+λ(x1+x2)(n1,n2)Tt−αsin⁡(x1+x2)(n1,n2)Tet. It is easy to check that the exact solution for this problem is   u(x,t)=t2(x12,x22)T,p(x,t)=sin⁡(x1+x2)et. We note that the boundary conditions used above are not pure Neumann conditions; instead, they are mixed Dirichlet–Neumann conditions. As pointed out in Remark 2.3 (c), the approach and methods of this article also apply to this case, the only change is to replace the test and trial space $${\mathbf{H}}^1_\perp({\it{\Omega}})$$ by $${\mathbf{H}}^1({\it{\Omega}})$$ with some appropriately built-in Dirichlet boundary condition in Definition 2.2. The goal of doing this test problem is to compute the order of the exact errors and to show that the theoretical error bounds proved in the previous section are sharp. Table 1 lists the computed $$L^\infty(0,T; L^2({\it{\Omega}}))$$ and $$L^2(0,T;H^1({\it{\Omega}}))$$-norm errors and the convergence rates with respect to $$h$$ at the terminal time $$T$$. In the test, $${\it{\Delta}} t=10^{-5}$$ is used so that the time error is negligible. Evidently, the spatial rates of convergence are consistent with that proved in the convergence theorem. Table 1. Spatial errors and convergence rates of Test 1    $$L^\infty(L^2)$$ error  $$L^\infty(L^2)$$ order  $$L^2(H^1)$$ error  $$L^2(H^1)$$ order  $$h=0.16$$  2.0789e-3     5.5045e-2     $$h=0.08$$  5.9674e-4  1.8006  2.9431e-2  0.9032  $$h=0.04$$  1.6227e-4  1.8787  1.5332e-2  0.9408  $$h=0.02$$  4.0971e-5  1.9857  7.6968e-3  0.9942     $$L^\infty(L^2)$$ error  $$L^\infty(L^2)$$ order  $$L^2(H^1)$$ error  $$L^2(H^1)$$ order  $$h=0.16$$  2.0789e-3     5.5045e-2     $$h=0.08$$  5.9674e-4  1.8006  2.9431e-2  0.9032  $$h=0.04$$  1.6227e-4  1.8787  1.5332e-2  0.9408  $$h=0.02$$  4.0971e-5  1.9857  7.6968e-3  0.9942  Figures 1 and 2 show, respectively, the surface plot of the computed pressure $$p$$ at the terminal time $$T$$ and the color plot of both the computed pressure $$p$$ and displacement $$\mathbf{u}$$ with mesh parameters $$h=0.02$$ and $${\it{\Delta}} t=10^{-5}$$. They coincide with the exact solution on the same space–time resolution. Fig. 1. View largeDownload slide Test 1: surface plot of the computed pressure $$p$$ at the terminal time $$T$$. Fig. 1. View largeDownload slide Test 1: surface plot of the computed pressure $$p$$ at the terminal time $$T$$. Fig. 2. View largeDownload slide Test 1: computed pressure $$p$$ (color plot) and displacement $$\mathbf{u}$$ (arrow plot) at $$T$$. Fig. 2. View largeDownload slide Test 1: computed pressure $$p$$ (color plot) and displacement $$\mathbf{u}$$ (arrow plot) at $$T$$. Test 2. In this test, we consider so-called Barry–Mercer’s problem, which is a Benchmark test problem for the poroelasticity model (2.26)–(2.28) (cf. Phillips & Wheeler, 2007b, 2009 and the references therein). Again, $${\it{\Omega}}=[0,1]\times [0,1]$$, but $$T=1$$. Barry–Mercer’s problem assumes no source, that is, $$\mathbf{f}\equiv 0$$ and $$\phi\equiv 0$$, and takes the following boundary conditions:   p =0on Γj×(0,T),j=1,3,4,p =p2on Γj×(0,T),j=2,u1 =0on  Γj×(0,T),j=1,3,u2 =0on  Γj×(0,T),j=2,4,σn−αpn =f1:=(0,αp)Ton ∂ΩT, where   p2(x1,t)={sin⁡twhen x∈[0.2,0.8)×(0,T),0others.  The boundary segments $${\it{\Gamma}}_j, j=1,2,3,4$$, which are defined in Test 1, and the above boundary conditions are depicted in Fig. 3. Also, the initial conditions for Barry–Mercer’s problem are $$\mathbf{u}(x,0) \equiv \mathbf{0}$$ and $$p(x,0)\equiv 0$$. We remark that Barry–Mercer’s problem has a unique solution that is given by an infinite series (cf. Phillips & Wheeler, 2009). Fig. 3. View largeDownload slide Test 2: boundary conditions. Fig. 3. View largeDownload slide Test 2: boundary conditions. Figures 4 and 5 display, respectively, the computed pressure $$p$$ (surface plot) and the computed displacement $$\mathbf{u}$$ (arrow plot). We note that the arrows near the boundary match very well with those on the boundary. Our numerical solution approximates the exact solution of Barry–Mercer’s problem very well and does not produce any oscillation in computed pressure. Fig. 4. View largeDownload slide Test 2: surface plot of the computed pressure $$p$$ at the terminal time $$T$$. Fig. 4. View largeDownload slide Test 2: surface plot of the computed pressure $$p$$ at the terminal time $$T$$. Fig. 5. View largeDownload slide Test 2: computed pressure $$p$$ (color plot) and displacement (arrow plot) at $$T$$. Fig. 5. View largeDownload slide Test 2: computed pressure $$p$$ (color plot) and displacement (arrow plot) at $$T$$. Test 3. This test problem is taken from Phillips & Wheeler (2009). Again, we consider problems (2.26)–(2.28) with $${\it{\Omega}}=[-0.5,0.5]\times [-0.5,0.5]$$. Let $${\it{\Gamma}}_j$$ be same as in Test 1 and $$c_0=0, E=10^5, \nu=0.49, \mu=33557, \lambda=1.6443\times10^6$$ and $$T=0.001$$. There is no source, that is, $$\mathbf{f}\equiv 0$$ and $$\phi\equiv 0$$. The boundary conditions are taken as   −κμf(∇p−ρfg)⋅n =0on ∂ΩT,u =0on Γ3×(0,T),σn−αpn =f1on Γj×(0,T),j=1,2,4, where $$\mathbf{f}_1=(f_1^1, f_1^2)$$ and   f11≡0on ∂ΩT,f12={0on Γj×(0,T),j=1,2,3,−1on Γ4×(0,T).  The computational domain $${\it{\Omega}}$$ and the above boundary conditions are depicted in Fig. 6. Also, the zero initial conditions are assigned for both $$\mathbf{u}$$ and $$p$$ in this test. Fig. 6. View largeDownload slide Test 3: boundary conditions. Fig. 6. View largeDownload slide Test 3: boundary conditions. Figures 7 and 8 display, respectively, the surface and color plot of the computed pressure, the arrow plot of the displacement vector and the deformation of the whole $${\it{\Omega}}$$. There is no oscillation in the computed pressure, and the arrows near the boundary match very well with arrows on the boundary, even when the Poisson ratio $$\nu=0.49$$, in which case the material is nearly incompressible. Note the finite element solution of the pressure deteriorates in the nearly incompressible case in Phillips & Wheeler (2009). Fig. 7. View largeDownload slide Test 3: computed pressure $$p$$: surface plot (left) and color plot (right). Fig. 7. View largeDownload slide Test 3: computed pressure $$p$$: surface plot (left) and color plot (right). Fig. 8. View largeDownload slide Test 3: arrow plot of the computed displacement (left) and deformation of $${\it{\Omega}}$$ (right). Fig. 8. View largeDownload slide Test 3: arrow plot of the computed displacement (left) and deformation of $${\it{\Omega}}$$ (right). We remark that the ‘locking phenomenon’ was observed in the simulation of Phillips & Wheeler (2009) at $$T=0.001$$ for this problem, namely, the computed pressure exhibits some oscillation at $$T=0.001$$. The reason for the locking phenomenon was explained as follows: when time step $${\it{\Delta}} t$$ is small, the displacement vector $$\mathbf{u}$$ is almost divergence free in the short time, while the numerical solution does not observe this nearly divergence free property, which results in the locking. However, at later times, the displacement vector is no longer divergence free, so no locking exists at later times. It is clear that our numerical solution does not exhibit the locking phenomenon at $$T=0.001$$. This is because our multiphysics reformulation weakly imposes the condition $$\mbox{div } \mathbf{u}=q$$, and hence $$\mathbf{u}$$ automatically becomes nearly divergence free when $$q\approx 0$$ for $$0<t<<1$$. Moreover, the pressure $$p$$ is not a primary variable anymore in our reformulation; instead, $$p$$ becomes a derivative variable and it is computed using the new primary variables $$\xi$$ and $$\eta$$. Therefore, our numerical methods are insensitive to the regularity of $$p$$. Funding NSF grants (DMS-1016173 and DMS-1318486 to X.F. and DMS-1016173 and DMS-1318486 to Y.L.); National Natural Science Foundation of China grant (#10901047 to Z.G.). References Bercovier J. & Pironneau O. ( 1979) Error estimates for finite element solution of the Stokes problem in the primitive variables. Numer. Math. , 33, 211– 224. Google Scholar CrossRef Search ADS   Biot M. ( 1955) Theory of elasticity and consolidation for a porous anisotropic media. J. Appl. Phys. , 26, 182– 185. Google Scholar CrossRef Search ADS   Brenner S. C. ( 1993) A nonconforming mixed multigrid method for the pure displacement problem in planar linear elasticity. SIAM J. Numer. Anal ., 30, 116– 135. Google Scholar CrossRef Search ADS   Brenner S. C. & Scott L. R. ( 2008) The Mathematical Theory of Finite Element Methods , 3rd edn. New York: Springer. Google Scholar CrossRef Search ADS   Brezzi F. & Fortin M. ( 1992) Mixed and Hybrid Finite Element Methods.  New York: Springer. Ciarlet P. G. ( 1978) The Finite Element Method for Elliptic Problems . Amsterdam: North-Holland. Coussy O. ( 2004) Poromechanics , Hoboken, NJ: Wiley & Sons. Dauge M. ( 1988) Elliptic Boundary Value Problems on Corner Domains . Lecture Notes in Math., vol. 1341. Berlin: Springer. Google Scholar CrossRef Search ADS   Dautray R. & Lions J. L. ( 1990) Mathematical Analysis and Numerical Methods for Science and Technology,  vol. 1. New York: Springer. Doi M. & Edwards S. F. ( 1986) The Theory of Polymer Dynamics . Oxford: Clarendon Press. Feng X., Ge Z. & Li Y. ( 2014) Multiphysics finite element methods for a poroelasticity model. arXiv:1411.7464 [matth.NA]. Feng X. & He Y. ( 2010) Fully discrete finite element approximations of a polymer gel model. SIAM J. Numer. Anal ., 48, 2186– 2217. Google Scholar CrossRef Search ADS   Girault V. & Raviart P. A. ( 1981) Finite Element Method for Navier-Stokes Equations: Theory and Algorithms . Berlin, Heidelberg, New York: Springer. Hamley I. ( 2007) Introduction to Soft Matter.  Hoboken, NJ: John Wiley & Sons. Lee J. J. ( 2014) Guaranteed locking-free finite element methods for Biot’s consolidation model in poroelasticity, arXiv:1403.7038v2 [math.NA]. Lee J. J. ( 2016) Robust error analysis of coupled mixed methods for Biot’s consolidation model. J. Sci. Comput ., 69, 610– 632. Google Scholar CrossRef Search ADS   Lee J. J., Mardal K. & Winther R. ( 2017) Parameter-robust discretization and preconditioning of Biot’s consolidation model. SIAM J. Sci. Comput. , 39, A1– A24. Google Scholar CrossRef Search ADS   Murad M. A. & Loula A. F. D. ( 1992) Improved accuracy in finite element analysis of Biot’s consolidation problem. Comput. Methods in Appl. Mech. and Engr ., 95, 359– 382. Google Scholar CrossRef Search ADS   Oyarzúa R. & Ruiz-Baier R. ( 2016) Locking-free finite element methods for poroelasticity. SIAM J. Numer. Anal. , 54, 2951– 2973. Google Scholar CrossRef Search ADS   Phillips P. J. & Wheeler M. F. ( 2007a) A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous in time case. Comput. Geosci ., 11, 131– 144. Google Scholar CrossRef Search ADS   Phillips P. J. & Wheeler M. F. ( 2007b) A coupling of mixed and continuous Galerkin finite element methods for poroelasticity II: the discrete in time case. Comput. Geosci ., 11, 145– 158. Google Scholar CrossRef Search ADS   Phillips P. J. & Wheeler M. F. ( 2009) Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach. Comput. Geosci ., 13, 1– 15. Google Scholar CrossRef Search ADS   Roberts J. E. & Thomas J. M. ( 1991) Mixed and hybrid methods. Handbook of Numerical Analysis  ( Ciarlet P. G. & Lions J. L. eds), vol. II. New York: North-Holland, pp. 523– 639. Google Scholar CrossRef Search ADS   Schowalter R. E. ( 2000) Diffusion in poro-elastic media. J. Math. Anal ., 251, 310– 340. Google Scholar CrossRef Search ADS   Schreyer-Bennethum L. ( 2007) Theory of flow and deformation of swelling porous materials at the macroscale. Computers and Geotechnics , 34, 267– 278. Google Scholar CrossRef Search ADS   Tanaka T. & Fillmore D. J. ( 1979) Kinetics of swelling of gels. J. Chem. Phys ., 70, 1214. Google Scholar CrossRef Search ADS   Temam R. ( 2000) Navier-Stokes Equations: Theory and Numerical Analysis . AMS Chelsea Publishing, Providence, RI: American Mathematical Society. Terzaghi K. ( 1943) Theoretical Soil Mechanics . New York: John Wiley & Sons. Google Scholar CrossRef Search ADS   Yamaue T. & Doi M. ( 2004) Swelling dynamics of constrained thin-plate under an external force. Phys. Rev. E  70, 011401. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Analysis of a multiphysics finite element method for a poroelasticity model

, Volume 38 (1) – Jan 1, 2018
30 pages

/lp/ou_press/analysis-of-a-multiphysics-finite-element-method-for-a-poroelasticity-VcEirFBDu9
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx003
Publisher site
See Article on Publisher Site

### Abstract

Abstract This article concerns finite element approximations of a quasi-static poroelasticity model in displacement–pressure formulation, which describes the dynamics of poro-elastic materials under an applied mechanical force on the boundary. To better describe the multiphysics process of deformation and diffusion for poroelastic materials, we first present a reformulation of the original model by introducing two pseudo-pressures. We then propose a time-stepping algorithm that decouples the reformulated partial differential equation (PDE) problem at each time step into two sub-problems: one of which is a generalized Stokes problem for the displacement vector field (of the solid network of the poroelastic material) along with one pseudo-pressure field, and the other is a diffusion problem for the other pseudo-pressure field (of the solvent of the material). To make this multiphysics approach feasible numerically, two critical issues must be resolved: the first one is the uniqueness of the generalized Stokes problem and the other is to find a good boundary condition for the diffusion equation, so that it also becomes uniquely solvable. To address the first issue, we discover certain conserved quantities for the PDE solution that provide ideal candidates for a needed boundary condition for the pseudo-pressure field. The solution to the second issue is to use the generalized Stokes problem to generate a boundary condition for the diffusion problem. A practical advantage of the time-stepping algorithm allows one to use any convergent Stokes solver (and its code) together with any convergent diffusion equation solver (and its code) to solve the poroelasticity model. In this article, the Taylor–Hood mixed finite element method combined with the $$P_1$$-conforming finite element method is used as an example to demonstrate the viability of the proposed multiphysics approach. It is proved that the solutions of the fully discrete finite element methods fulfill a discrete energy law, which mimics the differential energy law satisfied by the PDE solution and converges optimally in the energy norm. Moreover, it is showed that the proposed formulation also has a built-in mechanism to overcome so-called ‘locking phenomenon’ associated with the numerical approximations of the poroelasticity model. Numerical experiments are presented to show the performance of the proposed approach and methods, and to demonstrate the absence of ‘locking phenomenon’ in our numerical experiments. This article also presents a detailed PDE analysis for the poroelasticity model, especially it is proved that this model converges to the well-known Biot’s consolidation model from soil mechanics as the constrained specific storage coefficient tends to zero. As a result, the proposed approach and methods are robust under such a limit process. 1. Introduction A poroelastic material (or medium) is a fluid–solid interaction (FSI) system at pore scale, and poro-mechanics is a branch of continuum mechanics and acoustics that studies the behavior of fluid-saturated porous materials. If the solid is an elastic material, then the subject of the study is known as poroelasticity. Moreover, the elastic material may be governed by linear or nonlinear constitutive law, which then leads, respectively, to linear and nonlinear poroelasticity. Examples of poroelastic materials include soil, polymer gels and medicine pills. Poroelastic materials exhibit an important state of matter found in a wide variety of mechanical, biomedical and chemical systems (cf. Biot, 1955; Hamley, 2007; Terzaghi, 1943; Tanaka & Fillmore, 1979; Doi & Edwards, 1986; Coussy, 2004; Yamaue & Doi, 2004 and the references therein). They also possess some fascinating properties, in particular, they display thixotropy, which means that they become fluid when agitated, but resolidify when resting. In general, the behavior of a poroelastic material is described by a multiphysics FSI process at pore scale. However, for poroelastic materials to be considered in this article, one important physical phenomenon of the multiphysics process is not explicitly revealed in their mathematical model instead, it is hidden in the model. This is a different feature (and a source of extra difficulties) of poroelastic materials compared with standard (macroscopic) FSI systems. This article considers a general quasi-static model of linear poroelasticity, which is broad enough to contain the well-known Biot’s consolidation model from soil mechanics (cf. Murad & Loula, 1992; Schreyer-Bennethum, 2007) and the Doi’s model for polymer gels (cf. Yamaue & Doi, 2004; Feng & He, 2010). The quasi-static feature is due to the assumption that the acceleration of the solid (described by the second-order time derivative of the displacement vector field) is assumed to be negligible. For the Biot’s model, some locking-free formulations have recently been studied in Lee (2014, 2016), Oyarzúa & Ruiz-Baier (2016) and Lee et al. (2017). We refer the reader to Terzaghi (1943), Coussy (2004) and Phillips & Wheeler (2007a) for a derivation of the model and to Schowalter (2000) for its mathematical analysis. When the parameter $$c_0$$, called the constrained specific storage coefficient, vanishes in the model, it reduces to the Boit’s model and Doi’s model arising from two distinct applications. Their mathematical analysis can be found in Feng & He (2010), and their finite element numerical approximations based on two very different approaches were carried out in Murad & Loula (1992) and Feng & He (2010), respectively. In Phillips & Wheeler (2007a,b), the authors proposed and analysed a semidiscrete and a fully discrete mixed finite element method, which simultaneously approximate the pressure and its gradient along with the displacement vector field. Since the implicit Euler scheme is used for the time discretization, a combined linear system must be solved at each time step. It is observed in the numerical tests that the proposed fully discrete mixed finite method may exhibit a ‘locking phenomenon’ in the sense that the computed pressure oscillates and its accuracy deteriorates when a rapidly changed initial pressure is given. It is explained in Phillips & Wheeler (2009) that such a ‘locking phenomenon’ is caused by the difficulty of satisfying the nearly divergence-free condition of the displacement field for very small time $$t$$. As a comparison, we recall that the ‘volumetric locking’ is used in poromechanics to describe the phenomenon when a finite element solution deteriorates as the Poisson ratio approaches $$0.5$$, that is, when the material becomes incompressible or nearly incompressible. The goal of this article is to present a multiphysics approach for approximating the poroelasticity model. A key idea of this approach is to derive a multiphysics reformulation for the original model that clearly reveals the underlying multiple physics process (i.e., the deformation and diffusion) of the pore-scale FSI system. At end, two pseudo-pressures are introduced, one of them is shown to satisfy a diffusion equation, whereas the displacement vector field along with the other pseudo-pressure variable is shown to satisfy a generalized Stokes system. It should be noted that the original pressure is eliminated in the reformulation, and thus it is not approximated as a primary variable; instead, it is computed as a linear combination of the two pseudo-pressures. On the basis of this multiphysics reformulation, we propose a time-stepping algorithm that decouples (or couples) the reformulated PDE problem at each time step into two sub-problems, a generalized Stokes problem for the displacement vector field along with a pseudo-pressures and a diffusion problem for another pseudo-pressure field. To make this multiphysics approach feasible numerically, two critical issues must be resolved: the first one is the uniqueness of the generalized Stokes problem and the other is to find a good boundary condition for the diffusion equation so that it also becomes uniquely solvable. To overcome these difficulties, we discover certain conserved quantities for the PDE solution that can be imposed as needed boundary conditions for the subproblems. Moreover, we demonstrate that, regardless the choice of discretization methods, the proposed formulation has a built-in mechanism to overcome the ‘locking phenomenon’ associated with numerical approximations of the poroelasticity model. The remainder of this article is organized as follows. In Section 2, we present a complete PDE analysis of the poroelasticity model that emphasizes the energy law of the underlying model. Several conserved quantities are derived for the PDE solution. Moreover, it is proved that the poroelasticity model converges to the Biot’s consolidation model as the constrained specific storage coefficient $$c_0\to 0$$. In Section 3, we propose and analyse some fully discrete finite element methods based on the multiphysics reformulation. Both coupled and decoupled time stepping are considered and compared. The Taylor–Hood mixed finite element method combined with the $$P_1$$-conforming finite element method is chosen as an example for spatial discretization. It is proved that the solution of the fully discrete finite element methods fulfills a discrete energy law that mimics the differential energy law satisfied by the PDE solution. Optimal order error estimates in the energy norm are also established. Finally, in Section 4, several benchmark numerical experiments are provided to show the performance of the proposed approach and methods, and to demonstrate the absence of ‘locking phenomenon’ in our numerical experiments. 2. PDE model and its analysis 2.1 Preliminaries $${\it{\Omega}} \subset \mathbb{R}^d \,(d=1,2,3)$$ denotes a bounded polygonal domain with the boundary $${\partial}{\it{\Omega}}$$. The standard function space notation is adopted in this article, and their precise definitions can be found in Temam (1977), Ciarlet (1978) and Brenner & Scott (2008). In particular, $$(\cdot,\cdot)$$ and $$\langle \cdot,\cdot\rangle$$ denote, respectively, the standard $$L^2({\it{\Omega}})$$ and $$L^2({\partial}{\it{\Omega}})$$ inner products. For any Banach space $$B$$, we let $$\mathbf{B}=[B]^d$$ and use $$\mathbf{B}^\prime$$ to denote its dual space. In particular, we use $$(\cdot,\cdot)_{\small\rm dual}$$ to denote the dual product on $${\mathbf{H}}^1({\it{\Omega}})' \times {\mathbf{H}}^1({\it{\Omega}})$$, and $$\Vert\,{\cdot}\,\Vert_{L^p(B)}$$ is a shorthand notation for $$\Vert\,{\cdot}\,\Vert_{L^p((0,T);B)}$$. We also introduce the function spaces   L02(Ω):={q∈L2(Ω);(q,1)=0},X:=H1(Ω). It is well known (Temam, 1977) that the following so-called inf–sup condition holds in the space $${\mathbf{X}}\times L^2_0({\it{\Omega}})$$:   supv∈X(divv,φ)‖v‖H1(Ω)≥α0‖φ‖L2(Ω)∀φ∈L02(Ω),α0>0. (2.1) Let   RM:={r:=a+b×x;a,b,x∈Rd} denote the space of infinitesimal rigid motions. It is well known (Temam, 1977; Girault & Raviart, 1981; Brenner, 1993) that $${\mathbf{RM}}$$ is the kernel of the strain operator $${\varepsilon}$$, that is, $${\mathbf{r}}\in {\mathbf{RM}}$$ if and only if $${\varepsilon}({\mathbf{r}})=0$$. Hence, we have   ε(r)=0divr=0,∀r∈RM. (2.2) Let $${\mathbf{L}}^2_\bot({\partial}{\it{\Omega}})$$ and $${\mathbf{H}}^1_\bot({\it{\Omega}})$$ denote, respectively, the subspaces of $${\mathbf{L}}^2({\partial}{\it{\Omega}})$$ and $${\mathbf{H}}^1({\it{\Omega}}),$$ which are orthogonal to $${\mathbf{RM}}$$, that is,   H⊥1(Ω):={v∈H1(Ω);(v,r)=0∀r∈RM},L⊥2(∂Ω):={g∈L2(∂Ω);⟨g,r⟩=0∀r∈RM}. It is well known (Dautray & Lions, 1990) that there exists a constant $$c_1>0$$ such that   infr∈RM‖v+r‖L2(Ω)≤c1‖ε(v)‖L2(Ω)∀v∈H1(Ω). Hence, for each $${\mathbf{v}}\in {\mathbf{H}}^1_\bot({\it{\Omega}})$$ there holds   ‖v‖L2(Ω)=infr∈RM‖v+r‖L2(Ω)2−‖r‖L2(Ω)2≤c1‖ε(v)‖L2(Ω), (2.3) which, and the well-known Korn’s inequality (Dautray & Lions, 1990), yield that for some $$c_2>0$$  ‖v‖H1(Ω) ≤c2[‖v‖L2(Ω)+‖ε(v)‖L2(Ω)] ≤c2(1+c1)‖ε(v)‖L2(Ω)∀v∈H⊥1(Ω). (2.4) By Lemma 2.1 of Brenner (1993), we know that for any $$q\in L^2({\it{\Omega}})$$, there exists $${\mathbf{v}}\in {\mathbf{H}}^1_\bot({\it{\Omega}})$$ such that $$\text{div}\,\, {\mathbf{v}}=q$$ and $$\|{\mathbf{v}}\|_{H^1({\it{\Omega}})} \leq C\|q\|_{L^2({\it{\Omega}})}$$. An immediate consequence of this lemma is that there holds the following alternative version of the inf–sup condition:   supv∈H⊥1(Ω)(divv,φ)‖v‖H1(Ω)≥α1‖φ‖L2(Ω)∀φ∈L02(Ω),α1>0. (2.5) Throughout the article, we assume $${\it{\Omega}} \subset \mathbb{R}^d$$ is a bounded polygonal domain such that $${\it{\Delta}}: H^1_0({\it{\Omega}}) \cap H^2({\it{\Omega}}) \rightarrow L^2({\it{\Omega}})$$ is an isomorphism (cf. Girault & Raviart, 1981; Dauge, 1988). In addition, $$C$$ is used to denote a generic positive (pure) constant which may be different in different places. 2.2 PDE model and its multiphysics reformulation The quasi-static poroelasticity model to be studied in this article is given by Phillips & Wheeler (2007a)   −divσ(u)+α∇p =fin ΩT:=Ω×(0,T)⊂Rd×(0,T), (2.6)  (c0p+αdivu)t+divvf =ϕin ΩT, (2.7) where   σ(u) :=με(u)+λdivuI,ε(u):=12(∇u+∇uT), (2.8)  vf :=−Kμf(∇p−ρfg). (2.9) Here, $${\mathbf{u}}$$ denotes the displacement vector of the solid and $$p$$ denotes the pressure of the solvent. $${\mathbf{f}}$$ is the body force and $$\mathbf{g}$$ is the gravitational acceleration, which is assumed to be a constant vector. $$I$$ denotes the $$d\times d$$ identity matrix, and $${\varepsilon}({\mathbf{u}})$$ is known as the strain tensor. The parameters in the model are Lamé constants $${\lambda}$$ and $$\mu$$ (warning: for the sake of notation brevity later, we use $$\mu$$ instead of $$2\mu$$ in (2.8) as often seen in the literature, in other words, precisely, here $$\mu/2$$ is a Lamé constant); the permeability tensor $$K=K(x)$$, which is assumed to be symmetric and uniformly positive definite in the sense that there exist positive constants $$K_1$$ and $$K_2$$ such that $$K_1|\zeta|^2\leq K(x)\zeta\cdot \zeta \leq K_2 |\zeta|^2$$ for a.e. $$x\in{\it{\Omega}}$$ and any $$\zeta\in \mathbf{\mathbb{R}}^d$$; the solvent viscosity $$\mu_f$$, Biot–Willis constant $$\alpha$$, and the constrained specific storage coefficient $$c_0$$. In addition, $$\sigma({\mathbf{u}})$$ is called the (effective) stress tensor. $$\widehat{\sigma}({\mathbf{u}},p):=\sigma({\mathbf{u}})-\alpha p I$$ is the total stress tensor. $${\mathbf{v}}_f$$ is the volumetric solvent flux and (2.9) is the well-known Darcy’s law. We assume that $$\rho_f\not\equiv 0$$, which is a realistic assumption. To close the above system, suitable boundary and initial conditions must be prescribed. The following set of boundary and initial conditions will be considered in this article:   σ^(u,p)n=σ(u)n−αpn =f1on ∂ΩT:=∂Ω×(0,T), (2.10)  vf⋅n=−Kμf(∇p−ρfg)⋅n =ϕ1on ∂ΩT, (2.11)  u=u0,p =p0in Ω×{t=0}. (2.12) We note that in some engineering literature, the second Lamé constant $$\mu$$ is also called the shear modulus and denoted by $$G$$, and $$B:={\lambda} +\frac23 G$$ is called the bulk modulus. $${\lambda}, \mu$$ and $$B$$ are computed from the Young’s modulus$$E$$ and the Poisson ratio$$\nu$$ by the following formulas:   λ=Eν(1+ν)(1−2ν),μ=G=E2(1+ν),B=E3(1−2ν). Unlike the existing approaches in the literature (Murad & Loula, 1992; Phillips & Wheeler, 2007a), in this article we will not directly approximate the above original model; instead, we first derive a (multiphysics) reformulation for the model, we then approximate the reformulated model. This is a key idea of this article, and it will be seen in the later sections that this new approach is advantageous. To this end, we introduce new variables   q:=divu,η:=c0p+αq,ξ:=αp−λq. It is easy to check that   p=κ1ξ+κ2η,q=κ1η−κ3ξ, (2.13) where   κ1:=αα2+λc0,κ2:=λα2+λc0,κ3:=c0α2+λc0. (2.14) Then (2.6)–(2.9) can be written as   −μdivε(u)+∇ξ =fin ΩT, (2.15)  κ3ξ+divu =κ1ηin ΩT, (2.16)  ηt−1μfdiv[K(∇(κ1ξ+κ2η)−ρfg)] =ϕin ΩT, (2.17) where $$p$$ and $$q$$ are related to $$\xi$$ and $$\eta$$ through the algebraic equations in (2.13). The boundary and initial conditions (2.10)–(2.12) can be rewritten as   σ(u)n−α(κ1ξ+κ2η)n =f1on ∂ΩT:=∂Ω×(0,T), (2.18)  −Kμf(∇(κ1ξ+κ2η)−ρfg)⋅n =ϕ1on ∂ΩT, (2.19)  u=u0,p =p0in Ω×{t=0}. (2.20) It is now clear that $$({\mathbf{u}}, \xi)$$ satisfies a generalized Stokes problem for a given $$\eta$$, where $$\kappa_3\xi$$ in (2.16) acts as a penalty term and $$\eta$$ satisfies a diffusion problem for a given $$\xi$$. This new formulation reveals the underlying deformation and diffusion multiphysics process, which occurs in the poroelastic material. In particular, the diffusion part of the process is hidden in the original formulation, but is apparent in the new formulation. To make use of the above reformulation for computation, we need to address a crucial issue of the uniqueness for the generalized Stokes problem and the diffusion problem after they are decoupled. The difficulty will be overcome by discovering some invariant quantities for the solution of the PDE model and using them to impose some appropriate boundary conditions for both subproblems (cf. Lemma 2.9). 2.3 Analysis of the PDE model We start this section with a definition of weak solutions to problems (2.6)–(2.12). For convenience, we assume that $${\mathbf{f}}, {\mathbf{f}}_1,\phi$$ and $$\phi_1$$ all are independent of $$t$$ in the remaining of the article. We note that all the results of this article can be easily extended to the case of time-dependent source functions. Definition 2.1 Let $${\mathbf{u}}_0\in{\mathbf{H}}^1({\it{\Omega}}), {\mathbf{f}}\in{\mathbf{L}}^2({\it{\Omega}}), {\mathbf{f}}_1\in {\mathbf{L}}^2({\partial}{\it{\Omega}}), p_0\in L^2({\it{\Omega}}), \phi\in L^2({\it{\Omega}})$$ and $$\phi_1\in L^2({\partial}{\it{\Omega}})$$. Assume $$c_0>0$$ and $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$ for any $${\mathbf{v}}\in \mathbf{RM}$$. Given $$T > 0$$, a tuple $$({\mathbf{u}},p)$$ with   u∈L∞(0,T;H⊥1(Ω)),p∈L∞(0,T;L2(Ω))∩L2(0,T;H1(Ω)),pt,(divu)t∈L2(0,T;H1(Ω)′) is called a weak solution to (2.6)–(2.12), if it holds for almost every $$t \in [0,T]$$  μ(ε(u),ε(v))+λ(divu,divv)−α(p,divv) =(f,v)+⟨f1,v⟩∀v∈H1(Ω), (2.21)   ((c0p+αdivu)t,φ)dual+1μf(K(∇p−ρfg),∇φ) =(ϕ,φ)+⟨ϕ1,φ⟩|,∀φ∈H1(Ω), (2.22)  u(0)=u0,p(0)=p0. (2.23) Similarly, we can define weak solutions to problems (2.15)–(2.17) and (2.18)–(2.20). Definition 2.2 Let $${\mathbf{u}}_0\in {\mathbf{H}}^1({\it{\Omega}}), {\mathbf{f}} \in {\mathbf{L}}^2({\it{\Omega}}), {\mathbf{f}}_1 \in {\mathbf{L}}^2({\partial}{\it{\Omega}}), p_0\in L^2({\it{\Omega}}), \phi\in L^2({\it{\Omega}})$$ and $$\phi_1\in L^2({\partial}{\it{\Omega}})$$. Assume $$c_0>0$$ and $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$ for any $${\mathbf{v}}\in \mathbf{RM}$$. Given $$T > 0$$, a $$5$$-tuple $$({\mathbf{u}},\xi,\eta,p,q)$$ with   u∈L∞(0,T;H⊥1(Ω)),ξ∈L∞(0,T;L2(Ω)),η∈L∞(0,T;L2(Ω))∩H1(0,T;H1(Ω)′),q∈L∞(0,T;L2(Ω)),p∈L∞(0,T;L2(Ω))∩L2(0,T;H1(Ω)) is called a weak solution to (2.15)–(2.17) and (2.18)–(2.20) if it holds for almost every $$t \in [0,T]$$  μ(ε(u),ε(v))−(ξ,divv) =(f,v)+⟨f1,v⟩∀v∈H1(Ω), (2.24)  κ3(ξ,φ)+(divu,φ) =κ1(η,φ)∀φ∈L2(Ω), (2.25)  (ηt,ψ)dual+1μf(K(∇p −ρfg),∇ψ) =(ϕ,ψ)+⟨ϕ1,ψ⟩∀ψ∈H1(Ω), (2.26)  p:=κ1ξ+κ2η,q:=κ1η−κ3ξ, (2.27)  η(0)=η0: =c0p0+αq0, (2.28) where $$q_0:=\text{div}\,\, {\mathbf{u}}_0$$, $$u_0$$ and $$p_0$$ are same as in Definition 2.1. Remark 2.3 (a) We note that $$p$$ and $$q$$ are derivative variables in the system, after $${\mathbf{u}}$$, $$\xi$$ and $$\eta$$ are computed, $$p$$ and $$q$$ are simply updated by their algebraic expressions in (2.27). (b) Equation (2.26) implicitly imposes the following boundary condition for $$\eta$$:   κ2K∇η⋅n=Kρfg⋅n−κ1K∇ξ⋅n. (2.29) (c) By an Aubin–Lions Lemma (Temam, 1977), we conclude that $${\mathbf{u}}\in C^0([0,T];{\mathbf{L}}^2({\it{\Omega}}))$$, $$p\in C^0([0,T]; L^2({\it{\Omega}}))$$ and $$\eta\in C^0([0,T]; L^2({\it{\Omega}}))$$. Hence, the initial conditions in (2.23) and in (2.28) all make sense. (d) It should be pointed out that the only reason for introducing the space $${\mathbf{H}}_\perp^1({\it{\Omega}})$$ in the above two definitions is that the boundary condition (2.10) is a pure ‘Neumann condition’. If it is replaced by a pure Dirichlet condition or by a mixed Dirichlet–Neumann condition, there is no need to introduce this space. Neither do we need to impose the compatibility condition $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$$$\forall {\mathbf{v}}\in \mathbf{RM}$$. This fact will be used in our numerical experiments in Section 4. We also note that from the analysis point of view, the pure Neumann condition case is the most difficult case. Lemma 2.4 Every weak solution $$({\mathbf{u}},p)$$ of problems (2.21)–(2.23) satisfies the following energy law:   E(t)+1μf∫0t(K(∇p−ρfg),∇p)ds−∫0t(ϕ,p)ds −∫0t⟨ϕ1,p⟩,ds=E(0) (2.30) for all $$t\in [0,T]$$, where   E(t): =12[μ‖ε(u(t))‖L2(Ω)2+λ‖divu(t)‖L2(Ω)2+c0‖p(t)‖L2(Ω)2−2(f,u(t))−2⟨f1,u(t)⟩]. (2.31) Moreover,   ‖(c0p+αdivu)t‖L2(0.T;H1(Ω)′) ≤1μf‖K∇p−ρfg‖L2(ΩT)+‖ϕ‖L2(ΩT)+‖ϕ1‖L2(∂ΩT)<∞. (2.32) Likewise, weak solutions of (2.24)–(2.28) satisfy a similar energy law which is a rewritten version of (2.30) in the new variables. Lemma 2.5 Every weak solution $$({\mathbf{u}},\xi,\eta,p,q)$$ of problems (2.24)–(2.28) satisfies the following energy law:   J(t)+1μf∫0t(K(∇p−ρfg),∇p)ds−∫0t(ϕ,p)ds −∫0t⟨ϕ1,p⟩ds=J(0) (2.33) for all $$t\in [0,T]$$, where   J(t): =12[μ‖ε(u(t))‖L2(Ω)2+κ2‖η(t)‖L2(Ω)2+κ3‖ξ(t)‖L2(Ω)2−2(f,u(t))−2⟨f1,u(t)⟩]. (2.34) Moreover,   ‖ηt‖L2(0.T;H1(Ω)′) ≤1μf‖K∇p−ρfg‖L2(ΩT)+‖ϕ‖L2(ΩT)+‖ϕ1‖L2(∂ΩT)<∞. (2.35) We skip the proofs of the above two lemmas to save space and refer the interested reader to Feng et al. (xxxx) for their proofs. The above energy law immediately implies the following solution estimates. Lemma 2.6 There exists a positive constant $$C_1=C_1\bigl(\|{\mathbf{u}}_0\|_{H^1({\it{\Omega}})}, \|p_0\|_{L^2({\it{\Omega}})},$$$$\|{\mathbf{f}}\|_{L^2({\it{\Omega}})},\|{\mathbf{f}}_1\|_{L^2({\partial} {\it{\Omega}})},\|\phi\|_{L^2({\it{\Omega}})}, \|\phi_1\|_{L^2({\partial}{\it{\Omega}})} \bigr)$$ such that   μ‖ε(u)‖L∞(0,T;L2(Ω))+κ2‖η‖L∞(0,T;L2(Ω)) +κ3‖ξ‖L∞(0,T;L2(Ω))+K1μf‖∇p‖L2(0,T;L2(Ω))≤C1. (2.36)   ‖u‖L∞(0,T;L2(Ω))≤C1,‖p‖L∞(0,T;L2(Ω))≤C1(κ212+κ1κ3−12). (2.37)   ‖p‖L2(0,T;L2(Ω))≤C1,‖ξ‖L2(0,T;L2(Ω))≤C1κ1−1(1+κ212). (2.38) We note that (2.38) follows from (2.36), (2.3), the Poincaré inequality and (2.50) below, and the relation $$p=\kappa_1\xi +\kappa_2\eta$$. Furthermore, by exploiting the linearity of the PDE system, we have the following a priori estimates for the weak solution. Theorem 2.7 Suppose that $${\mathbf{u}}_0$$ and $$p_0$$ are sufficiently smooth, then there exist positive constants $$C_2=C_2\bigl(C_1,\|{\nabla} p_0\|_{L^2({\it{\Omega}})} \bigr)$$ and $$C_3=C_3\bigl(C_1,C_2, \|{\mathbf{u}}_0\|_{H^2({\it{\Omega}})},\|p_0\|_{H^2({\it{\Omega}})} \bigr)$$ such that there hold the following estimates for the solution to problems (2.15)–(2.17) and (2.10)–(2.12):   μ‖ε(ut)‖L2(0,T;L2(Ω))+κ2‖ηt‖L2(0,T;L2(Ω)) +κ3‖ξt‖L2(0,T;L2(Ω))+K1μf‖∇p‖L∞(0,T;L2(Ω))≤C2. (2.39)  μ‖ε(ut)‖L∞(0,T;L2(Ω))+κ2‖ηt‖L∞(0,T;L2(Ω)) +κ3‖ξt‖L∞(0,T;L2(Ω))+K1μf‖∇pt‖L2(0,T;L2(Ω))≤C3. (2.40)   ‖ηtt‖L2(H1(Ω)′)≤K2μfC3. (2.41) Proof. On noting that $${\mathbf{f}}, {\mathbf{f}}_1,\phi$$ and $$\phi_1$$ all are assumed to be independent of $$t$$, differentiating (2.24) and (2.25) with respect to $$t$$, taking $${\mathbf{v}}={\mathbf{u}}_t$$ and $$\varphi=\xi_t$$ in (2.24) and (2.25), respectively, and adding the resulting equations yield   μ‖ε(ut)‖L2(Ω)2=(qt,ξt)=κ1(ηt,ξt)−κ3‖ξt‖L2(Ω)2. (2.42) Setting $$\psi=p_t=\kappa_1\xi_t + \kappa_2\eta_t$$ in (2.26) gives   κ1(ηt,ξt)+κ2‖ηt‖L2(Ω)2+K2μfddt‖∇p−ρfg‖L2(Ω)2=ddt[(ϕ,p)+⟨ϕ1,p⟩]. (2.43) Adding (2.42) and (2.43) and integrating in $$t$$, we get for $$t\in[0,T]$$  K2μf‖∇p(t)−ρfg‖L2(Ω)2+∫0t[μ‖ε(ut)‖L2(Ω)2+κ2‖ηt‖L2(Ω)2+κ3‖ξt‖L2(Ω)2]ds =K2μf‖∇p0−ρfg‖L2(Ω)2+(ϕ,p(t)−p0)+⟨ϕ1,p(t)−p0⟩, which readily infers (2.39). To show (2.40), first differentiating (2.24) one time with respect to $$t$$ and setting $${\mathbf{v}}={\mathbf{u}}_{tt}$$, differentiating (2.25) twice with respect to $$t$$ and setting $$\varphi=\xi_t$$ and adding the resulting equations, we get   μ2ddt‖ε(ut)‖L2(Ω)2=(qtt,ξt)=κ1(ηtt,ξt)−κ32ddt‖ξt‖L2(Ω)2. (2.44) Secondly, differentiating (2.26) with respect $$t$$ one time and taking $$\psi=p_t=\kappa_1\xi_t + \kappa_2\eta_t$$ yield   κ1(ηtt,ξt)+κ22ddt‖ηt‖L2(Ω)2+Kμf‖∇pt‖L2(Ω)2=0. (2.45) Finally, adding the above two inequalities and integrating in $$t$$ give for $$t\in [0,T]$$  μ‖ε(ut(t))‖L2(Ω)2+κ2‖ηt(t)‖L2(Ω)2+κ3‖ξt(t)‖L2(Ω)2+2Kμf∫0t‖∇pt‖L2(Ω)2ds =μ‖ε(ut(0))‖L2(Ω)2+κ2‖ηt(0)‖L2(Ω)2+κ3‖ξt(0)‖L2(Ω)2. (2.46) Hence, (2.40) holds. Equation (2.41) follows immediately from the following inequality   (ηtt,ψ)=−1μf(K∇p,∇ψ)≤Kμf‖∇p‖L2(Ω)‖∇ψ‖L2(Ω)∀ψ∈H01(Ω), (2.40) and the definition of the $$H^{1}({\it{\Omega}})'$$-norm. The proof is complete. □ Remark 2.8 As expected, the above high order norm solution estimates require $$p_0\in H^1({\it{\Omega}}), {\mathbf{u}}_t(0)\in {\mathbf{L}}^2({\it{\Omega}}), \eta_t(0)\in L^2({\it{\Omega}})$$ and $$\xi_t(0)\in L^2({\it{\Omega}})$$. The values of $${\mathbf{u}}_t(0), \eta_t(0)$$ and $$\xi_t(0)$$ can be computed using the PDEs as follows. It follows from (2.17) that $$\eta_t(0)$$ satisfies   ηt(0)=ϕ+1μfdiv[K(∇p0−ρfg)]. Hence, $$\eta_t(0)\in L^2({\it{\Omega}})$$ provided that $$p_0\in H^2({\it{\Omega}})$$. To find $${\mathbf{u}}_t(0)$$ and $$\xi_t(0)$$, differentiating (2.15) and (2.16) with respect to $$t$$ and setting $$t=0$$, we get   −μdivε(ut(0))+∇ξt(0) =0in Ω,κ3ξt(0)+divut(0) =κ1ηt(0)in Ω. Hence, $${\mathbf{u}}_t(0)$$ and $$\xi_t(0)$$ can be determined by solving the above generalized Stokes problem. The next lemma shows that weak solutions of problems (2.24)–(2.28) preserve some ‘invariant’ quantities, it turns out that these ‘invariant’ quantities play a vital role in the construction of our time-splitting scheme to be introduced in the next section. Lemma 2.9 Every weak solution $$({\mathbf{u}},\xi,\eta,p, q)$$ to problems (2.24)–(2.28) satisfies the following relations:   Cη(t):=(η(⋅,t),1)=(η0,1)+[(ϕ,1)+⟨ϕ1,1⟩]tt≥0. (2.47)  Cξ(t):=(ξ(⋅,t),1)=1d+μκ3[μκ1Cη(t)−(f,x)−⟨f1,x⟩]. (2.48)  Cq(t):=(q(⋅,t),1)=κ1Cη(t)−κ3Cξ(t). (2.49)  Cp(t):=(p(⋅,t),1)=κ1Cξ(t)+κ2Cη(t). (2.50)  Cu(t):=⟨u(⋅,t)⋅n,1⟩=Cq(t). (2.51) We skip the proof of the above lemma again to save space and refer the interested reader to Feng et al. (xxxx) for its proof. Remark 2.10 We note that $$C_\eta, C_\xi, C_q$$ and $$C_p$$ all are (known) linear functions of $$t$$, and they become (known) constants when $$\phi\equiv 0$$ and $$\phi_1\equiv 0$$. With the help of the above lemmas, we can show the solvability of problems (2.6)–(2.12). Theorem 2.11 Let $${\mathbf{u}}_0\in{\mathbf{H}}^1({\it{\Omega}}), {\mathbf{f}}\in{\mathbf{L}}^2({\it{\Omega}}), {\mathbf{f}}_1\in {\mathbf{L}}^2({\partial}{\it{\Omega}}), p_0\in L^2({\it{\Omega}}), \phi\in L^2({\it{\Omega}})$$ and $$\phi_1\in L^2({\partial}{\it{\Omega}})$$. Suppose $$c_0>0$$ and $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$ for any $${\mathbf{v}}\in \mathbf{RM}$$. Then there exists a unique solution to problems (2.6)–(2.12) in the sense of Definition 2.1, likewise, there exists a unique solution to problems (2.15)–(2.17) and (2.18)–(2.20) in the sense of Definition 2.2. Proof. We only outline the main steps of the proof and leave the details to the interested reader. First, Because the PDE system is linear, the existence of weak solution can be proved by the standard Galerkin method and compactness argument (cf. Temam, 1977). We note that the energy laws established in Lemmas 2.4 and 2.5 guarantee the required uniform estimates for the Galerkin approximate solutions. Secondly, to show the uniqueness, suppose there are two sets of weak solutions, again by the linearity of the PDE system it is trivial to show that the difference of the solutions satisfy the same PDE system with zero initial and boundary data. The energy law immediately implies that the difference must be zero, and hence the uniqueness is verified. □ We conclude this section by establishing a convergence result for the solution of problems (2.15)–(2.17) and (2.18)–(2.20) when the constrained specific storage coefficient $$c_0$$ tends to $$0$$. Such a convergence result is useful and significant for the following two reasons. First, as mentioned earlier, the poroelasticity model studied in this article reduces to the well-known Biot’s consolidation model from soil mechanics (cf. Murad & Loula, 1992; Phillips & Wheeler, 2007a) and Doi’s model for polymer gels (cf. Yamaue & Doi, 2004; Feng & He, 2010). Secondly, it proves that the proposed approach and methods of this article are robust under such a limit process. Theorem 2.12 Let $${\mathbf{u}}_0\in{\mathbf{H}}^1({\it{\Omega}}), {\mathbf{f}}\in{\mathbf{L}}^2({\it{\Omega}}), {\mathbf{f}}_1\in {\mathbf{L}}^2({\partial}{\it{\Omega}}), p_0\in L^2({\it{\Omega}}), \phi\in L^2({\it{\Omega}})$$ and $$\phi_1\in L^2({\partial}{\it{\Omega}})$$. Suppose $$({\mathbf{f}},{\mathbf{v}})+\langle {\mathbf{f}}_1, {\mathbf{v}} \rangle =0$$ for any $${\mathbf{v}}\in \mathbf{RM}$$. Let $$({\mathbf{u}}_{c_0},\eta_{c_0},\xi_{c_0},p_{c_0},q_{c_0})$$ denote the unique weak solution to problems (2.15)–(2.17) and (2.18)–(2.20). Then there exists $$({\mathbf{u}}_*, \eta_*,\xi_*, p_*, q_*)\in {\mathbf{L}}^\infty(0,T;{\mathbf{H}}^1_\perp({\it{\Omega}}))\times L^\infty(0,T;L^2({\it{\Omega}}))\times L^2(0,T; L^2({\it{\Omega}}))\times L^2(0,T;H^1({\it{\Omega}})\times L^\infty(0,T;L^2({\it{\Omega}})))$$ such that $$({\mathbf{u}}_{c_0},\eta_{c_0},\xi_{c_0},p_{c_0},q_{c_0})$$ converges weakly to $$({\mathbf{u}}_*, \eta_*,\xi_*, p_*, q_*)$$ in the above product space as $$c_0\to 0$$. Proof. It follows immediately from (2.35)–(2.38) and Korn’s inequality that $${\mathbf{u}}_{c_0}$$ is uniformly bounded (in $$c_0$$) in $${\mathbf{L}}^\infty(0,T; {\mathbf{H}}^1_\perp({\it{\Omega}}))$$. $$\sqrt{\kappa_2} \eta_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^\infty(0,T; L^2({\it{\Omega}}))\cap H^1(0,T; H^{1}({\it{\Omega}})')$$. $$\sqrt{\kappa_3}\xi_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^\infty(0,T; L^2({\it{\Omega}}))$$. $$\xi_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^2({\it{\Omega}}_T)$$. $$p_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^2(0,T; H^1({\it{\Omega}}))$$. $$q_{c_0}$$ is uniformly bounded (in $$c_0$$) in $$L^\infty(0,T; L^2({\it{\Omega}}))$$. On noting that $$\lim_{c_0\to 0}\kappa_1=\frac{1}{\alpha}$$, $$\lim_{c_0\to 0}\kappa_2=\frac{{\lambda}}{\alpha^2}$$ and $$\lim_{c_0\to 0}\kappa_3=0$$, by the weak compactness of reflexive Banach spaces and Aubin–Lions Lemma (Dautray & Lions, 1990), we have that there exist $$({\mathbf{u}}_*, \eta_*,\xi_*, p_*, q_*)\in {\mathbf{L}}^\infty(0,T;{\mathbf{H}}^1_\perp({\it{\Omega}})) \times L^\infty(0,T;L^2({\it{\Omega}}))\times L^2(0,T; L^2({\it{\Omega}}))\times L^2(0,T;H^1({\it{\Omega}}))\times L^\infty(0,T;L^2({\it{\Omega}}))$$ and a subsequence of $$({\mathbf{u}}_{c_0},\eta_{c_0},\xi_{c_0},$$$$p_{c_0},q_{c_0})$$ (still denoted by the same notation) such that as $$c_0\to 0$$ (a subsequence of $$c_0$$, to be exact) $${\mathbf{u}}_{c_0}$$ converges to $${\mathbf{u}}_*$$ weak $$*$$ in $${\mathbf{L}}^\infty(0,T; {\mathbf{H}}^1_\perp({\it{\Omega}}))$$ and weakly in $${\mathbf{L}}^2(0,T; {\mathbf{H}}^1_\perp({\it{\Omega}}))$$. $$\sqrt{\kappa_2} \eta_{c_0}$$ converges to $$\frac{\sqrt{{\lambda}}}{\alpha} \eta_*$$ weak $$*$$ in $$L^\infty(0,T; L^2({\it{\Omega}}))$$ and strongly in $$L^2({\it{\Omega}}_T)$$. $$\kappa_3\xi_{c_0}$$ converges to $$0$$ strongly in $$L^2({\it{\Omega}}_T)$$. $$\xi_{c_0}$$ converges to $$\xi_*$$ weakly in $$L^2({\it{\Omega}}_T)$$. $$p_{c_0}$$ converges to $$p_*$$ weakly in $$L^2(0,T; H^1({\it{\Omega}}))$$. $$q_{c_0}$$ converges to $$p_*$$ weak $$*$$ in $$L^\infty(0,T; L^2({\it{\Omega}}))$$ and weakly in $$L^2({\it{\Omega}}_T)$$. Then setting $$c_0\to 0$$ in (2.24)–(2.28) yields (note that the dependence of the solution on $$c_0$$ is suppressed there)   2μ(ε(u∗),ε(v))−(ξ∗,divv) =(f,v)+⟨f1,v⟩∀v∈H1(Ω),(divu∗,φ) =1α(η∗,φ)∀φ∈L2(Ω),((η∗)t,ψ)dual+1μf(K(∇p∗ −ρfg),∇ψ) =(ϕ,ψ)+⟨ϕ1,ψ⟩∀ψ∈H1(Ω),p∗:=1αξ∗+λα2η∗,q∗:=1αη∗in L2(ΩT),u∗(0)=u0,η∗(0)=η0: =αq0. Equivalently,   2μ(ε(u∗),ε(v))−(ξ∗,divv) =(f,v)+⟨f1,v⟩∀v∈H1(Ω),(divu∗,φ) =(q∗,φ)∀φ∈L2(Ω),α((q∗)t,ψ)dual+1μf(K(∇(λα−1q∗+α−1ξ∗) −ρfg),∇ψ) =(ϕ,ψ)+⟨ϕ1,ψ⟩∀ψ∈H1(Ω),p∗ :=1α(ξ∗+λq∗)in L2(ΩT),q∗(0)=q0 :=divu0. Hence, $$({\mathbf{u}}_*, \eta_*,\xi_*,p_*,q_*)$$ is a weak solution of Biot’s consolidation model (cf. Yamaue & Doi, 2004; Feng & He, 2010). By the uniqueness of its solutions, we conclude that the whole sequence $$({\mathbf{u}}_{c_0},\eta_{c_0},\xi_{c_0},$$$$p_{c_0},q_{c_0})$$ converges to $$({\mathbf{u}}_*, \eta_*,\xi_*,p_*,q_*)$$ as $$c_0\to 0$$ in the above sense. The proof is complete. □ 3. Fully discrete finite element methods The goal of this section is to design and analyse some fully discrete finite element methods for the poroelasticity model based on the above new formulation. 3.1 Formulation of fully discrete finite element methods In this section, we consider the space–time discretization that combines the above splitting algorithm with appropriately chosen spatial discretization methods. To this end, we introduce some notation. Assume $${\it{\Omega}}\in\mathbb{R}^d (d=2, 3)$$ is a polygonal domain. Let $$\mathcal{T}_h$$ be a quasi-uniform triangulation or rectangular partition of $${\it{\Omega}}$$ with mesh size $$h$$ and $$\bar{{\it{\Omega}}}=\bigcup_{K\in\mathcal{T}_h}\bar{K}$$. Also, let $$({\mathbf{X}}_h, M_h)$$ be a stable mixed finite element pair, that is, $${\mathbf{X}}_h\subset {\mathbf{H}}^1({\it{\Omega}})$$ and $$M_h\subset L^2({\it{\Omega}})$$ satisfy the inf–sup condition   supvh∈Xh(divvh,φh)‖vh‖H1(Ω)≥β0‖φh‖L2(Ω)∀φh∈M0h:=Mh∩L02(Ω), β0>0. (3.1) A number of stable mixed finite element spaces $$({\mathbf{X}}_h, M_h)$$ have been known in the literature (Brezzi & Fortin, 1992). A well-known example is the following so-called Taylor–Hood element (cf. Bercovier & Pironneau, 1979; Brezzi & Fortin, 1992; Roberts & Thomas, 1991):   Xh ={vh∈C0(Ω¯);vh|K∈P2(K)∀K∈Th},Mh ={φh∈C0(Ω¯);φh|K∈P1(K)∀K∈Th}. In the next subsection, we shall only present the analysis for the Taylor–Hood element, we remark that the analysis can be extended to other stable mixed elements. However, piecewise constant space $$M_h$$ is not recommended because that would result in no rate of convergence for the approximation of the pressure $$p$$ (see Section 3.3). Finite element approximation space $$W_h$$ for $$\eta$$ variable can be chosen independently, any piecewise polynomial space is acceptable provided that $$W_h \supset M_h$$. Especially, $$W_h\subset L^2({\it{\Omega}})$$ can be chosen as a fully discontinuous piecewise polynomial space, although it is more convenient to choose $$W_h$$ to be a continuous (respectively, discontinuous) space if $$M_h$$ is a continuous (respectively, discontinuous) space. The most convenient choice is $$W_h =M_h$$, which will be adopted in the remainder of this article. Recall that $${\mathbf{RM}}$$ denotes the space of the infinitesimal rigid motions (see Section 2), evidently, $${\mathbf{RM}}\subset {\mathbf{X}}_h$$. We now introduce the $$L^2$$-projection $$\mathcal{P}_R$$ from $${\mathbf{L}}^2({\it{\Omega}})$$ to $${\mathbf{RM}}$$. For each $${\mathbf{v}}\in {\mathbf{L}}^2({\it{\Omega}})$$, $$\mathcal{P}_R{\mathbf{v}}_h\in {\mathbf{RM}}$$ is defined by   (PRvh,r)=(vh,r)∀r∈RM. Moreover, we define   Vh:=(I−PR)Xh={vh∈Xh;(vh,r)=0∀r∈RM}. (3.2) It is easy to check that $${\mathbf{X}}_h={\mathbf{V}}_h\bigoplus {\mathbf{RM}}$$. It was proved in Feng & He (2010) that there holds the following alternative version of the above inf–sup condition:   supvh∈Vh(divvh,φh)‖vh‖H1(Ω)≥β1‖φh‖L2(Ω)∀φh∈M0h,β1>0. (3.3) Finite element algorithm (FEA): (i) Compute $${\mathbf{u}}^0_h\in {\mathbf{V}}_h$$ and $$q^0_h\in W_h$$ by   uh0 =Rhu0,ph0=Qhp0,qh0=Qhq0 (q0=divu0),ηh0 =c0ph0+αqh0,ξh0=αph0−λqh0, where $$\mathcal{Q}_h$$ denotes the $$L^2$$-projection operator (see (3.23)). (ii) For $$n=0,1,2, \cdots$$, do the following two steps. Step 1: Solve for $$({\mathbf{u}}^{n+1}_h,\xi^{n+1}_h, \eta^{n+1}_h)\in {\mathbf{V}}_h\times M_h \times W_h$$ such that   μ(ε(uhn+1),ε(vh))−(ξhn+1,divvh)=(f,vh)+⟨f1,vh⟩∀vh∈Vh, (3.4)  κ3(ξhn+1,φh)+(divuhn+1,φh)=κ1(ηhn+θ,φh)∀φh∈Mh; (3.5)   (dtηhn+1,ψh)+1μf(K(∇(κ1ξhn+1+κ2ηhn+1) −ρfg,∇ψh)=(ϕ,ψh)+⟨ϕ1,ψh⟩∀ψh∈Wh, (3.6) where $$\theta=0$$ or $$1$$. Step 2: Update $$p^{n+1}_h$$ and $$q^{n+1}_h$$ by   phn+1 =κ1ξhn+1+κ2ηhn+θ,qhn+1=κ1ηhn+1−κ3ξhn+1. (3.7) Remark 3.1 (a) At each time step, problems (3.4)–(3.5) is a generalized Stokes problem with a mixed boundary condition for $$({\mathbf{u}}_h^{n+1},\xi_h^{n+1})$$. When $$\theta=0$$, the well posedness of the generalized Stokes problem follows from the fact that (3.4)–(3.5) is not a saddle point problem; instead, it results in a symmetric and positive definite linear system because $$\kappa_3>0$$. When $$\theta=1$$, the well posedness of (3.4)–(3.6) is an immediate consequence of the discrete energy law given in Lemma 3.3 below. (b) When $$\theta=0$$, Step 1 consists of two decoupled sub-problems that can be solved independently. On the other hand, when $$\theta=1$$, it is a coupled system, all three unknown functions must be solved together. 3.2 Stability analysis of fully discrete finite element methods The primary goal of this subsection is to derive a discrete energy law, which mimics the PDE energy law (2.30). It turns out that such a discrete energy law only holds if $$h$$ and $${\it{\Delta}} t$$ satisfy the mesh constraint $${\it{\Delta}} t=O(h^2)$$ when $$\theta=0$$, but for all $$h,{\it{\Delta}} t>0$$ when $$\theta=1$$. Before discussing the stability of (FEA), we first show that the numerical solution satisfies all side constraints which are fulfilled by the PDE solution. Lemma 3.2 Let $$\{({\mathbf{u}}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0}$$ be defined by the (FEA), then there hold   (ηhn,1) =Cη(tn)for n=0,1,2,⋯, (3.8)  (ξhn,1) =Cξ(tn−1+θ)for n=1−θ,1,2,⋯, (3.9)  ⟨uhn⋅n,1⟩ =Cu(tn−1+θ)for n=1−θ,1,2,⋯. (3.10) We skip the proof of the above lemma again to save space and refer the interested reader to Feng et al. (xxxx) for its proof. The next lemma establishes an identity that mimics the continuous energy law for the solution of (FEA). Lemma 3.3 Let $$\{({\mathbf{u}}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0}$$ be defined by (FEA), then there holds the following identity:   Jh,θℓ+Sh,θℓ=Jh,θ0for ℓ≥1,θ=0,1, (3.11) where   Jh,θℓ:=12[μ‖ε(uhℓ+1)‖L2(Ω)2+κ2‖ηhℓ+θ‖L2(Ω)2+κ3‖ξhℓ+1‖L2(Ω)2−2(f,uhℓ+1)−2⟨f1,uhℓ+1⟩],Sh,θℓ:=Δt∑n=1ℓ[μΔt2‖dtε(uhn+1)‖L2(Ω)2+1μf(K∇phn+1−Kρfg,∇phn+1) +κ2Δt2‖dtηhn+θ‖L2(Ω)2+κ3Δt2‖dtξhn+1‖L2(Ω)2−(ϕ,phn+1)−⟨ϕ1,phn+1⟩ −(1−θ)κ1Δtμf(Kdt∇ξhn+1,∇phn+1)].phn+1:=κ1ξhn+1+κ2ηhn+θ. Proof. Since the proof for the case $$\theta=1$$ is exactly the same as that of the PDE energy law, so we omit it and leave it to the interested reader to explore. Here, we only consider the case $$\theta=0$$. Based on (3.5), we can define $$\eta_h^{-1}$$ by   κ1(ηh−1,φh)=κ3(ξh0,φh)+(divuh0,φh). (3.12) Setting $${\mathbf{v}}_h=d_t{\mathbf{u}}_h^{n+1}$$ in (3.4), $$\varphi_h=\xi_h^{n+1}$$ in (3.5) and $$\psi_h= p_h^{n+1}$$ in (3.6) after lowing the super-index from $$n+1$$ to $$n$$ on the both sides of (3.6), we get   μ2dt‖ε(uhn+1)‖L2(Ω)2+μ2Δt‖dtε(uhn+1)‖L2(Ω)2 =dt(f,uhn+1)+dt⟨f1,uhn+1⟩+(ξhn+1,divdtuhn+1), (3.13)  κ3(dtξhn+1,ξhn+1)+(divdtuhn+1,ξhn+1)=κ1(dtηhn,ξhn+1), (3.14)   (dtηhn,phn+1)+1μf(K(∇(κ1ξhn+κ2ηhn)−ρfg),∇phn+1) =(ϕ,phn+1)+⟨ϕ1,phn+1⟩. (3.15) The first term on the left-hand side of (3.15) can be rewritten as   (dtηhn,phn+1) =(dtηhn,κ1ξhn+1+κ2ηhn) =κ1(dtηhn,ξhn+1)+κ2Δt2‖dtηhn‖L2(Ω)2+κ22dt‖ηhn‖L2(Ω)2. (3.16) Moreover,   1μf(K∇(κ1ξhn+κ2ηhn)−Kρfg,∇phn+1) =1μf(K∇phn+1−Kρfg,∇phn+1)−κ1Δtμf(Kdt∇ξhn+1,∇phn+1). (3.17)  κ3(dtξhn+1,ξhn+1)=κ32dt‖ξhn+1‖L2(Ω)2+κ3Δt2‖dtξhn+1‖L2(Ω)2. (3.18) Adding (3.13)–(3.15), using (3.16)–(3.18) and applying the summation operator $${\it{\Delta}} t\sum_{n=1}^{\ell}$$ to the both sides of the resulting equation yield the desired equality (3.11). The proof is complete. □ In the case $$\theta=1$$, (3.11) gives the desired solution estimates without any mesh constraint. On the other hand, when $$\theta=0$$, since the last term in the expression of $$S^\ell_{h,\theta}$$ does not have a fixed sign, and hence it needs to be controlled to ensure the positivity of $$S^\ell_{h,\theta}$$. Corollary 3.4 Let $$\{({\mathbf{u}}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0}$$ be defined by (FEA) with $$\theta=0$$, then there holds the following inequality:   Jh,0ℓ+S^h,0ℓ≤Jh,00for ℓ≥1, (3.19) provided that $${\it{\Delta}} t=O(h^2)$$. Here,   S^h,0ℓ:=Δt∑n=1ℓ[μΔt4‖dtε(uhn+1)‖L2(Ω)2+12μf‖∇phn+1‖L2(Ω)2−1μf(ρfg,∇phn+1)+κ2Δt2‖dtηhn‖L2(Ω)2+κ3Δt2‖dtξhn+1‖L2(Ω)2−(ϕ,phn+1)−⟨ϕ1,phn+1⟩]. Proof. By Schwarz inequality and inverse inequality (3.22), we get   κ1K1Δtμf(dt∇ξhn+1,∇phn+1) ≤κ12K22μf‖K12∇ξhn+1−K12∇ξhn‖L2(Ω)2 +12μfK2‖K12∇phn+1‖L2(Ω)2 ≤K22κ122μfh2‖ξhn+1−ξhn‖L2(Ω)2+12μf‖∇phn+1‖L2(Ω)2. (3.20) To bound the first term on the right-hand side of (3.20), we appeal to the inf–sup condition and get   ‖ξhn+1−ξhn‖L2 ≤1β1supvh∈Vh(divvh,ξhn+1−ξhn)‖∇vh‖L2(Ω) ≤μβ1supvh∈Vh(ε(un+1−un),ε(vh))‖∇vh‖L2(Ω) ≤μβ1Δt‖dtε(uhn+1)‖L2. (3.21) Substituting (3.21) into (3.20) and combining it with (3.11) imply (3.19) provided that $${\it{\Delta}} t \leq C(\mu_f\beta_1^2)(2\mu c_1^2\kappa_1^2)^{-1} h^2$$. The proof is complete. □ 3.3 Convergence analysis The goal of this section is to analyse the fully discrete finite element algorithm proposed in the previous subsection. Precisely, we shall derive optimal order error estimates for (FEA) in both $$L^\infty(0,T;L^2({\it{\Omega}}))$$ and $$L^2(0,T; H^1({\it{\Omega}}))$$-norm. To this end, we first list some facts, which are well known in the literature (Brezzi & Fortin, 1992; Brenner & Scott, 2008) about finite element functions. We first recall the following inverse inequality for polynomial functions Ciarlet (1978):   ‖∇φh‖L2(K)≤c1h−1‖φh‖L2(K)∀φh∈Pr(K),K∈Th. (3.22) For any $$\varphi\in L^2({\it{\Omega}})$$, we define its $$L^2$$-projection $$\mathcal{Q}_h: L^2\rightarrow W_h$$ as   (Qhφ,ψh)=(φ,ψh)ψh∈Wh. (3.23) It is well known that the projection operator $$\mathcal{Q}_h: L^2\rightarrow W_h$$ satisfies (cf. Brenner & Scott, 2008), for any $$\varphi\in H^s({\it{\Omega}}) (s\geq1)$$,   ‖Qhφ−φ‖L2(Ω)+h‖∇(Qhφ−φ)‖L2(Ω)≤Chℓ‖φ‖Hℓ(Ω)ℓ=min{2,s}. (3.24) We like to point out that when $$W_h\not\subset H^1({\it{\Omega}})$$, the second term on the left-hand side of (3.24) has to be replaced by the broken $$H^1$$-norm. Next, for any $$\varphi\in H^1({\it{\Omega}})$$, we define its elliptic projection $$\mathcal{S}_h\varphi$$ by   (K∇Shφ,∇φh) =(K∇φ,∇φh)∀φh∈Wh, (3.25)  (Shφ,1) =(φ,1). (3.26) It is well known that the projection operator $$\mathcal{S}_h: H^1({\it{\Omega}})\rightarrow W_h$$ satisfies (cf. Brenner & Scott, 2008), for any $$\varphi\in H^s({\it{\Omega}}) (s>1)$$,   ‖Shφ−φ‖L2(Ω)+h‖∇(Shφ−φ)‖L2(Ω)≤Chℓ‖φ‖Hℓ(Ω)ℓ=min{2,s}. (3.27) Finally, for any $${\mathbf{v}}\in {\mathbf{H}}^1_\perp({\it{\Omega}})$$, we define its elliptic projection $$\mathcal{R}_h{\mathbf{v}}$$ by   (ε(Rhv),ε(wh))=(ε(v),ε(wh))wh∈Vh. (3.28) It is easy to show that the projection $$\mathcal{R}_h{\mathbf{v}}$$ satisfies (cf. Brenner & Scott, 2008), for any $${\mathbf{v}}\in {\mathbf{H}}^1_\perp({\it{\Omega}})\cap {\mathbf{H}}^s({\it{\Omega}}) (s>1)$$,   ‖Rhv−v‖L2(Ω)+h‖∇(Rhv−v)‖L2(Ω)≤Chm‖v‖Hm(Ω)m=min{3,s}. (3.29) To derive error estimates, we introduce the following error notation   Eun :=u(tn)−uhn,Eξn:=ξ(tn)−ξhn,Eηn:=η(tn)−ηhn,Epn :=p(tn)−phn,Eqn:=q(tn)−qhn. It is easy to check that   Epn=κ1Eξn+κ2Eηn,Eqn=κ3Eξn+κ1Eηn. (3.30) Also, we denote   Eun =u(tn)−Rh(u(tn))+Rh(u(tn))−uhn:=Λun+Θun,Eξn =ξ(tn)−Sh(ξ(tn))+Sh(ξ(tn))−ξhn:=Λξn+Θξn,Eηn =η(tn)−Sh(η(tn))+Sh(η(tn))−ηhn:=Ληn+Θηn,Epn =p(tn)−Qh(p(tn))+Qh(p(tn))−phn:=Λpn+Θpn. Lemma 3.5 Let $$\{ ({\mathbf{u}}_h^n, \xi_h^n, \eta_h^n) \}_{n\geq0}$$ be generated by the (FEA) and $${\it{\Lambda}}_{{\mathbf{u}}}^n, {\it{\Theta}}_{{\mathbf{u}}}^n, {\it{\Lambda}}_{\xi }^n, {\it{\Theta}}_{\xi }^n, {\it{\Lambda}}_{\eta }^n$$ and $${\it{\Theta}}_{\eta }^n$$ be defined as above. Then there holds the following identity:   Ehℓ+Δt∑n=1ℓ[1μf(K∇Θ^pn+1−Kρfg,Θ^pn+1) +Δt2(μ‖dtε(Θun+1)‖L2(Ω)2+κ2‖dtΘηn+θ‖L2(Ω)2+κ3‖dtΘξn+1‖L2(Ω)2)] =Eh0+Δt∑n=1ℓ[(Λξn+1,divdtΘun+1)−(divdtΛun+1,Θξn+1)] +(Δt)2∑n=1ℓ(dt2ηh(tn+1),Θξn+1)+Δt∑n=1ℓ(Rhn+1,Θ^pn+1) +(1−θ)(Δt)2∑n=1ℓκ1μf(Kdt∇Θξn+1,∇Θ^pn+1), (3.31) where   Θ^pn+1 :=κ1∇Θξn+1+κ2∇Θηn+θ, (3.32)  Ehℓ :=12[μ‖ε(Θuℓ+1)‖L2(Ω)2+κ2‖Θηℓ+θ‖L2(Ω)2+κ3‖Θξℓ+1‖L2(Ω)2], (3.33)  Rhn+1 :=−1Δt∫tntn+1(s−tn)ηtt(s)ds. (3.34) Proof. Subtracting (3.4) from (2.24), (3.5) from (2.25), (3.6) from (2.26), respectively, we get the following error equations:   μ(ε(Eun+1),ε(vh))−(Eξn+1,divvh)=0∀vh∈Vh, (3.35)  κ3(Eξn+1,φh)+(divEun+1,φh)  =κ1(Eηn+θ,φh)+Δt(dtη(tn+1),φh)∀φh∈Mh, (3.36)   (dtEηn+1,ψh)+1μf(K∇(κ1Eξn+1+κ2Eηn+1)−Kρfg,∇ψh) =(Rhn+1,ψh)∀ψh∈Wh, (3.37)  Eu0=0,Eξ0=0,Eη−1=0. (3.38) Using the definition of the projection operators $$\mathcal{Q}_h, \mathcal{S}_h, \mathcal{R}_h$$, we have   μ(ε(Θun+1),ε(vh))−(Θξn+1,divvh)=(Λξn+1,divvh)∀vh∈Vh, (3.39)  κ3(Θξn+1,φh)+(divΘun+1,φh)=κ1(Θηn+θ,φh) −(divΛun+1,φh)+Δt(dtη(tn+1),φh)∀φh∈Mh, (3.40)   (dtΘηn+1,ψh)+1μf(K∇(κ1Θξn+1+κ2Θηn+1)−Kρfg,∇ψh) =(Rhn+1,ψh)∀ψh∈Wh, (3.41)  Eu0=0,Eξ0=0,Eη−1=0. (3.42) Equation (3.31) follows from setting $${\mathbf{v}}_h=d_t{\it{\Theta}}^{n+1}_{{\mathbf{u}}}$$ in (3.39), $$\varphi_h={\it{\Theta}}^{n+1}_\xi$$ (after applying the difference operator $$d_t$$ to the equation (3.40)), $$\psi_h=\hat{{\it{\Theta}}}_p^{n+1}=\kappa_1{\it{\Theta}}_\xi^{n+1} +\kappa_2{\it{\Theta}}_\eta^{n+\theta}$$ in (3.41), adding the resulting equations and applying the summation operator $${\it{\Delta}} t\sum^{\ell}_{n=1}$$ to both sides. □ Theorem 3.6 Let $$\{(u_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0}$$ be defined by (FEA), then there holds the error estimate for $$\ell \leq N$$   max0≤n≤ℓ[μ‖ε(Θun+1)‖L2(Ω)+κ2‖Θηn+θ‖L2(Ω)+κ3‖Θξn+1‖L2(Ω)] +[Δt∑n=0ℓ1μf‖Θ^pn+1‖L22]12≤C1(T)Δt+C2(T)h2, (3.43) provided that $${\it{\Delta}} t =O(h^2)$$ when $$\theta=0$$ and $${\it{\Delta}} t >0$$ when $$\theta=1$$. Here   C1(T) =C‖ηt‖L2((0,T);L2(Ω))2+C‖(η)tt‖L2((0,T);H1(Ω)′), (3.44)  C2(T) =C‖ξ‖L∞((0,T);H2(Ω))+C‖ξt‖L2((0,T);H2(Ω)) +C‖div(u)t‖L2((0,T);H2(Ω)). (3.45) Proof. To derive the above inequality, we need to bound each term on the right-hand side of (3.31). Using the fact $${\it{\Theta}}_{{\mathbf{u}}}^0=\mathbf{0}$$, $${\it{\Theta}}_\xi^0=0$$ and $${\it{\Theta}}_\eta^{-1}=0$$, we have   Ehℓ+Δt∑n=1ℓ[1μf(K∇Θ^pn+1−Kρfg,∇Θ^pn+1) +Δt2(μ‖dtε(Θun+1)‖L2(Ω)2+κ2‖dtΘηn+θ‖L2(Ω)2+κ3‖dtΘξn+1‖L2(Ω)2)] =(Δt)2∑n=1ℓ(dt2η(tn+1),Θξn+1)+Δt∑n=1ℓ(Rhn+1,Θ^pn+1)+μ2‖ε(Θu1)‖L2(Ω)2 +Δt∑n=1ℓ[(Λξn+1,divdtΘun+1)−(divdtΛun+1,Θξn+1)] +(1−θ)(Δt)2∑n=1ℓκ1μf(Kdt∇Θξn+1,∇Θ^pn+1). (3.46) We now estimate each term on the right-hand side of (3.46). To bound the first term on the right-hand side of (3.46), we first use the summation by parts formula and $$d_t\eta_h(t_0)=0$$ to get   ∑n=0ℓ(dt2η(tn+1),Θξn+1)=1Δt(dtη(tl+1),Θξl+1)−∑n=1ℓ(dtη(tn),dtΘξn+1). (3.47) Now, we bound each term on the right-hand side of (3.47) as follows:   1Δt(dtη(tℓ+1),Θξℓ+1)≤1Δt‖dtη(tℓ+1)‖L2(Ω)‖Θξℓ+1‖L2(Ω) ≤1Δt‖ηt‖L2((tℓ,tℓ+1);L2(Ω))⋅1β1supvh∈Vhμ(ε(Θuℓ+1),ε(vh))−(Λξℓ+1,divvh)‖∇vh‖L2(Ω) ≤Cμβ1Δt‖ηt‖L2((tℓ,tℓ+1);L2(Ω))[‖ε(Θuℓ+1)‖L2(Ω)+1μ‖Λξℓ+1‖L2(Ω)] ≤μ4(Δt)2‖ε(Θuℓ+1)‖L2(Ω)2+Cμβ12‖ηt‖L2((tℓ,tℓ+1);L2(Ω))2+Cβ12‖Λξℓ+1‖L2(Ω)2, (3.48)   ∑n=1ℓ(dtη(tn),dtΘξn+1)≤∑n=1ℓ‖dtη(tn)‖L2(Ω)‖dtΘξn+1‖L2(Ω) ≤∑n=1ℓ‖dtη(tn)‖L2(Ω)⋅1β1supvh∈Vhμ(dtε(Θun+1),ε(vh))−(dtΛξn+1,divvh)‖∇vh‖L2(Ω) ≤Cμβ1∑n=1ℓ‖dtη(tn)‖L2(Ω)[‖dtε(Θun+1)‖L2(Ω)+1μ‖dtΛξn+1‖L2(Ω)] ≤∑n=1ℓ[μ4‖dtε(Θun+1)‖L2(Ω)2+Cβ12‖dtΛξn+1‖L2(Ω)]+Cμβ12‖ηt‖L2((0,T);L2(Ω))2. (3.49) The second term on the right-hand side of (3.46) can be bounded as   |(Rhn+1,Θ^pn+1)| ≤‖Rhn+1‖H1(Ω)′‖∇Θ^pn+1‖L2(Ω) ≤K14μf‖∇Θ^pn+1‖L2(Ω)2+μfK1‖Rhn+1‖H1(Ω)′2 ≤K14μf‖∇Θ^pn+1‖L2(Ω)2+μfΔt3K1‖ηtt‖L2((tn,tn+1);H1(Ω)′)2, (3.50) where we have used the fact that   ‖Rhn+1‖H1(Ω)′2≤Δt3∫tntn+1‖ηtt‖H1(Ω)′2dt. The fourth term on the right-hand side of (3.46) can be bounded by   Δt∑n=1ℓ[(Λξn+1,divdtΘun+1)−(divdtΛun+1,Θξn+1)] =(Λξℓ+1,divΘuℓ+1)−Δt∑n=1ℓ[(dtΛξn+1,divΘun)+(divdtΛun+1,Θξn+1)]tag ≤12‖Λξℓ+1‖L2(Ω)2+12‖divΘuℓ+1‖L2(Ω)2+12Δt∑n=1ℓ[‖dtΛξn+1‖L2(Ω)2 +C‖ε(Θun+1)‖L2(Ω)2+‖divdtΛun+1‖L2(Ω)2+‖Θξn+1‖L2(Ω)2]. (3.51) When $$\theta=0$$ we also need to bound the last term on the right-hand side of (3.46), which is carried out below.    ∑n=1ℓ(dt∇Θξn+1,∇Θ^pn+1)≤∑n=1ℓ‖dtΘξn+1‖L2(Ω)‖∇Θ^pn+1‖L2(Ω) ≤∑n=1ℓsupvh∈Vhμ(dtε(Θun+1),ε(vh))−(dtΛξn+1,divvh)‖∇vh‖L2(Ω)‖∇Θ^pn+1‖L2(Ω) ≤∑n=1ℓ[μ2κ1Δth2β12‖dtε(Θun+1)‖L22+κ1Δth2β12‖dtΛξn+1‖L2+κ1−14Δt‖∇Θ^pn+1‖L2(Ω)2]. (3.52) Substituting (3.47)–(3.52) into (3.46) and rearranging terms, we get   μ‖ε(Θuℓ+1)‖L2(Ω)2+κ2‖Θηℓ+θ‖L2(Ω)2+κ3‖Θξℓ+1‖L2(Ω)2 +Δt∑n=1ℓ1μf‖K∇Θ^pn+1‖L2(Ω)2 ≤4μf(Δt)2K1‖ηtt‖L2((0,T);H1(Ω)′)2+4μ(Δt)2β12‖ηt‖L2((0,T);L2(Ω))2, +‖Λξℓ+1‖L2(Ω)2+Δt∑n=1ℓ‖dtΛξn+1‖L2(Ω)2+Δt∑n=1ℓ‖divdtΛun+1‖L2(Ω)2, (3.53) provided that $${\it{\Delta}} t\leq \mu_f\beta_1^2(4\mu\kappa_1^2K)^{-1} h^2$$ when $$\theta=0$$, but it holds for all $${\it{\Delta}} t>0$$ when $$\theta=1$$. Hence, (3.43) follows from using the approximation properties of the projection operators $$\mathcal{Q}_h, \mathcal{R}_h$$ and $$\mathcal{S}_h$$. The proof is complete. □ We conclude this section by stating the main theorem of the section. Theorem 3.7 The solution of the (FEA) satisfies the following error estimates:    max0≤n≤N[μ‖∇(u(tn)−uhn)‖L2(Ω)+κ2‖η(tn)−ηhn‖L2(Ω) +κ3‖ξ(tn)−ξhn‖L2(Ω)]≤C^1(T)Δt+C^2(T)h2. (3.54)   [Δt∑n=0N1μf‖∇(p(tn)−phn)‖L2(Ω)2]12≤C^1(T)Δt+C^2(T)h, (3.55) provided that $${\it{\Delta}} t=O(h^2)$$ when $$\theta=0$$ and $${\it{\Delta}} t>0$$ when $$\theta=1$$. Here,   C^1(T) :=C1(T),C^2(T) :=C2(T)+‖ξ‖L∞((0,T);H2(Ω))+‖η‖L∞((0,T);H2(Ω)) +‖∇u‖L∞((0,T);H2(Ω)). Proof. The above estimates follow immediately from an application of the triangle inequality on   u(tn)−uhn =Λun+Θun,ξ(tn)−ξhn=Λξn+Θξn,η(tn)−ηhn =Ληn+Θηn,p(tn)−phn=Λ^pn+Θ^pn and appealing to (3.24), (3.29) and Theorem 3.6. □ 4. Numerical experiments In this section, we shall present three two-dimensional numerical experiments to validate theoretical results for the proposed numerical methods, to numerically examine the performances of the approach and methods as well as to compare them with existing methods in the literature on two benchmark problems. One of these two problems was used to demonstrate the ‘locking phenomenon’ in Phillips & Wheeler (2009). Our numerical result shows that such a ‘locking phenomenon’ does not occur in our numerical methods, it confirms the fact that our approach and methods have a built-in mechanism to prevent the ‘locking phenomenon’. Test 1. Let $${\it{\Omega}}=[0,1]\times [0,1]$$, $${\it{\Gamma}}_1=\{(1,x_2); 0\leq x_2\leq 1\}$$, $${\it{\Gamma}}_2=\{(x_1,0); 0\leq x_1\leq 1\}$$, $${\it{\Gamma}}_3=\{(0,x_2); 0\leq x_2\leq 1\}$$, $${\it{\Gamma}}_4=\{(x_1,1); 0\leq x_1\leq 1\}$$ and $$T=0.001$$. We consider problems (2.26)–(2.28) with following source functions:   f =−(λ+μ)t(1,1)T+αcos⁡(x1+x2)et(1,1)T,ϕ =(c0+2κμf)sin⁡(x1+x2)et+α(x1+x2) and the following boundary and initial conditions:   p =sin⁡(x1+x2)eton ∂ΩT,u1 =12x12ton Γj×(0,T),j=1,3,u2 =12x22ton Γj×(0,T),j=,2,4,σn−αpn =f1on ∂ΩT,u(x,0)=0,p(x,0) =sin⁡(x1+x2)in Ω, where   f1(x,t)=μ(x1n1,x2n2)Tt+λ(x1+x2)(n1,n2)Tt−αsin⁡(x1+x2)(n1,n2)Tet. It is easy to check that the exact solution for this problem is   u(x,t)=t2(x12,x22)T,p(x,t)=sin⁡(x1+x2)et. We note that the boundary conditions used above are not pure Neumann conditions; instead, they are mixed Dirichlet–Neumann conditions. As pointed out in Remark 2.3 (c), the approach and methods of this article also apply to this case, the only change is to replace the test and trial space $${\mathbf{H}}^1_\perp({\it{\Omega}})$$ by $${\mathbf{H}}^1({\it{\Omega}})$$ with some appropriately built-in Dirichlet boundary condition in Definition 2.2. The goal of doing this test problem is to compute the order of the exact errors and to show that the theoretical error bounds proved in the previous section are sharp. Table 1 lists the computed $$L^\infty(0,T; L^2({\it{\Omega}}))$$ and $$L^2(0,T;H^1({\it{\Omega}}))$$-norm errors and the convergence rates with respect to $$h$$ at the terminal time $$T$$. In the test, $${\it{\Delta}} t=10^{-5}$$ is used so that the time error is negligible. Evidently, the spatial rates of convergence are consistent with that proved in the convergence theorem. Table 1. Spatial errors and convergence rates of Test 1    $$L^\infty(L^2)$$ error  $$L^\infty(L^2)$$ order  $$L^2(H^1)$$ error  $$L^2(H^1)$$ order  $$h=0.16$$  2.0789e-3     5.5045e-2     $$h=0.08$$  5.9674e-4  1.8006  2.9431e-2  0.9032  $$h=0.04$$  1.6227e-4  1.8787  1.5332e-2  0.9408  $$h=0.02$$  4.0971e-5  1.9857  7.6968e-3  0.9942     $$L^\infty(L^2)$$ error  $$L^\infty(L^2)$$ order  $$L^2(H^1)$$ error  $$L^2(H^1)$$ order  $$h=0.16$$  2.0789e-3     5.5045e-2     $$h=0.08$$  5.9674e-4  1.8006  2.9431e-2  0.9032  $$h=0.04$$  1.6227e-4  1.8787  1.5332e-2  0.9408  $$h=0.02$$  4.0971e-5  1.9857  7.6968e-3  0.9942  Figures 1 and 2 show, respectively, the surface plot of the computed pressure $$p$$ at the terminal time $$T$$ and the color plot of both the computed pressure $$p$$ and displacement $$\mathbf{u}$$ with mesh parameters $$h=0.02$$ and $${\it{\Delta}} t=10^{-5}$$. They coincide with the exact solution on the same space–time resolution. Fig. 1. View largeDownload slide Test 1: surface plot of the computed pressure $$p$$ at the terminal time $$T$$. Fig. 1. View largeDownload slide Test 1: surface plot of the computed pressure $$p$$ at the terminal time $$T$$. Fig. 2. View largeDownload slide Test 1: computed pressure $$p$$ (color plot) and displacement $$\mathbf{u}$$ (arrow plot) at $$T$$. Fig. 2. View largeDownload slide Test 1: computed pressure $$p$$ (color plot) and displacement $$\mathbf{u}$$ (arrow plot) at $$T$$. Test 2. In this test, we consider so-called Barry–Mercer’s problem, which is a Benchmark test problem for the poroelasticity model (2.26)–(2.28) (cf. Phillips & Wheeler, 2007b, 2009 and the references therein). Again, $${\it{\Omega}}=[0,1]\times [0,1]$$, but $$T=1$$. Barry–Mercer’s problem assumes no source, that is, $$\mathbf{f}\equiv 0$$ and $$\phi\equiv 0$$, and takes the following boundary conditions:   p =0on Γj×(0,T),j=1,3,4,p =p2on Γj×(0,T),j=2,u1 =0on  Γj×(0,T),j=1,3,u2 =0on  Γj×(0,T),j=2,4,σn−αpn =f1:=(0,αp)Ton ∂ΩT, where   p2(x1,t)={sin⁡twhen x∈[0.2,0.8)×(0,T),0others.  The boundary segments $${\it{\Gamma}}_j, j=1,2,3,4$$, which are defined in Test 1, and the above boundary conditions are depicted in Fig. 3. Also, the initial conditions for Barry–Mercer’s problem are $$\mathbf{u}(x,0) \equiv \mathbf{0}$$ and $$p(x,0)\equiv 0$$. We remark that Barry–Mercer’s problem has a unique solution that is given by an infinite series (cf. Phillips & Wheeler, 2009). Fig. 3. View largeDownload slide Test 2: boundary conditions. Fig. 3. View largeDownload slide Test 2: boundary conditions. Figures 4 and 5 display, respectively, the computed pressure $$p$$ (surface plot) and the computed displacement $$\mathbf{u}$$ (arrow plot). We note that the arrows near the boundary match very well with those on the boundary. Our numerical solution approximates the exact solution of Barry–Mercer’s problem very well and does not produce any oscillation in computed pressure. Fig. 4. View largeDownload slide Test 2: surface plot of the computed pressure $$p$$ at the terminal time $$T$$. Fig. 4. View largeDownload slide Test 2: surface plot of the computed pressure $$p$$ at the terminal time $$T$$. Fig. 5. View largeDownload slide Test 2: computed pressure $$p$$ (color plot) and displacement (arrow plot) at $$T$$. Fig. 5. View largeDownload slide Test 2: computed pressure $$p$$ (color plot) and displacement (arrow plot) at $$T$$. Test 3. This test problem is taken from Phillips & Wheeler (2009). Again, we consider problems (2.26)–(2.28) with $${\it{\Omega}}=[-0.5,0.5]\times [-0.5,0.5]$$. Let $${\it{\Gamma}}_j$$ be same as in Test 1 and $$c_0=0, E=10^5, \nu=0.49, \mu=33557, \lambda=1.6443\times10^6$$ and $$T=0.001$$. There is no source, that is, $$\mathbf{f}\equiv 0$$ and $$\phi\equiv 0$$. The boundary conditions are taken as   −κμf(∇p−ρfg)⋅n =0on ∂ΩT,u =0on Γ3×(0,T),σn−αpn =f1on Γj×(0,T),j=1,2,4, where $$\mathbf{f}_1=(f_1^1, f_1^2)$$ and   f11≡0on ∂ΩT,f12={0on Γj×(0,T),j=1,2,3,−1on Γ4×(0,T).  The computational domain $${\it{\Omega}}$$ and the above boundary conditions are depicted in Fig. 6. Also, the zero initial conditions are assigned for both $$\mathbf{u}$$ and $$p$$ in this test. Fig. 6. View largeDownload slide Test 3: boundary conditions. Fig. 6. View largeDownload slide Test 3: boundary conditions. Figures 7 and 8 display, respectively, the surface and color plot of the computed pressure, the arrow plot of the displacement vector and the deformation of the whole $${\it{\Omega}}$$. There is no oscillation in the computed pressure, and the arrows near the boundary match very well with arrows on the boundary, even when the Poisson ratio $$\nu=0.49$$, in which case the material is nearly incompressible. Note the finite element solution of the pressure deteriorates in the nearly incompressible case in Phillips & Wheeler (2009). Fig. 7. View largeDownload slide Test 3: computed pressure $$p$$: surface plot (left) and color plot (right). Fig. 7. View largeDownload slide Test 3: computed pressure $$p$$: surface plot (left) and color plot (right). Fig. 8. View largeDownload slide Test 3: arrow plot of the computed displacement (left) and deformation of $${\it{\Omega}}$$ (right). Fig. 8. View largeDownload slide Test 3: arrow plot of the computed displacement (left) and deformation of $${\it{\Omega}}$$ (right). We remark that the ‘locking phenomenon’ was observed in the simulation of Phillips & Wheeler (2009) at $$T=0.001$$ for this problem, namely, the computed pressure exhibits some oscillation at $$T=0.001$$. The reason for the locking phenomenon was explained as follows: when time step $${\it{\Delta}} t$$ is small, the displacement vector $$\mathbf{u}$$ is almost divergence free in the short time, while the numerical solution does not observe this nearly divergence free property, which results in the locking. However, at later times, the displacement vector is no longer divergence free, so no locking exists at later times. It is clear that our numerical solution does not exhibit the locking phenomenon at $$T=0.001$$. This is because our multiphysics reformulation weakly imposes the condition $$\mbox{div } \mathbf{u}=q$$, and hence $$\mathbf{u}$$ automatically becomes nearly divergence free when $$q\approx 0$$ for $$0<t<<1$$. Moreover, the pressure $$p$$ is not a primary variable anymore in our reformulation; instead, $$p$$ becomes a derivative variable and it is computed using the new primary variables $$\xi$$ and $$\eta$$. Therefore, our numerical methods are insensitive to the regularity of $$p$$. Funding NSF grants (DMS-1016173 and DMS-1318486 to X.F. and DMS-1016173 and DMS-1318486 to Y.L.); National Natural Science Foundation of China grant (#10901047 to Z.G.). References Bercovier J. & Pironneau O. ( 1979) Error estimates for finite element solution of the Stokes problem in the primitive variables. Numer. Math. , 33, 211– 224. Google Scholar CrossRef Search ADS   Biot M. ( 1955) Theory of elasticity and consolidation for a porous anisotropic media. J. Appl. Phys. , 26, 182– 185. Google Scholar CrossRef Search ADS   Brenner S. C. ( 1993) A nonconforming mixed multigrid method for the pure displacement problem in planar linear elasticity. SIAM J. Numer. Anal ., 30, 116– 135. Google Scholar CrossRef Search ADS   Brenner S. C. & Scott L. R. ( 2008) The Mathematical Theory of Finite Element Methods , 3rd edn. New York: Springer. Google Scholar CrossRef Search ADS   Brezzi F. & Fortin M. ( 1992) Mixed and Hybrid Finite Element Methods.  New York: Springer. Ciarlet P. G. ( 1978) The Finite Element Method for Elliptic Problems . Amsterdam: North-Holland. Coussy O. ( 2004) Poromechanics , Hoboken, NJ: Wiley & Sons. Dauge M. ( 1988) Elliptic Boundary Value Problems on Corner Domains . Lecture Notes in Math., vol. 1341. Berlin: Springer. Google Scholar CrossRef Search ADS   Dautray R. & Lions J. L. ( 1990) Mathematical Analysis and Numerical Methods for Science and Technology,  vol. 1. New York: Springer. Doi M. & Edwards S. F. ( 1986) The Theory of Polymer Dynamics . Oxford: Clarendon Press. Feng X., Ge Z. & Li Y. ( 2014) Multiphysics finite element methods for a poroelasticity model. arXiv:1411.7464 [matth.NA]. Feng X. & He Y. ( 2010) Fully discrete finite element approximations of a polymer gel model. SIAM J. Numer. Anal ., 48, 2186– 2217. Google Scholar CrossRef Search ADS   Girault V. & Raviart P. A. ( 1981) Finite Element Method for Navier-Stokes Equations: Theory and Algorithms . Berlin, Heidelberg, New York: Springer. Hamley I. ( 2007) Introduction to Soft Matter.  Hoboken, NJ: John Wiley & Sons. Lee J. J. ( 2014) Guaranteed locking-free finite element methods for Biot’s consolidation model in poroelasticity, arXiv:1403.7038v2 [math.NA]. Lee J. J. ( 2016) Robust error analysis of coupled mixed methods for Biot’s consolidation model. J. Sci. Comput ., 69, 610– 632. Google Scholar CrossRef Search ADS   Lee J. J., Mardal K. & Winther R. ( 2017) Parameter-robust discretization and preconditioning of Biot’s consolidation model. SIAM J. Sci. Comput. , 39, A1– A24. Google Scholar CrossRef Search ADS   Murad M. A. & Loula A. F. D. ( 1992) Improved accuracy in finite element analysis of Biot’s consolidation problem. Comput. Methods in Appl. Mech. and Engr ., 95, 359– 382. Google Scholar CrossRef Search ADS   Oyarzúa R. & Ruiz-Baier R. ( 2016) Locking-free finite element methods for poroelasticity. SIAM J. Numer. Anal. , 54, 2951– 2973. Google Scholar CrossRef Search ADS   Phillips P. J. & Wheeler M. F. ( 2007a) A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous in time case. Comput. Geosci ., 11, 131– 144. Google Scholar CrossRef Search ADS   Phillips P. J. & Wheeler M. F. ( 2007b) A coupling of mixed and continuous Galerkin finite element methods for poroelasticity II: the discrete in time case. Comput. Geosci ., 11, 145– 158. Google Scholar CrossRef Search ADS   Phillips P. J. & Wheeler M. F. ( 2009) Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach. Comput. Geosci ., 13, 1– 15. Google Scholar CrossRef Search ADS   Roberts J. E. & Thomas J. M. ( 1991) Mixed and hybrid methods. Handbook of Numerical Analysis  ( Ciarlet P. G. & Lions J. L. eds), vol. II. New York: North-Holland, pp. 523– 639. Google Scholar CrossRef Search ADS   Schowalter R. E. ( 2000) Diffusion in poro-elastic media. J. Math. Anal ., 251, 310– 340. Google Scholar CrossRef Search ADS   Schreyer-Bennethum L. ( 2007) Theory of flow and deformation of swelling porous materials at the macroscale. Computers and Geotechnics , 34, 267– 278. Google Scholar CrossRef Search ADS   Tanaka T. & Fillmore D. J. ( 1979) Kinetics of swelling of gels. J. Chem. Phys ., 70, 1214. Google Scholar CrossRef Search ADS   Temam R. ( 2000) Navier-Stokes Equations: Theory and Numerical Analysis . AMS Chelsea Publishing, Providence, RI: American Mathematical Society. Terzaghi K. ( 1943) Theoretical Soil Mechanics . New York: John Wiley & Sons. Google Scholar CrossRef Search ADS   Yamaue T. & Doi M. ( 2004) Swelling dynamics of constrained thin-plate under an external force. Phys. Rev. E  70, 011401. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations