Mathematical analysis of a penalisation strategy for incompressible elastodynamics

Abstract

This work addresses the mathematical analysis – by means of asymptotic analysis – of a penalisation strategy for the full discretisation of elastic wave propagation problems in quasi-incompressible media that has been recently developed by the authors. We provide a convergence analysis of the solution to the continuous version of the penalised problem towards its formal limit when the penalisation parameter tends to infinity. Moreover, as a fundamental intermediate step we provide an asymptotic analysis of the convergence of solutions to quasi-incompressible problems towards solutions to purely incompressible problems when the incompressibility parameter tends to infinity. Finally, we further detail the regularity assumptions required to guarantee that the mentioned convergence holds.

Keywords

Elastodynamics incompressibility asymptotic analysis

1. Introduction

Numerous finite element method (FEM) formulations have been specifically developed in the last decades to accurately solve the elasticity equations in pure and nearly-incompressible solids, due to the numerous applications in computational mechanics, f.e. for biological tissues. In this regard, we are interested in analysing the propagation of elastic waves in heterogenous, anisotropic and nearly-incompressible biological tissues (e.g. the heart), in the context of a recent biomedical imaging technique called “transient elastography”, that is raising a growing interest in clinical applications.

In this context we have recently introduced in [7] a novel formulation for wave propagation in incompressible solids. We recall that the proposed method is based on the extension to incompressible elasticity of existing numerical schemes suitable for viscous incompressible flows, due to the strong similarities between the respective governing equations. Indeed, when the bulk modulus λ of an elastic material is very large, this enforces the divergence of the displacement to be close to zero. Furthermore, at the limit $λ \to \infty$ , the pressure can be introduced as a Lagrange multiplier associated with the incompressibility constraint [4] and thus the similarity with unsteady Stokes equations.

For this reason, we have developed in [7] a new numerical scheme inspired by fractional-step algorithms and penalisation techniques (see [9,10,18–20,25,26]) for the resolution of incompressible elasticity. In particular, we propose a conservative time discretisation that treats implicitly only the terms corresponding to “informations” traveling at infinite velocity (i.e. the incompressibility constraint). Therefore, if efficient methods for explicit time-discretisation are used (e.g. Spectral Finite Elements with mass lumping), our algorithm requires at each time step the resolution of a scalar Poisson problem (that can be performed by efficient algorithms [8]) and few matrix-vector multiplications for the explicit methods.

In more detail, the proposed penalisation strategy – that can be formulated in the continuous setting – introduces a consistency error depending on the penalisation parameter $α^{- 1}$ . The objective of this work is the analysis of this consistency error at the continuous level. We point out that we are interested in the approximation of the quasi-incompressible problem. However, the penalisation strategy is constructed as an approximation of the pure incompressible problem. Therefore, in this work we also analyse the discrepancy between the pure and quasi-incompressible formulations (as the reader will see, the two analyses share some similarities). Our analysis is restricted to the continuous framework, but it can be seen as a preliminary step to the convergence analysis in a fully discrete setting.

The work is organised as follows. In Section 2, after introducing some preliminary results and notations, we give simple – but not necessarily standard – results for the formulations of the continuous elastodynamic problem in quasi-incompressible and pure incompressible media. Moreover, we recall the novel formulation that takes into account incompressibility by penalisation.

Section 3 is devoted to the derivation of the convergence estimates (in $H^{1}$ -norm in space) between the penalised formulation and the standard incompressible problem at the limit $α^{- 1} \to \infty$ , and between the quasi-incompressible formulation and the standard incompressible problem at the limit $λ \to \infty$ . By triangular inequality we retrieve the convergence estimate between the penalised and quasi-incompressible formulations.

In order to retrieve higher-order estimates, we also tackle the convergence estimates in $L^{2}$ -norm in space. To do so, we introduce in Section 4 the elasto-static operator associated with the elastodynamic problem and use intrinsic regularity properties of this operator, following a similar approach to the one presented by [17] for the Stokes problem to retrieve the aforementioned convergence estimates.

Section 5 is devoted to a further investigation of the regularity properties required to derive the obtained convergence estimates. In the Appendix that follows we provide a detailed proof on these aspects.

2. Setting of the problem

2.1. Preliminaries and notations

We first introduce several notations and recall some important inequalities that are fundamental to the analysis of our problem.

Hilbert spaces and dual spaces. In this work we denote by $Ω \subset R^{d}$ , with $d = 2$ or $d = 3$ , an open, connected and bounded domain with Lipschitz boundary. We introduce the following notations to define the Hilbert spaces for the elastic displacements $\begin{matrix} H : = {\underline{v} \in L^{2} {(Ω)}^{d}}, X : = H_{0}^{1} {(Ω)}^{d}, X^{'} = H^{- 1} {(Ω)}^{d}, \end{matrix}$ where $H$ is equipped with the usual $L^{2}$ -scalar product and $X$ is equipped with the usual $H^{1}$ -scalar product. For the sake of simplicity, we have considered homogenous Dirichlet conditions on the boundary of the propagation domain. We also need to consider divergence-free displacements. Hence, following [15], we introduce the subspaces of $X$ and $X^{'}$ respectively $\begin{matrix} V : = {\underline{v} \in X ∣ div \underline{v} = 0}, V^{0} : = {\underline{ϕ} \in X^{'} ∣ {⟨ \underline{ϕ}, \underline{v} ⟩}_{X^{'}, X} = 0, \forall \underline{v} \in V} . \end{matrix}$ We also introduce the space $D$ satisfying $V \subset D$ and defined as $\begin{matrix} D = {\underline{v} \in H (div; Ω) ∣ div \underline{v} = 0 in Ω, \underline{v} \cdot \underline{n} = 0 on \partial Ω} . \end{matrix}$ Note that $D$ is a Hilbert space when equipped with the scalar product of $H (div; Ω)$ . Pressure is a variable of interest and is sought in the spaces $\begin{matrix} L : = L_{0}^{2} (Ω) : = {q \in L^{2} (Ω) | \int_{Ω} q = 0}, M : = {q \in H^{1} (Ω) ∣ q \in L}, \end{matrix}$ where $L$ (respectively $M$ ) is equipped with the usual $L^{2}$ (respectively $H^{1}$ ) scalar product. As usual, we identify $L$ and $H$ with their dual spaces in what follows. Finally, we introduce the following subspace of $M$ , $\begin{matrix} N : = {q \in M ∣ Δ q \in L^{2} (Ω), \underline{\nabla} q \cdot \underline{n} \in L^{2} (\partial Ω)}, \end{matrix}$ where $\underline{n}$ denotes the unitary outward normal of the domain Ω.

Gradients. We recall the Corollary 2.4 of [15]: the operator $\underline{\nabla}$ is an isomorphism of $L$ onto $V^{0}$ . As a consequence, there exists $C > 0$ such that $\begin{matrix} (2.1) & ‖ q ‖_{L} ⩽ C ‖ \underline{\nabla} q ‖_{X^{'}}, \forall q \in L . \end{matrix}$

A trace inequality. Following the proof of Theorem 1.5.1.10 provided in [16], one can show (using also Poincaré-Friedrich’s inequality) that there exists $C > 0$ such that $\begin{matrix} (2.2) & ‖ q ‖_{L^{2} (\partial Ω)} ⩽ C ‖ q ‖_{L}^{\frac{1}{2}} ‖ \underline{\nabla} q ‖_{H}^{\frac{1}{2}}, \forall q \in M . \end{matrix}$

The Korn inequality. We recall here the well-known Korn inequality (see for instance Theorem 6.3-4 of [11]). There exists $C > 0$ such that $\begin{matrix} (2.3) & ‖ \underline{y} ‖_{X} ⩽ C {‖ \underline{\underline{ε}} (\underline{y}) ‖}_{L^{2} (Ω)}, \end{matrix}$ where $\underline{\underline{ε}} (\underline{y})$ is the linearised strain tensor associated with (or the symmetrised gradient of) the field $\underline{y}$ .

Space–time function spaces. Henceforth, for the sake of conciseness, we use the following notation: given a function space $A$ , we define $\begin{matrix} C^{k} (A) : = C^{k} ([0, T]; A), k \in {0, 1, 2, \dots}, \end{matrix}$ where $T > 0$ is a given final time of observation and we use the notation $\begin{matrix} ‖ u ‖_{L^{\infty} (A)} = sup_{t \in [0, T]} {‖ u (t) ‖}_{A} \end{matrix}$ as a space–time norm of a function belonging in $C^{0} (A)$ . Furthermore, we introduce the following Sobolev spaces, for all $k \in N^{*}$ : $\begin{matrix} W_{0}^{k} (0, T) : = {\underline{v} \in W^{k, 1} (0, T) ∣ \partial_{t}^{p} \underline{v} (t = 0) = \underline{0} \forall p < k} . \end{matrix}$ Note that this definition makes sense since, for all $k \in N^{*}$ , there exists an imbedding from $W_{0}^{k, 1} (0, T)$ to $C^{k - 1} ([0, T])$ [1]. Given a function space $A$ , we define the set of functions in $W_{0}^{k} (0, T)$ with values in $A$ , namely, $W_{0}^{k} (A) : = W_{0}^{k} (0, T; A)$ .

Uniform estimates. In what follows our objective is to obtain uniform estimates, with respect to the large parameter λ and the small parameter α, of solutions to elastodynamic problems. Therefore, up to Section 4 included, we use the following notation $\begin{matrix} a_{λ, α} ≲ b_{λ, α} \end{matrix}$ to denote that there exists a positive scalar C independent of λ and α such that the scalars $a_{λ, α}$ and $b_{λ, α}$ (that depend on λ and α) satisfy $\begin{matrix} a_{λ, α} ⩽ C b_{λ, α}, \end{matrix}$ for all λ sufficiently large and all α sufficiently small.

2.2. The quasi-incompressible elastodynamic equations

For the sake of simplicity, we assume that all quantities in the elastodynamic problem are non-dimensional (see [7] for a presentation of the non-dimensionalisation process). As a model problem, we consider the elastic wave propagation in a quasi-incompressible solid, which is described by the following partial differential equation system:

For $\underline{f}$ given and sufficiently regular, find ${\underline{y}}_{λ}$ sufficiently smooth, such that $\begin{matrix} (QI) & \{\begin{matrix} ρ \partial_{t}^{2} {\underline{y}}_{λ} - \underline{div} (\underline{\underline{C}} \underline{\underline{ε}} ({\underline{y}}_{λ})) - λ \underline{\nabla} div {\underline{y}}_{λ} = \underline{f} & in Ω_{T} : = Ω \times (0, T), \\ {\underline{y}}_{λ} (t = 0) = \underline{0}, \partial_{t} {\underline{y}}_{λ} (t = 0) = \underline{0} & in Ω, \end{matrix} \end{matrix}$ with $λ > > 1$ the bulk modulus, that is large due to quasi-incompressibility, $ρ (\underline{x})$ the positive density of the medium and $\underline{\underline{C}} (\underline{x})$ an elasticity fourth-order tensor, that is well defined almost everywhere in Ω and is symmetric, coercive and bounded, i.e. there exist two strictly positive scalars c, C such that $\begin{matrix} (2.4) & c | \underline{\underline{ε}} |^{2} ⩽ \underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} : \underline{\underline{ε}} ⩽ C | \underline{\underline{ε}} |^{2}, \forall \underline{\underline{ε}} \in R_{sym}^{d \times d}, a.e. in Ω, \end{matrix}$ where $R_{sym}^{d \times d}$ is the set of symmetric $d \times d$ second-order tensors and the symbol : denotes the tensor contraction product. For the sake of simplicity, we take $ρ \equiv 1$ in what follows. Note that (QI) can be deduced from the standard classical elastodynamic equation by assuming that the classical elasticity tensor $\tilde{\underline{\underline{C}}}$ has the following form $\begin{matrix} \tilde{\underline{\underline{C}}} (\underline{x}) \underline{\underline{ε}} = λ \underline{\underline{I}} \underline{\underline{ε}} : \underline{\underline{I}} + \underline{\underline{C}} (\underline{x}) \underline{\underline{ε}}, \forall \underline{\underline{ε}} \in R_{sym}^{d \times d}, a.e. in Ω, \end{matrix}$ where $\underline{\underline{I}}$ is the identity second-order tensor. Note that the tensor $\underline{\underline{C}}$ accounts for a general anisotropic medium however if homogeneous isotropic elastodynamics is considered, then one has $\underline{\underline{C}} (\underline{x}) = 2 μ \underline{\underline{I}}$ . Existence and uniqueness results for problem (QI) are well-known. As an illustration, it is proved in [21] (Theorem 9.1, see also Definition 9.1 for the definition of weak solution) that the following proposition holds:

Proposition 2.1.
Let $\underline{f} \in L^{1} (H)$ . Then, a weak solution ${\underline{y}}_{λ}$ to ( QI ) exists, it is unique and satisfies $\begin{matrix} (2.5) & {\underline{y}}_{λ} \in C^{1} (H) \cap C^{0} (X) . \end{matrix}$

As a consequence of Proposition 2.1, one can show (see [21], Theorem 9.2) the following energy estimates: $\begin{matrix} E (t) = \int_{0}^{t} \int_{Ω} \underline{f} \partial_{t} {\underline{y}}_{λ} d Ω, \end{matrix}$ with $\begin{matrix} E (t) = \frac{1}{2} ‖ \partial_{t} {\underline{y}}_{λ} ‖_{H}^{2} + \frac{1}{2} \int_{Ω} \underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{y}}_{λ}) : \underline{\underline{ε}} ({\underline{y}}_{λ}) d Ω + \frac{λ}{2} ‖ div {\underline{y}}_{λ} ‖_{H}^{2} . \end{matrix}$ It is standard to show, using the Grönwall lemma, $\begin{matrix} (2.6) & sup_{t \in [0, T]} E^{\frac{1}{2}} (t) ≲ \int_{0}^{T} ‖ \underline{f} ‖_{H} d s . \end{matrix}$ Therefore, we can assert by the Korn inequality that $\begin{matrix} (2.7) & ‖ \partial_{t} {\underline{y}}_{λ} ‖_{L^{\infty} (H)} + ‖ {\underline{y}}_{λ} ‖_{L^{\infty} (X)} + λ^{\frac{1}{2}} ‖ div {\underline{y}}_{λ} ‖_{L^{\infty} (L)} ≲ 1 . \end{matrix}$ Under further regularity assumptions on $\underline{f}$ in time, we can formally differentiate (QI) with respect to time and apply Proposition 2.1 again. As a consequence, we can state the following lemma.
Lemma 2.2.
Let $\underline{f} \in W_{0}^{k} (H)$ with $k ⩾ 1$ . Then, the solution ${\underline{y}}_{λ}$ to ( QI ) satisfies $\begin{matrix} {\underline{y}}_{λ} \in C^{k + 1} (H) \cap C^{k} (X) . \end{matrix}$ Moreover, $\begin{matrix} {‖ \partial_{t}^{k + 1} {\underline{y}}_{λ} ‖}_{L^{\infty} (H)} + {‖ \partial_{t}^{k} {\underline{y}}_{λ} ‖}_{L^{\infty} (X)} + λ^{\frac{1}{2}} {‖ div \partial_{t}^{k} {\underline{y}}_{λ} ‖}_{L^{\infty} (L)} ≲ 1 . \end{matrix}$

Assume now $\underline{f} \in W_{0}^{1} (H)$ . An alternative formulation to (QI) consists in introducing a scalar function $p_{λ}$ such that the couple $({\underline{y}}_{λ}, p_{λ}) \in [C^{2} (H) \cap C^{1} (X)] \times C^{1} (L)$ satisfies $\begin{matrix} (QIM) & \{\begin{matrix} \partial_{t}^{2} {\underline{y}}_{λ} - \underline{div} (\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{y}}_{λ})) - \underline{\nabla} p_{λ} = \underline{f} & in Ω_{T}, \\ div {\underline{y}}_{λ} = λ^{- 1} p_{λ} & in Ω_{T}, \\ {\underline{y}}_{λ} (t = 0) = \underline{0}, \partial_{t} {\underline{y}}_{λ} (t = 0) = \underline{0} & in Ω . \end{matrix} \end{matrix}$ Note that the mixed formulation (QIM) is obtained straightforwardly from (QI) and therefore existence and uniqueness results for (QIM) can be easily deduced. Theorem 2.3.
Let $\underline{f} \in W_{0}^{k} (H)$ with $k ⩾ 1$ . Then, the unique solution $({\underline{y}}_{λ}, p_{λ})$ to ( QIM ) satisfies $\begin{matrix} ({\underline{y}}_{λ}, p_{λ}) \in C^{k + 1} (H) \cap C^{k} (X) \times C^{k} (L) . \end{matrix}$ Moreover, $\begin{matrix} {‖ \partial_{t}^{k + 1} {\underline{y}}_{λ} ‖}_{L^{\infty} (H)} + {‖ \partial_{t}^{k} {\underline{y}}_{λ} ‖}_{L^{\infty} (X)} + λ^{- \frac{1}{2}} {‖ \partial_{t}^{k} p_{λ} ‖}_{L^{\infty} (L)} + {‖ \partial_{t}^{k - 1} p_{λ} ‖}_{L^{\infty} (L)} ≲ 1 . \end{matrix}$
Proof.
The proof can be divided into the following three steps:

Step 1. The uniform estimate (w.r.t. λ) on $\partial_{t}^{k + 1} {\underline{y}}_{λ} (t)$ in the $H$ -norm, $\partial_{t}^{k} {\underline{y}}_{λ} (t)$ in the $X$ -norm and $\partial_{t}^{k} p_{λ} (t)$ in the $L$ -norm are obtained applying Lemma 2.2.

Step 2. In order to obtain the final estimation one should observe that, using the first equation of (QIM), one can retrieve $\begin{matrix} ‖ \underline{\nabla} p_{λ} ‖_{X^{'}} ⩽ {‖ \partial_{t}^{2} {\underline{y}}_{λ} ‖}_{X^{'}} + {‖ \underline{div} \underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{y}}_{λ}) ‖}_{X^{'}} ≲ ({‖ \partial_{t}^{2} {\underline{y}}_{λ} ‖}_{H} + ‖ {\underline{y}}_{λ} ‖_{X}) . \end{matrix}$ From the preliminary result given in Section 2.1 (see (2.1)), we find the result of the Theorem for $k = 1$ using the estimations obtained in the first step.

Step 3. The estimation for higher values of k is obtained by successive differentiations of the first equation of (QIM). □

2.3. The incompressible limit

Since λ is large in the application that we consider, it is natural to approximate the solution to (QIM) by the solution obtained at the limit as $λ \to \infty$ . In more detail, if we define $(\underline{y}, p)$ and some corrector functions $({\underline{y}}_{1}, p_{1})$ and $({\underline{y}}_{2}, p_{2})$ such that $\begin{matrix} (2.8) & {\underline{y}}_{λ} = \underline{y} + λ^{- 1} {\underline{y}}_{1} + λ^{- 2} {\underline{y}}_{2} + \dots, p_{λ} = p + λ^{- 1} p_{1} + λ^{- 2} p_{2} + \dots, \end{matrix}$ then a standard formal asymptotic analysis procedure shows that $(\underline{y}, p)$ satisfies a pure incompressible problem, that reads:

Given $\underline{f} \in W_{0}^{1} (H)$ , find $(\underline{y}, p) \in [C^{2} (H) \cap C^{1} (V)] \times C^{0} (L)$ such that $\begin{matrix} (IM) & \{\begin{matrix} \partial_{t}^{2} \underline{y} - \underline{div} (\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} (\underline{y})) - \underline{\nabla} p = \underline{f} & in Ω_{T}, \\ div \underline{y} = 0 & in Ω_{T}, \\ \underline{y} (t = 0) = \underline{0}, \partial_{t} \underline{y} (t = 0) = \underline{0} & in Ω . \end{matrix} \end{matrix}$ Note that $\partial_{t}^{2} \underline{y}$ is continuous with values in $D$ (in particular it has zero divergence). Moreover, note that p is sought in $C^{0} (L)$ , although one could expect $C^{1} (L)$ . This can be explained by the estimation of Theorem 2.3 with $k = 1$ . Indeed, it can be seen that $p_{λ}$ is uniformly bounded with respect to λ, but we do not have such estimate for $\partial_{t} p_{λ}$ . Note that we do not need to add initial conditions on p – as it was the case of (QIM), due to the definition of $p_{λ}$ – since, at $t = 0$ , we have $\begin{matrix} \partial_{t}^{2} \underline{y} (t = 0) = \underline{\nabla} p (t = 0) in Ω, div \partial_{t}^{2} \underline{y} (t = 0) = 0 in Ω . \end{matrix}$ As a consequence, $\begin{matrix} \{\begin{matrix} - Δ p (t = 0) = 0 & in Ω \\ \underline{\nabla} p \cdot \underline{n} (t = 0) = 0 & on \partial Ω \end{matrix} ⟹ p (t = 0) = 0 . \end{matrix}$ Problem (QIM) can be seen as a penalised formulation of (IM). Existence of the solution to system (IM) is deduced from (QI) by taking the limit as $λ \to \infty$ , uniqueness is a consequence of stability estimates (see [6]). Therefore, we can state the following result.

Proposition 2.4.
Let $\underline{f} \in W_{0}^{k} (H)$ with $k ⩾ 1$ . Then, there exists a unique couple $(\underline{y}, p) \in [C^{k + 1} (H) \cap C^{k} (V)] \times C^{k - 1} (L)$ solution to ( IM ).

2.4. A penalised quasi-incompressible formulation

The last formulation under consideration is the one introduced in [7], called (QIP). It corresponds to an approximation by penalisation of the problem (IM), inspired by existing formulations for incompressible fluid dynamics [4]:

Given $\underline{f} \in W_{0}^{1} (H)$ , find the couple $({\underline{\tilde{y}}}_{α}, {\tilde{p}}_{α}) \in [C^{2} (H) \cap C^{1} (X)] \times C^{1} (M)$ such that $\begin{matrix} (QIP) & \{\begin{matrix} \partial_{t}^{2} {\underline{\tilde{y}}}_{α} - \underline{div} (\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{\tilde{y}}}_{α})) - \underline{\nabla} {\tilde{p}}_{α} = \underline{f} & in Ω_{T}, \\ div {\underline{\tilde{y}}}_{α} = - α Δ {\tilde{p}}_{α} & in Ω_{T}, \\ {\underline{\tilde{y}}}_{α} (t = 0) = \underline{0}, \partial_{t} {\underline{\tilde{y}}}_{α} (t = 0) = \underline{0} & in Ω . \end{matrix} \end{matrix}$ We note that (QIP) differs from (QIM) for the introduction of the Laplace operator in the second equation. As a consequence, the pressure ${\tilde{p}}_{α}$ must be sought in a more regular space, in order to give an appropriate meaning to the introduced Laplace operator. Hence, a boundary condition is required for the second equation on the pressure, that now has a trace (we recall that we use homogeneous Dirichlet conditions on the displacement). A common but arbitrary choice is to use homogenous Neumann boundary conditions $\begin{matrix} (2.9) & \underline{\nabla} {\tilde{p}}_{α} \cdot \underline{n} = 0 on \partial Ω \times (0, T) . \end{matrix}$ Consequently, for small values of α, (QIP) corresponds to another approximation of the pure incompressible mixed formulation (IM). Using the same arguments as before, i.e., writing $\begin{matrix} (2.10) & {\underline{\tilde{y}}}_{α} = \underline{y} + α {\tilde{\underline{y}}}_{1} + α^{2} {\tilde{\underline{y}}}_{2} + \dots, {\tilde{p}}_{α} = p + α {\tilde{p}}_{1} + α^{2} {\tilde{p}}_{2} + \dots, \end{matrix}$ one can see that $({\underline{\tilde{y}}}_{α}, {\tilde{p}}_{α})$ is a formal approximation of $(\underline{y}, p)$ , solution to the pure incompressible mixed formulation (IM). The interest of the formulation (QIP) is discussed in [7]. In a nutshell: after space–time discretisation, the space–time discrete pressure field can be computed, at each time step $Δ t$ , from the known displacement field by solving a scalar Laplace equation, and afterwards the updated displacement field is computed explicitly. This offers a real benefit compared to either fully explicit or fully implicit schemes. Note that the parameter α cannot be chosen arbitrarily small, due to a CFL-type stability condition. In practice, to obtain an efficient and accurate numerical scheme, α is chosen proportional to $Δ t^{2}$ .

Concerning the analysis of the formulation (QIP), if we eliminate ${\tilde{p}}_{α}$ in the first equation in (QIP), we retrieve the following variational formulation $\begin{matrix} (2.11) & {(\partial_{t}^{2} {\underline{\tilde{y}}}_{α}, \underline{w})}_{H} + a_{α} ({\underline{\tilde{y}}}_{α}, \underline{w}) = (\underline{f}, \underline{w}), \forall \underline{w} \in X, \end{matrix}$ with, for all $\underline{y}$ and $\underline{w}$ in $X$ , $\begin{matrix} a_{α} (\underline{y}, \underline{w}) : = {(\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} (\underline{y}), \underline{\underline{ε}} (\underline{w}))}_{H} + α^{- 1} b (\underline{y}, \underline{w}), b (\underline{y}, \underline{w}) : = \int_{Ω} (- Δ^{- 1} div \underline{y}) div \underline{w} d Ω, \end{matrix}$ where $- Δ^{- 1} : L \to M$ stands for the inverse Laplace operator on $M$ , equipped with a homogeneous Neumann boundary condition.

Proposition 2.5.
The bilinear form $b : X \times X \mapsto R$ is symmetric semi-definite positive and for all $\underline{y}$ and $\underline{w}$ in $X$ , $\begin{matrix} (2.12) & | b (\underline{y}, \underline{w}) | ⩽ ‖ \underline{y} ‖_{H} ‖ \underline{w} ‖_{H} . \end{matrix}$
Proof.
For any $\underline{y} \in X$ we introduce $q_{\underline{y}} \in N$ as the unique solution of $\begin{matrix} (2.13) & - Δ q_{\underline{y}} = div \underline{y} in Ω, \underline{\nabla} q_{\underline{y}} \cdot \underline{n} = 0 on \partial Ω . \end{matrix}$ We have $b (\underline{y}, \underline{w}) = {(q_{\underline{y}}, div \underline{w})}_{L} = - {(q_{\underline{y}}, Δ q_{\underline{w}})}_{L}$ , hence, using Green’s formula we obtain $\begin{matrix} (2.14) & b (\underline{y}, \underline{w}) = {(\underline{\nabla} q_{\underline{y}}, \underline{\nabla} q_{\underline{w}})}_{H} . \end{matrix}$ Using the fact that from (2.13) we have $‖ \underline{\nabla} q_{\underline{y}} ‖_{H} ⩽ ‖ \underline{y} ‖_{H}$ , we obtain (2.12) and finally (2.14) also shows the symmetric semi-definite positive property. □

From Proposition 2.5 one can see that $a_{α} : X \times X \mapsto R$ is a symmetric, continuous, coercive, bilinear form. Consequently, existence, uniqueness and regularity of the solution of System (QIP) can be retrieved using Theorem 9.1 of [21]. Moreover, using a similar approach to the one used to derive Theorem 2.3, one can retrieve the following theorem: Theorem 2.6.
Let $\underline{f} \in W_{0}^{k} (H)$ with $k ⩾ 1$ . Then, there exists a unique couple $({\underline{\tilde{y}}}_{α}, {\tilde{p}}_{α}) \in [C^{k + 1} (H) \cap C^{k} (X)] \times C^{k} (M)$ solution to ( QIP ). Moreover, $\begin{matrix} (2.15) & {‖ \partial_{t}^{k + 1} {\underline{\tilde{y}}}_{α} ‖}_{L^{\infty} (H)} + {‖ \partial_{t}^{k} {\underline{\tilde{y}}}_{α} ‖}_{L^{\infty} (X)} + α^{\frac{1}{2}} {‖ \underline{\nabla} \partial_{t}^{k} {\tilde{p}}_{α} ‖}_{L^{\infty} (H)} + {‖ \partial_{t}^{k - 1} {\tilde{p}}_{α} ‖}_{L^{\infty} (L)} ≲ 1 . \end{matrix}$
Remark 2.7.
The results presented in the next section also hold if the second equation of (QIP) is replaced by $\begin{matrix} div {\underline{\tilde{y}}}_{α} = - α \underline{div} (\underline{\underline{A}} \underline{\nabla} {\tilde{p}}_{α}) in Ω_{T}, with \underline{\underline{A}} \underline{\nabla} {\tilde{p}}_{α} \cdot \underline{n} = 0 on \partial Ω \times (0, T), \end{matrix}$ where $\underline{\underline{A}} \in C^{1} {(\overline{Ω})}^{d \times d}$ is a definite positive matrix field (uniformly coercive). As suggested in [8, Section 6], this has a practical interest when computing ${\tilde{p}}_{α}$ in a discrete setting.
Remark 2.8.
One can think that (QIP) is constructed as a singular perturbation of (IM) by the term $- α Δ$ . However, Eq. (2.11) shows that the underlying problem solved by the displacement field ${\underline{\tilde{y}}}_{α}$ is a penalised problem with penalisation parameter $α^{- 1}$ and where the divergence of ${\underline{\tilde{y}}}_{α}$ is penalised in a weak norm.

3. Convergence estimates

The objective of this section is to provide an error estimate between the formulations (QIM) and (QIP) in the energy norm. More precisely, we want to show that the difference ${\underline{y}}_{λ} - {\underline{\tilde{y}}}_{α}$ decreases with $λ^{- 1}$ and α in $H^{1}$ -norm in space To do so, we consider (IM) as the reference formulation, and we derive the error estimates associated with the approximation of (IM) by (QIM), i.e. we study $\underline{y} - {\underline{y}}_{λ}$ , or by (QIP), i.e. we study $\underline{y} - {\underline{\tilde{y}}}_{α}$ . Then, we infer the sought estimate using the triangular inequality. The obtained estimates involve non-integer exponents. Therefore, we introduce the non-negative rational numbers $δ (m)$ and $η (m)$ defined as $\begin{matrix} δ (0) : = 0, δ (m) : = \sum_{n = 1}^{m} \frac{1}{2^{n}}, δ (m) \underset{m \to + \infty}{⟶} 1 \end{matrix}$ and $\begin{matrix} η (0) : = 0, η (m) : = \sum_{n = 1}^{m} \frac{1}{4^{n}}, η (m) \underset{m \to + \infty}{⟶} \frac{1}{3} . \end{matrix}$

Remark 3.1.
One could choose to construct a penalised approximation of the displacement ${\underline{y}}_{λ}$ by directly introducing the unknown $({\underline{\tilde{y}}}_{λ, α}, {\tilde{p}}_{λ, α})$ satisfying $\begin{matrix} \partial_{t}^{2} {\underline{\tilde{y}}}_{λ, α} - \underline{div} (\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{\tilde{y}}}_{λ, α})) - \underline{\nabla} {\tilde{p}}_{λ, α} = \underline{f}, div {\underline{\tilde{y}}}_{λ, α} = λ^{- 1} {\tilde{p}}_{λ, α} - α Δ {\tilde{p}}_{λ, α}, \end{matrix}$ and then directly proving that ${\underline{y}}_{λ} - {\underline{\tilde{y}}}_{λ, α}$ goes to zero with $λ^{- 1}$ and α. However, the procedure that we present has the main advantage that, in practice, an exact value of λ is not required to construct the penalised problem. Moreover, the analysis of the difference $\underline{y} - {\underline{y}}_{λ}$ is interesting in order to understand the range of validity of the pure incompressible approximations.

3.1. Convergence estimate of the quasi-incompressible formulation

In order to perform the error analysis of the approximation of (QIM) and (QIP) by triangular inequality, we first compute the error between formulations (IM) and (QIM). Let us define the quantities $\begin{matrix} (3.1) & {\underline{e}}_{λ} = {\underline{y}}_{λ} - \underline{y}, q_{λ} = p_{λ} - p . \end{matrix}$ Let $\underline{f} \in W_{0}^{k} (H)$ . The couple $({\underline{e}}_{λ}, q_{λ}) \in [C^{k + 1} (H) \cap C^{k} (X)] \times C^{k - 1} (L)$ satisfies $\begin{matrix} (3.2) & \{\begin{matrix} \partial_{t}^{2} {\underline{e}}_{λ} - div \underline{\underline{C}} \underline{\underline{ε}} ({\underline{e}}_{λ}) - \underline{\nabla} q_{λ} = \underline{0} & in Ω_{T}, \\ div {\underline{e}}_{λ} = λ^{- 1} q_{λ} + λ^{- 1} p & in Ω_{T}, \\ {\underline{e}}_{λ} (t = 0) = \underline{0}, \partial_{t} {\underline{e}}_{λ} (t = 0) = \underline{0} & in Ω . \end{matrix} \end{matrix}$ To derive the sought convergence estimate, we provide the stability estimates associated with (3.1) in the theorem below.

Theorem 3.2.
Let $\underline{f} \in W_{0}^{ℓ + 2} (H)$ with $ℓ ⩾ 0$ . Then $\begin{matrix} (3.3) & {‖ \partial_{t}^{ℓ + 1} {\underline{e}}_{λ} ‖}_{L^{\infty} (H)} + {‖ \partial_{t}^{ℓ} {\underline{e}}_{λ} ‖}_{L^{\infty} (X)} + λ^{- \frac{1}{2}} {‖ \partial_{t}^{ℓ} q_{λ} ‖}_{L^{\infty} (L)} ≲ λ^{- \frac{1}{2}} . \end{matrix}$ Moreover, if $ℓ ⩾ 1$ we have, for all $0 ⩽ r ⩽ ℓ - 1$ , $\begin{matrix} (3.4) & {‖ \partial_{t}^{r + 1} {\underline{e}}_{λ} ‖}_{L^{\infty} (H)} + {‖ \partial_{t}^{r} {\underline{e}}_{λ} ‖}_{L^{\infty} (X)} ≲ λ^{- \frac{1}{2} - \frac{1}{2} δ (ℓ - r)}, \end{matrix}$ and $\begin{matrix} (3.5) & {‖ \partial_{t}^{r} q_{λ} ‖}_{L^{\infty} (L)} ≲ λ^{- \frac{1}{2} - \frac{1}{2} δ (ℓ - r - 1)} . \end{matrix}$
Proof.
We define the following energy functionals, for every integer $0 ⩽ r ⩽ ℓ$ , $\begin{matrix} E^{(r)} = \frac{1}{2} {‖ \partial_{t}^{r + 1} {\underline{e}}_{λ} ‖}_{H}^{2} + \frac{1}{2} \int_{Ω} \underline{\underline{C}} \underline{\underline{ε}} (\partial_{t}^{r} {\underline{e}}_{λ}) : \underline{\underline{ε}} (\partial_{t}^{r} {\underline{e}}_{λ}) d Ω + \frac{λ^{- 1}}{2} {‖ \partial_{t}^{r} q_{λ} ‖}_{L}^{2} \in C^{1} ([0, T]) . \end{matrix}$ Since we have considered vanishing initial data and both the source term and its derivatives vanish at the initial time, one can show that $\begin{matrix} sup_{t \in [0, T]} E^{(r - 1)} (t) ≲ sup_{t \in [0, T]} E^{(r)} (t) \end{matrix}$ for $1 ⩽ r ⩽ ℓ$ . Finally it is possible to retrieve, thanks to the Korn inequality, $\begin{matrix} {‖ \partial_{t}^{r + 1} {\underline{e}}_{λ} (t) ‖}_{H} + {‖ \partial_{t}^{r} {\underline{e}}_{λ} (t) ‖}_{X} ≲ E^{(r)} (t) . \end{matrix}$

Step 1. We prove (3.3) for the case $ℓ = 0$ . If we perform a scalar product of the first equation in (3.2) and $\partial_{t} {\underline{e}}_{λ}$ , we obtain, using the second equation in (3.2), $\begin{matrix} (3.6) & \frac{d}{d t} E^{(0)} = - λ^{- 1} {(\partial_{t} p, q_{λ})}_{L} ⩽ λ^{- \frac{1}{2}} ‖ \partial_{t} p ‖_{L} λ^{- \frac{1}{2}} ‖ q_{λ} ‖_{L} ≲ λ^{- \frac{1}{2}} ‖ \partial_{t} p ‖_{L} {(E^{(0)})}^{\frac{1}{2}} . \end{matrix}$ By the Grönwall lemma, it is possible to show that, for all $t \in [0, T]$ , $\begin{matrix} (3.7) & {(E^{(0)} (t))}^{\frac{1}{2}} ≲ λ^{- \frac{1}{2}} \int_{0}^{t} {‖ \partial_{t} p (s) ‖}_{L} d s . \end{matrix}$ Since $\underline{f} \in W_{0}^{2} (H)$ by assumption and due to Proposition 2.4, we can assert that $p \in C^{1} (L)$ . Consequently, from Eq. (3.7) we can prove that the energy behaves at worst proportionally to $λ^{- \frac{1}{2}}$ . Hence, using the Korn inequality, we recover the estimation (3.3) for $ℓ = 0$ . The estimations for higher values of ℓ are obtained in a similar way. More precisely, after differentiating with respect to time (3.1) and by similar computations as before, we obtain $\begin{matrix} (3.8) & \frac{d}{d t} E^{(ℓ)} ⩽ λ^{- \frac{1}{2}} {‖ \partial_{t}^{ℓ + 1} p ‖}_{L} λ^{- \frac{1}{2}} {‖ \partial_{t}^{ℓ} q_{λ} ‖}_{L} ≲ λ^{- \frac{1}{2}} {‖ \partial_{t}^{ℓ + 1} p ‖}_{L} {(E^{(ℓ)})}^{\frac{1}{2}} . \end{matrix}$ Since $f \in W_{0}^{ℓ + 2} (H)$ , we have $p \in C^{ℓ + 1} (L)$ . Hence, we obtain (3.3), i.e., $\begin{matrix} (3.9) & sup_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{2}} ≲ λ^{- \frac{1}{2}} . \end{matrix}$

Step 2. We assume now $ℓ ⩾ 1$ . Differentiating with respect to time the first equation of (3.2), we have, for every $t \in [0, T]$ , $\begin{array}{rcl} {‖ \underline{\nabla} \partial_{t}^{ℓ - 1} q_{λ} (t) ‖}_{X^{'}} & ⩽ & {‖ \partial_{t}^{ℓ + 1} {\underline{e}}_{λ} (t) ‖}_{X^{'}} + {‖ div \underline{\underline{C}} \underline{\underline{ε}} (\partial_{t}^{ℓ - 1} {\underline{e}}_{λ}) (t) ‖}_{X^{'}} \\ ≲ & {‖ \partial_{t}^{ℓ + 1} {\underline{e}}_{λ} (t) ‖}_{H} + {‖ \partial_{t}^{ℓ - 1} {\underline{e}}_{λ} (t) ‖}_{X} ≲ sup_{s \in [0, T]} {(E^{(ℓ)} (s))}^{\frac{1}{2}} . \end{array}$ This shows that, thanks to (2.1), $\begin{matrix} (3.10) & {‖ \partial_{t}^{ℓ - 1} q_{λ} ‖}_{L^{\infty} (L)} ≲ sup_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{2}} . \end{matrix}$ Moreover, since $\underline{f} \in W_{0}^{ℓ + 2} (H)$ , we have in particular $\underline{f} \in W_{0}^{ℓ + 1} (H)$ . Therefore, (3.8) holds by formally replacing ℓ by $ℓ - 1$ and we deduce that, using (3.10), $\begin{matrix} \frac{d}{d t} E^{(ℓ - 1)} (t) ≲ λ^{- 1} {‖ \partial_{t}^{ℓ - 1} q_{λ} ‖}_{L^{\infty} (L)} ≲ λ^{- 1} sup_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{2}}, \end{matrix}$ which implies $\begin{matrix} (3.11) & sup_{t \in [0, T]} E^{(ℓ - 1)} (t) ≲ λ^{- 1} sup_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{2}} . \end{matrix}$ Inequality (3.11) relates $E^{(ℓ - 1)}$ and ${(E^{(ℓ)})}^{\frac{1}{2}}$ which allows to derive (3.4) next.

Step 3. We now exploit the recursive relation that can be deduced from (3.9) and (3.11), namely, for $ℓ ⩾ 1$ and all $0 ⩽ r ⩽ ℓ - 1$ , $\begin{matrix} sup_{t \in [0, T]} E^{(ℓ)} (t) ≲ λ^{- 1} and sup_{t \in [0, T]} E^{(r)} (t) ≲ λ^{- 1} sup_{t \in [0, T]} {(E^{(r + 1)} (t))}^{\frac{1}{2}} . \end{matrix}$ Then, by induction one can show the estimate $\begin{matrix} (3.12) & sup_{t \in [0, T]} E^{(r)} (t) ≲ λ^{- 1 - δ (ℓ - r)} \end{matrix}$ from which Eq. (3.4) follows. Estimation (3.5) is obtained as a direct consequence of (3.10) and (3.12). □

3.2. Convergence estimate of the penalised quasi-incompressible formulation

We now compare formulations (IM) and (QIP). Let us consider that $\underline{f} \in W_{0}^{k} (H)$ and define the discrepancies $\begin{matrix} (3.13) & {\underline{\tilde{e}}}_{α} = {\underline{\tilde{y}}}_{α} - \underline{y}, {\tilde{q}}_{α} = {\tilde{p}}_{α} - p, \end{matrix}$ where $({\underline{\tilde{e}}}_{α}, {\tilde{q}}_{α}) \in C^{k + 1} (H) \cap C^{k} (X) \times C^{k - 1} (L)$ is solution to $\begin{matrix} (3.14) & \{\begin{matrix} \partial_{t}^{2} {\underline{\tilde{e}}}_{α} - div \underline{\underline{C}} \underline{\underline{ε}} ({\underline{\tilde{e}}}_{α}) - \underline{\nabla} {\tilde{q}}_{α} = \underline{0} & in Ω_{T}, \\ div {\underline{\tilde{e}}}_{α} = - α Δ {\tilde{q}}_{α} - α Δ p & in Ω_{T}, \\ {\underline{\tilde{e}}}_{α} (t = 0) = \underline{0}, \partial_{t} {\underline{\tilde{e}}}_{α} (t = 0) = \underline{0} & in Ω . \end{matrix} \end{matrix}$ As before, to derive the sought convergence estimate we provide the stability estimates associated with (3.13) in the theorem below.

Theorem 3.3.
Let $\underline{f} \in W_{0}^{ℓ + 2} (H)$ with $ℓ ⩾ 0$ and let $p \in C^{ℓ + 1} (M)$ . Then, $\begin{matrix} (3.15) & {‖ \partial_{t}^{ℓ + 1} {\underline{\tilde{e}}}_{α} ‖}_{L^{\infty} (H)} + {‖ \partial_{t}^{ℓ} {\underline{\tilde{e}}}_{α} ‖}_{L^{\infty} (X)} + α^{\frac{1}{2}} {‖ \partial_{t}^{ℓ} \underline{\nabla} {\tilde{q}}_{α} ‖}_{L^{\infty} (H)} ≲ α^{\frac{1}{2}} . \end{matrix}$ Moreover, if $ℓ ⩾ 1$ and $p \in C^{ℓ} (N)$ we have, for all $0 ⩽ r ⩽ ℓ - 1$ , $\begin{matrix} (3.16) & {‖ \partial_{t}^{r + 1} {\underline{\tilde{e}}}_{α} ‖}_{L^{\infty} (H)} + {‖ \partial_{t}^{r} {\underline{\tilde{e}}}_{α} ‖}_{L^{\infty} (X)} ≲ α^{\frac{1}{2} + \frac{1}{2} η (ℓ - r)} \end{matrix}$ and $\begin{matrix} (3.17) & {‖ \partial_{t}^{r} {\tilde{q}}_{α} ‖}_{L^{\infty} (L)} ≲ α^{\frac{1}{2} + \frac{1}{2} η (ℓ - r - 1)} . \end{matrix}$
Proof.
Following the proof of Theorem 3.2, we define a family of energy functionals associated with Eq. (3.14) that read $\begin{matrix} E^{(r)} = \frac{1}{2} {‖ \partial_{t}^{r + 1} {\underline{\tilde{e}}}_{α} ‖}_{H}^{2} + \frac{1}{2} \int_{Ω} \underline{\underline{C}} \underline{\underline{ε}} (\partial_{t}^{r} {\underline{\tilde{e}}}_{α}) : \underline{\underline{ε}} (\partial_{t}^{r} {\underline{\tilde{e}}}_{α}) d Ω + \frac{α}{2} {‖ \underline{\nabla} \partial_{t}^{r} {\tilde{q}}_{α} ‖}_{H}^{2} \in C^{1} ([0, T]), \end{matrix}$ for every integer $0 ⩽ r ⩽ ℓ$ . As it has already been noticed, we have $\begin{matrix} sup_{t \in [0, T]} E^{(r - 1)} (t) ≲ sup_{t \in [0, T]} E^{(r)} (t) and {‖ \partial_{t}^{r + 1} {\underline{\tilde{e}}}_{α} (t) ‖}_{H} + {‖ \partial_{t}^{r} {\underline{\tilde{e}}}_{α} (t) ‖}_{X} ≲ E^{(r)} (t), \end{matrix}$ due to our choice of initial conditions and thanks to the Korn inequality. Furthermore ${\tilde{q}}_{α} \in C^{ℓ + 1} (M)$ by assumption.

Step 1. By standard energy analysis one can show that $\begin{array}{rcl} \frac{d}{d t} E^{(0)} & = & - α {(\partial_{t} \underline{\nabla} p, \underline{\nabla} {\tilde{q}}_{α})}_{H} \\ (3.18) & ⩽ & α^{\frac{1}{2}} ‖ \partial_{t} \underline{\nabla} p ‖_{H} α^{\frac{1}{2}} ‖ \underline{\nabla} {\tilde{q}}_{α} ‖_{H} ≲ α^{\frac{1}{2}} ‖ \partial_{t} \underline{\nabla} p ‖_{H} {(E^{(0)})}^{\frac{1}{2}} . \end{array}$ This implies by the Grönwall lemma that, for all $t \in [0, T]$ , $\begin{matrix} (3.19) & {(E^{(0)} (t))}^{\frac{1}{2}} ≲ α^{\frac{1}{2}} \int_{0}^{t} ‖ \partial_{t} \underline{\nabla} p ‖_{H} d s, \end{matrix}$ which gives (3.15) for $ℓ = 0$ , since $p \in C^{1} (M)$ by assumption. The estimates for $ℓ > 0$ are deduced in a similar way using the energy identity $\begin{matrix} (3.20) & \frac{d}{d t} E^{(ℓ)} = - α {(\partial_{t}^{ℓ + 1} \underline{\nabla} p, \underline{\nabla} \partial_{t}^{ℓ} {\tilde{q}}_{α})}_{H} ≲ α^{\frac{1}{2}} {‖ \partial_{t}^{ℓ + 1} \underline{\nabla} p ‖}_{H} {(E^{(ℓ)})}^{\frac{1}{2}} \end{matrix}$ and, using the Grönwall lemma, we obtain ${sup}_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{2}} ≲ α^{\frac{1}{2}}$ and deduce (3.15).

Step 2. We assume now $ℓ ⩾ 1$ . Similarly to (3.10), it is possible to show that $\begin{matrix} (3.21) & {‖ \partial_{t}^{(ℓ - 1)} {\tilde{q}}_{α} ‖}_{L^{\infty} (L)} ≲ sup_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{2}} . \end{matrix}$ At this stage, it is no longer possible to mimic the proof of Theorem 3.2, since (3.20) involves high-order derivatives compared to (3.8). Therefore, we must use the assumption that $p \in C^{1} (N)$ and we have, using the Green formulas and the trace inequality (2.2), $\begin{matrix} (3.22) & \begin{matrix} | (\underline{\nabla} {\tilde{q}}_{α}, \partial_{t} \underline{\nabla} p) | & ⩽ ‖ {\tilde{q}}_{α} ‖_{L} ‖ Δ \partial_{t} p ‖_{L^{2} (Ω)} + ‖ {\tilde{q}}_{α} ‖_{L^{2} (\partial Ω)} ‖ \underline{\nabla} \partial_{t} p \cdot \underline{n} ‖_{L^{2} (\partial Ω)} \\ ≲ ‖ {\tilde{q}}_{α} ‖_{L} + ‖ {\tilde{q}}_{α} ‖_{L^{2} (\partial Ω)} ≲ ‖ {\tilde{q}}_{α} ‖_{L} + ‖ {\tilde{q}}_{α} ‖_{L}^{\frac{1}{2}} ‖ \underline{\nabla} {\tilde{q}}_{α} ‖_{H}^{\frac{1}{2}} . \end{matrix} \end{matrix}$ Hence, from the estimate above and (3.20), we deduce that if $p \in C^{ℓ} (N)$ then $\begin{matrix} sup_{t \in [0, T]} E^{(ℓ - 1)} (t) ≲ α {‖ \partial_{t}^{ℓ - 1} {\tilde{q}}_{α} ‖}_{L^{\infty} (L)}^{\frac{1}{2}} ({‖ \partial_{t}^{ℓ - 1} {\tilde{q}}_{α} ‖}_{L^{\infty} (L)}^{\frac{1}{2}} + {‖ \underline{\nabla} \partial_{t}^{ℓ - 1} {\tilde{q}}_{α} ‖}_{L^{\infty} (H)}^{\frac{1}{2}}) . \end{matrix}$ Using (3.15) we have $‖ \underline{\nabla} \partial_{t}^{ℓ - 1} {\tilde{q}}_{α} ‖_{L^{\infty} (H)} ≲ 1$ and, thanks to (3.21), $\begin{matrix} (3.23) & \begin{matrix} sup_{t \in [0, T]} E^{(ℓ - 1)} (t) & ≲ α sup_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{2}} + α sup_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{4}} \\ ≲ α sup_{t \in [0, T]} {(E^{(ℓ)} (t))}^{\frac{1}{4}} . \end{matrix} \end{matrix}$ The inequality that we have obtained relates $E^{(ℓ - 1)}$ and ${(E^{(ℓ)})}^{\frac{1}{4}}$ (this differs from what we have obtained in the proof of Theorem 3.2).

Step 3. We now derive (3.16) by exploiting the following recursive relation: for $ℓ ⩾ 1$ and all $0 ⩽ r ⩽ ℓ - 1$ $\begin{matrix} sup_{t \in [0, T]} E^{(ℓ)} (t)) ≲ α and sup_{t \in [0, T]} E^{(r)} (t) ≲ α sup_{t \in [0, T]} {(E^{(r + 1)} (t))}^{\frac{1}{4}} . \end{matrix}$ By induction one can show the estimate $\begin{matrix} (3.24) & sup_{t \in [0, T]} E^{(r)} (t) ≲ α^{1 + η (ℓ - r)}, \end{matrix}$ and (3.16) directly follows. The estimation (3.17) is obtained as a straightforward consequence of (3.16) and (3.21). □

3.3. Discussions

Using the triangular inequality, from Theorems 3.2 and 3.3 we can deduce a first estimate on the error between the formulations (QIM) and (QIP). In the case of minimal regularity, we have the following results.

Corollary 3.4.
Let $\underline{f} \in W_{0}^{2} (H)$ , and let $p \in C^{1} (M)$ . Let the couple $({\underline{y}}_{λ}, p_{λ})$ be the solution to problem ( QIM ) and $({\underline{\tilde{y}}}_{α}, {\tilde{p}}_{α})$ the solution to problem ( QIP ), then $\begin{matrix} (3.25) & ‖ \partial_{t} {\underline{y}}_{λ} - \partial_{t} {\underline{\tilde{y}}}_{α} ‖_{L^{\infty} (H)} + ‖ {\underline{y}}_{λ} - {\underline{\tilde{y}}}_{α} ‖_{L^{\infty} (X)} ≲ λ^{- \frac{1}{2}} + α^{\frac{1}{2}} . \end{matrix}$

If an arbitrarily high regularity in time is assumed, then one can deduce the following result.
Corollary 3.5.
Let $\underline{f} \in W_{0}^{k} (H)$ and $p \in C^{k - 1} (M) \cap C^{k - 2} (N)$ for every $k \in N^{}$ , then for all* $ε > 0$ we have $\begin{matrix} (3.26) & ‖ \partial_{t} {\underline{y}}_{λ} - \partial_{t} {\underline{\tilde{y}}}_{α} ‖_{L^{\infty} (H)} + ‖ {\underline{y}}_{λ} - {\underline{\tilde{y}}}_{α} ‖_{L^{\infty} (X)} ≲ λ^{- 1 + ε} + α^{\frac{2}{3} - ε} . \end{matrix}$

We can observe that the order of convergence in α is lower than the one in λ. There is a simple technical explanation for this. Inspecting Eq. (3.22) one can see that the boundary contribution $‖ {\tilde{q}}_{α} ‖_{L^{2} (\partial Ω)}$ deteriorates the estimation (3.23) (we obtain a power 1/4 instead of 1/2 as in (3.11)). This a typical phenomenon related to boundary layer effects. It is not surprising that this boundary layer appears, since we recall that the choice of boundary condition (2.9) was arbitrary. Finally note that the regularity assumptions on the pressure field p implicitly require some smoothness of the coefficients, source term and geometry of the domain. We refer to Remark 5.4 for more detail on this aspect.
4. Convergence estimates in $H$ -norm

Estimations of Theorems 3.2 and 3.3 involve a control in $H^{1}$ -norm in space. It is possible to retrieve a higher order of convergence with respect to λ and α in $L^{2}$ -norm in space. Hence, the objective of this section is to obtain these estimates. More precisely, in what follows we show that the error due to the approximation of (IM) by (QIM) and (QIP) is proportional to $λ^{- 1}$ and α, respectively. As previously, the convergence estimate obtained for the penalised problem requires some additional regularity assumptions that are studied in more detail in Section 5.

4.1. Definition of an elasto-static operator

As a preliminary step we introduce the elasto-static operator associated with the underlying elastodynamic problem. Following the approach proposed in [17], we define a continuous operator $S : H \to V$ such that $\begin{matrix} (4.1) & \underline{r} \in H \mapsto S \underline{r} \in V, \end{matrix}$ and the couple $(S \underline{r}, q_{r}) \in V \times L$ is the unique solution to $\begin{matrix} (4.2) & \{\begin{matrix} {(\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} (S \underline{r}), \underline{\underline{ε}} (\underline{w}))}_{H} - {(q_{r}, div \underline{w})}_{L} = {(\underline{r}, \underline{w})}_{H} & \forall \underline{w} \in X, \\ {(div S \underline{r}, z)}_{L} = 0 & \forall z \in L . \end{matrix} \end{matrix}$ Using the Korn inequality and the property $‖ q ‖_{L} ⩽ C ‖ \underline{\nabla} q ‖_{X^{'}}$ (see (2.1)), it is possible to show that there exists $C > 0$ such that, for all $\underline{r} \in H$ $\begin{matrix} (4.3) & ‖ S \underline{r} ‖_{X} + ‖ q_{r} ‖_{L} ⩽ C ‖ \underline{r} ‖_{H} . \end{matrix}$ Note that C only depends on the domain Ω and the constitutive law considered. It can be proven that $S$ is a positive self-adjoint operator for the scalar product in $H$ , hence we can define a norm $‖ \cdot ‖_{*}$ as $\begin{matrix} (4.4) & ‖ \underline{r} ‖_{*}^{2} : = {(\underline{r}, S \underline{r})}_{H}, \end{matrix}$ that is a weaker norm than the norm in $H$ . In what follows, we use the following definition.

Definition 4.1.
The operator $S$ satisfies the hidden regularity property if, for all $\underline{r} \in H$ , $\begin{matrix} \underline{div} \underline{\underline{C}} (x) \underline{\underline{ε}} (S \underline{r}) \in H, q_{r} \in M, \end{matrix}$ and there exists $C > 0$ such that $\begin{matrix} (4.5) & {‖ \underline{div} \underline{\underline{C}} (x) \underline{\underline{ε}} (S \underline{r}) ‖}_{H} + ‖ \underline{\nabla} q_{r} ‖_{H} ⩽ C ‖ \underline{r} ‖_{H} . \end{matrix}$

This property of $S$ is a fundamental aspect to study the convergence property of the penalised formulation. Moreover, it can be used as a sufficient condition in Corollary 3.4. Corollary 4.2.
Let $S$ satisfy the hidden regularity property and let $\underline{f} \in W_{0}^{k} (H)$ . Then, the solution p to problem ( IM ) belongs to $C^{k - 1} (M)$ and the estimation ( 3.25 ) of Corollary 3.4 holds.
Proof.
We observe that $(\underline{y}, p) = (S \underline{r}, q_{r})$ , with $\underline{r} = \underline{f} - \partial_{t}^{2} \underline{y} \in C^{k - 1} (H)$ , up to modifications of $\underline{f}$ on zero-measure sets. The result follows due to the assumption that $S$ satisfies the hidden regularity property. □

4.2. Convergence estimate of the quasi-incompressible formulation

For the rest of this section we use the same notation as in Section 3. First, we analyse again the approximation of (IM) by formulation (QIM) from a new perspective. Preliminarily, let us consider a Helmholtz decomposition of the solution ${\underline{y}}_{λ}$ to formulation (QIM) [24]. We can introduce a scalar field ${\hat{p}}_{λ}$ such that $\begin{matrix} (4.6) & {\underline{y}}_{λ} = - λ^{- 1} \underline{\nabla} {\hat{p}}_{λ} + {\underline{y}}_{λ}^{0} in Ω_{T}, \end{matrix}$ with ${\underline{y}}_{λ}^{0} (t) \in D$ (in particular its divergence and normal trace are zero) and ${\hat{p}}_{λ} (t) \in M$ is given by $\begin{matrix} (4.7) & \{\begin{matrix} - Δ {\hat{p}}_{λ} = p_{λ} & in Ω_{T}, \\ \underline{\nabla} {\hat{p}}_{λ} \cdot \underline{n} = 0 & in \partial Ω \times (0, T) . \end{matrix} \end{matrix}$ We introduce ${\underline{r}}_{λ} : = {\underline{y}}_{λ}^{0} - \underline{y}$ and we underline that, if $\underline{f} \in W_{0}^{k} (H)$ , then $\begin{matrix} (4.8) & ({\underline{r}}_{λ}, {\hat{p}}_{λ}) \in C^{k} (D) \times C^{k} (M), ({\underline{e}}_{λ}, q_{λ}) \in [C^{k + 1} (H) \cap C^{k} (X)] \times C^{k - 1} (L), \end{matrix}$ where the couple $({\underline{e}}_{λ}, q_{λ})$ is solution to (3.2) and also satisfies $\begin{matrix} (4.9) & \{\begin{matrix} \partial_{t}^{2} {\underline{e}}_{λ} - \underline{div} (\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{e}}_{λ})) - \underline{\nabla} q_{λ} = \underline{0} & in Ω_{T}, \\ - Δ {\hat{p}}_{λ} = p_{λ} & in Ω_{T}, \\ {\underline{e}}_{λ} = - λ^{- 1} \underline{\nabla} {\hat{p}}_{λ} + {\underline{r}}_{λ} & in Ω_{T} . \end{matrix} \end{matrix}$ We now provide the stability analysis of (4.9) and show the following result.

Theorem 4.3.
Let $\underline{f} \in W_{0}^{ℓ + 2} (H)$ , $ℓ ⩾ 0$ . Then $\begin{matrix} (4.10) & {‖ \partial_{t}^{ℓ} {\underline{e}}_{λ} ‖}_{L^{\infty} (H)} ≲ λ^{- 1} . \end{matrix}$
Proof.
We choose $ℓ = 0$ for the proof, but the general case can be easily deduced. From (4.8) one can see that ${\underline{e}}_{λ} \in C^{3} (H) \cap C^{2} (X)$ and $q_{λ} \in C^{1} (L)$ . A same property holds for $({\underline{y}}_{λ}, p_{λ})$ with in addition $p_{λ} \in C^{2} (L)$ (due to the second equation in (QIM)) and, therefore, $\begin{matrix} ({\underline{r}}_{λ}, {\hat{p}}_{λ}) \in C^{2} (H) \times C^{2} (L) . \end{matrix}$ We choose $S \partial_{t} {\underline{r}}_{λ} \in C^{1} (V) \subset C^{0} (X)$ as a test function (in the dual sense of $X$ ) in the first equation of (4.9). Since by definition $S {\underline{r}}_{λ} = 0$ along the boundary, we get $\begin{matrix} (4.11) & {(\partial_{t}^{2} {\underline{e}}_{λ}, S \partial_{t} {\underline{r}}_{λ})}_{H} + {(\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{e}}_{λ}), \underline{\underline{ε}} (S \partial_{t} {\underline{r}}_{λ}))}_{H} = {⟨ \underline{\nabla} q_{λ}, S \partial_{t} {\underline{r}}_{λ} ⟩}_{X^{'}, X} . \end{matrix}$ The right-hand side of Eq. (4.11) vanishes, due to the fact that $div S \partial_{t} {\underline{r}}_{λ} = 0$ inside the domain. From Eq. (4.2) we can write, for all $t \in [0, T]$ , $\begin{matrix} (4.12) & \begin{matrix} {(\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{e}}_{λ}), \underline{\underline{ε}} (S \partial_{t} {\underline{r}}_{λ}))}_{H} & = - {⟨ \underline{div} (\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} (S \partial_{t} {\underline{r}}_{λ})), {\underline{e}}_{λ} ⟩}_{X^{'}, X} \\ = - {⟨ \underline{\nabla} \partial_{t} q_{r_{λ}}, {\underline{e}}_{λ} ⟩}_{X^{'}, X} + {(\partial_{t} {\underline{r}}_{λ}, {\underline{e}}_{λ})}_{H}, \end{matrix} \end{matrix}$ where we have used the property ${⟨ \partial_{t} {\underline{r}}_{λ}, {\underline{e}}_{λ} ⟩}_{X^{'}, X} = {(\partial_{t} {\underline{r}}_{λ}, {\underline{e}}_{λ})}_{H}$ (see Theorem 1, Section 5.9 of [14]). We emphasize that the first term at the right-hand side of Eq. (4.12) is a duality pairing in $X$ , since $\partial_{t} q_{r_{λ}} (t) \in L$ by definition of $S$ . Using the third equation in system (4.9) and Eq. (4.12), we can rewrite Eq. (4.11) as $\begin{matrix} (4.13) & {(λ^{- 1} \partial_{t}^{2} \underline{\nabla} {\hat{p}}_{λ} + \partial_{t}^{2} {\underline{r}}_{λ}, S \partial_{t} {\underline{r}}_{λ})}_{H} - {⟨ \underline{\nabla} \partial_{t} q_{r_{λ}}, {\underline{e}}_{λ} ⟩}_{X^{'}, X} + {(\partial_{t} {\underline{r}}_{λ}, {\underline{e}}_{λ})}_{H} = 0 . \end{matrix}$ Using again the fact that by definition $div S \partial_{t} {\underline{r}}_{λ} = 0$ , Eq. (4.13) can be further simplified into $\begin{matrix} (4.14) & {(\partial_{t}^{2} {\underline{r}}_{λ}, S \partial_{t} {\underline{r}}_{λ})}_{H} + {(\partial_{t} {\underline{r}}_{λ}, {\underline{e}}_{λ})}_{H} = {⟨ \underline{\nabla} \partial_{t} q_{r_{λ}}, {\underline{e}}_{λ} ⟩}_{X^{'}, X} = - {(\partial_{t} q_{r_{λ}}, div {\underline{e}}_{λ})}_{H} . \end{matrix}$ We use again the second and third equation in (4.9) to infer that Eq. (4.14) is equivalent to (see (4.4) for the definition of $‖ \cdot ‖_{}^{2}$ ) $\begin{matrix} (4.15) & \frac{1}{2} \frac{d}{d t} (‖ \partial_{t} {\underline{r}}_{λ} ‖_{}^{2} + ‖ {\underline{r}}_{λ} ‖_{H}^{2}) = - λ^{- 1} {(\partial_{t} q_{r_{λ}}, p_{λ})}_{H} . \end{matrix}$ Integrating by parts in time and considering zero initial conditions, we retrieve $\begin{matrix} (4.16) & ‖ \partial_{t} {\underline{r}}_{λ} ‖_{}^{2} + ‖ {\underline{r}}_{λ} ‖_{H}^{2} = - 2 \int_{0}^{t} λ^{- 1} {(q_{r_{λ}}, \partial_{t} p_{λ})}_{L} d s + 2 λ^{- 1} {(q_{r_{λ}}, p_{λ})}_{L} . \end{matrix}$ Taking into account the Cauchy–Schwarz inequality and Eq. (4.3), we obtain $\begin{matrix} ‖ \partial_{t} {\underline{r}}_{λ} ‖_{}^{2} + ‖ {\underline{r}}_{λ} ‖_{H}^{2} & ≲ \int_{0}^{t} λ^{- 1} ‖ q_{r_{λ}} ‖_{L} ‖ \partial_{t} p_{λ} ‖_{L} d s + λ^{- 1} ‖ q_{r_{λ}} ‖_{L} ‖ p_{λ} ‖_{L} \\ ≲ \int_{0}^{t} λ^{- 1} ‖ {\underline{r}}_{λ} ‖_{H} ‖ \partial_{t} p_{λ} ‖_{L} d s + λ^{- 1} ‖ {\underline{r}}_{λ} ‖_{H} ‖ p_{λ} ‖_{L} . \end{matrix}$ Since $f \in W_{0}^{2} (H)$ , by application of Theorem 2.3 we have $\begin{matrix} (4.17) & ‖ p_{λ} ‖_{L^{\infty} (L)} + ‖ \partial_{t} p_{λ} ‖_{L^{\infty} (L)} ≲ 1, \end{matrix}$ and therefore, using standard estimation techniques (e.g. Young’s inequality), we obtain $\begin{matrix} (4.18) & ‖ {\underline{r}}_{λ} ‖_{L^{\infty} (H)} ≲ λ^{- 1} . \end{matrix}$ In order to retrieve an estimation on ${\underline{e}}_{λ}$ , we recall that ${\underline{e}}_{λ} = λ^{- 1} \underline{\nabla} {\hat{p}}_{λ} + {\underline{r}}_{λ}$ . Furthermore, since $- Δ {\hat{p}}_{λ} = p_{λ}$ , by the Poincaré-Wirtinger inequality and (4.17) we have $\begin{matrix} {‖ \underline{\nabla} {\hat{p}}_{λ} (t) ‖}_{L^{\infty} (H)} ≲ {‖ p_{λ} (t) ‖}_{L^{\infty} (L)} ≲ 1 . \end{matrix}$ Consequently, from Eq. (4.18) we obtain the following estimate on ${\underline{e}}_{λ}$ , $\begin{matrix} (4.19) & ‖ {\underline{e}}_{λ} ‖_{L^{\infty} (H)} ⩽ λ^{- 1} ‖ \underline{\nabla} {\hat{p}}_{λ} ‖_{L^{\infty} (H)} + ‖ {\underline{r}}_{λ} ‖_{L^{\infty} (H)} ≲ λ^{- 1}, \end{matrix}$ thus concluding the proof. □

4.3. Convergence estimate of the penalised quasi-incompressible formulation

We now consider the discrepancy between the solution ${\underline{\tilde{y}}}_{α}$ to problem (QIP) and $\underline{y}$ , solution to the pure incompressible mixed problem (IM). In an analogous way to Eq. (4.6), we can introduce the Helmholtz decomposition for the solution ${\underline{\tilde{y}}}_{α}$ to formulation (QIP). We obtain $\begin{matrix} (4.20) & {\underline{\tilde{y}}}_{α} = - \underline{\nabla} {\tilde{p}}_{α} + {\underline{\tilde{y}}}_{α}^{0}, in Ω_{T}, \end{matrix}$ with ${\underline{\tilde{y}}}_{α}^{0} (t) \in D$ and ${\tilde{p}}_{α}$ solution to (QIP). We also define ${\underline{\tilde{r}}}_{α} : = {\underline{\tilde{y}}}_{α}^{0} - \underline{y}$ . Note that, under the assumption $\underline{f} \in W_{0}^{k} (H)$ , we have $\begin{matrix} {\underline{\tilde{r}}}_{α} \in C^{k} (D) and ({\underline{\tilde{e}}}_{α}, {\tilde{q}}_{α}) \in [C^{k + 1} (H) \cap C^{k} (X)] \times C^{k - 1} (M), \end{matrix}$ where $({\underline{\tilde{e}}}_{α}, {\tilde{q}}_{α})$ is solution to (3.14) and also satisfies $\begin{matrix} (4.21) & \{\begin{matrix} \partial_{t}^{2} {\underline{\tilde{e}}}_{α} - \underline{div} (\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} ({\underline{\tilde{e}}}_{α})) - \underline{\nabla} {\tilde{q}}_{α} = \underline{0} & in Ω_{T}, \\ {\underline{\tilde{e}}}_{α} = - α \underline{\nabla} {\tilde{p}}_{α} + {\underline{\tilde{r}}}_{α} & in Ω_{T} . \end{matrix} \end{matrix}$ The stability analysis of (4.21) leads to the following result.

Theorem 4.4.
Let $\underline{f} \in W_{0}^{ℓ + 2} (H)$ , $ℓ ⩾ 0$ , and let $S$ satisfy the hidden regularity property. Then, ${\underline{\tilde{e}}}_{α}$ defined by Eq. ( 3.13 ) satisfies $\begin{matrix} (4.22) & {‖ \partial_{t}^{ℓ} {\underline{\tilde{e}}}_{α} ‖}_{L^{\infty} (H)} ≲ α . \end{matrix}$
Proof.
The proof of Theorem 4.4 follows the main arguments of the proof of Theorem 4.3. Therefore, we only provide the main ideas of the proof and point out where the additional assumption that $S$ satisfies the hidden regularity property is required. We choose $ℓ = 0$ and we introduce the couple $(S {\underline{\tilde{r}}}_{α}, {\tilde{q}}_{r_{α}}) \in C^{2} (X) \times C^{2} (L)$ solution to $\begin{matrix} (4.23) & \{\begin{matrix} - \underline{div} (\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} (S {\underline{\tilde{r}}}_{α})) + \underline{\nabla} {\tilde{q}}_{r_{α}} = {\underline{\tilde{r}}}_{α} & in Ω_{T}, \\ div S {\underline{\tilde{r}}}_{α} = 0 & in Ω_{T} . \end{matrix} \end{matrix}$ Note however that, by Definition 4.1 and assuming that $S$ satisfies the hidden regularity property, ${\tilde{q}}_{r_{α}} \in C^{2} (M)$ and $\begin{matrix} ‖ {\tilde{q}}_{r_{α}} ‖_{L^{\infty} (M)} ≲ ‖ {\underline{\tilde{r}}}_{α} ‖_{L^{\infty} (H)} . \end{matrix}$ First, we retrieve an analogous estimation to Eq. (4.16) that is suitable for this formulation. To do so, we test the first equation in Eq. (4.21) with $S \partial_{t} {\underline{\tilde{r}}}_{α}$ . After similar considerations to those outlined in the previous proof, we derive (see (4.4) for the definition of $‖ \cdot ‖_{}^{2}$ ) $\begin{matrix} (4.24) & \frac{1}{2} \frac{d}{d t} (‖ \partial_{t} {\underline{r}}_{λ} ‖_{}^{2} + ‖ {\underline{r}}_{λ} ‖_{H}^{2}) = {⟨ \underline{\nabla} \partial_{t} {\tilde{q}}_{r_{α}}, {\underline{\tilde{e}}}_{α} ⟩}_{X^{'}, X} = {(\underline{\nabla} \partial_{t} {\tilde{q}}_{r_{α}}, {\underline{\tilde{e}}}_{α})}_{H} . \end{matrix}$ By definition of ${\underline{\tilde{e}}}_{α}$ we deduce $\begin{matrix} (4.25) & ‖ \partial_{t} {\underline{\tilde{r}}}_{α} ‖_{}^{2} + ‖ {\underline{\tilde{r}}}_{α} ‖_{H}^{2} = 2 \int_{0}^{t} α {(\underline{\nabla} {\tilde{q}}_{r_{α}}, \underline{\nabla} \partial_{t} {\tilde{p}}_{α})}_{H} d s - 2 α {(\underline{\nabla} {\tilde{q}}_{r_{α}}, \underline{\nabla} {\tilde{p}}_{α})}_{H} . \end{matrix}$ Using the Cauchy–Schwarz inequality and the property that $S$ satisfies the hidden regularity property, we obtain $\begin{matrix} (4.26) & ‖ \partial_{t} {\underline{\tilde{r}}}_{α} ‖_{}^{2} + ‖ {\underline{\tilde{r}}}_{α} ‖_{H}^{2} ≲ α \int_{0}^{t} ‖ {\underline{\tilde{r}}}_{α} ‖_{H} ‖ \partial_{t} \underline{\nabla} {\tilde{p}}_{α} ‖_{H} d s + α ‖ {\underline{\tilde{r}}}_{α} ‖_{H} ‖ \underline{\nabla} {\tilde{p}}_{α} ‖_{H} . \end{matrix}$ Since $ℓ = 0$ , from Theorem 3.3 we retrieve that $\begin{matrix} ‖ \partial_{t} \underline{\nabla} {\tilde{p}}_{α} ‖_{L^{\infty} (H)} + ‖ \underline{\nabla} {\tilde{p}}_{α} ‖_{L^{\infty} (H)} ≲ 1 . \end{matrix}$ By a similar reasoning to the previous proof, we get $\begin{matrix} ‖ {\underline{\tilde{r}}}_{α} ‖_{L^{\infty} (H)} ≲ α ⟹ ‖ {\underline{\tilde{e}}}_{α} ‖_{L^{\infty} (H)} ⩽ α ‖ \underline{\nabla} {\tilde{p}}_{α} ‖_{L^{\infty} (H)} + ‖ {\underline{\tilde{r}}}_{α} ‖_{L^{\infty} (H)} ≲ α, \end{matrix}$ that is the result of the theorem. □

4.4. Discussions and conclusion

Using the results provided in Theorems 4.3 and 4.4 and by triangular inequality, we are finally able to retrieve a higher-order estimate in the norm on $H$ of the error performed if we approximate (QIM) by (QIP).

Corollary 4.5.
Let $\underline{f} \in W_{0}^{2} (H)$ and let $S$ satisfy the hidden regularity property. Then, $\begin{matrix} (4.27) & ‖ {\underline{y}}_{λ} - {\underline{\tilde{y}}}_{α} ‖_{L^{\infty} (H)} ≲ λ^{- 1} + α . \end{matrix}$

As already mentioned, in a fully discrete context, α should be chosen proportional to $Δ t^{2}$ in order to obtain an efficient and accurate numerical scheme. Assuming $λ^{- 1} < < α$ and in light of the estimates (4.22) and (4.27), we expect the fully discrete numerical scheme of [7] to perform with an error in $Δ t^{2}$ compared to the discretisation of the quasi-incompressible formulation and/or the purely incompressible formulation. This is exactly what was observed in the numerical convergence experiments of [7].
5. Further insights on the hidden regularity property

There is room for improvement of Corollaries 4.2 and 4.5 since an a priori regularising property of the elasto-static operator $S$ is assumed. The objective of this section is to derive sufficient conditions on the properties of the elastic tensor $\underline{\underline{C}}$ and the domain geometry so that it can be proved that $S$ satisfies the hidden regularity property. Concerning the description of the domain, we need to introduce the following standard definition: a bounded Lipschitz domain is said to be of class $C^{k, 1}$ if its boundary can be parametrised by Lipschitz-continuous local maps for which the first k derivatives are Lipschitz-continuous. We refer the reader to the Definition 1.2.1.1 of [16] for a complete definition (or the proof in Appendix).

5.1. Isotropic elasticity

If the medium considered is isotropic, then the elasticity tensor $\underline{\underline{C}}$ is given by $\begin{matrix} (5.1) & {\underline{\underline{C}}}_{iso} (\underline{x}) \underline{\underline{ε}} = \tilde{λ} (\underline{x}) (\underline{\underline{ε}} : \underline{\underline{I}}) \underline{\underline{I}} + 2 μ (\underline{x}) \underline{\underline{ε}} \end{matrix}$ where $\underline{\underline{I}}$ is the second-order identity tensor, $μ (\underline{x})$ is the shear modulus and $\tilde{λ} (\underline{x})$ corresponds to small variations of the high bulk modulus λ.

Theorem 5.1.
If Ω is either a convex polygonal in $R^{2}$ or a domain of class $C^{1, 1}$ in $R^{d}$ , and if $\underline{\underline{C}} (\underline{x}) \equiv {\underline{\underline{C}}}_{iso} (\underline{x})$ with $\begin{matrix} (5.2) & (\tilde{λ}, μ) \in W^{1, \infty} {(Ω)}^{2}, \tilde{λ} (\underline{x}) + 2 μ (\underline{x}) ⩾ c, μ (\underline{x}) ⩾ c \end{matrix}$ for some constant $c > 0$ , then $S$ satisfies the hidden regularity property.
Proof.
We follow the approach proposed in [5]. The main idea is to reduce our problem to the analysis of a Stokes problem and use well-known $H^{2}$ -regularity results. In the case of isotropic elasticity we have, since $div S \underline{r} = 0$ , $\begin{matrix} {\underline{\underline{C}}}_{iso} (\underline{x}) \underline{\underline{ε}} (S \underline{r}) = 2 μ (\underline{x}) \underline{\underline{ε}} (S \underline{r}) . \end{matrix}$ Therefore, $(S \underline{r}, q_{r}) \in X \times L$ is solution to $\begin{matrix} (5.3) & \{\begin{matrix} - \underline{div} \underline{\underline{ε}} (S \underline{r}) + \underline{\nabla} \frac{q_{r}}{μ} = \frac{1}{μ} (\underline{r} + \underline{\underline{ε}} (S \underline{r}) \underline{\nabla} μ - \frac{q_{r}}{μ} \underline{\nabla} μ) = : \underline{s} & in Ω, \\ div S \underline{r} = 0 & in Ω . \end{matrix} \end{matrix}$ We have the algebraic relation $\underline{div} \underline{\underline{ε}} (S \underline{r}) = Δ S \underline{r} + \underline{\nabla} div S \underline{r}$ , hence Problem (5.3) is equivalent to: find $(S \underline{r}, p_{r}) \in X \times L$ solution to $\begin{matrix} - Δ S \underline{r} + \underline{\nabla} p_{r} = \underline{s}, div S \underline{r} = 0, in Ω . \end{matrix}$ Due to the results of [22] or [3] (Theorem 4.6), if the assumptions of the theorem hold, this problem has a unique solution and there exists $C > 0$ such that, for all $\underline{r} \in H$ , $\begin{matrix} ‖ S \underline{r} ‖_{H^{2} (Ω)} + ‖ p_{r} ‖_{M} ⩽ C ‖ \underline{s} ‖_{H} . \end{matrix}$ Finally, thanks to the definition of $\underline{s}$ and (4.3), it can be observed that $‖ \underline{s} ‖_{H} ⩽ C ‖ \underline{r} ‖_{H}$ . Combining the property $q_{r} = μ p_{r} + μ c_{q}$ for some constant $c_{q}$ and the regularity property of (5.2), we obtain inequality (4.5). □
Remark 5.2.
In (IM) we have $div \underline{y} = \underline{\underline{ε}} (\underline{y}) : \underline{\underline{I}} = 0$ . Hence, the value of $\tilde{λ} (\underline{x})$ can be chosen arbitrarily without changing the solution $(\underline{y}, p)$ . Therefore, in the definition of the penalised formulation (QIP), the value of the parameter can also be chosen arbitrarily (e.g. $\tilde{λ} (\underline{x}) \equiv 0$ ) and Corollaries 3.4, 3.5 and 4.5 still hold.

5.2. Anisotropic elasticity

Anisotropic elasticity tensors can be equivalently represented by a set of $d^{2}$ matrices of dimension $d \times d$ , that we denote ${{\underline{\underline{A}}}_{k ℓ} (\underline{x})}_{k ℓ}$ , such that $\begin{matrix} (5.4) & {(\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} (\underline{v}), \underline{\underline{ε}} (\underline{w}))}_{H} = \sum_{k = 1}^{d} \sum_{ℓ = 1}^{d} {({\underline{\underline{A}}}_{k ℓ} (\underline{x}) \underline{\nabla} v_{ℓ}, \underline{\nabla} w_{k})}_{L}, \end{matrix}$ with $\underline{v} = (v_{1}, \dots, v_{d})$ and $\underline{w} = (w_{1}, \dots, w_{d})$ . We refer to [12] for further detail on this notation. We have not found in the literature any results – similar to Theorem 5.1 – that deal with incompressible elasticity in the general case of heterogeneous anisotropic materials. Nevertheless, one can cite numerous authors for the analysis of the general elliptic system (see in particular [2,13,16]), as well as the specific treatment of Stokes problem (see [3, Proposition 4.3] and reference therein). Using the techniques introduced in [16] together with the presentation proposed in [13] we prove in Appendix the theorem below. It states that $S$ satisfies the hidden regularity property if the parameters and the domain parametrisation are smooth enough.

Theorem 5.3.
Assume that the domain Ω is of class $C^{1, 1}$ and $\underline{\underline{C}}$ satisfies ( 2.4 ) as well as ${A_{k ℓ}} \subset W^{1, \infty} {(Ω)}^{d \times d}$ , then $S$ satisfies the hidden regularity property.
Remark 5.4.
By following the same ideas of the proof given in the Appendix, it is reasonable to conjecture that if the assumptions of the theorem above hold and, in addition,
Ω is of class $C^{2, 1}$ ;

${A_{k ℓ}} \subset W^{2, \infty} {(Ω)}^{d \times d}$ ;

$\underline{f} \in W_{0}^{k} (H^{1} {(Ω)}^{d})$ ;
then p belongs to $C^{k - 1} (N)$ . This enables to verify the assumptions of Theorem 3.3 and Corollary 3.5.

5.3. Application to fibered media

One application of this work is the modelling of wave propagation in the heart or in muscles in general. Such fibered media are suitably described by transverse isotropic elasticity tensors of the form $\begin{matrix} {\underline{\underline{C}}}_{Tiso} (\underline{x}) \underline{\underline{ε}} = \tilde{λ} (\underline{x}) (\underline{\underline{ε}} : \underline{\underline{I}}) \underline{\underline{I}} + μ (\underline{x}) \underline{\underline{ε}} + (\underline{\underline{T}} (\underline{x}) : \underline{\underline{ε}}) \underline{\underline{T}} (\underline{x}), \underline{\underline{T}} (\underline{x}) = \underline{τ} (\underline{x}) \otimes \underline{τ} (\underline{x}), \end{matrix}$ where $\underline{τ} (\underline{x})$ is the fibre direction. Obviously, for any $τ (\underline{x}) \in R^{d}$ the coercivity property (2.4) holds if it holds for $τ (\underline{x}) \equiv \underline{0}$ . By a straightforward application of Theorem 5.3 we have the following result.

Corollary 5.5.

Let the domain Ω be of class $C^{1, 1}$ , let $\underline{\underline{C}} (\underline{x}) \equiv {\underline{\underline{C}}}_{Tiso}$ and ( 5.2 ) hold as well as $\begin{matrix} \underline{τ} \in W^{1, \infty} {(Ω)}^{3} . \end{matrix}$ Then, $S$ satisfies the hidden regularity property.

Footnotes

Proof of Theorem 5.3

This appendix is dedicated to the proof of Theorem 5.3. Although most of the material presented here is rather classical, we have not found in the literature the sought estimate (i.e. Eq. (4.5)) and therefore, for the sake of completeness, we prove here the result. We consider the coercive bilinear form $\begin{matrix} (A.1) & a (\underline{v}, \underline{w}) = {(\underline{\underline{C}} (\underline{x}) \underline{\underline{ε}} (\underline{v}), \underline{\underline{ε}} (\underline{w}))}_{H} = \sum_{k = 1}^{d} \sum_{ℓ = 1}^{d} {({\underline{\underline{A}}}_{k ℓ} (\underline{x}) \underline{\nabla} v_{ℓ}, \underline{\nabla} w_{k})}_{H} . \end{matrix}$ We will prove an a priori stronger result than the mentioned hidden regularity: our objective is to show that the unique solution $(\underline{v}, q) \equiv (S \underline{r}, q_{r}) \in X \times L$ to $\begin{matrix} (A.2) & \{\begin{matrix} a (\underline{v}, \underline{w}) - {(q, div \underline{w})}_{L} = {(\underline{r}, \underline{w})}_{H} & \forall \underline{w} \in X, \\ {(div \underline{v}, z)}_{L} = 0 & \forall z \in L, \end{matrix} \end{matrix}$ with $\underline{r} \in H$ that satisfies $\begin{matrix} (A.3) & ‖ \underline{v} ‖_{X} + ‖ q ‖_{L} ⩽ C ‖ \underline{r} ‖_{H}, \end{matrix}$ for some $C > 0$ depending only on Ω and the tensor $\underline{\underline{C}}$ , also satisfies $(\underline{v}, q) \in H^{2} (Ω) \times M$ with norms uniformly bounded by $‖ \underline{r} ‖_{H}$ . We will proceed by local arguments, the main difficulty being to show that both the $H^{2}$ -regularity and a local estimate of the form (A.3) hold in neighbourhoods of any point of the boundary. We tackle this difficulty for the case $d = 2$ for the sake of clarity. The case $d = 3$ does not provide any additional difficulty. Before giving the details of the proof, we first state some important preliminary definitions.

Preliminary: The strongly elliptic condition. For symmetric elasticity tensors it can be proved (see [23, Proposition 3.10 Chap. 4]) that the coercivity property (2.4) implies the weaker condition of strong ellipticity.

Condition (A.4) is the standard assumption that one can find in the literature on regularity of elliptic systems. In general, it allows to use the Gårding inequality (see e.g. [23, Chapter 6, Proposition 1.5] to show that the bilinear form $a (\cdot, \cdot)$ is elliptic, although in our case the Korn inequality is sufficient. Finally, Eq. (A.4) is also used to show that second-order normal derivatives close to the boundary have the appropriate regularity (see Step 5 below).

Preliminary: Boundary parametrisation. We now give a more detailed description of the boundary pertaining to domains of class $C^{1, 1}$ (see again Definition 1.2.1.1 of [16]). For every ${\underline{x}}_{0} \in \partial Ω$ there exists a neighbourhood V of ${\underline{x}}_{0}$ and a new orthogonal coordinate system $\hat{\underline{x}}$ defined by $\begin{matrix} \hat{\underline{x}} (\underline{x}) = \underline{\underline{R}} (\underline{x} - {\underline{x}}_{R}), with {\underline{\underline{R}}}^{- 1} = {\underline{\underline{R}}}^{t}, \end{matrix}$ such that V is a box in this new coordinate system $\begin{matrix} V = {\hat{\underline{x}} ∣ \hat{\underline{x}} \in (- a_{1}, a_{1}) \times (- a_{2}, a_{2})}, \hat{\underline{x}} (\underline{x}) = \underline{\underline{R}} (\underline{x} - {\underline{x}}_{R}), with {\underline{\underline{R}}}^{- 1} = {\underline{\underline{R}}}^{t}, \end{matrix}$ and there exists a Lipschitz continuous function $φ : [- a_{1}, a_{1}] \to [- \frac{a_{2}}{2}, \frac{a_{2}}{2}]$ with Lipschitz first derivatives such that $\begin{matrix} V_{Ω} = Ω \cap V = {\hat{\underline{x}} \in V ∣ {\hat{x}}_{2} < φ ({\hat{x}}_{1})} and V_{Γ} = \partial Ω \cap V = {\hat{\underline{x}} \in V ∣ {\hat{x}}_{2} = φ ({\hat{x}}_{1})} . \end{matrix}$ Then, we denote by Ψ and Φ the mappings $\begin{matrix} Ψ (\hat{\underline{x}}) = {({\hat{x}}_{1}, {\hat{x}}_{2} - φ ({\hat{x}}_{1}))}^{t} in V, Φ (\hat{\underline{y}}) = {({\hat{y}}_{1}, {\hat{y}}_{2} + φ ({\hat{y}}_{1}))}^{t} in U, \end{matrix}$ where we have defined the domain $U = Ψ (V)$ and denoted by $\hat{\underline{y}}$ any point in U. We have $Ψ \in C^{1, 1} (V)$ and $Φ = Ψ^{- 1} \in C^{1, 1} (U)$ . The deformation gradient of the mapping Φ is given by $\begin{matrix} \underline{\underline{F}} (\hat{\underline{y}}) = \underline{\underline{\nabla}} Φ (\hat{\underline{y}}) = (\begin{matrix} 1 & 0 \\ φ^{'} ({\hat{y}}_{1}) & 1 \end{matrix}) \in W^{1, \infty} {(U)}^{2 \times 2} with det \underline{\underline{F}} = 1 . \end{matrix}$ Note that Φ is a volume-preserving mapping and in particular, if $q \in L_{0}^{2} (V)$ (i.e. the set of $L^{2}$ functions of V with zero average), then $\hat{q} (\cdot) = q ({\underline{\underline{R}}}^{t} Φ (\cdot) + {\underline{x}}_{R}) \in L_{0}^{2} (U)$ . We refer the reader to Fig. 1 for an illustrative example of the boundary parametrisation.

Preliminary: Definition of embedded subdomains. We define $\begin{matrix} U_{Ω} = Ψ (V_{Ω}), U_{Γ} = Ψ (V_{Γ}) and U_{s} = {| \hat{y} | < s, {\hat{y}}_{2} < 0} \end{matrix}$ as the negative semi-disk of radius s and we denote by $V_{s}$ the set $V_{s} = Φ (U_{s})$ . Let r and δ be positive scalars such that $\begin{matrix} Ψ ({\underline{x}}_{0}) \in U_{r} and U_{r + 3 δ} \subset U_{Ω} . \end{matrix}$ We denote the space ${\tilde{L}}^{1} (U_{s}) \subset L^{1} (U_{Ω})$ for $s ⩽ r + 3 δ$ as the space of functions of $L^{1} (U_{s})$ that are extended by 0 to $U_{Ω}$ . We define the smooth indicator function ${\hat{χ}}_{s} \in C^{1, 1} (U_{Ω})$ as $\begin{matrix} {\hat{χ}}_{s} = \{\begin{matrix} 1 & in U_{s}, \\ {(\frac{| \hat{y} | - s}{δ} - 1)}^{2} (2 \frac{| \hat{y} | - s}{δ} + 1) & s < | \hat{y} | < s + δ, \\ 0 & otherwise, \end{matrix} \end{matrix}$ and we define $χ_{s} \in C^{1, 1} (V_{Ω})$ by $χ_{s} = {\hat{χ}}_{s} \circ Ψ$ . The function $χ \equiv χ_{r}$ is used in the localisation process that is explained below (Step 1 of the proof).

Preliminary: Difference quotients. Now assume that $s ⩽ r + 2 δ$ . Following the standard strategy (see e.g. [16]) we define, for $| h | < δ$ , the translation operator $τ_{h}$ , $\begin{matrix} τ_{h} : & {\tilde{L}}^{1} (U_{s}) ⟶ {\tilde{L}}^{1} (U_{s + δ}), \\ w (\hat{\underline{y}}) ⟶ τ_{h} w (\hat{\underline{y}}) = \{\begin{matrix} w (\hat{\underline{y}} + h {\hat{\underline{e}}}_{1}) & \hat{\underline{y}} + h {\hat{\underline{e}}}_{1} \in U_{s}, \\ 0 & otherwise . \end{matrix} \end{matrix}$ One can easily see that $τ_{- h} τ_{h}$ is the identity operator in ${\tilde{L}}^{1} (U_{s})$ for $s ⩽ r + δ$ . Moreover, for all w and v in $L^{2} (U_{Ω}) \cap {\tilde{L}}^{1} (U_{s})$ we have the property $\begin{matrix} (A.5) & {(w, τ_{h} v)}_{L^{2} (U_{Ω})} = {(τ_{- h} w, v)}_{L^{2} (U_{Ω})} . \end{matrix}$ The operator $τ_{h}$ commutes with derivatives and for all $w \in H_{0}^{1} (U_{Ω}) \cap {\tilde{L}}^{1} (U_{s})$ we have $τ_{h} w \in H_{0}^{1} (U_{Ω}) \cap {\tilde{L}}^{1} (U_{s + δ})$ . For $2 ⩽ p ⩽ \infty$ , we have $\begin{matrix} (A.6) & {‖ \frac{τ_{h} - 1}{h} w ‖}_{L^{p} (U_{Ω})} ⩽ ‖ w ‖_{W^{1, p} (U_{Ω})}, \forall w \in W^{1, p} (U_{Ω}) \cap {\tilde{L}}^{1} (U_{s}) . \end{matrix}$ Finally, if $w \in L^{2} (U_{Ω}) \cap {\tilde{L}}^{1} (U_{s})$ and, for h sufficiently small, $‖ (τ_{h} - 1) w 3 ‖_{L^{2} (U_{Ω})} ⩽ C h$ for some $C > 0$ independent of h, then $\begin{matrix} (A.7) & \frac{\partial w}{\partial {\hat{y}}_{1}} \in L^{2} (U_{Ω}) and {‖ \frac{\partial w}{\partial {\hat{y}}_{1}} ‖}_{L^{2} (U_{Ω})} ⩽ C \end{matrix}$ with the same constant C. The definition of $τ_{h}$ is naturally extended to vector fields by applying the operator component-wise.

In what follows we will use the notation ≲ to compare quantities up to a constant that depends only on the domain Ω as well as on the elasticity parameters, the neighbourhood V, the radius r and the increment δ, but it is independent of h and $\underline{r}$ (the source term of (A.2)). We mimic [13] and separate the proof in several steps that we briefly describe below.

Step 1: Localisation. By algebraic manipulations we derive the equations for $(χ \underline{v}, χ q)$ , where χ is the smooth indicator function supported in some neighbourhood of $V \equiv V ({\underline{x}}_{0})$ for ${\underline{x}}_{0} \in \partial Ω$ as described in the preliminary considerations.

Step 2: Local maps. Using the mapping Φ, the equations obtained in Step 1 are written in the domain $U_{Ω}$ by a simple change of variables. The transformed unknowns ${\hat{\underline{v}}}_{χ}$ and ${\hat{q}}_{χ}$ satisfy then a local problem of mixed type. Note that this is different from the results provided in reference works on elliptic problems (e.g. [13]) where the local problems are elliptic.

Step 3: Local estimate of the Lagrange multiplier difference quotients. Using estimation techniques for mixed problem we deduce the following estimate of the Lagrange multiplier ${\hat{q}}_{χ}$ : $\begin{matrix} (A.8) & {‖ \frac{τ_{h} - 1}{h} {\hat{q}}_{χ} ‖}_{L^{2} (U_{Ω})} ≲ ‖ \underline{r} ‖_{H} + {‖ \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} ‖}_{H^{1} {(U_{Ω})}^{2}} . \end{matrix}$

Step 4: Local estimate of the tangential derivatives. Continuing the analysis of the local problem we deduce a local estimate of ${\hat{\underline{v}}}_{χ}$ from which we deduce a local estimate of ${\hat{q}}_{χ}$ : $\begin{matrix} (A.9) & {‖ \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} ‖}_{H_{0}^{1} {(U_{Ω})}^{2}} ≲ ‖ \underline{r} ‖_{H} \underset{(A.8)}{\Rightarrow} {‖ \frac{τ_{h} - 1}{h} {\hat{q}}_{χ} ‖}_{L^{2} (U_{Ω})} ≲ ‖ \underline{r} ‖_{H} . \end{matrix}$ Letting h go to 0 we use the properties of difference quotients to show that tangential derivatives of $({\hat{\underline{v}}}_{χ}, {\hat{q}}_{χ})$ (i.e. derivatives with respect to ${\hat{y}}_{1}$ ) belong to $L^{2} (U_{Ω})$ and are bounded by $‖ \underline{r} ‖_{H}$ .

Step 5: Local estimate of the normal derivatives. Interpreting the local mixed problem in the distributional sense, after algebraic manipulation involving the strongly elliptic condition, we derive that the couple $({\hat{\underline{v}}}_{χ}, {\hat{q}}_{χ})$ has normal derivatives in $L^{2} (U_{Ω})$ and $\begin{matrix} (A.10) & ‖ {\hat{\underline{v}}}_{χ} ‖_{H^{2} {(U_{Ω})}^{2}} + ‖ \underline{\nabla} {\hat{q}}_{χ} ‖_{L^{2} {(U_{Ω})}^{2}} ≲ ‖ \underline{r} ‖_{H} . \end{matrix}$

Step 6: Global estimate. First, thanks to the smoothness of the mapping Φ, we deduce an adequate local estimate on $V_{Ω}$ from (A.10). Then, we argue that the presented approach can be extended to any point ${\underline{x}}_{0}$ in the closure of Ω – and corresponding neighbourhood $V_{Ω} ({\underline{x}}_{0})$ – and that there exists a finite number of neighbourhoods $V_{Ω} ({\underline{x}}_{0})$ that constitutes a cover of Ω. Therefore, we can show that $\underline{v} \in H^{2} {(Ω)}^{2}$ and $q \in M$ as well as $\begin{matrix} (A.11) & ‖ \underline{v} ‖_{H^{2} {(Ω)}^{2}} + ‖ \underline{\nabla} q ‖_{L^{2} {(Ω)}^{2}} ≲ ‖ \underline{r} ‖_{H}, \end{matrix}$ thus concluding the proof.

We now give more details on each step of the proof.

Step 1: Localisation. Using the variational formulation (A.2), by direct algebraic computations one can show that $\begin{matrix} (A.12) & \{\begin{matrix} a (χ \underline{v}, \underline{w}) - {(χ q, div \underline{w})}_{L} = {({\underline{f}}_{χ, 0}, \underline{w})}_{H} - \sum_{k = 1}^{2} {({\underline{f}}_{χ, k}, \underline{\nabla} w_{k})}_{H} & \forall \underline{w} \in X, \\ {(div χ \underline{v}, z)}_{L} = {(g_{χ}, z)}_{L} & \forall z \in L, \end{matrix} \end{matrix}$ where $g_{χ} = \underline{\nabla} χ \cdot \underline{v} \in M$ and the ${\underline{f}}_{χ, k} \in H$ for $k \in {0, 1, 2}$ are given by $\begin{matrix} {\underline{f}}_{χ, 0} = χ \underline{r} + q \underline{\nabla} χ + \sum_{k = 1}^{2} \sum_{ℓ = 1}^{2} ({\underline{\underline{A}}}_{k ℓ} \underline{\nabla} v_{ℓ} \cdot \underline{\nabla} χ) {\underline{e}}_{k}, {\underline{f}}_{χ, k} = \sum_{ℓ = 1}^{2} v_{ℓ} {\underline{\underline{A}}}_{k ℓ} \underline{\nabla} χ, k \in {1, 2}, \end{matrix}$ where ${\underline{e}}_{k}$ is the k-th element of the canonical basis in $R^{2}$ . Observe that, because of the inclusion ${A_{k ℓ}} \subset W^{1, \infty} {(Ω)}^{2 \times 2}$ , the ${\underline{f}}_{χ, k}$ for $k \in {1, 2}$ belong to $H^{1} {(Ω)}^{2}$ and therefore we can simplify (A.12) by writing $\begin{matrix} {({\underline{f}}_{χ, 0}, \underline{w})}_{H} - \sum_{k = 1}^{2} {({\underline{f}}_{χ, k}, \underline{\nabla} w_{k})}_{H} = {({\underline{f}}_{χ, 0}, \underline{w})}_{H} + \sum_{k = 1}^{2} {(div {\underline{f}}_{χ, k}, w_{k})}_{H} = {({\underline{f}}_{χ}, \underline{w})}_{H} \end{matrix}$ with ${\underline{f}}_{χ} \in H$ and $‖ {\underline{f}}_{χ} ‖_{H} ≲ ‖ \underline{r} ‖_{H}$ . Finally, observe that the support of $g_{χ}$ and ${\underline{f}}_{χ}$ is included in $V_{r + δ}$ .

Step 2: Local maps. We denote by $X_{V}$ (respectively $L_{V}$ ) the subspace of functions of $X$ (respectively $L$ ) that vanish outside the domain $V_{Ω}$ . We choose $\underline{w} \in X_{V}$ and $z \in L_{V}$ and define $\begin{matrix} (A.13) & \hat{\underline{w}} (\cdot) = \underline{w} ({\underline{\underline{R}}}^{t} Φ (\cdot) + {\underline{x}}_{R}) \in H_{0}^{1} {(U_{Ω})}^{2}, \hat{z} (\cdot) = z ({\underline{\underline{R}}}^{t} Φ ((\cdot)) + {\underline{x}}_{R}) \in L_{0}^{2} (U_{Ω}), \end{matrix}$ and we use the notation $({\hat{\underline{v}}}_{χ}, {\hat{q}}_{χ}, {\underline{\hat{f}}}_{χ}, {\hat{g}}_{χ})$ for the functions – compactly supported in $U_{r + δ}$ – obtained using the same transformation as above, applied to $(χ \underline{v}, χ q, {\underline{f}}_{χ}, g_{χ})$ respectively. Then, using standard calculus it is possible to prove that for all $\underline{w} \in X_{V}$ and $z \in L_{V}$ $\begin{matrix} (A.14) & \{\begin{matrix} \hat{a} ({\hat{\underline{v}}}_{χ}, \hat{\underline{w}}) - \hat{b} ({\hat{q}}_{χ}, \hat{\underline{w}}) = {({\underline{\hat{f}}}_{χ}, \hat{\underline{w}})}_{L^{2} {(U_{Ω})}^{2}}, \\ \hat{b} (\hat{z}, {\hat{\underline{v}}}_{χ}) = {({\hat{g}}_{χ}, \hat{z})}_{L^{2} (U_{Ω})}, \end{matrix} \end{matrix}$ with $\begin{matrix} \hat{a} (\hat{\underline{v}}, \hat{\underline{w}}) = \sum_{k = 1}^{2} \sum_{ℓ = 1}^{2} \int_{U_{Ω}} {\underline{\underline{\hat{A}}}}_{k ℓ} \underline{\nabla} {\hat{v}}_{ℓ} \cdot \underline{\nabla} {\hat{w}}_{k} d \hat{\underline{y}}, \hat{b} (\hat{z}, \hat{\underline{w}}) = {(\hat{z}, div {\underline{\underline{F}}}^{- 1} \underline{\underline{R}} \hat{\underline{w}})}_{L^{2} (U_{Ω})} \end{matrix}$ and $\begin{matrix} {\underline{\underline{\hat{A}}}}_{k ℓ} (\hat{\underline{y}}) = {\underline{\underline{F}}}^{- 1} (\hat{\underline{y}}) \underline{\underline{R}} (\hat{\underline{y}}) {\underline{\underline{A}}}_{k ℓ} ({\underline{\underline{R}}}^{t} Φ (\hat{\underline{y}}) + {\underline{x}}_{R}) {\underline{\underline{R}}}^{t} (\hat{\underline{y}}) {\underline{\underline{F}}}^{- t} (\hat{\underline{y}}) . \end{matrix}$ Since for any $\hat{\underline{w}} \in H_{0}^{1} {(U_{Ω})}^{2}$ (respectively $\hat{z} \in L_{0}^{2} (U_{Ω})$ ) there exists $\underline{w} \in X_{V}$ (respectively $z \in L_{V}$ ) such that (A.13) holds, Eq. (A.14) is true for all test functions $(\hat{\underline{w}}, \hat{z}) \in H_{0}^{1} {(U_{Ω})}^{2} \times L_{0}^{2} (U_{Ω})$ .

Step 3: Local estimate of the Lagrange multiplier difference quotients. We choose as a test function in (A.14) $\begin{matrix} \hat{\underline{w}} = \frac{τ_{- h} - 1}{h} \tilde{\underline{w}}, \end{matrix}$ where $\tilde{\underline{w}} \in H_{0}^{1} {(U_{Ω})}^{2} \cap {\tilde{L}}^{1} (U_{r + 2 δ})$ . Using the Eq. (A.5) and by classic algebraic manipulation (see also the similar computations in [16, Lemma 2.2.2.1]), it is possible to show that $\begin{matrix} (A.15) & \hat{a} ({\hat{\underline{v}}}_{χ}, \frac{τ_{- h} - 1}{h} \tilde{\underline{w}}) = \hat{a} (\frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ}, \tilde{\underline{w}}) + \tilde{a} ({\hat{\underline{v}}}_{χ}, \tilde{\underline{w}}), \end{matrix}$ with $\begin{matrix} \tilde{a} ({\hat{\underline{v}}}_{χ}, \tilde{\underline{w}}) = \sum_{k = 1}^{2} \sum_{ℓ = 1}^{2} \int_{U_{Ω}} (\frac{τ_{h} - 1}{h} ({\underline{\underline{\hat{A}}}}_{k ℓ} {\hat{χ}}_{2 r})) \underline{\nabla} {\hat{v}}_{χ, ℓ} \cdot \underline{\nabla} {\tilde{w}}_{k} d \hat{\underline{y}} . \end{matrix}$ Using property (A.3), the fact that the ${\underline{\underline{\hat{A}}}}_{k ℓ}$ have uniformly bounded derivatives and Eq. (A.6), one can show that $\begin{matrix} | \tilde{a} ({\hat{\underline{v}}}_{χ}, \tilde{\underline{w}}) | ≲ ‖ \underline{r} ‖_{H} ‖ \underline{\underline{\nabla}} \tilde{\underline{w}} ‖_{L^{2} {(U_{Ω})}^{2 \times 2}} . \end{matrix}$ Similarly we have $\begin{matrix} (A.16) & b ({\hat{q}}_{χ}, \frac{τ_{- h} - 1}{h} \tilde{\underline{w}}) = b (\frac{τ_{h} - 1}{h} {\hat{q}}_{χ}, \tilde{\underline{w}}) + \tilde{b} ({\hat{q}}_{χ}, \tilde{\underline{w}}), \end{matrix}$ with $\begin{matrix} \tilde{b} ({\hat{q}}_{χ}, \tilde{\underline{w}}) = \int_{U_{Ω}} {\hat{q}}_{χ} div [(\frac{τ_{- h} - 1}{h} ({\underline{\underline{F}}}^{- 1} \underline{\underline{R}} {\hat{χ}}_{2 r})) \tilde{\underline{w}}] d \hat{\underline{y}} . \end{matrix}$ Thanks to Piola’s identity, $\underline{div} {\underline{\underline{R}}}^{t} {\underline{\underline{F}}}^{- t} = \underline{0}$ (see [14, Section 8.1]), one can check that $\begin{matrix} | \tilde{b} ({\hat{q}}_{χ}, \tilde{\underline{w}}) | = | \int_{U_{Ω}} {\hat{q}}_{χ} (\frac{τ_{- h} - 1}{h} ({\underline{\underline{F}}}^{- 1} \underline{\underline{R}} {\hat{χ}}_{2 r})) : \underline{\underline{\nabla}} \tilde{\underline{w}} d \hat{\underline{y}} | ≲ ‖ r ‖_{H} ‖ \underline{\underline{\nabla}} \tilde{\underline{w}} ‖_{L^{2} {(U_{Ω})}^{2 \times 2}} . \end{matrix}$ Using (A.14), (A.15) and (A.16) we find $\begin{array}{l} b (\frac{τ_{h} - 1}{h} {\hat{q}}_{χ}, \tilde{\underline{w}}) \\ (A.17) & = - {({\underline{\hat{f}}}_{χ}, \frac{τ_{- h} - 1}{h} \tilde{\underline{w}})}_{L^{2} {(U_{Ω})}^{2}} - \tilde{b} ({\hat{q}}_{χ}, \tilde{\underline{w}}) + \hat{a} (\frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ}, \tilde{\underline{w}}) + \tilde{a} ({\hat{\underline{v}}}_{χ}, \tilde{\underline{w}}) . \end{array}$ Observe that $(τ_{h} - 1) {\hat{q}}_{χ} \in L_{0}^{2} (U_{r + 2 δ})$ . Hence, from [15, Corollary 2.4] there exists $\tilde{\underline{v}} \in H_{0}^{1} {(U_{r + 2 δ})}^{2}$ such that $\begin{matrix} (A.18) & div \tilde{\underline{v}} = \frac{τ_{h} - 1}{h} {\hat{q}}_{χ}, ‖ \tilde{\underline{v}} ‖_{H^{1} {(U_{Ω})}^{2}} ≲ {‖ \frac{τ_{h} - 1}{h} {\hat{q}}_{χ} ‖}_{L^{2} (U_{Ω})} . \end{matrix}$ We choose as a test function in (A.17) the function $\tilde{\underline{w}} = \underline{\underline{F}} \underline{\underline{R}} \tilde{\underline{v}}$ in $U_{r + 2 δ}$ that is extended by 0 over $U_{Ω}$ , so we have $\tilde{\underline{w}} \in H_{0}^{1} {(U_{Ω})}^{2} \cap {\tilde{L}}^{1} (U_{r + 2 δ})$ . Since $\underline{\underline{F}}$ has bounded derivatives, using (A.18) we have $\begin{matrix} (A.19) & ‖ \tilde{\underline{w}} ‖_{H^{1} {(U_{Ω})}^{2}} ≲ ‖ \tilde{\underline{v}} ‖_{H^{1} {(U_{Ω})}^{2}} ≲ {‖ \frac{τ_{h} - 1}{h} {\hat{q}}_{χ} ‖}_{L^{2} (U_{Ω})} . \end{matrix}$ With this choice of test function we deduce from (A.17) $\begin{array}{rcl} b (\frac{τ_{h} - 1}{h} {\hat{q}}_{χ}, \tilde{\underline{w}}) = {‖ \frac{τ_{h} - 1}{h} {\hat{q}}_{χ} ‖}_{L^{2} (U_{Ω})}^{2} & ≲ & ‖ {\underline{\hat{f}}}_{χ} ‖_{L^{2} {(U_{Ω})}^{2}} {‖ \frac{τ_{- h} - 1}{h} \tilde{\underline{w}} ‖}_{L^{2} {(U_{Ω})}^{2}} \\ + (‖ r ‖_{H} + {‖ \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} ‖}_{H^{1} {(U_{Ω})}^{2}}) ‖ \tilde{\underline{w}} ‖_{H^{1} {(U_{Ω})}^{2}} . \end{array}$ Using (A.6), the expression of ${\underline{f}}_{χ}$ and (A.19), we obtain the estimate (A.8).

Step 4: Local estimate of the tangential derivatives. We now choose $\begin{matrix} \hat{z} = \frac{τ_{- h} - 1}{h} \frac{τ_{h} - 1}{h} {\hat{q}}_{χ} \in L^{2} (U_{Ω}) \cap {\tilde{L}}^{1} (U_{r + 2 δ}) \end{matrix}$ as a test function in the second equation of (A.14). Thanks to (A.16) we obtain $\begin{matrix} (A.20) & \hat{b} (\frac{τ_{h} - 1}{h} {\hat{q}}_{χ}, \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ}) + \tilde{b} ({\hat{q}}_{χ}, \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ}) = {(\frac{τ_{h} - 1}{h} {\hat{g}}_{χ}, \frac{τ_{h} - 1}{h} {\hat{q}}_{χ})}_{L^{2} (U_{Ω})} . \end{matrix}$ Moreover, it is possible to choose $\begin{matrix} \tilde{\underline{w}} = \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} \in L^{2} (U_{Ω}) \cap {\tilde{L}}^{1} (U_{3 r}) \end{matrix}$ as a test function in (A.17). Using (A.20), we get $\begin{array}{l} \hat{a} (\frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ}, \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ}) \\ = {({\underline{\hat{f}}}_{χ}, \frac{τ_{- h} - 1}{h} \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ})}_{L^{2} {(U_{Ω})}^{2}} + {(\frac{τ_{h} - 1}{h} {\hat{g}}_{χ}, \frac{τ_{h} - 1}{h} {\hat{q}}_{χ})}_{L^{2} (U_{Ω})} - \tilde{a} ({\hat{\underline{v}}}_{χ}, \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ}) . \end{array}$ Thanks to Proposition A.1, Gårding’s inequality can be used. We obtain the estimate $\begin{array}{rcl} {‖ \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} ‖}_{H_{0}^{1} {(U_{Ω})}^{2}}^{2} & ≲ & ‖ {\underline{\hat{f}}}_{χ} ‖_{L^{2} {(U_{Ω})}^{2}} {‖ \frac{τ_{- h} - 1}{h} \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} ‖}_{L^{2} {(U_{Ω})}^{2}} \\ + ‖ \frac{τ_{h} - 1}{h} {\hat{g}}_{χ} ‖ (‖ \underline{r} ‖_{H} + {‖ \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} ‖}_{H^{1} {(U_{Ω})}^{2}}) \\ + ‖ {\hat{\underline{v}}}_{χ} ‖_{H_{0}^{1} {(U_{Ω})}^{2}} {‖ \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} ‖}_{H_{0}^{1} {(U_{Ω})}^{2}} + {‖ \frac{τ_{h} - 1}{h} {\hat{\underline{v}}}_{χ} ‖}_{L^{2} {(U_{Ω})}^{2}}^{2} . \end{array}$ Using (A.3) and (A.6), the estimation of ${\underline{\hat{f}}}_{χ}$ , the property that ${\hat{g}}_{χ}$ belongs to $H^{1} (U_{Ω})$ and Young’s inequality, one can show that (A.9) holds. This implies, at the limit as h goes to 0 and thanks to (A.7), $\begin{matrix} \frac{\partial {\hat{\underline{v}}}_{χ}}{\partial {\hat{y}}_{1}} \in H_{0}^{1} {(U_{Ω})}^{2}, \frac{\partial {\hat{q}}_{χ}}{\partial {\hat{y}}_{1}} \in L^{2} (U_{Ω}) and {‖ \frac{\partial {\hat{\underline{v}}}_{χ}}{\partial {\hat{y}}_{1}} ‖}_{H_{0}^{1} {(U_{Ω})}^{2}} + {‖ \frac{\partial {\hat{q}}_{χ}}{\partial {\hat{y}}_{1}} ‖}_{L^{2} (U_{Ω})} ≲ ‖ \underline{r} ‖_{H} . \end{matrix}$

Step 5: Local Estimate of the normal derivatives. We define the matrix field $\underline{\underline{\hat{A}}}$ and the vector field $\underline{\hat{b}}$ as follows: $\begin{matrix} \underline{\underline{\hat{A}}} {\underline{e}}_{k} \cdot {\underline{e}}_{ℓ} = {\underline{\underline{\hat{A}}}}_{k ℓ} {\underline{e}}_{2} \cdot {\underline{e}}_{2}, \underline{\hat{b}} \cdot {\underline{e}}_{ℓ} = {\underline{\underline{R}}}^{t} {\underline{\underline{F}}}^{- t} {\underline{e}}_{2} \cdot {\underline{e}}_{ℓ} . \end{matrix}$ Interpreting (A.14) in the distributional sense one can show that $\begin{matrix} (\begin{matrix} \underline{\underline{\hat{A}}} & \underline{\hat{b}} \\ {\underline{\hat{b}}}^{t} & 0 \end{matrix}) \frac{\partial}{\partial {\hat{y}}_{2}} (\begin{matrix} \frac{\partial {\hat{\underline{v}}}_{χ}}{\partial {\hat{y}}_{2}} \\ {\hat{q}}_{χ} \end{matrix}) = (\begin{matrix} {\underline{\tilde{f}}}_{χ} \\ {\tilde{g}}_{χ} \end{matrix}) with ‖ {\underline{\tilde{f}}}_{χ} ‖_{L^{2} {(U_{Ω})}^{2}} + ‖ {\tilde{g}}_{χ} ‖_{L^{2} (U_{Ω})} ≲ ‖ \underline{r} ‖_{H} . \end{matrix}$ Since the strong elliptic condition is satisfied, $\underline{\underline{\hat{A}}}$ is invertible (with an inverse uniformly bounded in $\hat{\underline{y}}$ ). Consequently, one can see by a Schur complement technique that we have $\begin{matrix} (A.21) & {\underline{\hat{b}}}^{t} {\underline{\underline{\hat{A}}}}^{- 1} \underline{\hat{b}} \frac{\partial {\hat{q}}_{χ}}{\partial {\hat{y}}_{2}} = {\underline{\hat{b}}}^{t} {\underline{\underline{\hat{A}}}}^{- 1} {\underline{\tilde{f}}}_{χ} - {\tilde{g}}_{χ} . \end{matrix}$ Moreover, $1 ≲ | \underline{\hat{b}} |^{2} ≲ | {\underline{\hat{b}}}^{t} {\underline{\underline{\hat{A}}}}^{- 1} \underline{\hat{b}} |$ almost everywhere in $U_{Ω}$ . Therefore one can show, together with the results of Step 4, that $q_{χ}$ belongs to $H^{1} (U_{Ω})$ and $\begin{matrix} ‖ \underline{\nabla} {\hat{q}}_{χ} ‖_{L^{2} {(U_{Ω})}^{2}} ≲ ‖ \underline{r} ‖_{H} . \end{matrix}$ Using again (A.21), one can show that ${\hat{\underline{v}}}_{χ} \in H^{2} {(U_{Ω})}^{2}$ and (A.10) can be deduced.

Step 6: Global estimate. We have proved that $\begin{matrix} ‖ \underline{v} \circ Φ ‖_{H^{2} {(U_{r})}^{2}} + {‖ \underline{\nabla} (q \circ Φ) ‖}_{L^{2} {(U_{r})}^{2}} ≲ ‖ \underline{r} ‖_{H} . \end{matrix}$ This implies, since $Φ \in C^{1, 1} (U)$ (in particular it has bounded second-order derivatives), $\begin{matrix} (A.22) & ‖ \underline{v} ‖_{H^{2} (V_{r})} + ‖ \underline{\nabla} q ‖_{L^{2} (V_{r})} ≲ ‖ \underline{r} ‖_{H} . \end{matrix}$ Repeating Steps 1–5 for any ${\underline{x}}_{0} \in Ω$ , it is possible to show that (A.22) holds for a neighbourhood $V ({\underline{x}}_{0})$ of any ${\underline{x}}_{0}$ in the closure of Ω. Since the domain Ω is bounded in $R^{2}$ , there exists a finite cover of Ω such that (A.22) holds and therefore one can deduce the global estimate (A.11).

References

R.A.

Adams and

J.J.

Fournier, Sobolev Spaces, Vol. 140, Elsevier, 2003.

Agmon,

Douglis and

Nirenberg, Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions II, Communications on Pure and Applied Mathematics 17(1) (1964), 35–92. doi:10.1002/cpa.3160170104.

Amrouche and

Girault, Decomposition of vector spaces and application to the Stokes problem in arbitrary dimension, Czechoslovak Mathematical Journal 44(1) (1994), 109–140.

Boffi,

Brezzi,

Fortin et al., Mixed Finite Element Methods and Applications, Vol. 44, Springer, 2013.

S.C.

Brenner and

L.-Y.

Sung, Linear finite element methods for planar linear elasticity, Mathematics of Computation 59(200) (1992), 321–338. doi:10.1090/S0025-5718-1992-1140646-2.

Caforio, Mathematical modelling and numerical simulation of elastic wave propagation in soft tissues with application to cardiac elastography, Ph.D. thesis, Université Paris-Saclay, 2019.

Caforio and

Imperiale, A conservative penalisation strategy for the semi-implicit time discretisation of the incompressible elastodynamics equation, Advanced Modeling and Simulation in Engineering Sciences 5(1) (2018), 30. doi:10.1186/s40323-018-0121-8.

Caforio and

Imperiale, A high-order spectral element fast Fourier transform for the Poisson equation, SIAM Journal on Scientific Computing 41(5) (2019), A2747–A2771.

A.J.

Chorin, Numerical solution of the Navier–Stokes equations, Mathematics of Computation 22(104) (1968), 745–762. doi:10.1090/S0025-5718-1968-0242392-2.

10.

A.J.

Chorin, On the convergence of discrete approximations to the Navier–Stokes equations, Mathematics of Computation 23(106) (1969), 341–353. doi:10.1090/S0025-5718-1969-0242393-5.

11.

P.G.

Ciarlet, Mathematical Elasticity. Vol. I, Studies in Mathematics and Its Applications, Vol. 20, 1988.

12.

Cohen and

Fauqueux, Mixed spectral finite elements for the linear elasticity system in unbounded domains, SIAM Journal on Scientific Computing 26(3) (2005), 864–884. doi:10.1137/S1064827502407457.

13.

Costabel and

Dauge, Corner singularities and analytic regularity for linear elliptic systems. Part I: Smooth domains, hal-00453934v2, 2010.

14.

L.C.

Evans, Partial Differential Equations, American Mathematical Society, 2010.

15.

Girault and

P.-A.

Raviart, Finite Element Approximation of the Navier–Stokes Equations, Springer-Verlag, 1979, Tech. rep.

16.

Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman Advanced Publishing Programm, 1985.

17.

J.-L.

Guermond, Un résultat de convergence d’ordre deux en temps pour l’approximation des équations de Navier–Stokes par une technique de projection incrémentale, ESAIM: Mathematical Modelling and Numerical Analysis 33(1) (1999), 169–189. doi:10.1051/m2an:1999101.

18.

J.-L.

Guermond,

Minev and

Shen, An overview of projection methods for incompressible flows, Computer Methods in Applied Mechanics and Engineering 195(44–47) (2006), 6011–6045. doi:10.1016/j.cma.2005.10.010.

19.

J.-L.

Guermond and

Quartapelle, On the approximation of the unsteady Navier–Stokes equations by finite element projection methods, Numerische Mathematik 80(2) (1998), 207–238. doi:10.1007/s002110050366.

20.

J.-L.

Guermond and

Quartapelle, On stability and convergence of projection methods based on pressure Poisson equation, International Journal for Numerical Methods in Fluids 26(9) (1998), 1039–1053. doi:10.1002/(SICI)1097-0363(19980515)26:9<1039::AID-FLD675>3.0.CO;2-U.

21.

Joly, The mathematical model for elastic wave propagation, in: Effective Computational Methods for Wave Propagation, Numerical Insights, CRC Press, 2008.

22.

Kellogg and

Osborn, A regularity result for the Stokes problem in a convex polygon, Journal of Functional Analysis 21(4) (1976), 397–431. doi:10.1016/0022-1236(76)90035-5.

23.

J.E.

Marsden and

Hughes, Mathematical Foundations of Elasticity, Dover Publications, 1994.

24.

Monk, Finite Element Methods for Maxwell’s Equations, Oxford University Press, 2003.

25.

Temam, Une méthode d’approximation de la solution des équations de Navier–Stokes, Bull. Soc. Math. France 98(4) (1968), 115–152. doi:10.24033/bsmf.1662.

26.

Temam, Navier–Stokes Equations: Theory and Numerical Analysis, Vol. 343, American Mathematical Soc., 2001.

Mathematical analysis of a penalisation strategy for incompressible elastodynamics

Abstract

Keywords

1. Introduction

2. Setting of the problem

2.1. Preliminaries and notations

2.2. The quasi-incompressible elastodynamic equations

Proposition 2.4. Let f _ ∈ W 0 k ( H ) with k ⩾ 1 . Then, there exists a unique couple ( y _ , p ) ∈ [ C k + 1 ( H ) ∩ C k ( V ) ] × C k − 1 ( L ) solution to ( IM ). 2.4. A penalised quasi-incompressible formulation

4.1. Definition of an elasto-static operator

5.1. Isotropic elasticity

Footnotes

Proof of Theorem 5.3

References

Proposition 2.4.
Let $\underline{f} \in W_{0}^{k} (H)$ with $k ⩾ 1$ . Then, there exists a unique couple $(\underline{y}, p) \in [C^{k + 1} (H) \cap C^{k} (V)] \times C^{k - 1} (L)$ solution to ( IM ).

2.4. A penalised quasi-incompressible formulation