A conservative penalisation strategy for the semi-implicit time discretisation of the incompressible elastodynamics equation

The principal aim of this work is to provide an adapted numerical scheme for the approximation of elastic wave propagation in incompressible solids. We rely on high-order conforming finite element with mass lumping for space discretisation and implicit/explicit, second-order, energy-preserving time discretisation. The time step restriction only depends on the shear wave velocity and at each time step a Poisson problem must be solved to account for the incompressibility constraint that is imposed by penalisation techniques.


Introduction
Since large applications in computational mechanics concern nearly and pure incompressible elasticity (living tissues, biomaterials), a great effort has been made in the last 40 years to provide accurate finite element method (FEM) formulations for the approximation of elasticity in incompressible solids. However, the majority of the works proposed in the literature only deals with static computations. The main contributions to date can be mainly divided into two categories: pure displacement methods and mixed methods.
On the one hand, displacement-based FEM can provide accurate solutions of quasi or pure incompressible elasticity problems; nonetheless, the space resolution necessary to provide an accurate solution is far greater than the one required for a compressible material [1,2]. Indeed, these methods can suffer of undesirable limitations, such as ill-conditioning of the stiffness matrix, spurious or incorrect pressures and numerical locking (severe stiffening near the incompressible limit) [2], especially if low-order shape functions are adopted, due to the enforcement of the incompressibility constraint-i.e. the requirement that the displacement field is divergence-free. Locking is due to the fact that, in case of incompressible materials, volumetric strains approach zero, while the pressure field is of the order of the boundary traction, therefore it cannot be computed from strains, but it must be calculated directly from the equilibrium equations [2,3]. Several methods have been proposed to improve accuracy of displacement-based methods. Among them, we cite the reduced/selective integration method proposed in [4], the B-bar method [5] and the F-bar method [6]. These methods circumvent volumetric locking by reducing the number of discrete incompressibility constraints enforced at the quadrature points [7]. We also cite a method adapted to nearly-incompressible materials with volumetric energy penalty function in the framework of the continuum-based absolute nodal coordinate formulation (ANCF) [8].
On the other hand, mixed finite element methods [9] have proven effective and even necessary to obtain accurate results in the resolution of incompressible fluid flows and incompressible elasticity. In these methods, the constrained problem is rewritten in form of an unconstrained saddle-point problem, due to the introduction of a second variable (namely, the pressure). However, not all mixed methods are stable. In fact, the convergence properties of this formulation are governed by stability considerations, involving ellipticity requirements and the famous Ladyzenskaya-Babǔska-Brezzi (LBB) inf-sup condition [10,11]. For example, equal-order interpolation both for displacement and pressure field does not satisfy the LBB condition for classical mixed FEM [9]. If this stability condition is not satisfied, severe unphysical oscillations in the pressure field can appear, named "spurious pressure modes". Stabilised methods have been proposed to overcome the limitations of classical mixed FE formulations in the field of incompressible fluid dynamics (see [12][13][14] and references therein for throughout reviews on the subject). Brezzi and collaborators [15] proposed to extend the equation accounting for incompressibility in Stokes flows by adding a laplacian of the pressure field. Other methods are based on the addition of artificial high-order differential terms to the discrete continuity equation, in order to let the formulation satisfy the LBB condition (for more details, see [16]). We cite for example the Streamline Upwind Petrov-Galerkin method [17][18][19]. Similar stabilised methods have been recently extended to the context of linear elasticity [20][21][22][23][24].
All the aforementioned methods can be straightforwardly extended to the discretisation of dynamic equations using implicit time discretisation (e.g. Implicit Euler scheme or implicit Newmark schemes). However, at each time step, the resolution of the resulting linear system is required for the computation of the velocity field (or most often the displacement field in elasticity). We highlight that, in linear elasticity, this could be done in practice by performing a factorisation of the matrices to invert, since they are constant in time. However, for large scale problems, it is not possible to store any factorisation of the matrices and it is even difficult to store preconditioners. A popular approach to increase the efficiency of dynamic solvers was first proposed in computational fluid dynamics in the late 1960s [25][26][27][28], and it is called fractional-step projection. This family of methods aims at accurately solving the equation governing viscous incompressible fluids by performing a time-discretisation in which viscosity and incompressibility are treated in two separate steps. In more detail, the first half-step (Burgers step) corresponds to an elliptic Boundary Value Problem (BVP) for an intermediate velocity, accounting for viscosity diffusion and advection. The second half-step (projection step) represents an inviscid problem where the end-of-step, divergence-free velocity is computed, along with pressure distribution. This step essentially consists in solving a Poisson problem. In this way, at each time step two decoupled elliptic equations are solved, and this is very advantageous for large scale simulations [29][30][31].
Less effort has been made to develop efficient methods for the treatment of the incompressibility constraint in elastodynamics. Since the underlying physics is wave propagation, fully explicit methods seem to be good candidates to obtain efficient schemes (this is possible for nearly-incompressible media). We cite in this context a recent method for quasi-incompressible elastodynamics [32] with linear finite elements, extended in [33] to the purely incompressible case. Furthermore, a recent velocity/stress formulation for the simulation of linear elastodynamics in nearly-incompressible solids with weakly enforced boundary conditions was proposed by Scovazzi and collaborators [34]. However, the stability constraint (CFL condition) is drastically limited by the enforcement of incompressibility. In this regard, the efficiency of the fractional-step projection algorithm mentioned above, along with the similarities with viscous incompressible fluids, suggests the possibility to adopt the main ideas of this method to design an efficient scheme for incompressible elastodynamics.
A first method that integrates a fractional time-step method for Lagrangian formulations of elastodynamics problems has been proposed by Lahiri and collaborators [35]. In this paper, the authors use variational integrators that take advantage of the Hamilton variational principle to construct a discrete approximation of the integral of the Lagrangian over a given interval, and they adopt linear finite element discretisation in space.
We have developed a numerical scheme that carefully takes into account the intrinsic properties of the wave equation. Namely, we construct a conservative time discretisation, and treat implicitly only the terms corresponding to "informations" travelling at infinite velocity (i.e. the incompressibility constraint) by solving a Poisson problem. The fully discrete scheme that we propose to efficiently solve the incompressible elastodynamic problem is where A h denotes a stiffness operator (independent of the compressibility parameters), B T h and B h correspond to discrete gradient and divergence operators, and C h C T h represents a discrete laplacian operator. The coefficient α corresponds to a penalisation parameter and (ỹ n α,h ,p n α,h , f n h ) denote respectively the displacement field, the pressure and the source term.
Therefore, if effective methods for explicit time-discretisation are used, our algorithm requires at each time step the resolution of a scalar Poisson problem (that can be done by several, efficient algorithms) and few matrix-vector multiplications for the explicit methods.
Note that the stability condition for Scheme (1) reads for some operator norm · that we define later. Contrarily to the standard results one could expect (i.e. Δt 2 4 A h −1 , see [36]), the stability condition is slightly more constraining due to the factor (4α − 1)/4α. Our approach is a strategy to approximate elastic wave propagation in quasiincompressible media. In order to do that, we introduce a first good approximation: the pure incompressible formulation. Then, we construct a penalised formulation to approximate the pure incompressible problem. Our procedure is justified by arguments of stability, computational cost and numerical convergence of the resulting schemes.
Finally, the work we propose is closely related to the Selective Mass Scaling Method (SMS) developed in [37] for wave propagation in quasi-incompressible materials. In this method, the mass matrix is modified: artificial mass is added, hence obtaining very good stability properties at the cost of computing a mass matrix that cannot be lumped anymore. The problem is then implicit for the displacement field, and thus it not possible to use fast solvers (like FFT solvers) for the inversion of a scalar Poisson problem.
The article is organised as follows. In "Continuous framework" section we provide two standard formulations of the continuous elastodynamic problem for quasi-incompressible and pure incompressible media along with a novel formulation for the treatment of incompressibility by penalisation. Furthermore, we derive the variational formulation associated with each problem. "Space discretisation" section deals with the abstract framework for space discretisation of the incompressible elastodynamics equation by Spectral Element Method. In "Time discretisation" section we provide the time discretisations for each formulation by finite difference. Then, a stability analysis based on energy considerations is performed for each scheme in "Stability analysis" section and pros and cons are discussed. Numerical results, including convergence curves to the solution of the incompressible elastodynamics equations for different choices of materials, are shown in "Two-dimensional numerical convergence results" section. Finally, a three-dimensional illustration in a more realistic test case for elastography imaging is shown in "A threedimensional test case" section.

The equation of elastodynamics
Given a domain Ω ⊂ R d smooth enough, with d = 2 or d = 3, we introduce the following notations to define Hilbert spaces for the elastic displacements For the sake of simplicity, we consider homogenous Dirichlet condition on the boundary of the propagation domain. We also need to consider divergence-free displacements. Hence, we introduce the following subspace of X Pressure is a variable of interest, and is sought in the spaces As usual, we identify L and H with their dual in what follows. Moreover, for the sake of conciseness, we define, given a function space A on Ω, where T > 0 is a given final time of observation. Furthermore, we introduce Ω T := [0, T ] × Ω. Our aim is to analyse the propagation of elastic waves in heterogenous, anisotropic, incompressible solids and we consider as a reference model the following partial differential equation (PDE) system: For f given and sufficiently regular, find y with λ ∈ R + the bulk modulus, that is assumed to be large, due to nearly-incompressibility, ρ(x) the strictly positive density of the medium and C(x) the elasticity tensor which is symmetric, coercive and bounded, i.e. there exist two positive scalars c, C such that

Non-dimensionalisation
Since we are going to consider a limit process where the bulk modulus tends to infinity, a first step is to non-dimensionalise our system of equations. To do so, we introduce a typical length scale L of the domain Ω, a typical observation time τ and a shear modulus μ, and we define a non-dimensionalised displacement as follows: where Ω is a rescaled domain and T = T /τ is of the order of unity. Note that t and x refer now to non-dimensionalised variables. We also introduce the non-dimensionalised quantitieŝ Then, the equation of elastodynamics can be recast as We also introduce here the pressure field p λ associated with the displacement field y λ . Its non-dimensionalised counterpart is given bŷ For the sake of simplicity, we drop the notation· throughout the rest of the paper.

The mixed and penalised formulations
Existence and uniqueness results for problem (QI) are well-known (see [36]). An alternative, equivalent formulation to (QI) is obtained by introducing the scalar function Since λ is large, it is natural to approximate the solution of (QIM) by the solution obtained at the limit as λ goes to infinity. More precisely, if one defines (y, p) and some corrector functions (y 1 , p 1 ) and (y 2 , p 2 ) such that then a standard asymptotic analysis procedure shows that (y, p) satisfies a pure incompressible problem. This formulation reads Observe that the function p in (IM) acts as a Lagrange multiplier to enforce incompressibility. If we assume that (IM) is the standard problem to solve, then Eq. (QIM) can be seen as an approximate penalised formulation of (IM). Existence, uniqueness and regularity results of problem (IM) can also be induced from (QI) by a limit process as λ → ∞.
We now introduce another formulation by penalisation of the problem (IM), inspired from existing formulations for the Stokes problem [14]. It reads First, note that (QIP) is different from (QIM) due to the introduction of the laplacian operator in the second equation. Furthermore, one can observe that (QIP) differs from the other formulations by several important but subtle aspects. First, the pressurep α is sought in a more regular space, namely C 0 (M), that is mandatory to give an appropriate meaning to the introduced Laplace operator. Second, system (QIP) is not a closed set of equations. Indeed, a boundary condition is required for the second equation for the pressure, that now has a trace (we recall that we use homogeneous Dirichlet conditions on the displacement). A standard choice is to use homogenous Neumann boundary conditions Then, using the same arguments as before, i.e., writing one can see that (ỹ α ,p α ) can be approximated by (y, p), solution of the pure incompressible mixed formulation (IM). Reciprocally, (ỹ α ,p α ) represents another approximation of the pure incompressible mixed formulation (IM). Moreover, eliminating the unknownp α , it is possible to rewrite (QIP) as where Δ −1 : L → M stands for the inverse Laplace operator equipped with a homogeneous Neumann boundary condition. It is possible to prove that the operator −div Cε(·) − α −1 ∇(−Δ) −1 div(·) defines a self-adjoint coercive bilinear form. Consequently, existence and uniqueness of the solution can be retrieved from (6) by standard theory.
Remark 1 Note that the choice (4) is arbitrary. In the context of the Stokes equations it was observed (see [31] and reference therein) that it results in a boundary layer that deteriorates the approximation of the gradient of the solution. Correcting terms can be introduced in specific cases (see again [31]) but their analysis is more difficult.

Weak formulation of the continuous PDE
Let us first consider the weak formulation associated with (IM). Given f sufficiently regular, it reads: where we have defined the bilinear forms m : Note that the bilinear form m(·, ·) is symmetric and positive. Furthermore, due to Eq. (2) and the Korn inequality, the bilinear form a(·, ·) is symmetric and coercive in X. We can write (7) as a set of equations written in X × L for all t ∈ [0, T ]. To do so, we introduce the linear continuous operators M : and we define the divergence operator B : X → L and its transpose B T : L → X such that, ∀ (y, q) ∈ X × L, where (·, ·) L denotes the standard scalar product of L 2 (Ω). Note that the operator M can be continuously extended to an operator from H to H and, if ρ ≡ 1, then it is the identity operator. Finally, the variational formulation can be equivalently rewritten and Analogously, the variational formulation associated with system (QIM) reads, equivalently, and where I is the identity operator from L to L. Concerning the problem (QIP) with the boundary condition (4), we propose a variational formulation that reads It is not straightforward to write equations in dual spaces from the variational formulation (10). This is due to the fact that we have changed the functional space in which the pressure is sought. To do so, we introduce the divergence operator C : H → M and its transpose where (·, ·) H denotes the standard scalar product on L 2 (Ω) d . Then, by identification of the operator B as an operator from X to M (instead of an operator from X to L), one can show that (10) Note that the operator C corresponds to an extension of the operator B. Inversely, the gradient operator B T represents the extension of C T in a larger space. This distinction in the notation is not relevant in the continuous framework, but it will be fundamental at the discrete level.

Space discretisation
Let us now consider a regular finite-dimensional space X h ⊂ X for the discretisation of the displacement field and M h ⊂ M for the discretisation of the pressure field. These spaces are obtained by finite element approximation of X and M, respectively. Furthermore, inspired by [29], in order to discretise accordingly the variational formulation (10) that we propose, we introduce a third finite-dimensional space denoted Y h ⊂ H that should satisfy, for the sake of simplicity, where (·, ·) Y h stands for the approximation of the scalar product in H that is defined using quadrature formulae (it is a symmetric coercive and continuous bilinear form in Y h for the norm in H). In a more general way, we assume that quadrature formulae do not corrupt symmetry and positivity properties of the bilinear forms. We introduce the discrete divergence operator C h : Y h → M h and the discrete gradient operator h corresponds to the operator C T ≡ −∇ applied to functions in M h , considering the resulting functions in the larger space Y h . Then, we define another discrete divergence operator B h : X h → M h and another discrete gradient operator One can observe that, for all ( Note also that the following commutative diagrams (taken from [29]) hold where the subscript h in the notation of the bilinear forms stands for the use of quadrature rule in the computation of integrals. The finite element approximation of (IMv) reads and where f h (t) ∈ X h denotes some approximation of f . Analogously to the pure incompressible mixed formulation, we can retrieve the space discretisation associated with (QIMv). It reads Find (y and where (·, ·) L h stands for the approximation of the scalar product in L by quadrature formulae. Finally, we are able to give the space discretisation associated with the novel formulation we propose (QIPv). It reads for any reasonable choice of finite element spaces and quadrature rule. This is obviously also true for I h and M h .
Let us insist on the importance of the introduction of the space Y h . First, it is related to the definition of the quadrature formulae in H in the definition of (12). Second, even if exact integration is performed, the introduction of the space Y h allows us to take into account the fact that, in general, the gradient of functions in M h does not belong to X h . Indeed, if this was the case, then C h and B h would be the same operator and the penalisation strategy would be useless in terms of computational efficiency (see the discussion in "Fully discrete schemes" section).
In the numerical results we present in this work we use high-order Spectral Finite Element, as in [38] and [39]. Since we consider a simple geometry for our purposes, we construct a quasi-uniform triangulation of Ω composed of quadrangles or hexahedra whereK is the unit square or the unit cube and ∀ i ∈ {1, 2, . . . , N }, F i denotes the invertible transformation of the reference elementK to the deformed element K i . Then, we define where Q r is the set of polynomials with degree r 1 in each variable of space. To obtain mass-lumping (meaning that M h can be inverted trivially) one must choose the quadrature points in the computation of m h (·, ·) at the same location as the interpolation points (see [40]). Sufficient accuracy is obtained if the interpolation and quadrature points correspond in the reference elements to the Gauss-Lobatto points of same order (r in our case, that gives (r + 1) d quadrature/interpolation points in the reference elements). For the computation of a h (·, ·) we use the same quadrature rule as for the computation of m h (·, ·) (which gives sufficient accuracy on non-distorted mesh, see [41]). Then, we choose Note that the major difference here is that we choose a lower-order finite element space, but M h is still constructed using continuous finite element. It is mentioned in [9] that this choice is compatible with X h in the sense that a discrete inf-sup condition is satisfied (the importance of this condition is further detailed below). In addition, we assume that (·, ·) L h is computed using Gauss-Lobatto quadrature points of order r − 1. Hence, mass-lumping is achieved and the operator I h is easily invertible. Finally, we define and (·, ·) Y h is computed using the Gauss-Lobatto quadrature points of order r − 1. Note that we have the inclusion ∇M h ⊂ Y h only if all the F i are affine. We believe that this is only a technical limitation. However, numerical results are provided only in that case.

Time discretisation
This section deals with the time discretisation of the semi-discrete formulations obtained by FE approximation in space. We consider only finite difference schemes that are centred, in order to preserve energy conservation at the discrete level. In what follows, the fully discrete schemes for the standard formulations (IM) and (QIM) are provided. Moreover, we propose the fully discrete scheme for the novel formulation (QIP).

Fully discrete schemes
Let us consider a time interval [0, T ], with T > 0, and define the partition t n = n Δt, with n ∈ {0, 1, . . . , N }, and Δ t = T /N . The fully discrete scheme corresponding to (IM) for n ∈ {0, 1, . . . , N } is constructed based on a simple second-order finite difference scheme, namely a leapfrog scheme.
then, y n+1 h is given by One can see that the pressure field is an intermediate unknown that acts as a Lagrange multiplier to enforce the constraint B h y n+1 h = 0. Note that we use the notation p n h by analogy with the quasi-incompressible schemes that we present in what follows. Moreover, the system (16) is well-posed if B T h is injective and B h surjective. This corresponds to verify that the LBB condition is satisfied, i.e. there exists a constant c > 0 such that where, if c does not depend on the discretisation parameter, one recovers an optimal convergence behaviour.
Note that for the incompressible, linear Stokes problem it is standard to show [42][43][44][45] that standard second-order discrete space-time discretisation can be achieved as soon as the LBB condition holds. We assume that similar results hold for the elastodynamic problem (IMh). Finally, note that Eq. (16) can be solved at each time step by iterative algorithms and, since the underlying problem is symmetric, the Conjugate Gradient method is a good candidate. However, it is important to highlight the fact that the opera- Remark 2 If non-zero initial displacement y λ (t = 0) = y 0 and/or initial velocity ∂ t y λ (t = 0) = v 0 are considered, or if the source term f h is not 0 at time t = 0, then, to preserve the expected second-order consistency, the computation of the first two iterates is performed as follows: where p 0 h is computed from (16) with n = 0, and (y h,0 , v h,0 ) belong to X h × X h and correspond to an approximation of (y 0 , v 0 ).
Observe that the Scheme (QIMnh) is fully explicit, due to the use of the mass-lumping technique (we recall that M h and I h are easily invertible). However, we show in the following section that, due to stability considerations, the maximum time step allowed is considerably reduced by the fact that the pressure term is treated explicitly. Observe that the first two iterates can be computed using (18) Note that here, for consistency reasons, we have rescaled the penalisation parameter by Δt 2 and assume α is independent of Δt. This choice should guarantee the second-order consistency in time that is expected from the leapfrog time discretisation. Note that we can directly rewrite the second equation in (QIPnh) as Consequently, this step is equivalent to solving a discrete Poisson problem for the pressure at each time step, with homogeneous Neumann boundary conditions on the boundary ∂Ω. One of the main advantages of this formulation is that the Poisson problem is very standard and its solution can be retrieved by fast solvers. In addition, observe that the first two iterates can be computed using (18) (20) is small.

Remark 4
The SMS method, first introduced in [46], consists in adding an inertial term. In the specific case of quasi-incompressible materials, it is suggested in [37] to add the volumetric contribution of the stiffness operator. In our framework, starting from (QIMnh), after eliminating the pressure term p n λ,h , by SMS strategy one could obtain the Scheme where β is some well-chosen parameter. Note that the operators I h and B h embed the definition of adequate quadrature formulae to avoid numerical locking effects (typically, reduced integration must be used The scheme above is an explicit-implicit hybrid scheme (see [47]) which has a stability condition depending only on (M h , A h ) (i.e. the operators related to the shear wave propagation only) as soon as θ 1/4. Although very good stability properties can be proved, at each time step one has to find the solution of a linear system for the displacement field and, thus, it not possible to use fast solvers for the inversion of a scalar Poisson problem. This last point motivated the introduction of the Scheme (QIPnh).

Stability analysis
The aim of this section is to find uniform estimates of the discrete energy of the different schemes, i.e.
where the constant C depends on the final time step T = N Δt and on the data of the continuous problem, but is independent of Δt and h. However, if explicit schemes are employed for time discretisation, the time step is limited by a stability condition depending on h. We refer to [36] for further reading. For the sake of simplicity, we consider a time t n such that the source term has vanished (i.e. f h (t) = 0, ∀t t n ). Then, the energy of the continuous problem is constant in time. In order to retrieve a discrete counterpart of this energy, we consider for every formulation, as a test function, the centred discrete approximation of the time derivative of the displacement at time t n .

Stability of scheme (IMnh)
Let us first consider formulation (IMnh). By scalar product in the first equation in ( We define the discrete energy at time n + 1 2 as Then, after some computations, using the symmetry properties of the operators M h and A h , we obtain from (25) the discrete conservation property It now remains to prove that E n+ 1 2 h is positive.

Proposition 1 A sufficient condition for the stability of scheme (IMnh) is
with Proof The proof is very standard (see [36]). Consequently, we only prove the positivity of the energy. We provide here some details for the sake of completeness. Let us first consider the definition of E is positive by definition, we can easily retrieve the lower bound for the energy Hence, the energy is positive if Finally, Eq. (31) can be rewritten as oncluding the proof.
Remark 5 Note that Proposition 1 introduces an abstract CFL condition. Moreover, we expect the following estimation with c s a positive constant, depending on the elasticity tensor driving the shear wave propagation. Consequently, we obtain the sufficient stability condition Therefore, the time step is not affected by the pressure wave propagation, that is travelling at an "infinite" velocity at the incompressible limit.

Stability of the scheme (QIMnh)
By similar reasoning, we can retrieve en energy estimation for the formulation (QIMnh

Proposition 2 A sufficient condition for the stability of scheme (QIMnh) is
with Remark 6 Note that, by definition, By similar reasoning to Eq. (33), we can now introduce a constant c p , related to the maximum generalised eigenvalue of the operator Therefore, we can assert Consequently, the stability condition (36) imposes a significant restriction on the time step.
Note that, because of the non-dimensionalisation, we expect c s to be close to the unity, whereas, c p is given by the velocity ratio between pressure and shear waves. For soft tissues this ratio is around 10 3 . This makes in practice the Scheme (QIMnh) not efficient and justifies our need to formulate more adequate methods for the limit-incompressibleproblem.

Stability of the scheme (QIPnh)
Finally, let us analyse the stability estimates related to the novel formulation (QIPnh). By analogous reasoning, we get from (20) where we have defined Then, the discrete energy at time n + 1 Before providing a stability estimate for (40), we introduce the following lemma.
Then P h is a projection and P h L(Y h ) 1.
Proof The proof consists in demonstrating that P 2 h = P h . That follows easily by the definition of P h .
We are now able to assert a stability condition for Scheme (QIPnh).

Proposition 3 A sufficient condition for the stability of scheme (QIPnh) is
and α > 1 4 sup Proof Again we prove the positivity of the energy E is positive, we can obtain the estimation Hence, we need to satisfy Finally, if (42) holds, then Eq. (45) can be rewritten as concluding the proof.  (42) to obtain the result of the corollary.

Remark 7
We observe that the CFL condition (47) is well-defined for α > 1/(4 ρ). Moreover, it is slightly worse than condition (28), since in the allowed range for α. Nevertheless, this condition depends only on the space discretisation and on the tensor C and the density ρ. Therefore, it is still very advantageous with respect to condition (36).
As a final comment, note that if the Scheme (QIPnh) is stable, i.e. (22) holds and Δt is sufficiently small, then with C independent of Δt. Therefore, denotingỹ As a consequence, there exists another constant C independent of Δt such that, by definition of Q h , we have sup n∈{0,1,...,N } This shows that B hỹ n+ 1 2 α,h goes to 0 with Δt, but for a weak norm that involves the inverse Laplace operator.

Two-dimensional numerical convergence results
In order to perform the numerical validation of the properties of scheme (QIPnh), we consider as a model problem the elastic wave propagation in a 2D medium, and we take into account different constitutive laws. In the interest of clarity, we provide the physical and numerical parameters used for the simulations with their original units of measure. For the sake of simplicity, we assume constant density ρ = 1050 kg m −1 and homogeneous Dirichlet boundary conditions on all the boundaries of the domain. Space discretisation is performed by high-order Spectral Finite Elements (of order 7 for the displacement, 6 for the pressure). The computational grid is a 1 m 2 square composed of N uniform elements of size h = 1/N in each direction. Concerning time discretisation, we adopt the discretisation introduced in scheme (QIPnh) with penalisation coefficient α = 1/3 ρ and we choose the time step as where we set ε = 0.2 to account for the fact that the expression above is an approximation of the CFL condition given by Eq. (42) when considering that ρ y h 2 Y h = m h (y h , y h ) and accounting for the approximation of the norm in (49) by a power iteration algorithm.
Note that the value of this parameter is not tight: ε = 0 also gives stable results. For each iteration of this scheme, the scalar pressure field (19) is computed by means of an in-house fast solver for the Poisson problem, based on Higher-Order Fourier Transform. Note that even when the medium is heterogeneous, the problem for the pressure wave remains homogeneous, due to incompressibility.
Henceforth, we present a space/time convergence analysis. To do so, we test scheme (QIPnh) against the pure incompressible scheme (IMnh) for different values of N in [10,40]. For each time step of each simulation, the pressure field is evaluated in Eq. (16) by means of the Conjugate Gradient method, with euclidian norm of the residue lower than 1e − 14 and maximum number of iterations N iter = 2000. The source term considered is a standard Gaussian profile in space multiplied by a profile in time corresponding to the first derivative of the Gaussian function. In greater detail, it reads with centre (x F , y F ) = (0.5, 0.5) m, covariance σ 2 s = 0.0005, σ 2 t = 0.0005 and mean t pulse = 0.6 s. We consider several constitutive laws, and we present four convergence curves for each example, that illustrate the error in L 2 (Ω T ), C 0 (L 2 (Ω)), L 2 (H 1 (Ω)) and C 0 (H 1 (Ω)) norms, respectively.

Homogeneous isotropic material
As a first example, we study a homogenous isotropic medium. A standard constitutive law for this type of medium is C ε(y) := μ ε(y). (51) with shear modulus μ = 80 kPa. Note that from Eq. (51) we retrieve in (QI) the standard elastodynamic problem for isotropic law. The convergence curves in Figs. 1 and 2 confirm that second-order convergence is preserved in L 2 (Ω T ). Note that it is slightly degraded in and L ∞ in time and L 2 norm in space. Furthermore, the error does not vary [or it varies marginally in L ∞ (L 2 (Ω))] if we consider a larger t end such that the wave has reached the boundaries of the computational domain. If we consider H 1 norm in space, the error  is slightly degraded with respect to the previous norms. Moreover, when we consider a larger time interval (so that we take into account the effects at the boundaries), the error increases. Figure 3 illustrates three snapshots at t = 1.0 s, t = 2. s and t = 2.65 s, of the displacement. The mesh is composed by N = 40 elements in each direction. Furthermore, Fig. 4 depicts the time evolution of the displacement in three locations, indicated in Fig. 3. Figures 5 and 6 are related to the acceleration. We do not plot the results concerning the velocity, since they do not provide relevant information. Note that the time evolution profiles in x and y direction of the location denoted P 1 correspond to the profiles in y and x direction of the location denoted P 2 , that is the symmetric to P 1 , due to isotropy.
Concerning the pressure, Fig. 7 indicates that this term corresponds to a correction mostly at the boundaries of the medium. In addition, symmetric points correspond to the same pressure time evolution, as depicted in Fig. 8.

Heterogeneous transversely isotropic material
In the perspective of an application to elastography imaging of biological tissues, we need to analyse more complex constitutive laws. For example, striated muscle tissue can be modelled as a transversely isotropic medium, i.e. there exists, at every point, a privileged direction represented by the unit vector τ 1 , related to the collagen fibre. A simple transversely isotropic law, inspired by [48] is with μ = 80 kPa, η = 3400 kPa. In what follows, we present the results for a transversely isotropic medium in which the fibre direction τ 1 varies linearly along the direction y. This configuration is inspired by the structure of myocardial tissue. In fact, it has been experimentally validated that muscle fibres are arranged in laminar structures, denoted sheets, of three to four muscle fibres in the thickness, that are oriented transversely to the heart wall [49,50]. Moreover, the fibre orientation in human left ventricle myocardium changes linearly throughout the wall thickness, from − 60 • close to the epicardium to + 60 • near the endocardium [51]. In order to model such a material, we consider an orientation of τ 1 varying linearly from − 60 • at the bottom of the domain (y = 0) to + 60 • at the top of the domain (y = 1). The resulting medium is highly heterogeneous. Figures 9 and 10 show that the error worsens moderately in transversely isotropic media w.r.t. isotropic media. In addition, the order of convergence remains two for L 2 norm in space, and it is slightly degraded in L 2 (H 1 (Ω)) and in L ∞ (H 1 (Ω)). Figure 11 illustrates three snapshots at t = 0.5 s, t = 0.75 s and t = 1.1 s of the displacement field, respectively. The mesh is composed by N = 40 elements in each direction. We do not provide the results concerning the velocity and the acceleration, since they not provide additional information. Note that the pressure field gives a significant contribution also in the interior of the domain, as it is illustrated in Fig. 12, due to anisotropy.

A three-dimensional test case
The main application we have in mind is the propagation of elastic waves in nearlyincompressible biological tissues in the context of ultrasound transient elastography [52,53]. Therefore, by way of illustration, we propose here a three-dimensional test case, inspired by elastography imaging of the myocardial tissue. Note that, due to dissipation, shear waves are fully attenuated in a few wavelength distance, i.e. some millimetres from the focal point of the ultrasonic beam [54]. As a consequence, we can consider a small region of interest and a simple geometry. However, there is a need for a fine space discretisation, compatible with the wavelength of shear waves. In this example we consider a transversely isotropic domain [defined by Eq. (52)], with fibre direction τ 1 oriented in the xy-plane and varying linearly along the direction z from −60 • w.r.t. the x-axis at the bottom of the domain (z = 0) to + 60 • at the top of the domain (z = 0.02 m). The geometry of the domain is a parallelepiped of length 0.04 m in x and y directions, and 0.02 m in z direction. This configuration represents an approximation of a region of interest in the left myocardium, selected in the anterior wall at the middle ventricular level. Density is set equal to ρ = 1050 kg m −1 and homogeneous Neumann boundary conditions are imposed on all the boundaries of the domain. High-order Spectral Finite Elements (of order 7 for the displacement, 6 for the pressure) are used for the space discretisation. The computational grid is composed of 24 uniform elements of size h = 1/24 in each direction, for a total of 14,480,427 degrees of freedom (DOF) for the displacement field, and 3,048,625 DOF for the pressure field. We adopt the time discretisation introduced in scheme (QIPnh) with penalisation coefficient α = 1/3ρ and we choose the time step as in Eq. (49), with ε = 0.2, as in "Two-dimensional numerical convergence results" section. For each time step of the scheme, the Poisson problem for the scalar pressure field (19) is computed by means of our in-house fast solver based on Higher-Order Fourier Transform, as in "Two-dimensional numerical convergence results" section. In this way, the time step used in the explicit time discretisation is constrained by the shear wave velocity only. Figure 13 depicts four snapshots of the solution at different time steps.
Note that we have also considered an improved approximation by post-processing of the problem, by combination of the solutions associated with α 0 = 1/3 ρ and α 1 = 2/5 ρ, respectively, considering N = 24 elements in each direction. Since we have not noticed a qualitative change in the behaviour of the solution, we are confident that the observed solution is a good approximation of the pure incompressible problem. Note also that, due to the resolution of the Poisson problem by our in-house fast solver, we have drastically reduced the computational cost of the scheme. In particular, the resolution of 1600 time steps of the problem on a 12-cores workstation (cores at 2.7 GHz and 64 GB of RAM at 1867 MHz), considering higher-order extension of the scheme (i.e. computation of the two solutions corresponding to α 0 and α 1 ) takes around 4 h.

Conclusions
In this article we have outlined a new numerical scheme that is suitable for approximating elastic wave propagation in incompressible media, based on a mixed explicit/implicit time discretisation. We have demonstrated that the time step of the resulting algorithm is only constrained by the shear wave velocity, and it is not limited by the enforcement of incompressibility (as it is the case for fully-explicit time discretisation). Furthermore, if effective methods are adopted for explicit time-discretisation, our algorithm only entails at each time step one resolution of a scalar Poisson problem-that can be performed by various, efficient algorithms-and few matrix-vector multiplications for the explicit part of the scheme. Moreover, we have presented a two-dimensional numerical test case, to demonstrate the favourable convergence properties of the scheme, and one three-dimensional example, to illustrate a realistic application to elastography imaging. A theoretical demonstration of second-order convergence in time of the scheme is out of the scope of this article and will be provided in future work.
Authors' contributions SI proposed the subject; FC and SI developed the method and implemented the scheme on an in-house finite element code (previously developed by SI); FC performed the finite element simulations and the numerical analysis under the technical supervision of SI; FC is the main author of the paper; critical revision was provided by SI. Both authors read and approved the final manuscript. 1 Inria, Université Paris-Saclay, Palaiseau, France, 2 LMS, Ecole Polytechnique, Université Paris-Saclay, CNRS, Palaiseau, France.