A high-order discretization of nonlinear poroelasticity

TL;DR

This paper presents a novel high-order discretization method for nonlinear poroelasticity that is stable, flexible on complex meshes, and efficient, enabling accurate simulations of fluid flow in deformable porous media.

Contribution

It introduces a nonconforming high-order discretization combining Hybrid High-Order and discontinuous Galerkin methods for nonlinear poroelasticity, supporting arbitrary orders and complex geometries.

Findings

The method is stable and accurate in 2D and 3D.

It handles complex meshes including polyhedral elements.

It allows for reduced computational cost through static condensation.

Abstract

In this work we construct and analyze a nonconforming high-order discretization method for the quasi-static single-phase nonlinear poroelasticity problem describing Darcean flow in a deformable porous medium saturated by a slightly compressible fluid. The nonlinear elasticity operator is discretized using a Hybrid High-Order method, while the Darcy operator relies on a Symmetric Weighted Interior Penalty discontinuous Galerkin scheme. The method is valid in two and three space dimensions, delivers an inf-sup stable discretization on general meshes including polyhedral elements and nonmatching interfaces, supports arbitrary approximation orders, and has a reduced cost thanks to the possibility of statically condensing a large subset of the unknowns for linearized versions of the problem. Moreover, the proposed construction can handle both nonzero and vanishing specific storage coefficients.

Tables

Table 1. Table 1: Convergence results on the Cartesian and Voronoi meshes for k = 1 𝑘 1 k=1 . OCV stands for order of convergence.

$h$	${\left({\sum }_{n=1}^N\tau {\Vert {\overrightarrow{\underset{¯}{u}}}_h^n-{\underset{¯}{I}}_h^k{\overline{\overrightarrow{u}}}^n\Vert }_{\text{matr}\epsilon ,h}^2\right)}^{\frac{1}{2}}$	OCV	${\left({\sum }_{n=1}^N\tau {\Vert p_h^n-{\pi }_h^k{\overline{p}}^n\Vert }_{\Omega }^2\right)}^{\frac{1}{2}}$	OCV
Cartesian mesh family
$6.25\cdot {10}^{-2}$	$3.10\cdot {10}^{-2}$	—	$0.39$	—
$3.12\cdot {10}^{-2}$	$8.52\cdot {10}^{-3}$	$1.86$	$9.65\cdot {10}^{-2}$	$2.00$
$1.56\cdot {10}^{-2}$	$2.22\cdot {10}^{-3}$	$1.94$	$2.44\cdot {10}^{-2}$	$1.98$
$7.81\cdot {10}^{-3}$	$5.61\cdot {10}^{-4}$	$1.99$	$6.18\cdot {10}^{-3}$	$1.99$
$3.91\cdot {10}^{-3}$	$1.41\cdot {10}^{-4}$	$2.00$	$1.56\cdot {10}^{-3}$	$1.99$
Voronoi mesh family
$6.50\cdot {10}^{-2}$	$3.28\cdot {10}^{-2}$	—	$0.27$	—
$3.15\cdot {10}^{-2}$	$8.48\cdot {10}^{-3}$	$1.87$	$6.58\cdot {10}^{-2}$	$1.96$
$1.61\cdot {10}^{-2}$	$2.20\cdot {10}^{-3}$	$2.01$	$1.63\cdot {10}^{-2}$	$2.08$
$9.09\cdot {10}^{-3}$	$5.72\cdot {10}^{-4}$	$2.36$	$4.24\cdot {10}^{-3}$	$2.36$
$4.26\cdot {10}^{-3}$	$1.42\cdot {10}^{-4}$	$1.83$	$1.05\cdot {10}^{-3}$	$1.84$

Equations

$$-\vec{\nabla }\cdot \text{matr}\sigma (\cdot ,{\vec{\nabla }}_{\mathrm{s}}\vec{u})+\alpha \vec{\nabla }$$

$$C_0{\mathrm{d}}_{t}p+\alpha {\mathrm{d}}_t(\vec{\nabla }\cdot \vec{u})-\vec{\nabla }\cdot (\text{matr}\kappa (\cdot )\vec{\nabla }p)$$

$$\underline{\kappa }|\vec{\xi }{|}^2\le \text{matr}\kappa (\vec{x})\vec{\xi }\cdot \vec{\xi }\le \overline{\kappa }|\vec{\xi }{|}^2.$$

$$\vec{u}$$

$$(\text{matr}\kappa (\cdot )\vec{\nabla }p)\cdot \vec{n}$$

$$C_{0}p(\cdot ,0)+\vec{\nabla }\cdot \vec{u}(\cdot ,0)={\varphi }^0(\cdot ).$$

$${\int }_{\Omega }{\varphi }^0=0,{\int }_{\Omega }g(\cdot ,t)=0,\text{ and }{\int }_{\Omega }p(\cdot ,t)=0\forall t\in (0,t_{\mathrm{F}}).$$

$$\Vert \phi {\Vert }_{C^m(V)}:-\underset{0\le i\le m}{\max }\underset{t\in [0,t_{\mathrm{F}}]}{\max }\Vert {\mathrm{d}}_t^i\phi (t){\Vert }_V.$$

$$(\phi ,\psi {)}_{H^m(V)}=\sum _{j=0}^m{\int }_0^{t_{\mathrm{F}}}({\mathrm{d}}_t^j\phi (t),{\mathrm{d}}_t^j\psi (t){)}_V\mathrm{d}t\forall \phi ,\psi \in H^m(V).$$

$$|\text{matr}\sigma (\vec{x},\text{matr}\tau ){|}_{d\times d}$$

$$\text{matr}\sigma (\vec{x},\text{matr}\tau ):\text{matr}\tau$$

$$\left(\text{matr}\sigma (\vec{x},\text{matr}\tau )-\text{matr}\sigma (\vec{x},\text{matr}\eta )\right):\left(\text{matr}\tau -\text{matr}\eta \right)$$

$$\text{matr}\sigma (\cdot ,\text{matr}\tau )=\lambda (\cdot )\mathrm{tr}(\text{matr}\tau ){\vec{I}}_d+2\mu (\cdot )\text{matr}\tau ,$$

$$\text{matr}\sigma (\cdot ,\text{matr}\tau )=\tilde{\lambda }(\cdot ,\mathrm{dev}(\text{matr}\tau ))\mathrm{tr}(\text{matr}\tau ){\vec{I}}_d+2\tilde{\mu }(\cdot ,\mathrm{dev}(\text{matr}\tau ))\text{matr}\tau ,$$

$$\text{matr}\sigma (\cdot ,\text{matr}\tau )=(1-D(\cdot ,\text{matr}\tau ))\mathfrak{C}(\cdot )\text{matr}\tau .$$

$$\underline{C}|\text{matr}\tau {|}_{d\times d}^2\le \mathfrak{C}(\vec{x})\text{matr}\tau :\text{matr}\tau \le \overline{C}|\text{matr}\tau {|}_{d\times d}^2\forall \text{matr}\tau \in {\mathbb{R}}^{d\times d},\forall \vec{x}\in \Omega .$$

$$\vec{U}:-{\text{H}}_0^1(\Omega )\text{and}P:-\{\begin{array}{ll}{H}^1(\Omega )& \text{if C0>0,}\\ H^1(\Omega )\cap L_0^2(\Omega )& \text{if C0=0,}\end{array}$$

$${\int }_0^{t_{\mathrm{F}}}a(\vec{u}(t),\vec{v})\phi (t)\mathrm{d}t+{\int }_0^{t_{\mathrm{F}}}b(\vec{v},p(t))\phi (t)\mathrm{d}t$$

$${\int }_0^{t_{\mathrm{F}}}\left[b(\vec{u}(t),q)-C_0(p(t),q{)}_{\Omega }\right]{\mathrm{d}}_t\phi (t)\mathrm{d}t+{\int }_0^{t_{\mathrm{F}}}c(p(t),q)\phi (t)\mathrm{d}t$$

$$(C_{0}p(0)+\vec{\nabla }\cdot \vec{u}(0),q{)}_{\Omega }$$

$$a(\vec{v},\vec{w}):-(\text{matr}\sigma (\cdot ,{\vec{\nabla }}_{\mathrm{s}}\vec{v}),{\vec{\nabla }}_{\mathrm{s}}\vec{w}{)}_{\Omega },b(\vec{v},q):--(\vec{\nabla }\cdot \vec{v},q{)}_{\Omega },c(q,r):-(\text{matr}\kappa (\cdot )\vec{\nabla }r,\vec{\nabla }q{)}_{\Omega }.$$

$${\mathrm{d}}_t\left[C_0(p,q{)}_{\Omega }-b(\vec{u},q)\right]+c(p,q)=(g,q{)}_{\Omega }\forall q\in P\text{ in }{L}^2((0,t_{\mathrm{F}})).$$

$${\mathrm{d}}_t\left(C_0{\int }_{\Omega }p(\cdot ,t)\right)={\int }_{\Omega }g(\cdot ,t)\text{ in }{L}^2((0,t_{\mathrm{F}})).$$

$$H^1({\mathcal{T}}_{\tau }):-\left\{\phi \in L^2((0,t_{\mathrm{F}}));{\phi }_{|(t^n,t^{n-1})}\in H^1((t^{n-1},t^n))\text{ for all }1\le n\le N\right\}.$$

$${\phi }^n:-\underset{t\to (t^n{)}^{-}}{\lim }\phi (t)\in V.$$

$${\overline{\psi }}^n:-{\tau }^{-1}{\int }_{t^{n-1}}^{t^n}\psi (t)\mathrm{d}t\in V,$$

$${\delta }_t{\phi }^n:-\frac{{\phi }^n-{\phi }^{n-1}}{\tau }\in V$$

$${\psi }^n-{\overline{\psi }}^n=\psi (t^n)-\frac{1}{\tau }{\int }_{t^{n-1}}^{t^n}\left(\psi (t^n)-{\int }_s^{t^n}{\mathrm{d}}_t\psi (t)\mathrm{d}t\right)\mathrm{d}s=\frac{1}{\tau }{\int }_{t^{n-1}}^{t^n}{\int }_s^{t^n}{\mathrm{d}}_t\psi (t)\mathrm{d}t\mathrm{d}s\le {\int }_{t^{n-1}}^{t^n}|{\mathrm{d}}_t\psi (t)|\mathrm{d}t.$$

$$\Vert {\phi }^n-{\overline{\phi }}^n{\Vert }_{\Omega }^2\le {\int }_{\Omega }{\left({\int }_{t^{n-1}}^{t^n}|{\mathrm{d}}_t\phi (\vec{x},t)|\mathrm{d}t\right)}^2\mathrm{d}\vec{x}\le \tau {\int }_{t^{n-1}}^{t^n}\Vert {\mathrm{d}}_t\phi (t){\Vert }_{\Omega }^2\mathrm{d}t\le \tau \Vert \phi {\Vert }_{H^1((t^{n-1},t^n);L^2(\Omega ))}^2,$$

$$\sum _{n=1}^N\tau \Vert {\phi }^n-{\overline{\phi }}^n{\Vert }_{\Omega }^2\le {\tau }^2\sum _{n=1}^N\Vert \phi {\Vert }_{H^1((t^{n-1},t^n);L^2(\Omega ))}^2-:{\tau }^2\Vert \phi {\Vert }_{H^1({\mathcal{T}}_{\tau };L^2(\Omega ))}^2.$$

$${\int }_X({\pi }_X^{k}v-v)w=0\forall w\in P^k(X).$$

$$|v-{\pi }_T^{k}v{|}_{H^m(T)}\le C_{\mathrm{ap}}{h}_T^{l-m}|v{|}_{H^l(T)}$$

$$h_T^{\frac{1}{2}}|v-{\pi }_T^{k}v{|}_{H^m({\mathcal{F}}_T)}\le C_{\mathrm{ap}}{h}_T^{l-m}|v{|}_{H^l(T)}$$

$$\underline{\vec{U}}^{k}_{h}\coloneq\left\{\underline{\vec{v}}_{h}=\big{(}(\vec{v}_{T})_{T\in\mathcal{T}_{h}},(\vec{v}_{F})_{F\in\mathcal{F}_{h}}\big{)}\;;\;\vec{v}_{T}\in\textbf{P}^{k}(T)\mbox{ for all }T\in\mathcal{T}_{h}\mbox{ and }\vec{v}_{F}\in\textbf{P}^{k}(F)\mbox{ for all }F\in\mathcal{F}_{h}\right\}.$$

$$({\vec{v}}_h{)}_{|T}={\vec{v}}_T\forall T\in {\mathcal{T}}_h.$$

Full text

Key words. Nonlinear poroelasticity, nonlinear Biot problem, Korn’s inequality, Hybrid High-Order methods, discontinuous Galerkin methods, polyhedral meshes

AMS subject classification. 65N08, 65N30, 76S05

1 Introduction

In this paper we analyze a Hybrid High-Order (HHO) discretization method for nonlinear poroelastic models. We overstep a previous work [6] devoted to the linear Biot model [4, 39] by incorporating more general, possibly nonlinear stress-strain constitutive laws [10]. The model is valid under the assumptions of small deformations of the rock matrix, small variations of the porosity, and small relative variations of the fluid density. The interest of the poroelastic models considered here is particularly manifest in geosciences applications [27, 28, 31], where fluid flows in geological subsurface, modeled as a porous media, induce a deformation of the rock matrix. The challenge is then to design a discretization method able to (i) treat a complex geometry with polyhedral meshes and nonconforming interfaces, (ii) handle possible heterogeneities of the poromechanical parameters and nonlinearities of the stress-strain relation, and (iii) deal with the numerical instabilities encountered in this type of coupled problem.

Let $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, denote a bounded connected polyhedral domain with boundary $\partial\Omega$ and outward normal $\vec{n}$. Without loss of generality, we assume that the domain is scaled so that its diameter is equal to 1. For a given finite time $t_{\rm F}>0$, volumetric load $\vec{f}$, fluid source $g$, we consider the nonlinear poroelasticity problem that consists in finding a vector-valued displacement field $\vec{u}$ and a scalar-valued pore pressure field $p$ solution of

[TABLE]

When discretizing the poroelasticity system (1), the main challenges are to ensure stability and convergence under mild assumptions on the nonlinear stress-strain relation and on the permeability field, and to prevent localized pressure oscillations arising in the case of poorly permeable, quasi-incompressible porous media. Since the latter issue is in part related to the saddle point structure in the coupled equations for $C_{0}=0$ and small $\underline{\kappa}$, the discrete spaces for the displacement and the pressure should satisfy an inf-sup condition. Indeed, as observed in [32, 26, 34] in the context of finite element discretizations of the linear poroelasticity problem, the inf-sup condition yields an $L^{2}$-estimate of the discrete pressure independent of $\underline{\kappa}^{-1}$, and allows one to prove the convergence of the approximate pressure towards the continuous pressure also in the incompressible case $C_{0}=0$. We notice, however, that the problem of spurious pressure oscillations is actually more involved than a simple saddle-point coupling issue. For instance, it has been recently pointed out in [35] that, even for discretization methods leading to an inf-sup stable discretization of the Stokes problem in the steady case, pressure oscillations can arise owing to a lack of monotonicity of the discrete operator. The robustness with respect to spurious oscillations has been numerically observed in [6, Section 6.2] for a HHO–dG discretization of the linear poroelasticity model.

In this work, we present and analyze a nonconforming space discretization of problem (1) where the nonlinear elasticity operator is discretized using the HHO method of [9] (c.f. also [19, 15]), while the Darcy operator relies on the Symmetric Weighted Interior Penalty (SWIP) method of [20]. The proposed method has several assets:

(i) it is valid in two and three space dimensions;

(ii) it delivers an inf-sup stable discretization on general spatial meshes including, e.g., polyhedral elements and nonmatching interfaces;

(iii) it allows one to increase the space approximation order to accelerate convergence in the presence of (locally) regular solutions.

Compared to the method proposed in [6] for the linear poroelasticity problem, there are two main differences in the design. First, for a given polynomial degree $k\geq 1$, the symmetric gradient reconstruction sits in the full space of tensor-valued polynomials of total degree $\leq k$, as opposed to symmetric gradients of vector-valued polynomials of total degree $\leq(k+1)$. Following [16, 9], this modification is required to obtain optimal convergence rates when considering nonlinear stress-strain laws. Second, the right-hand side of the discrete problem of Section 4.4 is obtained by taking the average in time of the loading force $\vec{f}$ and fluid source $g$ over a time step instead of their value at the end of the time step. This modification allows us to prove stability and optimal error estimates under significantly weaker time regularity assumptions on data (cf. Remarks 13 and 19). Finally, in Section 3.4 we give a new simple proof of a discrete counterpart of Korn’s inequality on HHO spaces, not requiring particular geometrical assumptions on the mesh. The interest of these results goes beyond the specific application considered here.

The material is organized as follows. In Section 2 we present the assumptions on the stress-strain law and the variational formulation of the nonlinear poroelasticity problem. In Section 3 we define the space and time meshes and the discrete spaces for the displacement and the pressure fields. In Section 4 we define the discrete counterparts of the elasticity, Darcy, and hydromechanical coupling operators and formulate the discrete problem. In Section 5 we prove the well-posedness of the scheme by deriving an a priori estimate on the discrete solution that holds also when the specific storage coefficient vanishes. The convergence analysis of the method is carried out in Section 6. Finally, Section 7 contains numerical tests to asses the performance of the method.

2 Continuous setting

In this section we introduce the notation for function spaces, formulate the assumptions on the stress-strain law, and derive a weak formulation of problem (1).

2.1 Notation for function spaces

Let $X\subset\overline{\Omega}$. Spaces of functions, vector fields, and tensor fields defined over $X$ are respectively denoted by italic capital, boldface Roman capital, and special Roman capital letters. The subscript “s” appended to a special Roman capital letter denotes a space of symmetric tensor fields. Thus, for example, $L^{2}(X),\textbf{L}^{2}(X)$, and $\mathbb{L}^{2}_{\rm s}(X)$ respectively denote the spaces of square integrable functions, vector fields, and symmetric tensor fields over $X$. For any measured set $X$ and any $m\in\mathbb{Z}$, we denote by $H^{m}(X)$ the usual Sobolev space of functions that have weak partial derivatives of order up to $m$ in $L^{2}(X)$, with the convention that $H^{0}(X)\coloneq L^{2}(X)$, while $C^{m}(X)$ and $C_{\rm c}^{\infty}(X)$ denote, respectively, the usual spaces of $m$-times continuously differentiable functions and infinitely continuously differentiable functions with compact support on $X$. We denote by $(\cdot,\cdot)_{X}$ and $(\cdot,\cdot)_{m,X}$ the usual scalar products in $L^{2}(X)$ and $H^{m}(X)$ respectively, and by $\|{\cdot}\|_{X}$ and $\|{\cdot}\|_{m,X}$ the induced norms.

For a vector space $V$ with scalar product $(\cdot,\cdot)_{V}$, the space $C^{m}(V)\coloneq C^{m}([0,t_{\rm F}];V)$ is spanned by $V$-valued functions that are $m$-times continuously differentiable in the time interval $[0,t_{\rm F}]$. The space $C^{m}(V)$ is a Banach space when equipped with the norm

[TABLE]

Similarly, the Hilbert space $H^{m}(V)\coloneq H^{m}((0,t_{\rm F});V)$ is spanned by $V$-valued functions of the time interval, and the norm $\|{\cdot}\|_{H^{m}(V)}$ is induced by the scalar product

[TABLE]

2.2 Stress-strain law

The following assumptions on the stress-strain relation are required to obtain a well-posed weak formulation of the nonlinear poroelasticity problem.

Assumption 1** (Stress-strain relation).**

We assume that the stress function $\matr{\sigma}:\Omega\times\mathbb{R}^{d\times d}_{\rm s}\to\mathbb{R}^{d\times d}_{\rm s}$ is a Carathéodory function, i.e., $\matr{\sigma}(\vec{x},\cdot)$ is continuous on $\mathbb{R}^{d\times d}_{\rm s}$ for almost every $\vec{x}\in\Omega$ and $\matr{\sigma}(\cdot,\matr{\tau})$ is measurable on $\Omega$ for all $\matr{\tau}\in\mathbb{R}^{d\times d}_{\rm s}$. Moreover, there exist real numbers $C_{\rm gr},C_{\rm cv}\in(0,+\infty)$ such that, for a.e. $\vec{x}\in\Omega$ and all $\matr{\tau},\matr{\eta}\in\mathbb{R}^{d\times d}_{\rm s}$, the following conditions hold:

[TABLE]

Above, we have introduced the Frobenius product such that, for all $\matr{\tau},\matr{\eta}\in\mathbb{R}^{d\times d}$, $\matr{\tau}:\matr{\eta}\coloneq\sum_{1\leq i,j\leq d}\tau_{ij}\eta_{ij}$ with corresponding matrix norm such that, for all $\matr{\tau}\in\mathbb{R}^{d\times d}$, $|\matr{\tau}|_{d\times d}\coloneq(\matr{\tau}:\matr{\tau})^{\frac{1}{2}}$.

Three meaningful examples for the stress-strain relation $\matr{\sigma}:\Omega\times\mathbb{R}^{d\times d}_{\rm s}\to\mathbb{R}^{d\times d}_{\rm s}$ in (1a) are:

• •

The (possibly heterogeneous) linear elasticity model given by the usual Hooke’s law

[TABLE]

where $\mu:\Omega\to[\mu_{*},\mu^{*}]$, with $0<\mu_{*}\leq\mu^{*}<+\infty$, and $\lambda:\Omega\to\mathbb{R}_{+}$ are the Lamé parameters.

• •

The nonlinear Hencky–Mises model of [33, 24] corresponding to the mechanical behavior law

[TABLE]

with nonlinear Lamé scalar functions $\tilde{\mu}:\Omega\times\mathbb{R}_{+}\to[\mu_{*},\mu^{*}]$ and $\tilde{\lambda}:\Omega\times\mathbb{R}_{+}\to\mathbb{R}_{+}$ depending on the deviatoric part $\operatorname{dev}(\matr{\tau})\coloneq\operatorname{tr}(\matr{\tau}^{2})-\frac{1}{d}\operatorname{tr}(\matr{\tau})^{2}$ of the strain.

• •

The isotropic reversible hyperelastic damage model [29], for which the stress-strain relation reads

[TABLE]

where $D:\Omega\times\mathbb{R}^{d\times d}_{\rm s}\to[0,1]$ is the scalar damage function and $\mathfrak{C}:\Omega\to\mathbb{R}^{d^{4}}$ is a fourth-order symmetric and uniformly elliptic tensor field, namely, for some strictly positive constants $\underline{C}$ and $\overline{C}$, it holds

[TABLE]

Being linear, the Cauchy stress tensor in (1c) clearly satisfies the previous assumptions. Moreover, under some mild requirements (cf. [11, 23]) on the nonlinear Lamé scalar functions $\tilde{\mu}$ and $\tilde{\lambda}$ in (1d) and on the damage function $D$ in (1e), it can be proven that also the Hencky–Mises model and the isotropic reversible damage model satisfy Assumption 1.

2.3 Weak formulation

At each time $t\in[0,t_{\rm F}]$, the natural functional spaces for the displacement $\vec{u}(t)$ and pore pressure $p(t)$ taking into account the boundary condition (1c) and the zero average constraint (1f) are, respectively,

[TABLE]

with $\textbf{H}_{0}^{1}(\Omega)\coloneq\left\{\vec{v}\in\textbf{H}^{1}(\Omega)\;;\;{\vec{v}}_{|\partial\Omega}=\vec{0}\right\}$ and $L^{2}_{0}(\Omega)\coloneq\left\{q\in L^{2}(\Omega)\;;\;\int_{\Omega}q=0\right\}$. We consider the following weak formulation of problem (1): For a loading term $\vec{f}\in L^{2}(\textbf{L}^{2}(\Omega))$, a fluid source $g\in L^{2}(L^{2}(\Omega))$, and an initial datum $\phi^{0}\in L^{2}(\Omega)$ that verify (1f) if $C_{0}>0$, find $\vec{u}\in L^{2}(\vec{U})$ and $p\in L^{2}(P)$ such that, for all $\vec{v}\in\vec{U}$, all $q\in P$, and all $\varphi\in C_{\rm c}^{\infty}((0,t_{\rm F}))$

[TABLE]

where we have defined the nonlinear function $a:\vec{U}\times\vec{U}\to\mathbb{R}$ and the bilinear forms $b:\vec{U}\times P\to\mathbb{R}$ and $c:P\times P\to\mathbb{R}$ such that, for all $\vec{v},\vec{w}\in\vec{U}$ and all $q,r\in P$,

[TABLE]

The first term in (1fa) is well defined thanks to the growth assumption (1ba). Moreover, owing to (1bb) together with Korn’s first inequality and Poincaré’s inequality, $a(\cdot,\cdot)$ and $c(\cdot,\cdot)$ are coercive on $\vec{U}$ and $P$, respectively. The strict monotonicity assumption (1bc) guarantees the uniqueness of the weak solution.

Remark 2* (Regularity of the fluid content and of the pore pressure).*

Using an integration by parts in time in (1fb), it is inferred that

[TABLE]

Therefore, defining the fluid content $\phi\coloneq C_{0}p+\vec{\nabla}\cdot\vec{u}$, we have that $t\mapsto(\phi(t),q)_{\Omega}\in H^{1}{((0,t_{\rm F}))}\subset C^{0}([0,t_{\rm F}])$ for all $q\in P$, and, as a result, (1fc) makes sense. Moreover, in the case $C_{0}>0$, taking $q=1$ in (1g) and owing to the definition of the bilinear form $c$ and the homogeneous Dirichlet condition (1c), we infer that

[TABLE]

Thus, $t\mapsto\int_{\Omega}p(\cdot,t)\in H^{1}{((0,t_{\rm F}))}$, namely the average of the pore pressure over $\Omega$ is a continuous function in $[0,t_{\rm F}]$.

3 Discrete setting

In this section we define the space and time meshes, recall the definition and properties of $L^{2}$-orthogonal projectors on local and broken polynomial spaces, and introduce the discrete spaces for the displacement and the pressure.

3.1 Space mesh

We consider here polygonal or polyhedral meshes corresponding to couples $\mathcal{M}_{h}\coloneq(\mathcal{T}_{h},\mathcal{F}_{h})$, where $\mathcal{T}_{h}$ is a finite collection of polygonal elements such that $h\coloneq\max_{T\in\mathcal{T}_{h}}h_{T}>0$ with $h_{T}$ denoting the diameter of $T$, while $\mathcal{F}_{h}$ is a finite collection of hyperplanar faces. It is assumed henceforth that the mesh $\mathcal{M}_{h}$ matches the geometrical requirements detailed in [22, Definition 7.2]; see also [21, Section 2]. To avoid dealing with jumps of the permeability coefficient inside elements, we additionally assume that $\mathcal{M}_{h}$ is compliant with the partition $P_{\Omega}$ on which $\matr{\kappa}$ is piecewise constant meaning that, for every $T\in\mathcal{T}_{h}$, there exists a unique subdomain $\omega\in P_{\Omega}$ such that $T\subset\omega$. For every mesh element $T\in\mathcal{T}_{h}$, we denote by $\mathcal{F}_{T}$ the subset of $\mathcal{F}_{h}$ containing the faces that lie on the boundary $\partial T$ of $T$. For each face $F\in\mathcal{F}_{T}$, $\vec{n}_{TF}$ is the (constant) unit normal vector to $F$ pointing out of $T$. Boundary faces lying on $\partial\Omega$ and internal faces contained in $\Omega$ are collected in the sets $\mathcal{F}_{h}^{{\rm b}}$ and $\mathcal{F}_{h}^{{\rm i}}$, respectively.

Our focus is on the so-called $h$-convergence analysis, so we consider a sequence of refined meshes that is regular in the sense of [21, Definition 3]. The mesh regularity assumption implies, in particular, that the diameter $h_{T}$ of a mesh element $T\in\mathcal{T}_{h}$ is uniformly comparable to the diameter $h_{F}$ of each face $F\in\mathcal{F}_{T}$, and that the number of faces in $\mathcal{F}_{T}$ is bounded above by an integer $N_{\partial}$ independent of $h$. We additionally assume that, for each mesh in the sequence, all the elements are star-shaped with respect to every point of a ball of radius uniformly comparable to the diameter of the element. This assumption is required to use the results of [8, Appendix A].

3.2 Time mesh

We subdivide $(0,t_{\rm F})$ into $N\in\mathbb{N}^{*}$ uniform subintervals, and introduce the timestep $\tau\coloneq t_{\rm F}/N$ and the discrete times $t^{n}\coloneq n\tau$ for all $0\leq n\leq N$. We define the space of piecewise $H^{1}$ functions on $(0,t_{\rm F})$ by

[TABLE]

Since each $\psi\in H^{1}((t^{n-1},t^{n}))$ has an absolutely continuous representative in $[t^{n-1},t^{n}]$, we can identify $\varphi\in H^{1}(\mathcal{T}_{\tau})$ with a left continuous function in $(0,t_{\rm F})$. Therefore, for any vector space $V$ and any $\varphi\in H^{1}(\mathcal{T}_{\tau};V)$, we set $\varphi^{0}\coloneq\varphi(0)$ and, for all $1\leq n\leq N$,

[TABLE]

If $\varphi\in C^{0}(V)$, this simply amounts to setting $\varphi^{n}\coloneq\varphi(t^{n})$. For all $n\geq 1$ and $\psi\in L^{1}(V)$, we define the time average of $\psi$ in $(t^{n-1},t^{n})$ as

[TABLE]

with the convention that $\overline{\psi}^{0}=0\in V$. We also let, for all $(\varphi^{i})_{0\leq i\leq N}\in V^{N+1}$ and all $1\leq n\leq N$,

[TABLE]

denote the backward approximation of the first derivative of $\varphi$ at time $t^{n}$.

We note a preliminary result that will be used in the convergence analysis of Section 6. Let $\psi\in H^{1}(\mathcal{T}_{\tau})$ and $1\leq n\leq N$. Identifying $\psi_{|(t^{n-1},t^{n})}\in H^{1}((t^{n-1},t^{n}))$ with its absolutely continuous representative in $(t^{n-1},t^{n}]$, one can use the fundamental theorem of calculus to infer that

[TABLE]

Thus, applying the previous result together with the Jensen inequality yields, for all $\varphi\in H^{1}(\mathcal{T}_{\tau};L^{2}(\Omega))$,

[TABLE]

where $\varphi(\vec{x},t)$ is a shorthand notation for $(\varphi(t))(\vec{x})$. As a result of (1i), we get

[TABLE]

3.3 $L^{2}$-orthogonal projectors on local and broken polynomial spaces

For $X\subset\overline{\Omega}$ and $k\in\mathbb{N}$, we denote by $P^{k}(X)$ the space spanned by the restriction to $X$ of scalar-valued, $d$-variate polynomials of total degree $k$. The $L^{2}$-projector $\pi_{X}^{k}:L^{1}(X)\to P^{k}(X)$ is defined such that, for all $v\in L^{1}(X)$,

[TABLE]

As a projector, $\pi_{X}^{k}$ is linear and idempotent so that, for all $v\in P^{k}(X)$, $\pi_{X}^{k}v=v$. When dealing with the vector-valued polynomial space $\textbf{P}^{k}(X)$ or with the tensor-valued polynomial space $\mathbb{P}^{k}(X)$, we use the boldface notation $\vec{\pi}_{X}^{k}$ for the corresponding $L^{2}$-orthogonal projectors acting component-wise. At the global level, we denote by $P^{k}(\mathcal{T}_{h})$, $\textbf{P}^{k}(\mathcal{T}_{h})$, and $\mathbb{P}^{k}(\mathcal{T}_{h})$, respectively, the spaces of scalar-valued, vector-valued, and tensor-valued broken polynomial functions on $\mathcal{T}_{h}$ of total degree $\leq k$, and by $\pi_{h}^{k}$ and $\vec{\pi}_{h}^{k}$ the $L^{2}$-projectors on $P^{k}(\mathcal{T}_{h})$ and $\textbf{P}^{k}(\mathcal{T}_{h})$, respectively. The following optimal approximation properties for the $L^{2}$-projector $\pi_{X}^{k}$ follow from [14, Lemmas 3.4 and 3.6]: There exists a strictly positive real number $C_{\rm ap}$ independent of $h$ such that, for all $T\in\mathcal{T}_{h}$, all $l\in\{0,\ldots,k+1\}$, all $m\in\{0,\ldots,l\}$, and all $v\in H^{l}(T)$,

[TABLE]

and, if $l\geq 1$ and $m\leq l-1$,

[TABLE]

where $|{\cdot}|_{H^{m}(\mathcal{F}_{T})}$ is the broken Sobolev seminorm on $\mathcal{F}_{T}$.

3.4 Discrete spaces

In this section we define the discrete spaces upon which the HHO method corresponding to a polynomial degree $k\geq 1$ is built.

3.4.1 Displacement

The discrete unknowns for the displacement are collected in the space

[TABLE]

For any $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h}$, we denote by $\vec{v}_{h}\in\textbf{P}^{k}(\mathcal{T}_{h})$ the broken polynomial vector field obtained patching element-based unknowns, so that

[TABLE]

The discrete unknowns corresponding to a function $\vec{v}\in\textbf{H}^{1}(\Omega)$ are obtained by means of the interpolator $\underline{\vec{I}}^{k}_{h}:\textbf{H}^{1}(\Omega)\to\underline{\vec{U}}^{k}_{h}$ such that

[TABLE]

For all $T\in\mathcal{T}_{h}$, we denote by $\underline{\vec{U}}^{k}_{T}$ and $\underline{\vec{I}}^{k}_{T}$ the restrictions to $T$ of $\underline{\vec{U}}^{k}_{h}$ and $\underline{\vec{I}}^{k}_{h}$, respectively, and, for any $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h}$, we let \underline{\vec{v}}_{T}\coloneq\big{(}\vec{v}_{T},(\vec{v}_{F})_{F\in\mathcal{F}_{T}}\big{)} collect the local discrete unknowns attached to $T$. At each time step, the displacement is sought in the following subspace of $\underline{\vec{U}}^{k}_{h}$ that strongly accounts for the homogeneous Dirichlet condition (1c):

[TABLE]

We next prove a discrete version of Korn’s first inequality on $\underline{\vec{U}}^{k}_{h,\mathrm{D}}$ that will play a key role in the analysis. To this purpose, we endow the space $\underline{\vec{U}}^{k}_{h}$ with the discrete strain seminorm $\|{\cdot}\|_{\matr{\varepsilon},h}$ defined, for all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h}$, such that

[TABLE]

We will also need the following continuous trace inequality, whose proof follows the arguments of [18, Lemma 1.49] (where a slightly different notion of mesh faces is considered): There exists a strictly positive real number $C_{\rm tr}$, independent of $h$ but possibly depending on the mesh regularity parameter, such that, for all $T\in\mathcal{T}_{h}$, all $\vec{v}_{T}\in\textbf{H}^{1}(T)$, and all $F\in\mathcal{F}_{T}$,

[TABLE]

Proposition 3** (Discrete Korn’s first inequality).**

There is a real number $C_{\rm K}>0$, only depending on $\Omega$, $d$, and the mesh regularity parameter, such that, for all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$,

[TABLE]

Remark 4* (Strain norm).*

An immediate consequence of (1p) is that the map $\|{\cdot}\|_{\matr{\varepsilon},h}$ is a norm in $\underline{\vec{U}}^{k}_{h,\mathrm{D}}$.

Proof.

We start by noticing that it holds, for all $\matr{\alpha}\in\mathbb{R}^{d\times d}_{\rm s}$ and all $\matr{\beta}\in\mathbb{R}^{d\times d}$, denoting by $\matr{\beta}_{\rm s}\coloneq\frac{1}{2}(\matr{\beta}+\matr{\beta}\trans)$ the symmetric part of $\matr{\beta}$,

[TABLE]

Let $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$. Since the divergence operator $\vec{\nabla}\cdot:\mathbb{H}^{1}_{\rm s}(\Omega)\to\textbf{L}^{2}(\Omega)$ is onto (c.f. [7, Section 9.1.1] and [1, Theorem 3.2]), there exists $\matr{\tau}_{\vec{v}_{h}}\in\mathbb{H}^{1}_{\rm s}(\Omega)$ such that $\vec{\nabla}\cdot\matr{\tau}_{\vec{v}_{h}}=\vec{v}_{h}$ and $\|\matr{\tau}_{\vec{v}_{h}}\|_{1,\Omega}\leq C_{\rm{sj}}\|\vec{v}_{h}\|_{\Omega}$, with $C_{\rm{sj}}>0$ independent of $h$. It follows that

[TABLE]

where, to pass to the second line, we have integrated by parts element by element and used the fact that $\matr{\tau}_{\vec{v}_{h}}$ has continuous normal traces across interfaces and that boundary unknowns are set to zero in order to insert $\vec{v}_{F}$ into the boundary term, while, to pass to the third line, we have used (1q) with $\matr{\alpha}=\matr{\tau}_{\vec{v}_{h}}$ and $\matr{\beta}=\vec{\nabla}\vec{v}_{T}$. Applying a Cauchy–Schwarz inequality on the integrals over the element, a generalized Hölder inequality with exponents $(2,2,+\infty)$ on the integrals over faces, using the fact that $\|\vec{n}_{TF}\|_{\textbf{L}^{\infty}(F)}\leq 1$, and invoking a discrete Cauchy–Schwarz inequality on the sum over $T\in\mathcal{T}_{h}$, we infer that

[TABLE]

where, to pass to the third line, we have estimated $\|\matr{\tau}_{\vec{v}_{h}}\|_{F}$ using the continuous trace inequality (1o) and used the fact that $h_{F}\leq h_{T}\leq{\rm diam}(\Omega)=1$ for any $T\in\mathcal{T}_{h}$ and $F\in\mathcal{F}_{T}$. Thus, invoking the boundedness of the divergence operator, we get

[TABLE]

which yields the conclusion with $C_{\rm K}=\sqrt{2N_{\partial}}C_{\rm sj}C_{\rm tr}$. ∎

3.4.2 Pore pressure

At each time step, the discrete pore pressure is sought in the space

[TABLE]

For any internal face $F\in\mathcal{F}_{h}^{{\rm i}}$, we denote by $T_{F,1},T_{F,2}\in\mathcal{T}_{h}$ the two mesh elements that share $F$, that is to say $F\subset\partial T_{F,1}\cap\partial T_{F,2}$ and $T_{F,1}\neq T_{F,2}$ (the ordering of the elements is arbitrary but fixed), and we set

[TABLE]

For all $q_{h}\in P_{h}^{k}$, we denote by $q_{T}$ the restriction of $q_{h}$ to an element $T\in\mathcal{T}_{h}$ and we define the discrete seminorm

[TABLE]

The fact that in (1s) boundary terms only appear on internal faces reflects the homogeneous Neumann boundary condition (1d).

Using the surjectivity of the divergence operator $\vec{\nabla}\cdot:\textbf{H}_{0}^{1}(\Omega)\to L^{2}_{0}(\Omega)$ and proceeding as in the proof of the discrete Korn inequality (1p), a discrete Poincaré–Wirtinger inequality in $P^{k}(\mathcal{T}_{h})$ is readily inferred, namely one has the existence of $C_{\rm P}>0$, only depending on $\Omega$, $d$, and the mesh regularity parameter such that, for all $q_{h}\in P^{k}(\mathcal{T}_{h})$,

[TABLE]

This result ensures, in particular, that the seminorm $\|{\cdot}\|_{\matr{\kappa},h}$ defined in (1s) is a norm on $P^{k}_{0}(\mathcal{T}_{h})$. For a proof of more general Sobolev inequalities on broken polynomial spaces, we refer the reader to [17] and [18, Section 5.1.2].

4 Discretization

In this section we define the discrete counterparts of the elasticity, hydro-mechanical coupling, and Darcy operators, and formulate the HHO–dG scheme for problem (1f).

4.1 Nonlinear elasticity operator

The discretization of the nonlinear elasticity operator closely follows [9]. We define the local symmetric gradient reconstruction $\vec{G}^{k}_{{\rm s},T}:\underline{\vec{U}}^{k}_{T}\to\mathbb{P}_{\rm s}^{k}(T)$ such that, for a given \underline{\vec{v}}_{T}=\big{(}\vec{v}_{T},(\vec{v}_{F})_{F\in\mathcal{F}_{T}}\big{)}\in\underline{\vec{U}}^{k}_{T}, $\vec{G}^{k}_{{\rm s},T}\underline{\vec{v}}_{T}\in\mathbb{P}_{\rm s}^{k}(T)$ solves

[TABLE]

Existence and uniqueness of $\vec{G}^{k}_{{\rm s},T}\underline{\vec{v}}_{T}$ follow from the Riesz representation theorem in $\mathbb{P}_{\rm s}^{k}(T)$ for the $L^{2}(T)^{d\times d}$-inner product. This definition is motivated by the following property.

Proposition 5** (Commuting property for the local symmetric gradient reconstruction).**

For all $\vec{v}\in\textbf{H}^{1}(T)$, it holds that

[TABLE]

Remark 6* (Approximation properties of the local symmetric gradient reconstruction).*

The commuting property (1u) combined with (1k) shows that $\vec{G}^{k}_{{\rm s},T}\underline{\vec{I}}^{k}_{T}\vec{v}$ optimally approximates $\vec{\nabla}_{\rm s}\vec{v}$ in $\mathbb{P}_{\rm s}^{k}(T)$.

Proof.

For all $\matr{\tau}\in\mathbb{P}^{k}(T)$, we can write

[TABLE]

where we have used (1q) with $\matr{\alpha}=\vec{G}^{k}_{{\rm s},T}\underline{\vec{I}}^{k}_{T}\vec{v}$ and $\matr{\beta}=\matr{\tau}$ in the first line, the definition (1t) of the local symmetric gradient with $\underline{\vec{v}}_{T}=\underline{\vec{I}}^{k}_{T}\vec{v}$ in the second line, and definition (1j) after observing that $\vec{\nabla}\cdot\matr{\tau}_{\rm s}\in\textbf{P}^{k-1}(T)\subset\textbf{P}^{k}(T)$ and $\matr{\tau}_{\rm s}\vec{n}_{TF}\in\textbf{P}^{k}(F)$ for all $F\in\mathcal{F}_{T}$ to remove the $L^{2}$-orthogonal projectors in the third line. In the fourth line, we have used an integration by parts, then invoked (1q) first with $\matr{\alpha}=\matr{\tau}_{\rm s}$ and $\matr{\beta}=\vec{\nabla}\vec{v}$, then with $\matr{\alpha}=\vec{\nabla}_{\rm s}\vec{v}$ and $\matr{\beta}=\matr{\tau}$, and we have used the definition (1j) of $\vec{\pi}_{T}^{k}$ to conclude. ∎

From $\vec{G}^{k}_{{\rm s},T}$, we define the local displacement reconstruction operator $\vec{r}^{k+1}_{T}:\underline{\vec{U}}^{k}_{T}\to\textbf{P}^{k+1}(T)$ such that, for all $\underline{\vec{v}}_{T}\in\underline{\vec{U}}^{k}_{T}$,

[TABLE]

and, denoting by $\partial_{i}$ the partial derivative with respect to the $i$th space variable, if $d=2$,

[TABLE]

while, if $d=3$,

[TABLE]

Optimal approximation properties for $\vec{r}^{k+1}_{T}\underline{\vec{I}}^{k}_{T}$ have been recently proved in [8, Appendix A] generalizing the ones of [19, Lemma 2]. The optimal approximation properties of $\vec{r}^{k+1}_{T}\underline{\vec{I}}^{k}_{T}$ are required to infer (1bf) below.

The discretization of the nonlinear elasticity operator is realized by the function $a_{h}:\underline{\vec{U}}^{k}_{h}\times\underline{\vec{U}}^{k}_{h}\to\mathbb{R}$ such that, for all $\underline{\vec{w}}_{h},\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h}$,

[TABLE]

where $\gamma>0$ denotes a user-dependent parameter and we penalize in a least-square sense the face-based residual $\vec{\varDelta}_{TF}^{k}:\underline{\vec{U}}^{k}_{T}\to\textbf{P}^{k}(F)$ such that, for all $T\in\mathcal{T}_{h}$, all $\underline{\vec{v}}_{T}\in\underline{\vec{U}}^{k}_{T}$, and all $F\in\mathcal{F}_{T}$,

[TABLE]

This definition ensures that $\vec{\varDelta}_{TF}^{k}$ vanishes whenever its argument is of the form $\underline{\vec{I}}^{k}_{T}\vec{w}$ with $\vec{w}\in\textbf{P}^{k+1}(T)$, a crucial property to obtain high-order error estimates (cf. [6, Theorem 12]). For further use, we note the following seminorm equivalence, which can be proved using the arguments of [19, Lemma 4]: For all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h}$,

[TABLE]

where $C_{\rm{eq}}>0$ is independent of $h$, and the discrete strain seminorm $\|{\cdot}\|_{\matr{\varepsilon},h}$ is defined by (1n). By (1bb), this implies the coercivity of $a_{h}$.

Remark 7* (Choice of the stabilization parameter).*

The constants $C_{\rm gr},C_{\rm cv}$ appearing in (1b) satisfy $C_{\rm cv}^{2}\leq C_{\rm gr}$. Indeed, owing to (1bb), the Cauchy–Schwarz inequality, and (1ba), it holds for all $\matr{\tau}\in\mathbb{R}^{d\times d}_{\rm s}$,

[TABLE]

Thus, we choose the stabilization parameter $\gamma$ in (1v) such that

[TABLE]

For the linear elasticity model (1c), we have $C_{\rm gr}=2\mu+d\lambda$ and $C_{\rm cv}=\sqrt{2\mu}$, so that a natural choice for the stabilization parameter is $\gamma=2\mu$.

4.2 Hydro-mechanical coupling

The hydro-mechanical coupling is realized by means of the bilinear form $b_{h}$ on $\underline{\vec{U}}^{k}_{h}\times P^{k}(\mathcal{T}_{h})$ such that, for all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h}$ and all $q_{h}\in P^{k}(\mathcal{T}_{h})$,

[TABLE]

where $q_{T}\coloneq q_{h|T}$ for all $T\in\mathcal{T}_{h}$. It can be checked using Cauchy–Schwarz inequalities together with the definition (1n) of the strain seminorm and discrete trace inequalities that there exists $C_{{\rm bd}}>0$ independent of $h$ such that

[TABLE]

Additionally, using the strongly enforced boundary condition in $\underline{\vec{U}}^{k}_{h,\mathrm{D}}$, it can be proved that

[TABLE]

Finally, we note the following lemma stating that the hybrid interpolator $\underline{\vec{I}}^{k}_{h}:\textbf{H}^{1}(\Omega)\to\underline{\vec{U}}^{k}_{h}$ is a Fortin operator.

Lemma 8** (Fortin operator).**

For all $\vec{v}\in\textbf{H}^{1}(\Omega)$ and all $q_{h}\in P^{k}(\mathcal{T}_{h})$, the interpolator $\underline{\vec{I}}^{k}_{h}$ satisfies

[TABLE]

where the strictly positive real number $C_{\rm st}$ is independent of $h$.

Proof.

(i) Proof of (1aba). Recalling the definitions (1n) of the discrete strain seminorm and (1m) of the global interpolator, we can write

[TABLE]

with $C_{\rm st}>0$ independent of $h$. To pass to the second line, we have used a triangle inequality after inserting $\pm\vec{\nabla}_{\rm s}\vec{v}$ into the first term, and we have used the linearity, idempotency, and boundedness of $\vec{\pi}_{F}^{k}$ to write $\|\vec{\pi}_{F}^{k}\vec{v}-\vec{\pi}_{T}^{k}\vec{v}\|_{F}=\|\vec{\pi}_{F}^{k}(\vec{v}-\vec{\pi}_{T}^{k}\vec{v})\|_{F}\leq\|\vec{v}-\vec{\pi}_{T}^{k}\vec{v}\|_{F}$. To conclude, we have used (1k) with $l=m=1$ and (1l) with $l=1$ and $m=0$ to bound the first and third term inside the summation.

(ii) Proof of (1abb). Recalling the definitions (1z) of $b_{h}(\cdot,\cdot)$ and (1m) of the global interpolator, we can write letting, for the sake of brevity, $q_{T}\coloneq q_{h|T}$ for all $T\in\mathcal{T}_{h}$, for all $q_{h}\in P^{k}(\mathcal{T}_{h})$,

[TABLE]

where we have used definition (1j) after observing that $\vec{\nabla}q_{T}\in\textbf{P}^{k-1}(T)\subset\textbf{P}^{k}(T)$ and $q_{T|F}\vec{n}_{TF}\in\textbf{P}^{k}(F)$ to remove the $L^{2}$-orthogonal projectors in the second line, and integration by parts over $T\in\mathcal{T}_{h}$ to conclude. ∎

As a result of the previous Lemma, one has the following inf-sup condition, cf. [10] for the proof, which follows the classical Fortin argument (see, e.g., [7, Section 8.4] for further details).

Proposition 9**.**

There is a strictly positive real number $\beta$ independent of $h$ such that, for all $q_{h}\in P^{k}_{0}(\mathcal{T}_{h})$,

[TABLE]

4.3 Darcy operator

The discretization of the Darcy operator is based on the Symmetric Weighted Interior Penalty method of [20], cf. also [18, Section 4.5]. For all $F\in\mathcal{F}_{h}^{{\rm i}}$ and all $q_{h}\in P^{k}(\mathcal{T}_{h})$, we define the jump and weighted average operators such that

[TABLE]

with $T_{F,1},T_{F,2}\in\mathcal{T}_{h}$, $T_{F,1}\neq T_{F,2}$, such that $F\subset\partial T_{F,1}\cap\partial T_{F,2}$ and $\kappa_{F,1},\kappa_{F,2}$ defined in (1r). The bilinear form $c_{h}$ on $P^{k}(\mathcal{T}_{h})\times P^{k}(\mathcal{T}_{h})$ is defined such that, for all $q_{h},r_{h}\in P^{k}(\mathcal{T}_{h})$,

[TABLE]

where we have introduced the broken gradient operator $\vec{\nabla}_{h}$ on $\mathcal{T}_{h}$ and we have denoted by $\varsigma>\underline{\varsigma}>0$ a user-defined penalty parameter chosen large enough to ensure the coercivity of $c_{h}$ (the proof is similar to [18, Lemma 4.51]):

[TABLE]

Since, under this condition, $c_{h}$ is a symmetric positive definite bilinear form on the broken polynomial space $P_{h}^{k}$, we can define an associated norm by setting $\|\cdot\|_{c,h}\coloneq c_{h}(\cdot,\cdot)^{\frac{1}{2}}$.

The following consistency result can be proved adapting the arguments of [18, Chapter 4] to homogeneous Neumann boundary conditions and will be instrumental for the analysis. We define the functional spaces P_{*}\coloneq\left\{r\in H^{1}(\Omega)\cap H^{2}({P_{\Omega}})\;;\;\text{\matr{\kappa}\vec{\nabla}r\cdot\vec{n}=0$on$\partial\Omega}\right\} and set $P_{*h}^{k}\coloneq P_{*}+P_{h}^{k}$. Extending the bilinear form $c_{h}$ to $P_{*h}^{k}\times P_{*h}^{k}$, it is inferred that, for all $r\in P_{*}$,

[TABLE]

4.4 Discrete problem

For all $1\leq n\leq N$, the discrete solution $(\underline{\vec{u}}_{h}^{n},p_{h}^{n})\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}\times P_{h}^{k}$ at time $t^{n}$ is such that, for all $(\underline{\vec{v}}_{h},q_{h})\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}\times P^{k}(\mathcal{T}_{h})$,

[TABLE]

Remark 10* (Initial condition).*

We observe that the initial displacement $\underline{\vec{u}}_{h}^{0}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$ and pressure $p_{h}^{0}\in P_{h}^{k}$ in (1aec) are not explicitly required to initialize the scheme. However, assuming that $\vec{f}\in C^{0}(\textbf{L}^{2}(\Omega))$, so that (1a) makes sense also for $t=0$, it is possible to compute the initial discrete fields $(\underline{\vec{u}}_{h}^{0},p_{h}^{0})$ by solving

[TABLE]

In the limit case $C_{0}=0$ the previous equations corresponds to a well-posed HHO discretization of a steady nonlinear Stokes-like problem. If $C_{0}>0$ we can take $q_{h}$ in (1afb) such that, for all $T\in\mathcal{T}_{h}$, $(q_{h})_{|T}=C_{0}^{-1}\operatorname{tr}(\vec{G}^{k}_{{\rm s},T}\underline{\vec{v}}_{T})$ and sum the resulting equation to (1afa). Owing to definitions (1t) and (1z), we obtain

[TABLE]

where the nonlinear function $\tilde{a}_{h}$ is defined as $a_{h}$ in (1v) but replacing the stress-strain law $\matr{\sigma}$ with

[TABLE]

According to [9, Theorem 7], the nonlinear elasticity problem (1ag) admits a solution. Once the initial displacement $\underline{\vec{u}}_{h}^{0}$ is computed, we set $p_{h}^{0}$ such that, for all $T\in\mathcal{T}_{h}$, $(p_{h}^{0})_{|T}=C_{0}^{-1}(\pi_{T}^{k}\phi^{0}-\operatorname{tr}(\vec{G}^{k}_{{\rm s},T}\underline{\vec{u}}_{T}^{0}))$.

Remark 11* (Time discretization).*

The modified backward Euler scheme obtained by taking time averages instead of pointwise evaluation of the right-hand sides in (1ae) can be interpreted as a low-order discontinuous Galerkin time-stepping method, cf. [36, 38].

Notice that other time discretizations could be used, but we have decided to focus on the backward Euler scheme to keep the proofs as simple as possible. From the practical point of view, at each time step $n$, the discrete nonlinear system (1ae) can be solved by the Newton method using as initial guess the solution at step $(n-1)$. The size of the linear system to be solved at each Newton iteration can be reduced by statically condensing a large part of the unknowns as described in [6, Section 5].

5 Stability and well-posedness

In this section we study the stability of problem (1ae) and prove its well-posedness. We start with an a priori estimate on the discrete solution not requiring conditions on the time step $\tau$ and robust with respect to vanishing storage coefficients and small permeability.

Proposition 12** (A priori estimate).**

Denote by $(\underline{\vec{u}}_{h}^{n},p_{h}^{n})_{1\leq n\leq N}$ the solution to (1ae). Under Assumption 1 on the stress-strain relation and the regularity on the data $\vec{f}$, $g$, and $\phi^{0}$ assumed in Section 2.3, it holds

[TABLE]

where $C>0$ denotes a real number independent of $h$, $\tau$, the physical parameters $C_{0}$ and $\matr{\kappa}$, and the final time $t_{\rm F}$. In (1aj), we have defined $s_{h}^{N}\coloneq\sum_{n=1}^{N}\tau p_{h}^{n}$ and we have adopted the convention that $C_{0}^{-1}\|\pi_{\Omega}^{0}{g}\|_{L^{2}(L^{2}(\Omega))}^{2}=0$ and $C_{0}^{-1}\|\pi_{\Omega}^{0}{\phi^{0}}\|_{\Omega}^{2}=0$ if $C_{0}=0$.

Remark 13*.*

In order to prove the a priori bound (1aj), no additional time regularity assumption on the loading term $\vec{f}$ and the mass source $g$ are needed, whereas the stability estimate of [6, Lemma 7], valid for linear stress-strain relation, requires $\vec{f}\in C^{1}(\textbf{L}^{2}(\Omega))$ and $g\in C^{0}(L^{2}(\Omega))$. On the other hand, Proposition 12 gives an estimate of the discrete displacement and pressure in the $L^{2}$-norm in time, while [6, Lemma 7] ensures a control in the $L^{\infty}$-norm in time. However, under additional requirements on the stress-strain law (for instance Assumption 16) and $H^{1}$-regularity in time of $\vec{f}$, a stronger version of (1aj) can be inferred, including in particular an estimate in the $L^{\infty}$-norm in time.

Proof.

(i) Estimate of $\|p_{h}^{n}-\pi_{\Omega}^{0}p_{h}^{n}\|_{\Omega}$. The growth property of the stress-strain function (1ba) together with the Cauchy–Schwarz inequality, assumption (1y) on the stabilization parameter, and the second inequality in (1w) yield, for all $1\leq n\leq N$ and all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$,

[TABLE]

Using the inf-sup condition (1ac), (1aa), and the mechanical equilibrium equation (1aea), we get, for any $1\leq n\leq N$,

[TABLE]

Therefore, owing to the discrete Korn inequality (1p) and to (1ak), we infer from the previous bound that

[TABLE]

(ii) Energy balance. For all $1\leq n\leq N$, summing (1aeb) at times $1\leq i\leq n$, taking $q_{h}=\tau^{2}p_{h}^{n}$ as a test function, and recalling the discrete initial condition (1aec) yields

[TABLE]

Moreover, using the linearity of $c_{h}$ and the formula $2x(x-y)=x^{2}+(x-y)^{2}-y^{2}$, the third term in the left-hand side of (1am) can be rewritten as

[TABLE]

where we have set $s_{h}^{0}\coloneq 0$, $s_{h}^{n}\coloneq\sum_{i=1}^{n}\tau p_{h}^{i}$ for any $1\leq n\leq N$, and observed that $p_{h}^{n}=\delta_{t}s_{h}^{n}$. Therefore, summing (1am) and (1aea) at discrete time $n$ with $\underline{\vec{v}}_{h}=\tau\underline{\vec{u}}_{h}^{n}$, leads to

[TABLE]

Summing the previous relation for $1\leq n\leq N$, telescoping out the appropriate summands, and using the coercivity property (1bb), assumption (1y), and the first inequality in (1w), we get

[TABLE]

with the notation $G^{n}\coloneq\sum_{i=1}^{n}\tau\overline{g}^{i}=\int_{0}^{t^{n}}g(t){\rm d}t$. We denote by $\mathcal{R}$ the right-hand side of (1an) and proceed to find a suitable upper bound.

(iii) Upper bound for $\mathcal{R}$. For the first term in the right-hand side of (1an), using the Cauchy–Schwarz, discrete Korn (1p), and Young inequalities, we obtain

[TABLE]

where we have used the Jensen inequality to infer that

[TABLE]

We estimate the second term in the right-hand side of (1an) by splitting it into two contributions as follows:

[TABLE]

where we have used its definition (1j) to move $\pi_{\Omega}^{0}$ from $p_{h}^{n}$ to $(G^{n}+\phi^{0})$ in the second term. Owing to the Cauchy–Schwarz, triangle, Jensen, and Young inequalities, and using (1al), we have

[TABLE]

Owing to the compatibility condition (1f) and the linearity of the $L^{2}$-projector, $\mathfrak{T}_{2}=0$ if $C_{0}=0$. Otherwise, using again the Cauchy–Schwarz, triangle, Jensen, and Young inequalities, leads to

[TABLE]

Finally, from (1ao), (1ap), (1aq), it follows that

[TABLE]

(iv) Conclusion. Passing the first two terms in the right-hand side of (1at) to the left-hand side of (1an) and multiplying both sides by a factor $4$, we obtain

[TABLE]

where we have denoted by $\mathcal{C}$ the last four summands in the right-hand side of (1at). In order to conclude we apply again (1al) to obtain a bound of the $L^{2}$-norm of the discrete pressure independent of the storage coefficient $C_{0}$. Indeed, owing to (1al) and $C_{\rm cv}^{2}\leq C_{\rm gr}$ (see Remark 7), it is inferred that

[TABLE]

Summing the previous relation to (1au) yields

[TABLE]

Thus, multiplying both sides of the previous relation by $\max\{C_{\rm eq}^{2}C_{\rm cv}^{-2},\,2C_{\rm gr}^{2}C_{\rm eq}^{6}C_{\rm cv}^{-2},\,1\}$ gives (1aj). ∎

Remark 14* (A priori bound for $C_{0}=0$).*

When $C_{0}=0$, the a priori bound (1aj) reads

[TABLE]

The conventions $C_{0}^{-1}\|\pi_{\Omega}^{0}g\|_{L^{2}(L^{2}(\Omega))}^{2}=0$ and $C_{0}^{-1}\|\pi_{\Omega}^{0}{\phi^{0}}\|_{\Omega}^{2}=0$ if $C_{0}=0$ are justified since the term $\mathfrak{T}_{2}$ in point (3) of the previous proof vanishes in this case thanks to the compatibility condition (1f).

We next proceed to discuss the existence and uniqueness of the discrete solutions. The proof of the following theorem hinges on the arguments of [13, Theorem 3.3].

Theorem 15** (Existence and uniqueness).**

Let Assumption 1 hold and let $(\mathcal{M}_{h})_{h\in{\cal H}}$ be a regular mesh sequence. Then, for all $h\in{\cal H}$ and all $N\in\mathbb{N}^{*}$, there exists a unique solution $(\underline{\vec{u}}_{h}^{n},p_{h}^{n})_{1\leq n\leq N}\in(\underline{\vec{U}}^{k}_{h,\mathrm{D}}\times P_{h}^{k})^{N}$ to (1ae).

Proof.

We define the linear stress-strain function ${\matr{\sigma}}^{\rm lin}:\mathbb{R}^{d\times d}_{\rm s}\to\mathbb{R}^{d\times d}_{\rm s}$ such that, for all $\matr{\tau}\in\mathbb{R}^{d\times d}_{\rm s}$,

[TABLE]

where $C_{\rm cv}$ is the coercivity constant of $\matr{\sigma}$ (see (1bb)), and we denote by $a_{h}^{\rm lin}$ the bilinear form obtained by replacing $\matr{\sigma}$ with ${\matr{\sigma}}^{\rm lin}$ in (1v). We consider the following auxiliary linear problem: For all $1\leq n\leq N$, find $(\underline{\vec{y}}_{h}^{n},p_{h}^{n})\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}\times P_{h}^{k}$ such that

[TABLE]

with initial condition as in (1aec). Since the previous system is linear and square and its solution satisfies the a priori estimate of Proposition 12, it is readily inferred that problem (1ay) admits a unique solution.

Now we observe that, thanks to the norm equivalence (1w), $a_{h}^{\rm lin}(\cdot,\cdot)$ is a scalar product on $\underline{\vec{U}}^{k}_{h,\mathrm{D}}$, and we define the mapping $\underline{\vec{\Phi}}_{h}:\underline{\vec{U}}^{k}_{h,\mathrm{D}}\to\underline{\vec{U}}^{k}_{h,\mathrm{D}}$ such that, for all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$,

[TABLE]

We want to show that $\underline{\vec{\Phi}}_{h}$ is an isomorphism. Let $\underline{\vec{v}}_{h},\underline{\vec{z}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$ be such that $\underline{\vec{\Phi}}_{h}(\underline{\vec{v}}_{h})=\underline{\vec{\Phi}}_{h}(\underline{\vec{z}}_{h})$. If $\underline{\vec{v}}_{h}\neq\underline{\vec{z}}_{h}$, owing to the norm equivalence (1w) and the fact that $\|\cdot\|_{\matr{\varepsilon},h}$ is a norm on $\underline{\vec{U}}^{k}_{h,\mathrm{D}}$, there is at least one $T\in\mathcal{T}_{h}$ such that $\vec{G}^{k}_{{\rm s},T}\underline{\vec{v}}_{T}\neq\vec{G}^{k}_{{\rm s},T}\underline{\vec{z}}_{T}$ or $\vec{\varDelta}_{TF}^{k}\underline{\vec{v}}_{T}\neq\vec{\varDelta}_{TF}^{k}\underline{\vec{z}}_{T}$ for some $F\in\mathcal{F}_{T}$. In both cases, owing to the definition of $a_{h}$ and the strict monotonicity assumption (1bc), it holds

[TABLE]

Thus, we infer by contradiction that $\underline{\vec{v}}_{h}=\underline{\vec{z}}_{h}$ and, as a result, $\underline{\vec{\Phi}}_{h}$ is injective. In order to prove that $\underline{\vec{\Phi}}_{h}$ is also onto, we recall the following result: If $(E,(\cdot,\cdot)_{E})$ is a Euclidean space and $\Psi:E\to E$ is a continuous map such that $\frac{(\Psi(x),x)_{E}}{\|x\|_{E}}\to+\infty$ as $\|x\|_{E}\to+\infty$, then $\Psi$ is surjective. Since $(\underline{\vec{U}}^{k}_{h,\mathrm{D}},\,a_{h}^{\rm lin}(\cdot,\cdot))$ is a Euclidean space and the coercivity (1bb) of $\matr{\sigma}$ together with the definition of $\matr{\sigma}^{\rm lin}$ yield $a_{h}^{\rm lin}(\underline{\vec{\Phi}}_{h}(\underline{\vec{v}}_{h}),\underline{\vec{v}}_{h})\geq a_{h}^{\rm lin}(\underline{\vec{v}}_{h},\underline{\vec{v}}_{h})$ for all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$, we deduce that $\underline{\vec{\Phi}}_{h}$ is an isomorphism. Let, for all $1\leq n\leq N$, $(\underline{\vec{y}}_{h}^{n},p_{h}^{n})\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}\times P_{h}^{k}$ be the solution to problem (1ay). By the surjectivity and injectivity of $\underline{\vec{\Phi}}_{h}$, for all $1\leq n\leq N$, there exists a unique $\underline{\vec{u}}_{h}^{n}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$ such that $\underline{\vec{\Phi}}_{h}(\underline{\vec{u}}_{h}^{n})=\underline{\vec{y}}_{h}^{n}$. By definition of $\underline{\vec{\Phi}}_{h}$ and $(\underline{\vec{y}}_{h}^{n})_{1\leq n\leq N}$, $(\underline{\vec{u}}_{h}^{n},p_{h}^{n})_{1\leq n\leq N}$ is therefore the unique solution of the discrete problem (1ae). ∎

6 Convergence analysis

In this section we study the convergence of problem (1ae) and prove optimal error estimates under the following additional assumptions on the stress-strain function $\matr{\sigma}$.

Assumption 16** (Stress-strain relation II).**

There exist real numbers $C_{\rm lp},C_{\rm mn}\in(0,+\infty)$ such that, for a.e. $\vec{x}\in\Omega$, and all $\matr{\tau},\matr{\eta}\in\mathbb{R}^{d\times d}_{\rm s}$,

[TABLE]

Remark 17* (Lipschitz continuity and strong monotonocity).*

It is readily seen, by taking $\matr{\eta}=\matr{0}$ in (1az), that Lipschitz continuity and strong monotonicity imply respectively the growth and coercivity properties of Assumption 1. Therefore, recalling (1x), it is inferred that the constants appearing in (1ba), (1bb), (1aza), and (1azb) satisfy

[TABLE]

It was proved in [2, Lemma 4.1] that the stress-strain relation for the Hencky–Mises model is strongly monotone and Lipschitz-continuous. Also the isotropic damage model satisfies Assumption 16 if the damage function in (1e) is, for instance, such that

[TABLE]

In order to prove a convergence rate of $(k+1)$ in space for both the displacement and pressure errors, we assume from this point on that the permeability tensor field $\matr{\kappa}$ is constant on $\Omega$, and that the following elliptic regularity holds (which is the case, e.g., when $\Omega$ is convex [25, 30]): There is a real number $C_{\rm el}>0$ only depending on $\Omega$ such that, for all $\psi\in L^{2}_{0}(\Omega)$, the unique function $\zeta\in P$ solution of the homogeneous Neumann problem

[TABLE]

is such that

[TABLE]

Let $(\underline{\vec{u}}_{h}^{n},p_{h}^{n})_{1\leq n\leq N}$ be the solution to (1ae). We consider, for all $1\leq n\leq N$, the discrete error components defined as

[TABLE]

where the global elliptic projection $\widehat{p}_{h}^{n}\in P_{h}^{k}$ is defined as the solution to

[TABLE]

Before proving the convergence of the scheme, we recall two preliminary approximation results for the projector $\underline{\vec{I}}^{k}_{h}$ and the projection $\widehat{p}_{h}^{n}$ that have been proved in [9, Theorem 16] and [6, Lemma 11], respectively. There is a strictly positive constant $C_{\rm pj}$ depending only on $\Omega$, $k$, and the mesh regularity parameter, such that,

• •

Assuming (1az) and $\vec{u}\in L^{2}(\vec{U}\cap\textbf{H}^{k+2}(\mathcal{T}_{h}))$ with $\matr{\sigma}(\cdot,\vec{\nabla}_{\rm s}\vec{u})\in L^{2}(\mathbb{H}^{k+1}_{\rm s}(\mathcal{T}_{h}))$, for a.e. $t\in(0,t_{\rm F})$ and all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$, it holds

[TABLE]

• •

Assuming the elliptic regularity (1bb) and $\overline{p}^{n}\in P\cap H^{k+1}(\mathcal{T}_{h})$, for all $1\leq n\leq N$, it holds

[TABLE]

Now we have all the ingredients to estimate the discrete errors defined in (1bc).

Theorem 18** (Error estimate).**

Let $(\vec{u},p)$ denote the unique solution to (1f), for which we assume

[TABLE]

with $\phi=C_{0}p+\vec{\nabla}\cdot\vec{u}$. If $C_{0}>0$, we further assume $\pi_{\Omega}^{0}p\in H^{1}(\mathcal{T}_{\tau};P^{0}(\Omega))\eqcolon H^{1}(\mathcal{T}_{\tau})$. Then, under Assumption 16 and the elliptic regularity (1bb), it holds

[TABLE]

where $C$ is a strictly positive constant independent of $h$, $\tau$, $C_{0}$, $\matr{\kappa}$, and $t_{\rm F}$, and, for the sake of brevity, we have defined $z_{h}^{N}\coloneq\sum_{n=1}^{N}\tau\epsilon_{h}^{n}$ and introduced the bounded quantities

[TABLE]

Remark 19* (Time regularity).*

In order to prove the previous error estimate, we only require the displacement $\vec{u}$ and the fluid content $\phi$ solving problem (1f) to be piecewise $H^{1}$-regular in $(0,t_{\rm F})$, whereas [6, Theorem 12] is established under the much stronger regularity $\vec{u}\in C^{2}(\textbf{H}^{1}(\Omega))$ and, if $C_{0}>0$, $p\in C^{2}(L^{2}(\Omega))$. Moreover, the assumptions $\phi\in H^{1}(\mathcal{T}_{\tau};L^{2}(\Omega))$ and, if $C_{0}>0$, $\pi_{\Omega}^{0}p\in H^{1}(\mathcal{T}_{\tau})$ are consistent with the time regularity results observed in Remark 2.

Proof.

(i) Estimate of $\|\underline{\vec{e}}_{h}^{n}\|_{\matr{\varepsilon},h}^{2}$. First we observe that, owing to (1a) and the definition of $b_{h}$ given in (1z), for all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$ and all $1\leq n\leq N$, we have

[TABLE]

where the residual linear form $\mathcal{R}^{n}:\underline{\vec{U}}^{k}_{h,\mathrm{D}}\to\mathbb{R}$ is defined such that, for all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$,

[TABLE]

Using the norm equivalence (1w), the strong monotonicity (1azb) of $\matr{\sigma}$ along with assumption (1y) on the stabilization parameter, the discrete mechanical equilibrium (1aea), and (1bi), yields

[TABLE]

Thus, owing to the previous relation and defining the dual norm

[TABLE]

we have that

[TABLE]

where the conclusion follows from Young’s inequality. Hence, rearranging, we arrive at

[TABLE]

(ii) Estimate of $C_{0}\|\epsilon_{h}^{n}\|_{\Omega}^{2}$. Using (1b), the fact that $(\vec{\nabla}\cdot\vec{u}(t),1)_{\Omega}=0$ to insert $\pi_{\Omega}^{0}q_{h}$, and the consistency property (1ad), we infer that, for all $q_{h}\in P_{h}^{k}$ and all $1\leq i\leq N$,

[TABLE]

Therefore, using the discrete mass conservation equation (1aeb), (1aa), the Fortin property (1abb) , the definition of the elliptic projection $\widehat{p}_{h}^{n}$, and (1bl), we obtain

[TABLE]

where, in order to pass to the last line, we have inserted $\pm\overline{p}^{i}$ into the first term inside brackets in the third line, we have defined, according to (1e), $\phi^{i}\coloneq C_{0}p^{i}+\vec{\nabla}\cdot\vec{u}^{i}$ for all $0\leq i\leq N$, and we have used the definition of the global $L^{2}$-projector $\pi_{\Omega}^{0}$. Moreover, setting $\widehat{p}_{h}^{0}\coloneq 0$, it follows from the initial condition (1aec), the boundary condition (1c), and (1aa) that

[TABLE]

For all $1\leq n\leq N$, summing (1bm) for $1\leq i\leq n$ with the choice $q_{h}=\tau\epsilon_{h}^{n}$, using (1bn), and proceeding as in the second step of the proof of Proposition 12, leads to

[TABLE]

where $z_{h}^{n}\coloneq\sum_{i=1}^{n}\tau\epsilon_{h}^{i}$ if $n\geq 1$ and $z_{h}^{0}\coloneq 0$. We bound the first term in the right-hand side of (1bq) applying the Cauchy–Schwarz and Young inequalities followed by the approximation result (1bg), yielding

[TABLE]

where, in order to pass to the second line, we have used the Cauchy–Schwarz inequality and adopted the notation $|\cdot|_{L^{2}((t^{n-1},t^{n});H^{m}(\mathcal{T}_{h}))}\coloneq\||\cdot|_{H^{m}(\mathcal{T}_{h})}\|_{L^{2}((t^{n-1},t^{n}))}$, for any $m\in\mathbb{N}$. We estimate the second and third terms using the Cauchy–Schwarz and the Young inequalities together with the time approximation result (1i) as follows:

[TABLE]

with $\eta$ denoting a positive real number that will be fixed later on in the proof. The relation obtained by plugging (1br) and (1bs) into (1bq) reads

[TABLE]

(iii) Estimate of $\|\epsilon_{h}^{n}-\pi_{\Omega}^{0}\epsilon_{h}^{n}\|_{\Omega}^{2}$. We proceed as in the first step of the proof of Proposition 12. Using the inf-sup condition (1ac), (1aa) followed by the definition (1bc) of the pressure error, the linearity of $b_{h}$, the mechanical equilibrium equation (1aea), and (1bi), we get, for all $1\leq n\leq N$,

[TABLE]

Moreover, the Lipschitz continuity of the stress-strain function (1aza), the Cauchy–Schwarz inequality, assumption (1y) on the stabilization parameter $\gamma$ together with (1ba), and the second inequality in (1w), lead to

[TABLE]

Therefore, plugging the previous bound into the last line of (1bw), yields

[TABLE]

Squaring and rearranging the previous relation and recalling that, owing to (1ba), $C_{\rm mn}^{2}\leq C_{\rm lp}$, it is inferred that

[TABLE]

(iv) Estimate of the dual norm of the residual. We split the residual linear form $\mathcal{R}^{n}$ defined in (1bj) into three contributions $\mathcal{R}^{n}\coloneq\mathcal{R}_{1}^{n}+\mathcal{R}_{2}^{n}+\mathcal{R}_{3}^{n}$, defined, for all $\underline{\vec{v}}_{h}\in\underline{\vec{U}}^{k}_{h,\mathrm{D}}$ and all $1\leq n\leq N$, such that

[TABLE]

The first contribution can be bounded proceeding as in (1bx), then using the stability property (1aba) of the interpolator $\underline{\vec{I}}^{k}_{h}$, the Cauchy–Schwarz inequality, and a Poincaré–Wirtinger inequality on the time interval $(t^{n-1},t^{n})$. By doing so, we get

[TABLE]

We estimate the residual linear form $\mathcal{R}_{2}^{n}$ defined in (1bzb) using the consistency property (1bf) and the Cauchy–Schwarz inequality, obtaining

[TABLE]

Finally, the third term in (1bz) can be bounded integrating by parts element-wise and using the Cauchy–Schwarz inequality, the trace inequality (1o), the consistency property (1bg), and again the Cauchy–Schwarz inequality, namely

[TABLE]

Therefore, combining (1ca), (1cb), and (1cc), it is inferred that

[TABLE]

with

[TABLE]

(v) Conclusion. Adding (1bk) to (1bv) with $\eta=4\beta^{2}C_{\rm mn}^{-2}C_{\rm lp}^{2}C_{\rm eq}^{6}$, using (1by) and (1cd), summing the resulting equation over $1\leq n\leq N$, and multiplying both sides by $2\tau$, we obtain

[TABLE]

with $\widetilde{C}\coloneq\max\left\{1,\,C_{\rm pj}^{2},\,2C_{\rm pj}^{2}C_{\rm mn}^{-2}(2C_{\rm eq}^{2}+C_{\rm eq}^{-6}),\,2C_{\rm lp}^{2}C_{\rm st}^{2}C_{\rm ap}^{2}(2C_{\rm eq}^{6}+C_{\rm eq}^{-2}),\,4\beta^{2}C_{\rm lp}^{2}C_{\rm eq}^{6}C_{\rm mn}^{-2}\right\}$ and

[TABLE]

Finally, multiplying both sides of (1ce) by $2(C_{\rm pj}^{-2}+1)\widetilde{C}$ yields (1bh) with $C=2(C_{\rm pj}^{-2}+1)\widetilde{C}^{2}$. ∎

7 Numerical results

We consider a regular exact solution in order to assess the convergence of the method. Specifically, we solve problem (1) in the square domain $\Omega=(0,1)^{2}$ with $t_{\rm F}=1$ and physical parameters $C_{0}=0$ and $\matr{\kappa}=\vec{I}_{d}$. As nonlinear constitutive law we take the Hencky–Mises relation given by

[TABLE]

It can be checked that the previous stress-strain relation satisfies Assumption 1. The exact displacement $\vec{u}$ and exact pressure $p$ are given by

[TABLE]

The volumetric load $\vec{f}$, the source term $g$, and the boundary conditions are inferred from the exact solution. We consider the Cartesian and Voronoi mesh families depicted in Figure 1 and polynomial degree $k=1$. The time step $\tau$ on the coarsest mesh is taken to be $0.2/2^{k+1}$ for every choice of the polynomial degree $k$, and it decreases with the mesh size $h$ according to the theoretical convergence rates, thus, $\tau_{l}/\tau_{l+1}=2^{k}h_{l}/h_{l+1}$. Table 1 displays convergence results for the two mesh families. The error measures are $\left(\sum_{n=1}^{N}\tau\|\underline{\vec{u}}_{h}^{n}-\underline{I}^{k}_{h}\overline{\vec{u}}^{n}\|_{\matr{\varepsilon},h}^{2}\right)^{\frac{1}{2}}$ for the displacement and $\left(\sum_{n=1}^{N}\tau\|p_{h}^{n}-\pi_{h}^{k}\overline{p}^{n}\|_{\Omega}^{2}\right)^{\frac{1}{2}}$ for the pressure. In all cases, the orders of convergence are in agreement with the theoretical predictions. In particular, it is observed that the optimal convergence rates stated in Theorem 18 are reached.

Acknowledgements

This work was partially funded by the Bureau de Recherches Géologiques et Minières. The work of M. Botti was additionally partially supported by Labex NUMEV (ANR-10-LABX-20) ref. 2014-2-006. The work of D. A. Di Pietro was additionally partially supported by project HHOMM (ANR-15-CE40-0005).