# Fully explicit numerical scheme for linearized wave propagation in nearly-incompressible soft hyperelastic solids

Giulia Merlini¹,², Jean-Marc Allain¹,², Sebastien Imperiale²,¹

¹LMS, Ecole Polytechnique, CNRS – Institut Polytechnique de Paris; ²Inria, Team-MΞDISIM, Inria-Saclay Ile de France 91120 Palaiseau

For correspondence:
sebastien.imperiale@inria.fr

## Abstract

The numerical approximation of wave propagation problems in nearly or pure incompressible solids faces several challenges such as locking and stability constraints. In this work we propose a stabilized Leapfrog scheme based on the use of Chebyshev polynomials to relax the stability condition, which is strongly limited by the enforcement of incompressibility. The scheme is fully explicit, second order accurate and energy-preserving. For the space discretization we use a mixed formulation with high-order spectral elements and mass-lumping. A strategy is proposed for an efficient and accurate computation of the pressure contribution with a new definition of the discrete Grad-div operator. Finally, we consider linear wave propagation problems in nearly-incompressible hyperelastic solids subject to static preload.

## Contents

1 Introduction 2

2 Modelling elastic wave propagation in nearly incompressible non-linear solid 4

2.1 Kinematics and the non-linear elastodynamic problem 4

2.2 Elastic wave propagation in quasi-static loaded solid 6

2.2.1 An example of strain energy function for nearly incompressible material 7

3 A fully explicit numerical scheme 8

3.1 Linearization at rest condition 8

3.2 Space discretization 10

3.2.1 Semi-discrete formulation 10

3.2.2 Spectral finite elements 11

3.2.3 Definition of divergence operator 12

3.2.4 A specific finite element 13

3.3 Time discretization 13

3.3.1 The leap-frog method 13

3.3.2 Stabilized leapfrog method 13

3.3.3 Leapfrog Chebyshev Methods 15

3.4 Efficiency estimation the LFC scheme 17

4 Numerical Results 19

4.1 Non-locking mixed formulation 19

4.2 Leapfrog Chebyshev Methods 21
<table><tbody><tr><td>4.2.1</td><td>Efficiency</td><td>21</td></tr><tr><td>4.2.2</td><td>Space-time Convergence analysis</td><td>22</td></tr><tr><td>5</td><td>Numerical application: elastic wave propagation in an incompressible hyper-elastic stretched cube</td><td>23</td></tr><tr><td>6</td><td>Conclusion</td><td>27</td></tr><tr><td>A</td><td>Derivation of elastic wave propagation in quasi-static loaded nearly-incompressible solid</td><td>30</td></tr><tr><td>A.1</td><td>Elastic wave propagation in quasi-static loaded solid</td><td>32</td></tr><tr><td>A.1.1</td><td>Wave propagation problem in the reference configuration</td><td>32</td></tr><tr><td>A.1.2</td><td>Wave propagation problem in the apparent configuration</td><td>34</td></tr></tbody></table>

# 1 Introduction

The accurate and efficient simulation of wave propagation in nearly incompressible solids represents a significant challenge in numerical modeling. Classically, the finite element approximation of solids relies on the displacement-based method, in which it is assumed that the displacement of a body is completely specified by the nodal point displacement [1]. However, the enforcement of incompressibility highly affects the accuracy of the solution, such that a finer space discretization will be needed, compared to the compressible case. This is due to the fact that volumetric strains are computed with a lower accuracy with respect to the displacements and are close to zero in the incompressible limit. Any error related to the computation of the volumetric strain will result in a large error in the computation of the stress (pressure) due to the multiplication by a large bulk modulus, and then a wrong prediction of the displacement obtained from the principle of virtual work [2]. With a higher spatial resolution, we may circumvent errors in the computation of the displacement field, but we will still have spurious or incorrect pressure. In addition, it is known [1, 3, 4, 5] that for nearly-incompressible solids, this method suffers from ill-conditioning of the stiffness matrix. In general, we may refer to formulation whose predictive capability is dependent from the bulk modulus as locking formulations [2, 6]. Nevertheless, efforts have been done to overcome these limitations, since most of large-scale structural analysis programs are based on displacement-based method [7]. We cite in particular the B-bar method for finite-deformation problems [8] and F-bar method for hyperelastic problems, which rely on the reduction of volumetric constraints through projection.

Mixed methods form another class of methods for finite element approximations, that have proven to be effective to avoid locking behavior in the framework of elasticity with nearly or pure incompressible materials. These methods are based on the introduction of the pressure as an additional unknown. The problem becomes now a saddle-point problem where the key aspect for the finite-element approximation is the introduction of an auxiliary finite element space for the computation of the pressure field. These methods are proven to be stable under the well-known LBB inf-sup condition derived by Brezzi and Babuska, together with the ellipticity condition [9, 10, 11], for which the choice of the finite element spaces is fundamental to guarantee stability and convergence of the method. When the stability is not given for a certain choice of the finite element spaces, it is possible to apply stabilization techniques such as cross-grid triangulation or with the introduction of a stabilization in order to relax the incompressibility constraint. We refer the reader to [12] for a review of stabilization methods [13] such as Galerkin-Least Square [14], Streamline Upwind Petrov Galerkin [15] and bubble functions [16, 17].

In addition, the simulation of dynamic problems must account for an appropriate discretization strategy in time. For elastic wave propagation problem it is generally prohibitive to use implicit schemes, as they require the inversion of the stiffness matrix at each time-step, which is too large to be stored even when considering factorization methods. Explicit is well suited given the phenomena that we are observing. In particular, the leapfrog (LF) scheme based on centered finite difference of second order time differential equation benefits of several properties:
second-order accuracy, energy conservation for linear problems and time reversibility [18]. However, the stability of the scheme is given under the classical CFL condition for the maximal time-step, which is highly constraining when considering nearly incompressibility. As the bulk modulus increases, the time-step goes to zero, with a corresponding increase of number of iterations and of computation time.

A possible alternative is to use fractional step projection methods, in which the idea is to split the original problem in two simpler ones, such that at each time step we compute implicitly a first half-step considering one sub-problem, and then we correct the obtained solution with the second half of the problem that we compute explicitly [19]. We cite in particular [5], which has been a starting point for this work in order to treat linear elastic wave propagation problems in nearly or pure incompressible solids. The conservative scheme proposed by Caforio et al. [5] is based on the relaxation of the volumetric component of the displacement solution, related to the incompressible constraint, by solving at each iteration in time a Poisson problem. Then, the explicit leapfrog time-discretization method is used to solve the second part of the problem to account for the "shear" part of the problem. Solving the Poisson problem must then involve the use of fast solvers, like FFT solvers.

In this work, we focus on the development of a fully explicit scheme, such that only the mass matrix will be inverted together with mass-lumping approximation [20]. Recently, [21] have proposed a new class of methods, named Leapfrog-Chebyshev (LFC) scheme for the relaxation of the CFL condition. They applied Chebyshev polynomials to the linear component of the problem, while the non-linear part is computed with the LF scheme. These schemes have proven to be effective together with local-time stepping methods (LFC-LTS) [22], which are developed to deal with stability constraint related to local mesh refinement. Here, we would like to apply the same idea, but for the incompressibility constraint. We propose, indeed, a novel explicit formulation to treat linear wave propagation in nearly-incompressible problems, such as follows

$$ \underline{u}_h^{n+1} - 2\underline{u}_h^n + \underline{u}_h^n + \Delta t^2 M_h^{-1} A_h^s \underline{u}_h^n + \mathcal{R}_{m+1}^{\nu} (\Delta t^2 \lambda M_h^{-1} \bar{A}_h^p) \underline{u}_h^n = \Delta t^2 M_h^{-1} \underline{f}^n $$

where $M_h$ is the lumped mass operator, $A_h^s$ corresponds to the shear-stiffness operator (independent of volumetric parameters), $\bar{A}_h^p$ is a discretized Grad-div operator multiplied by the bulk modulus $\lambda$ taken positive and very large, and $\mathcal{R}_{m+1}^{\nu}$ is the *damped* Chebyshev polynomial of order $m+1$, to be later presented. As explained in [21], it should be interpreted as a multi-rate scheme, such that we evaluate the shear part once, while the volumetric part is iterated as many times as the order $m+1$ of the Chebyshev polynomial. The scheme is stable under the following condition

$$ \Delta t^2 \le \frac{1}{\rho(M_h^{-1} A_h^s)} \quad \text{and} \quad m+1 \ge \sqrt{\lambda} \frac{\sqrt{\rho(\Delta t^2 M_h^{-1} \bar{A}_h^p)} e^{\nu/4}}{2} $$

where we denote $\rho(\cdot)$ as the spectral radius of the operator. We underline that the new time-step is now only dependent on the shear operator, while the incompressibility constraint determines the minimum order of the Chebyshev polynomial to guarantee stability.

In other words, the *objective* of this method is the same pursued by previous works [23, 5] in the framework of incompressible (or nearly incompressible) elastodynamics: have an efficient numerical scheme while maintaining the accuracy of the solution. In [23], the proposed approach relies on selective mass scaling (SMS) [24], which results in an artificial mass matrix that cannot be lumped anymore, such that its inversion becomes quite expensive for explicit time-domain computations. The previously mentioned work of [5] requires the use of FFT solvers, which are well suited for domains with a simple geometry, while here we would like to get rid of this constraint.

This method has been developed in the perspective of reproducing dynamic elastographic measurements of biological tissues, for the analysis of the material properties and the detection of pathologies. These materials are precisely characterized by a high water content that makes them almost incompressible, and they exhibit hyperelastic (nonlinear) behavior. We must also consider that in a physiological framework, these tissues are generally subjected to
certain loading conditions, which then generate a state of pre-stress and pre-deformation of the tissue that affects wave propagation during elastographic measurements. In the following work we have derived a novel formulation for linear wave propagation in nearly incompressible preloaded solids, for which we separate the Grad-div operator (volumetric term) from the stiffness operator. The isolation of the term related to the enforcement of incompressibility allows us to tackle the critical aspects related to accuracy, stability and efficiency that we have mentioned before with the presented Chebyshev Leapfrog scheme.

## 2 Modelling elastic wave propagation in nearly incompressible non-linear solid

In this section we introduce the continuous formulation of wave propagation problem in the framework of hyperelastic nearly-incompressible solids. In particular, we consider a small amplitude elastic wave propagating in the solid subjected to a quasi-static load with large deformations. We follow the procedure of linearization as described in [25, 26, 27, 28], with the primary aim of presenting the problem and the related notation. The second objective is to present the weak formulation for linear wave propagation problems in nearly-incompressible preloaded solids, for which the volumetric term is separated from the stiffness operator and is written as a Grad-div operator. The derivation of the formulation is detailed in Appendix A

### 2.1 Kinematics and the non-linear elastodynamic problem

We define $\Omega(t)$ as the domain occupied by a deformable solid at time $t$. We assume that its boundary $\Gamma(t)$ is sufficiently smooth and denote $\underline{n}(t)$ the unitary outward normal. For the formulation of the non-linear problem we use a Lagrangian description, which consists in representing the position of points in $\Omega(t)$ with respect to a reference configuration $\hat{\Omega}$. All the material points in $\hat{\Omega}$ are described by the coordinate system $\xi$ and we can retrieve the current (at time $t$) position of points in $\Omega(t)$ with the bijective deformation map $\underline{\phi}(\xi, t)$:

$$ \underline{\phi} : (\xi, t) \to \underline{x}(\xi, t) = \underline{\phi}(\xi, t), \quad \underline{\phi}(\hat{\Omega}) = \Omega(t). $$

From the definition of the bijective map, we define the displacement field

$$ \underline{u}(\xi, t) = \underline{\phi}(\xi, t) - \underline{\xi} = \underline{x}(\xi, t) - \underline{\xi} $$

and the deformation gradient with its determinant

$$ \underline{F}(\underline{u}(\xi, t)) = \underline{\nabla}_{\underline{\xi}} \underline{\phi} = \underline{1} + \underline{\nabla}_{\underline{\xi}} \underline{u} \quad \text{and} \quad J(\underline{u}) = \det \underline{F}(\underline{u}). $$

The volume element $d\hat{\Omega}$ in the current configuration is then given by the quantity $J(\underline{u})d\hat{\Omega}$, while the current area element $ds$ is equal to $J|\underline{F}^{-T} \cdot \underline{n}|d\hat{s}$. We suggest the reader to look into [29] for details. Before introducing the relevant mechanical quantities to describe the deformation of the solid, we underline that for second order tensors we use the dot-product and double dot-product defined as follows

$$ \underline{A} \underline{B} = A_{ik} B_{kj} \underline{e}_i \otimes \underline{e}_j, \quad \underline{A} : \underline{B} = \text{tr}(A^T B) = A_{ij} B_{ij} $$

We then introduce the right Cauchy-Green deformation tensor $\underline{C}(\underline{u}) = \underline{F}(\underline{u})^T \underline{F}(\underline{u})$, the Green-Lagrange strain tensor $\underline{\epsilon}$ and its linearized expression $\underline{\epsilon}$

$$ \underline{\epsilon}(\underline{u}) = \frac{1}{2}(\underline{C} - \underline{1}) = \frac{1}{2}(\underline{\nabla}_{\underline{\xi}}\underline{u} + \underline{\nabla}_{\underline{\xi}}\underline{u}^T + \underline{\nabla}_{\underline{\xi}}\underline{u}^T\underline{\nabla}_{\underline{\xi}}\underline{u}), \quad \underline{\epsilon}(\underline{u}) = \frac{1}{2}(\underline{\nabla}_{\underline{\xi}}\underline{u} + \underline{\nabla}_{\underline{\xi}}\underline{u}^T). $$

### Fundamental law of dynamics

Before presenting the functional spaces in which we will seek the displacement solutions, we define the boundary $\Gamma(t) = \Gamma_D \cup \Gamma_N(t)$, where $\Gamma_N(t)$ is the part of the boundary on which we have a free surface condition, while on the rest of the boundary we consider Dirichlet's condition,
meaning that the displacement vanishes. Then, the functional spaces of the displacement field are defined as follows

$$ \mathcal{V}(\hat{\Omega}) := \{\hat{\underline{u}} \in H^1(\hat{\Omega})^d \mid \hat{\underline{u}} = 0 \text{ on } \Gamma_D\}, \quad \mathcal{V}(\Omega) := \{\underline{u} \in H^1(\Omega)^d \mid \underline{u} = 0 \text{ on } \Gamma_D\} $$

with dimension $d = 2$ or $d = 3$, and

$$ \mathcal{W}(\hat{\Omega}) \subset \mathcal{V}(\hat{\Omega}), \quad \mathcal{W}(\Omega) \subset \mathcal{V}(\Omega) $$

with $\mathcal{W}$ as the space of admissible displacements for the non-linear problem. We consider that the solid is subjected to a force $f$ per unit mass. Given the observation time $T$, we write the standard weak formulation associated to the wave propagation problem defined in the current configuration $\Omega$:

For any $t \in [0, T]$, find $\underline{u}(t) \in \mathcal{W}(\Omega)$ such that

$$ \int_{\Omega} \varrho(\partial_{tt}\underline{u} \circ \underline{\phi}^{-1} - \underline{f}) \cdot \underline{v} \, d\Omega + \int_{\Omega} \underline{\sigma} : \underline{\nabla}_{\underline{x}} \underline{v} \, d\Omega = 0 \quad \forall \underline{v} \in \mathcal{V}(\Omega), \quad (1) $$

where $\varrho$ is the density in the current configuration and $\underline{\sigma}$ is the Cauchy stress tensor defined in $\Omega(t)$. In classical continuum mechanics, this problem is named *principle of virtual work* and it is formally equivalent to the *equilibrium equation* with boundary condition of traction in the current configuration $\Omega(t)$. In this work we consider wave propagation phenomena in materials with non-linear behavior, namely hyperelastic materials, so that a constitutive law must be introduced to relate internal stress and deformation of the solid due to external loads.

We then introduce the Second Piola Kirchhoff tensor $\underline{\Sigma}$, which corresponds to the symmetric stress tensor defined in the reference configuration, related to the Cauchy stress tensor by the following relation:

$$ \underline{\sigma} \circ \underline{\phi} = J^{-1} \underline{F} \cdot \underline{\Sigma} \cdot \underline{F}^T. \quad (2) $$

In the framework of hyperelastic materials, the state of stress within the solid depends on the *hyperelastic potential* $\mathcal{W}^e = \mathcal{W}^e(\underline{C})$ and, in the case of *isotropic* material, we have

$$ \mathcal{W}^e = \mathcal{W}^e(I_1, I_2, I_3) \quad \text{s.t.} \quad \underline{\Sigma} = \frac{\partial \mathcal{W}^e}{\partial \underline{\underline{\underline{\underline{C}}}}} (I_1, I_2, I_3) = 2 \sum_{i=1}^{3} \frac{\partial \mathcal{W}^e}{\partial I_i} \frac{\partial I_i}{\partial \underline{\underline{\underline{\underline{C}}}}} $$

where $I_1, I_2, I_3$ are the eigenvalues of $\underline{\underline{\underline{C}}}(\underline{u})$ known as *invariants*

$$ I_1 = \text{tr} \underline{\underline{\underline{C}}}, \quad I_2 = \text{tr}(\text{cot} \underline{\underline{\underline{C}}}), \quad I_3 = \det \underline{\underline{\underline{C}}} = J^2. $$

As we are dealing with materials characterized by a non-linear and nearly-incompressible response, we actually define the *hyperelastic potential* with the following splitting

$$ \mathcal{W}^e(\underline{\underline{\underline{C}}}) = \mathcal{W}^d(J_1, J_2) + \mathcal{W}^v(J_3), \quad (3) $$

where $J_1, J_2, J_3$ are the so-called *reduced invariants* defined by

$$ J_1 = I_1 I_3^{-1/3}, \quad J_2 = I_2 I_3^{-2/3}, \quad J_3 = I_3^{1/2} = J. \quad (4) $$

Reduced invariants allow to split the potential in a component related to *deviatoric* ($\mathcal{W}^d$) and *volumetric* ($\mathcal{W}^v$) deformations. The nearly-incompressibility is taken into account with $\mathcal{W}^v$ that constrains volume changes. This term can be interpreted as a dilatation strain energy function that depends only on $J$, for which a small variation in volume is penalized by a large *bulk modulus*. The following linearized formulation is valid for any material whose response can be described with Eq. (3), and extendable to account for a non-linear anisotropic behavior.
## 2.2 Elastic wave propagation in quasi-static loaded solid

We now consider that the load applied to the solid is defined as follows

$$ \hat{f}(\xi, t) = \hat{f}_0(\xi) + \delta\tilde{f}(\xi, t) \qquad (5) $$

The solid undergoes large deformation due to the application of the quasi-static load $\hat{f}_0(\xi)$, and it is subjected to a transient small-amplitude perturbation $\delta\tilde{f}(\xi, t)$ with $\delta$ as a small parameter, that generates the wave propagation. Solutions can be written as

$$ \underline{u}(\xi, t) = \underline{u}_0(\xi) + \delta\tilde{u}(\xi, t) + \mathcal{O}(\delta^2), \qquad (6) $$

where $\underline{u}_0$ is the governing solution of the quasi-static non-linear mechanics and $\delta\tilde{u}$ is the small-amplitude wave propagating in the loaded configuration.

From these assumptions, we can start the linearization process. In Appendix A we have detailed the derivation of the governing equations for hyperelastic nearly-incompressible solids. Here, we limit to present the principal steps to follow and the final result for the wave propagation problem. Writing the weak formulation (1) with respect to the reference configuration $\hat{\Omega}$ assumed to be known, we can use Eqs. (??) in order to linearize the mechanical quantities associated. We then separate the non-linear terms from the first-order terms, while we neglect higher order terms, i.e., terms multiplied by $\delta^p$ with $p \le 2$. We then obtain:

* a non-linear static problem associated to large static deformation $\underline{u}_0$,

* a wave propagation problem linearized around $\underline{u}_0$, with solution $\tilde{u}$,

both written in $\hat{\Omega}$. Assuming that we are able to solve the first non-linear static problem and compute $\underline{u}_0$, we can introduce the corresponding quasi-static loaded configuration. We refer to it as the *apparent* configuration $\Omega_0$, whose set of points $\underline{x}_0$ can be retrieved with the bijective deformation map $\underline{\phi}_0(\xi)$:

$$ \underline{\phi}_0 : \xi \rightarrow \underline{x}_0(\xi) = \underline{\phi}_0(\xi) $$

$$ \underline{\phi}_0(\xi) = \underline{u}_0(\xi) + \xi $$

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Reference $\hat{\Omega}$, apparent $\Omega_0$ and deformed $\Omega$ configurations related by the respective deformation maps. The blue triangle represents a small amplitude transient external excitation.</figcaption></figure>

As we are interested in rewriting the linearized formulation with respect to the apparent configuration, we also introduce the functional space $\mathcal{V}(\Omega_0) := \{\underline{v} \in H^1(\Omega_0)^d\}$ where we seek the linear solution. We then introduce the subscript $\cdot_{\underline{x}_0}$ to indicate quantities defined with respect to the apparent configuration $\Omega_0$.
We obtain (see Appendix A for details) that the formulation of the elastic wave propagation problem in the apparent configuration reads:

**Wave equation in apparent configuration Ω₀:** For any $t \in [0, T]$, find $\tilde{u}(t) \in V(\Omega_0)$ such that

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega_0} \varrho_0 \partial_{tt} \tilde{u} \cdot \underline{v} \, d\Omega_0 + \int_{\Omega_0} \underline{\varepsilon}_{\underline{x}_0} (\tilde{u}) : \tilde{\mathbf{C}}_0 : \underline{\varepsilon}_{\underline{x}_0} (\underline{v}) \, d\Omega_0 \\[1.5em] \displaystyle + \int_{\Omega_0} \underline{\nabla}_{\underline{x}_0} \tilde{u} \underline{\sigma}_{\underline{z}_0} : \underline{\nabla}_{\underline{x}_0} \underline{v} \, d\Omega_0 \\[1.5em] \displaystyle + \int_{\Omega_0} \lambda_0^* \mathrm{div}_{\underline{x}_0} (\tilde{u}) \mathrm{div}_{\underline{x}_0} (\underline{v}) \, d\Omega_0 = \int_{\Omega_0} \varrho_0 \tilde{f} \cdot \underline{v} \, d\Omega_0 \end{array} \right. \quad \forall \underline{v} \in V(\Omega_0) \quad (7) $$

with initial condition

$$ \tilde{u}(0) = 0 \quad \partial_t \tilde{u}(0) = 0 $$

We can then distinguish three stiffness terms that affect the propagation of the wave:

* the material stiffness given by the 4th-order tensor $\tilde{\mathbf{C}}_0$, that accounts for the linearized material behavior around a loaded configuration

$$ \int_{\Omega_0} \underline{\varepsilon}_{\underline{x}_0} (\tilde{u}) : \tilde{\mathbf{C}}_0 : \underline{\varepsilon}_{\underline{x}_0} (\underline{v}) \, d\Omega_0 $$

* the geometric stiffness depending on the prestress with $\underline{\sigma}_{\underline{z}_0}$ resulting from the non-linear static problem,

$$ \int_{\Omega_0} \underline{\nabla}_{\underline{x}_0} \tilde{u} \underline{\sigma}_{\underline{z}_0} : \underline{\nabla}_{\underline{x}_0} \underline{v} \, d\Omega_0 $$

* the volumetric stiffness that corresponds to the Grad-div operator

$$ \int_{\Omega_0} \lambda_0^* \mathrm{div}_{\underline{x}_0} (\tilde{u}) \mathrm{div}_{\underline{x}_0} (\underline{v}) \, d\Omega_0 $$

where $\lambda_0^*$ only depends on the chosen energy density function $\mathcal{W}^v$ to penalize volumetric deformations and $\det \underline{F}(\underline{u}_0)$.

As we mentioned at the beginning, rewriting the formulation of the problem in the apparent configuration $\Omega_0$ is a fundamental step to isolate the term related to the volumetric constraint, which depends on the non-linear solution only through the parameter $\lambda_0^*$. Nevertheless, this term will be proportional to the bulk modulus, which in the nearly-incompressible regime is much larger than all other material parameters. In order to penalize volumetric deformations, it is then equivalent to substitute $\lambda_0^*$ with a constant $\lambda$ taken sufficiently large. We usually consider three orders of magnitude larger than the shear parameter. In addition, the formulation shows that the volumetric component of the propagating wave can be computed with a Grad-div operator only, which is particularly advantageous in terms of computational cost. These two latter remarks are key points for the numerical strategy presented in the rest of the paper, developed for overcoming limitations related to the nearly-incompressible constraint.

### 2.2.1 An example of strain energy function for nearly incompressible material
Let us now consider an example of penalization term for nearly-incompressible solids, as the one proposed in [30]

$$ \mathcal{W}^v(J) = k(J - 1) - k \log(J) $$

where $k \in \mathbb{R}^+$ is the bulk modulus taken very large with respect to other material parameters. The stress related to the volumetric deformation will be defined

$$ \underline{\Sigma}^v = -pJ\underline{\underline{C}}^{-1} \quad \text{with} \quad p = k \frac{1 - J}{J} $$
By linearization, we obtain (from (67) and (69), see Appendix A) that

$$p_0 = k \frac{1 - J_0}{J_0} \quad \text{and} \quad \lambda^* = k,$$

such that the shear term will be characterized by the following material stiffness

$$\underline{\underline{\varepsilon}}_{\underline{x}_0}(\tilde{u}) : \tilde{C}_0 : \underline{\underline{\varepsilon}}_{\underline{x}_0}(\underline{v}) := J_0^{-1} \underline{\underline{F}}_0^T \underline{\underline{\varepsilon}}_{\underline{x}_0}(\tilde{u}) \underline{\underline{F}}_0 : \mathbf{C}_0^d : \underline{\underline{F}}_0^T \underline{\underline{\varepsilon}}_{\underline{x}_0}(\underline{v}) \underline{\underline{F}}_0 + 2 p_0 \underline{\underline{\varepsilon}}_{\underline{x}_0}(\tilde{u}) : \underline{\underline{\varepsilon}}_{\underline{x}_0}(\underline{v})$$

with

$$\underline{\underline{F}}_0 = \underline{\underline{F}}(\underline{u}_0), \quad J_0 = \det \underline{\underline{F}}_0 \quad \text{and} \quad \mathbf{C}_0^d = 4 \frac{\partial^2 \mathcal{W}^d}{\partial \underline{\underline{C}}^2} \bigg|_{\underline{u}_0}.$$

Instead, the term related to the pressure propagation simply reads

$$\int_{\Omega_0} k \operatorname{div}_{\underline{x}_0}(\tilde{u}) \operatorname{div}_{\underline{x}_0}(\underline{v}) d\Omega_0.$$

We remark that the internal pressure $p_0$ could be a source of instability, since it is negative for non-linear problems characterized by volumetric increase $J_0 > 1$. Furthermore, it has been seen that wave propagation problems obtained by linearization may suffer from loss of coercivity, depending on the chosen elastic potential and the static solution $\underline{u}_0$. To avoid a possible buckling issue, a check of the smallest eigenvalues of the stiffness operator should be performed. Further discussions concerning these aspects can be found in [28].

For the purpose of our problem, we will strongly penalize any volumetric deformation due to nearly-incompressibility $J_0 \sim 1$, hence we can assume the total tangent stiffness term to be coercive.

## 3 A fully explicit numerical scheme

The scope of this section is to propose a numerical scheme for the computation of the elastic wave propagation problem in a nearly-incompressible linearized solid. We will present the scheme considering the solid at rest, in order to simplify the presentation and highlight the main aspects of the scheme. We remark that this assumption doesn't affect the applicability of the scheme, since the operator related to the enforcement of incompressibility (the grad-div operator) is not affected by the preload, beside the parameter $\lambda_0^*$ which – although it may vary in space in general – is homogeneous for the specific case presented Section 2.2.1. In practice, the parameter $\lambda_0^*$ acts as a penalization parameter and can be approximated by a constant value. Therefore, in what follows we assume that

$$\lambda^* = \lambda \quad \text{with} \quad \nabla \lambda = 0.$$

Given the linearity of the problem, it is equivalent to write equation with Lagrangian or Eulerian formalism, such that we will write the wave propagation problem in the configuration $\Omega$. Once we have introduced the linear continuous operators associated to our problem, we will introduce a strategy for space discretization focused on the accuracy and the efficiency of the numerical approximation. Finally, a fully explicit method for the time-discretization is given. The method is based on a classical LF scheme together with the use of Chebyshev polynomials (as proposed in [21]), for the relaxation of the CFL condition in the novel application of nearly-incompressible solids.

### 3.1 Linearization at rest condition

For the presentation of the numerical scheme, we consider the case of a hyperelastic solid, linearized around at rest condition, such that the displacement field (6) reduces to $\underline{u}(\underline{x}, t) = \tilde{u}(\underline{x}, t)$.

We rewrite the formulation (7) as a set of equations in $V'$, the dual space of the Hilbert space

$$V = V(\Omega) = \{\underline{v} \in H^1(\Omega)^d \mid \underline{v} = 0 \text{ on } \Gamma_D\}.$$
Note that the scalar product on $\mathcal{V}$ is the standard $H^1(\Omega)$ scalar product. Let us define the bilinear form $m: \mathcal{V} \times \mathcal{V} \to \mathbb{R}$ and related linear continuous operator $M: \mathcal{V} \to \mathcal{V}'$ as

$$ \forall (\underline{u}, \underline{v}) \in \mathcal{V} \times \mathcal{V} \quad m(\underline{u}, \underline{v}) := \langle M\underline{u}, \underline{v} \rangle_{\mathcal{V}', \mathcal{V}} := \int_{\Omega} \underline{\rho}\underline{u} \cdot \underline{v} \, d\Omega $$

defined in $L^2(\Omega)$ space. Notice that this bilinear form is symmetric and coercive in $L^2(\Omega)^3$.

The bilinear form $a: \mathcal{V} \times \mathcal{V} \to \mathbb{R}$ is rewritten as the sum of two bilinear forms, one contributing to the shear-waves propagation $a_s: \mathcal{V} \times \mathcal{V} \to \mathbb{R}$ and another one contributing to pressure wave propagation $a_p: \mathcal{V} \times \mathcal{V} \to \mathbb{R}$. To these bilinear forms we associate linear continuous operators $A^s$ and $A^p: \mathcal{V} \to \mathcal{V}'$ defined by, for all $(\underline{u}, \underline{v}) \in \mathcal{V} \times \mathcal{V}$,

$$ a_s(\underline{u}, \underline{v}) := \langle A^s \underline{u}, \underline{v} \rangle_{\mathcal{V}', \mathcal{V}} := \int_{\Omega} \underline{\varepsilon}(\underline{u}) : \mathbf{C}^d : \underline{\varepsilon}(\underline{v}) \, d\Omega, $$

$$ a_p(\underline{u}, \underline{v}) := \langle A^p \underline{u}, \underline{v} \rangle_{\mathcal{V}', \mathcal{V}} := \int_{\Omega} \mathrm{div}(\underline{u}) \mathrm{div}(\underline{v}) \, d\Omega, $$

where $\mathbf{C}^d$ is a symmetric 4th-order tensor (obtained from the linearization of the deviatoric component of the stress, see Appendix A). The bilinear form $a(\cdot, \cdot) = a_s(\cdot, \cdot) + a_p(\cdot, \cdot)$ is continuous and weakly coercive i.e., $\exists c, C \in \mathbb{R}^+$ such that

$$ c||v||_{\mathcal{V}}^2 \le ||v||_{L^2(\Omega)}^2 + |a(v, v)| \quad \text{and} \quad |a(u, v)| \le C||u||_{\mathcal{V}}||v||_{\mathcal{V}} $$

while the symmetry is easily verified by looking at definition of $a_s(\cdot, \cdot)$ and $a_p(\cdot, \cdot)$. Under these assumptions, the wave propagation problem of (7) linearized around at rest condition reads: For any $t \in [0, T]$, find $\underline{u}(t) \in \mathcal{V}$ such that

$$ M \partial_{tt} \underline{u}(t) + A^s \underline{u}(t) + \lambda A^p \underline{u}(t) = \underline{f}(t) \quad \text{in } \mathcal{V}' \quad \text{and} \quad \underline{u}(0) = \partial_t \underline{u}(0) = 0 \quad \text{in } \Omega. \quad (\text{NI}) $$

Existence results of a unique solution for this problem are given for instance in [31].

Let us now present an equivalent formulation with the introduction of the pressure $p := \lambda \mathrm{div} \, \underline{u}$ as a second variable of the problem. We denote with $\mathcal{L}$ the function space of the scalar pressure field

$$ \mathcal{L} := \{q \in L^2(\Omega)\}, $$

and we define the bilinear form $d: \mathcal{L} \times \mathcal{V} \to \mathbb{R}$, with the related divergence operator $D: \mathcal{V} \to \mathcal{L}$ and its adjoint gradient operator $D^*: \mathcal{L} \to \mathcal{V}'$, such that $\forall (\underline{u}, q) \in \mathcal{V} \times \mathcal{L}$

$$ d(q, \underline{u}) := \langle D^* q, \underline{u} \rangle_{\mathcal{V}', \mathcal{V}} := (q, D \underline{u})_{\mathcal{L}} := \int_{\Omega} q \mathrm{div}(\underline{u}) \, d\Omega. $$

Note that the scalar product on $\mathcal{L}$ is the standard $L^2(\Omega)$ scalar product. Then, the mixed formulation reads: For any $t \in [0, T]$, find $(\underline{u}(t), p(t)) \in \mathcal{V} \times \mathcal{L}$ such that

$$ \left\{ \begin{array}{ll} M \partial_{tt} \underline{u}(t) + A^s \underline{u}(t) + D^* p(t) = \underline{f}(t) & \text{in } \mathcal{V}' \\ D \underline{u}(t) = \lambda^{-1} M^p p(t) & \text{in } \mathcal{L} \end{array} \right. \quad \text{and} \quad \underline{u}(0) = \partial_t \underline{u}(0) = 0 \quad \text{in } \Omega \quad (\text{NIM}) $$

with $M^p$ is actually the identity map from $\mathcal{L} \to \mathcal{L}$. As we are dealing with a nearly-incompressible solid, it is interesting to notice that this latter formulation corresponds to a sort of penalization method for the resolution of the pure incompressible problem, where the pressure $p$ would be a Lagrange multiplier [5].

Comparing (NI) and (NIM) one can observe that

$$ A^p = D^*(M^p)^{-1}D. \quad (8) $$

This observation has a discrete counterpart as shown in the Section below.
## 3.2 Space discretization

### 3.2.1 Semi-discrete formulation

The classical method to solve finite-element approximation problem is known as displacement-based, which relies on the assumption that nodal displacements completely define the displacement of the body and are indeed the unknowns of the problem [2, 1]. This nodal displacements belong to finite-dimensional space $V_h \subset V$ that is a space of piece-wise smooth and globally continuous functions. The space $V_h$ is equipped with the scalar product

$$ (\underline{u}_h, \underline{v}_h)_H := \int_{\Omega} \underline{u}_h \cdot \underline{v}_h \, d\Omega. $$

In order to derive the semi-discretized problem, we define the linear operators $M_h : V_h \to V_h$, $A_h^s : V_h \to V_h$ and $A_h^p : V_h \to V_h$ such that $\forall (\underline{u}_h, \underline{v}_h) \in V_h \times V_h$

$$ (M_h \underline{u}_h, \underline{v}_h)_H = \oint_{\Omega} \varrho \underline{u}_h \cdot \underline{v}_h \, d\Omega \quad (A_h^s \underline{u}_h, \underline{v}_h)_H = \int_{\Omega} \mathbf{C}^d \underline{\varepsilon}(\underline{u}_h) : \underline{\varepsilon}(\underline{v}_h) \, d\Omega $$

$$ (A_h^p \underline{u}_h, \underline{v}_h)_H = \int_{\Omega} \operatorname{div} \underline{u}_h \operatorname{div} \underline{v}_h \, d\Omega $$

where the notation $\oint$ stands for the use of a quadrature formula for the computation of the integrals. The adequate quadrature formula to use depends on the choice of the finite element space $V_h$ and is discussed later. The semi-discretized problem of (NI) reads: For any $t \in [0, T]$, find $\underline{u}_h(t) \in V_h$ such that

$$ M_h \partial_{tt} \underline{u}_h(t) + A_h^s \underline{u}_h(t) + \lambda A_h^p \underline{u}_h(t) = \underline{f}_h(t) \quad \text{in } V_h \quad (\text{NId}) $$

with null initial conditions and $f_h$ an approximation of $f$ in $V_h$. However, it is well-known in finite element methods that the approximation of nearly-incompressible problems with a displacement based method suffers from several limitations, such as need of highly-refined mesh, numerical locking (especially with low-order elements), ill-conditioning of the stiffness matrix and spurious pressure [6, 2, 1].

The classical technique to overcome these limitations is the mixed method, which relies on the interpolation of the pressure field in a separate finite element space, denoted here as $\mathcal{L}_h \subset \mathcal{L}$ equipped with the scalar product $(\cdot, \cdot)_\mathcal{L}$. Let us then consider the space discretization of the equivalent continuous problem defined with the mixed formulation (NIM). We introduce the discrete divergence operator $D_h : V_h \to \mathcal{L}_h$ and its adjoint gradient operator $D_h^* : \mathcal{L}_h \to V_h$ such that $\forall (\underline{u}_h, p_h) \in V_h \times \mathcal{L}_h$

$$ (D_h \underline{u}_h, p_h)_\mathcal{L} = (\underline{u}_h, D_h^* p_h)_\mathcal{H} \approx \int_{\Omega} p_h \operatorname{div} \underline{u}_h \, d\Omega. \quad (9) $$

Note that we are very vague at this point on the actual definition of the operator $D_h$. A precise description is given later in the manuscript. The spatial discretization of problem (NIM) reads: For any time $t \in [0, T]$, find $(\underline{u}_h(t), p_h(t)) \in V_h \times \mathcal{L}_h$ such that

$$ \left\{ \begin{array}{ll} M_h \partial_{tt} \underline{u}_h(t) + A_h^s \underline{u}_h(t) + D_h^* p_h(t) = \underline{f}_h(t) & \text{in } V_h \\ D_h \underline{u}_h(t) = \lambda^{-1} M_h^p p_h(t) & \text{in } \mathcal{L}_h \end{array} \right. \quad (10) $$

with null initial data and where $M_h^p$ is the linear operator $\mathcal{L}_h \to \mathcal{L}_h$ given by

$$ \forall (p_h, q_h) \in \mathcal{L}_h \times \mathcal{L}_h, \quad (M_h^p p_h, q_h)_\mathcal{L} = \oint_{\Omega} p_h q_h \, d\Omega, $$

where – once again – a quadrature formula, adapted to the space $\mathcal{L}_h$ is used to compute the integral. By means of Schur complement method, we first compute the pressure field

$$ p_h = \lambda(M_h^p)^{-1} D_h \underline{u}_h, $$
and we obtain the semi-discretized scheme

$$M_h \partial_{tt} \underline{u}_h + A_h^s \underline{u}_h + \lambda \bar{A}_h^p \underline{u}_h = \underline{f}_h \quad \text{with} \quad \bar{A}_h^p := D_h^*(M_h^P)^{-1} D_h. \quad \text{(NIMd)}$$

First, observe that $\bar{A}_h^p$ is defined using a discrete counterpart of the formula (8) for defining $A^p$. However – in general – the operator $\bar{A}_h^p$ differs from $A_h^p$ unless some specific conditions are met. For instance, the two operators coincide when integral are computed exactly and div $V_h \subset \mathcal{L}_h$.

### 3.2.2 Spectral finite elements

Let us now assume our spatial domain $\Omega$ to be a geometry, such that we can construct a quasi-uniform regular partition $\mathcal{T}_h$ of $\Omega$, made of quadrangles ($d=2$) or hexahedra ($d=3$). Each deformed element $K \in \mathcal{T}_h$ is obtained with an invertible transformation $\phi_K$ applied to the reference element $\hat{K} = [0, 1]^d$.

This setting allows the use of high-order spectral finite element, for which we can apply the mass-lumping approach. The method [20] consists in obtaining a diagonal mass matrix – which is trivial to invert – by choosing a quadrature rule which has quadrature points at the locations of the interpolation points of a Lagrangian basis of $V_h$. The lagrangian basis is – as usual for finite elements – constructed locally on the reference element. On the reference element the basis is constructed using Gauss-Lobatto points. For the sake of simplicity, in the forthcoming numerical tests we will always consider a mass density $\rho = 1$. However, we could have also considered a smooth space-dependent density while still making use of mass lumping. The density would be evaluated on each quadrature point of each elements. Discontinuous density could even be considered if the discontinuities are captured by the elements and by a careful evaluation of the density at discontinuities. We refer the reader to [20, 31] for further details.

Furthermore, we recall that mixed variational problems are stable only for certain choices of finite element spaces that satisfy the inf-sup condition from the so-called Ladyzenskaya-Babuska-Brezzi theorem [10, 9, 6]. Unstable mixed methods are characterized by the appearance of nonphysical oscillations of the pressure field. The proof of the inf-sup condition is quite complex and out of the scope of this work. Hence, we rely on the analysis held in [32], in which a numerical analysis of convergences properties is presented for various choice of finite element spaces for $V_h$ and $\mathcal{L}_h$.

We consider the $Q_n - Q_{n-2}^{disc}$ method, proposed by Maday et al. [33]. More precisely

$$V_h = \{\varphi \in C^0(\Omega)^d \mid \varphi|_K \circ \phi_K \in Q_n(\hat{K})^d\}, \quad (11)$$

$$\mathcal{L}_h = \{\psi \in L^2(\Omega) \mid \psi|_K \circ \phi_K \in Q_{n-2}(\hat{K})\}. \quad (12)$$

This choice benefits of higher efficiency for the computation of the pressure field, given the lower order for the pressure field and the discontinuity of elements, for which the pressure contribution is computed locally. We remark that in 3D problems, the lower order of the shape functions for the pressure field results in the loss of one order of accuracy, that is comparable with the loss of one order already happening on non-cartesian mesh due to mass-lumping [34].

Note that, with the discussed choices, the operator $M_h$ and $M_h^P$ are computed using respectively a quadrature formula on the set of $(n+1)^d$ and $(n-1)^d$ Gauss-Lobatto points on the reference element.

**Remark.** A first alternative would be the choice $Q_n - Q_{n-1}$ for which $\mathcal{L}_h$ is chosen as

$$\mathcal{L}_h = \{\psi \in C^0(\Omega) \mid \psi|_K \circ \phi_K \in Q_{n-1}(\hat{K})\}. \quad (13)$$

This choice is analyzed in [35], which is an extension of the widely used Taylor-Hood elements. They are characterized by optimal convergence properties, but become more costly when computing local contributions. A second alternative is the choice $Q_n - P_{n-1}^{disc}$, for which

$$\mathcal{L}_h = \{\psi \in L^2(\Omega) \mid \psi|_K \circ \phi_K \in P_{n-1}(\hat{K})\}. \quad (14)$$
The choice has been discarded since it does not allow important optimization procedure related to the use of finite elements constructed on $\hat{K} = [0, 1]^d$ by cartesian product (such as the spectral finite elements).

### 3.2.3 Definition of divergence operator

We now detail the definition of the discrete operator $D_h$ and its adjoint $D_h^*$, developed to increase the efficiency for the computation of the volumetric contribution. When considering the computation of the pressure term at the discrete level, we can write

$$ 
\begin{aligned}
\int_{\Omega} q_h \operatorname{div} \underline{w}_h \, d\Omega &= \sum_{K \in \mathcal{T}_h} \int_K q_h \operatorname{div} \underline{w}_h \, dK \\
&= \sum_{K \in \mathcal{T}_h} \int_{\hat{K}} \hat{q}_h \operatorname{div}(J_K \underline{\nabla} \phi_K^{-1} \hat{w}_h) \, d\hat{K} \\
&= \sum_{K \in \mathcal{T}_h} \int_{\hat{K}} \hat{q}_h \operatorname{cov}(\underline{\nabla} \phi_K) : \underline{\nabla} \hat{w}_h \, d\hat{K}
\end{aligned}
\quad (15) $$

where

$$ J_K = \det \underline{\nabla} \phi_K \quad \text{and} \quad \operatorname{cov}(\underline{\nabla} \phi_K) = J_K \underline{\nabla} \phi_K^{-T}. $$

Then, we have the displacement $\hat{w}_h$ and pressure $\hat{q}_h$ w.r.t the reference element $\hat{K}$. The computation performed on the element $\hat{K}$ is a classical procedure in finite-element analysis as it allows simplifying and reducing computations. The last line of (15) can be rewritten with the application of the *Green's formula*, and we obtain

$$ \int_{\Omega} q_h \operatorname{div} \underline{w}_h \, d\Omega = - \sum_{K \in \mathcal{T}_h} \int_{\hat{K}} \operatorname{cov}(\underline{\nabla} \phi_K) \underline{\nabla} \hat{q}_h \cdot \hat{w}_h \, d\hat{\Omega}_K + \sum_{K \in \mathcal{T}_h} \int_{\partial \hat{K}} \hat{q}_h \hat{w}_h \cdot \operatorname{cov}(\underline{\nabla} \phi_K) \hat{n} \, d\hat{s} \quad (16) $$

On this reformulation we base the definition of discrete operators $D_h$ and $D_h^*$, given by the application of a quadrature formula (Gauss Lobatto), to compute the integrals, so that we have

$$ 
\begin{aligned}
(D_h \underline{w}_h, q_h)_\mathcal{L} &:= (\underline{w}_h, D_h^* q_h)_\mathcal{H} \\
&:= - \sum_{K \in \mathcal{T}_h} \oint_{\hat{K}} \operatorname{cov}(\underline{\nabla} \phi_K) \underline{\nabla} \hat{q}_h \cdot \hat{w}_h \, d\hat{\Omega}_K + \sum_{K \in \mathcal{T}_h} \oint_{\partial \hat{K}} \hat{q}_h \hat{w}_h \cdot \operatorname{cov}(\underline{\nabla} \phi_K) \hat{n} \, d\hat{s} \quad (17)
\end{aligned}
$$

For the computation of the volumetric term we use a quadrature based on $(n+1)^d$ Gauss-Lobatto points. Moreover, we use a quadrature based on $(n+1)^{d-1}$ Gauss-Lobatto points on each facet of the reference element (4 when $d=2$, 6 when $d=3$).

An alternative would be to introduce $\tilde{D}_h$ and $\tilde{D}_h^*$ by

$$ (\tilde{D}_h \underline{w}_h, q_h)_\mathcal{L} := (\underline{w}_h, \tilde{D}_h^* q_h)_\mathcal{H} := \sum_{K \in \mathcal{T}_h} \oint_{\hat{K}} \hat{q}_h \operatorname{cov}(\underline{\nabla} \phi_K) : \underline{\nabla} \hat{w}_h \, d\hat{K} \quad (18) $$

using a quadrature based on $(n+1)^d$ Gauss-Lobatto points.

We now compare the computational cost associated to the application of $\tilde{D}_h$ or $D_h$. The application of the operator $\tilde{D}_h$ implies the computation of 9 derivatives for the gradient of the displacement field, plus the interpolation of the pressure field into $Q_n$. The use of the proposed operator $D_h$, instead, requires the computation of 3 derivatives of the $Q_{n-2}$ pressure field, 3 interpolations for the resulting gradient, and the computation of the boundary term that has a lower complexity. Without going into details our computations show that asymptotically, using $D_h$ over $\tilde{D}_h$ offers a reduction of 40% in terms of computational cost. Nevertheless, higher efficiency is expected for some specific choice of $n$ as we explain in the following section.
### 3.2.4 The $Q_4 - Q_2^{disc}$ element

In practice to further increase the efficiency of the numerical scheme, we use for the order of $V_h$ and $L_h$ the couple spaces such that the DoFs of the pressure element correspond to a subset of the DoFs of the displacement element. By taking $n = 4$, we have the highest element order with the mentioned feature, reducing the number of operations for the interpolation from $L_h$ to $V_h$.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Representation of spectral elements choice with related DoFs distribution. From the left, $Q_4(\hat{K})$ element of displacement space $V_h$ and $Q_2^{disc}(\hat{K})$ element for pressure space $L_h$.</figcaption></figure>

The computation is performed within a software (OndoMathX) developed within the research group, and the interpolation of the pressure field is specifically optimized in order to further speed up the computations. In our implementation, no matrix is assembled since all the computations are done locally, which allow to strongly reduce the CPU storage and cost. We measure numerically a speed-up of 70% for the computation of the pressure contribution compared to the use of (18).

## 3.3 Time discretization

### 3.3.1 The leap-frog method

We now develop a fully-discretized scheme associated to the mixed formulation (NIMd). Given the time interval $[0, T]$, with $T > 0$, we define the partition $t^n = n\Delta t$, with $n \in \{0, 1, \dots, N\}$ and the time-step defined as $\Delta t = T/N$, we can introduce a time-discretization based on Leap-Frog method

$$M_h \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + (A_h^s + \lambda \bar{A}_h^p) u_h^n = \bar{f}_h^n \quad \text{in } V_h \quad (LF)$$

The LF method is characterized by second order accuracy and energy conservation. It is fully explicit, which makes it useful to deal with large scale problem such as elastodynamic problems. However, the scheme is stable under the well-known CFL condition [31]

$$0 \le \Delta t^2 A_h^s + \Delta t^2 \lambda \bar{A}_h^p \le 4M_h \quad \Leftrightarrow \quad \Delta t^2 \le \frac{4}{\rho(M_h^{-1} A_h^s + \lambda M_h^{-1} \bar{A}_h^p)}, \quad (19)$$

where here $\rho(\cdot)$ now denotes the spectral radius of an operator. The enforcement of nearly-incompressibility corresponds to a large value of the bulk modulus $\lambda$. As a result, the maximal time-step to guarantee the stability of the scheme is strongly reduced, the number of iterations $N$ to reach the end of the simulation at time $T$ will proportionally increase.

For a physical interpretation, we could consider that with a large bulk modulus, the speed of pressure waves is much larger than shear waves (respectively $1500\,ms^{-1}$ and $10\,ms^{-1}$ in biological tissues) and our explicit scheme should have a time step adapted to the fastest wave propagation phenomena. Then, the pressure wave strongly limits the time-step and affects the efficiency of the computations.

### 3.3.2 Stabilized leapfrog method

As mentioned, the pressure operator $\lambda \bar{A}_h^p$ limits the time-step $\Delta t$, but at the same time it has a low computational cost, thanks to the proposed operators $D_h$ and $D_h^*$. On the other side, the shear operator $A^s$ has a lower spectral radius, around two order of magnitudes smaller
compared to $\lambda \bar{A}_h^p$ in soft tissues, but it has a higher computational cost as it involves the computation of the gradient of displacements and accounts for all the complexities of the material (typically heterogeneity and anisotropy). In this framework, we would like to separate the computation of shear and pressure wave operator. The shear wave operator is computed at a "slow scale", i.e., with a larger time step, whereas pressure wave operator is computed as if a time-step was defined at a "fast-scale". This general idea is the basis of the Chebyshev Leapfrog method [21], such that at each time-step the shear part is iterated only once and the volumetric operator is iterated $m + 1$ times. We remark that this double scale time-step is just to give an interpretation of the application of the polynomials in the form of an iterative approach, although we are not actually defining a time-step associated to the volumetric component. To our knowledge, this is the first time that this approach has been applied in the context of incompressible elastodynamics.

In order to use the strategy suggested above, the pressure operator is approximated with a polynomial $\mathcal{P}_m(x)$ of order $m$, given the condition $\mathcal{P}_m(0) = 1$ to be consistent with the leap-frog scheme,

$$M_h \frac{u_h^{n+1} - 2u_h^n + \underline{u}_h^{n-1}}{\Delta t^2} + A_h^s \underline{u}_h^n + \mathcal{P}_m(\Delta t^2 \lambda \bar{A}_h^p M_h^{-1}) \lambda \bar{A}_h^p \underline{u}_h^n = \underline{f}_h^n. \quad (20)$$

We rewrite the stability condition (19) for our modified scheme

$$0 \le \Delta t^2 A_h^s + \Delta t^2 \mathcal{P}_m(\Delta t^2 \lambda \bar{A}_h^p M_h^{-1}) \lambda \bar{A}_h^p \le 4M_h \quad (21)$$

from which we derive the following theorem for a polynomial $\mathcal{P}_m$ yet unspecified

**Theorem 3.1.** A sufficient condition for the stability of the scheme (20) is

$$\Delta t^2 \le \frac{1}{\rho(M_h^{-1} A_h^s)} \quad (22)$$

and

$$\mathcal{P}_m(x) x \in [0, 3] \quad \forall x \in [0, \rho(\Delta t^2 \lambda M_h^{-1} \bar{A}_h^p)] \quad (23)$$

*Proof.* From the stability condition of the modified scheme (21), we derive a sufficient condition

$$\left\{ \begin{array}{l} \Delta t^2 A_h^s \le M_h \\ 0 \le \Delta t^2 \mathcal{P}_m(\Delta t^2 \lambda \bar{A}_h^p M_h^{-1}) \lambda \bar{A}_h^p \le 3M_h. \end{array} \right. \quad (24)$$

The first condition (22) is then trivially obtained. Given that the pressure operator is self-adjoint and positive, we can then apply the Spectral Mapping theorem [36] to reduce the analysis to a spectral analysis. More precisely, for condition (23), we see that the second equation of (24) can be written in terms of the spectrum $\sigma$

$$\sigma \left( \Delta t^2 M_h^{-1} \mathcal{P}_m (\Delta t^2 \lambda \bar{A}_h^p M_h^{-1}) \lambda \bar{A}_h^p \right) \in [0, 3]$$

which is equivalent to

$$\sigma \left( \mathcal{P}_m (\Delta t^2 \lambda M_h^{-1} \bar{A}_h^p) \Delta t^2 \lambda M_h^{-1} \bar{A}_h^p \right) \in [0, 3]$$

Hence, the spectrum is in $[0, 3]$ if the polynomial $\mathcal{P}_m$ is such that

$$\mathcal{P}_m(x) x \in [0, 3] \quad \forall x \in [0, \rho(\Delta t^2 \lambda M_h^{-1} \bar{A}_h^p)] \quad (25)$$

concluding the proof.

<figure><img src="image_3.png" /></figure>

Notice that the implementation of the polynomial $\mathcal{P}_m$ should be obtained with an algorithm, for example like the Horner's rule or recurrence rules, so that the computation of the volumetric contribution of the displacement would correspond to the application of the argument of the polynomial for $m + 1$ times.
### 3.3.3 Leapfrog Chebyshev Methods

First note that – following (22) – the time step $\Delta t$ is chosen independently of the choice of $\mathcal{P}_m$. We set

$$ \Delta t = (1 - \epsilon) \frac{1}{\sqrt{\rho(M_h^{-1} A_h^s)}} \quad \text{with} \quad \epsilon \ll 1. \qquad (26) $$

For the choice of the polynomial, we use the so-called stabilized or damped Chebyshev Polynomials $\mathcal{P}_m^\nu$, proposed by [21], and obtained from scaling and shifting of Chebyshev polynomials of first kind $\mathcal{T}_m(x)$:

$$ \mathcal{P}_m^\nu(x) = \frac{2}{x} \left[ 1 - \frac{\mathcal{T}_{m+1}\left(\delta_m^\nu - \frac{x}{\omega_m^\nu}\right)}{\mathcal{T}_{m+1}\left(\delta_m^\nu\right)} \right] \qquad (27) $$

with parameters $\delta_m^\nu$ and $\omega_m^\nu$

$$ \delta_m^\nu = 1 + \frac{\nu}{(m+1)^2}, \quad \omega_m^\nu = 2 \frac{\mathcal{T}'_{m+1}(\delta_m^\nu)}{\mathcal{T}_m(\delta_m^\nu)}, $$

and first order derivative $\mathcal{T}'_{m+1}(x) = (m+1)\mathcal{U}_m(x)$ where $\mathcal{U}_m$ is the Chebyshev of second kind. We see that consistency with respect to the original scheme is guaranteed, since

$$ \mathcal{P}_m^\nu(x) = 2 \left[ 1 - \frac{\mathcal{T}_{m+1}\left(\delta_m^\nu - \frac{x}{\omega_m^\nu}\right)}{\mathcal{T}_{m+1}\left(\delta_m^\nu\right)} \right] = x + \mathcal{O}(x^2) \qquad (28) $$

and the Leap-Frog-Chebyshev scheme reads

$$ M_h \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_h^s u_h^n + \mathcal{P}_m^\nu (\Delta t^2 \lambda \bar{A}_h^p M_h^{-1}) \lambda \bar{A}_h^p u_h^n = f_h^n \qquad (\text{LFC}) $$

Chebyshev polynomials are sequences of polynomial functions that can be computed with recurrence relations and benefit of optimal properties for the uniform distribution of its roots within the interval $[-1, 1]$. The formulation (27) allows to change the interval in which Chebyshev polynomials oscillates, with the following properties derived in [21, 22].

**Property 3.1.** For $\nu \in [0, 1]$:

$$ \sup_{x \in [0, \alpha_m]} |x \mathcal{P}_m^\nu(x)| < 4 - \frac{2\nu}{1 + \nu} \quad \text{with} \quad \alpha_m = \left(2 + \frac{\nu}{(m+1)^2}\right) \omega_m^\nu \qquad (29) $$

and $\forall \beta \le \alpha_m$

$$ \inf_{x \in [\beta, \alpha_m]} |x \mathcal{P}_m^\nu(x)| > \min \left( \frac{\beta}{1 + \nu}, \frac{2\nu}{1 + \nu} \right) \qquad (30) $$

This choice of polynomials allows to have the largest stability region (in $x$) in which the Chebyshev $x\mathcal{P}_m^\nu(x)$ is comprised between the inferior and the superior limit, as given by Property 3.1.

We now present the constructed solution in order to respect the sufficient condition of Theorem 3.1. There are two parameters to choose: the parameter $\nu$ that defines the inferior and superior limit of $\mathcal{P}_m^\nu$, and the polynomial order $m$ that depends on the spectral radius of the pressure operator through the parameter $\alpha_m$.

By taking $\nu = 1$, we define the superior limit of the polynomial within the stability region $[0, \alpha_m]$.

$$ \sup_{x \in [0, \alpha_m]} |x \mathcal{P}_m^1(x)| < 3 \qquad (31) $$

which partially corresponds to the second condition (23) of the stability condition of theorem 3.1.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** On the left: Chebyshev Polynomials of first kind with order m from 0 to 4. On the right: Stabilized Chebyshev polynomials $x\mathcal{P}_m^\nu$ with order m from 1 to 4. Dashed vertical lines correspond to $\alpha_m$ which guarantees that for $x \in [0, \alpha_m]$ we have $x\mathcal{P}_m^\nu < 4 - \frac{2\nu}{1+\nu}$. In the figure, $\nu = 1$.</figcaption></figure>

Still, we remark that the stabilized Chebyshev polynomial is characterized by a superior limit only within the stability region $x \in [0, \alpha_m]$, for which we have to ensure that

$$ \sigma(\Delta t^2 \lambda M_h^{-1} \bar{A}_h^p) \in [0, \alpha_m] \qquad (32) $$

in order to guarantee condition (23). Since the positivity of the spectrum is given, we only need to respect

$$ \rho(\Delta t^2 \lambda M_h^{-1} \bar{A}_h^p) \le \alpha_m \qquad (33) $$

From an estimation of $\alpha_m$, we see that for $\nu > 0.7$ we have

$$ 4(m+1)^2 e^{-\nu/2} \le \alpha_m $$

By substitution, we derive the minimum order of the damped Chebyshev polynomial

$$ m+1 \ge \sqrt{\lambda} \frac{\sqrt{\rho(\Delta t^2 M_h^{-1} \bar{A}_h^p)} e^{1/4}}{2} . \qquad (34) $$

To sum up, the stability of the scheme (LFC) is given by the following proposition

**Proposition 3.2.** A sufficient condition to guarantee the stability of the Leapfrog Chebyshev scheme (LFC) is given by

Maximum time step

$$ \Delta t \le (1 - \epsilon) \frac{1}{\sqrt{\rho(M_h^{-1} A_h^s)}} \quad \text{with } \epsilon \ll 1 \qquad (35) $$

Minimum order of polynomial $x\mathcal{P}_m^1(x)$

$$ m+1 = \left\lceil \sqrt{\lambda} \frac{\sqrt{\rho(\Delta t^2 M_h^{-1} \bar{A}_h^p)} e^{1/4}}{2} \right\rceil \qquad (36) $$

Note that in Proposition 3.2 the operator $[\cdot]$ outputs the smallest larger integer. With the introduction of polynomial $\mathcal{R}_{m+1}(x) := x\mathcal{P}_m^1(x)$, the fully-discretized scheme reads

$$ \underline{u}_h^{n+1} - 2\underline{u}_h^n + \underline{u}_h^n + \Delta t^2 M_h^{-1} A_h^s \underline{u}_h^n + \mathcal{R}_{m+1}(\Delta t^2 \lambda M_h^{-1} \bar{A}_h^p) \underline{u}_h^n = \Delta t^2 M_h^{-1} \underline{f}^n \qquad (\text{LFC}) $$

In conclusion, we have developed a stable explicit scheme for linear elastic wave propagation in incompressible media. The stability of the scheme is based on the definition of a time step
$\Delta t$ from the shear operator only and on the order of a stabilized Chebyshev polynomials. One way to interpret this is to say that the large time step $\Delta t$ is chosen using the shear operator only, while the pressure terms are computed with a small time-step $\Delta t_m = \Delta t/(m + 1)$. As we will see in the Results section, the efficiency of the scheme w.r.t. the LF method relies on two points: First, the low computational cost of the pressure operator compared to the shear operator ; Second, the value of $\lambda$, the higher $\lambda$ is, higher will be the potential gain given by this strategy.

The implemented algorithm for computing the solution to (LFC) is given Algorithm 1. It is based on the recurrence relations of the Chebyshev polynomials. An auxiliary variable $\underline{z}_l^n$ is introduced, with $l = 0, ..., m + 1$, which corresponds to the contribution of the pressure operator at the fast scale within each time step $\Delta t^n$.

<caption>**Algorithm 1** Algorithm for the (LFC) scheme</caption>

$\underline{z}_0^n \leftarrow \underline{u}_h^n$
$\underline{z}_1^n \leftarrow \delta_m^1 \underline{u}_h^n - \frac{\Delta t^2 \lambda}{\omega_m^1} M_h^{-1} \bar{A}_h^p \underline{u}_h^n$
**for** $l = 2 : m + 1$ **do**
&nbsp;&nbsp;&nbsp;&nbsp; $\underline{z}_l^n \leftarrow 2 \left( \delta_m^1 I - \frac{\lambda \Delta t^2}{\omega_m^1} M_h^{-1} \bar{A}_h^p \right) \underline{z}_{l-1}^n - \underline{z}_{l-2}$
**end for**
$\mathcal{R}_{m+1}(\Delta t^2 \lambda M^{-1} \bar{A}_h^p) \underline{u}_h^n \leftarrow 2 \underline{u}_h^n - \frac{2 \underline{z}_{m+1}^n}{\mathcal{T}_{m+1}(\delta_m^1)}$
$\underline{u}_h^{n+1} \leftarrow \Delta t^2 M_h^{-1} \underline{f}_n - \underline{u}_h^{n-1} + 2 \underline{u}_h^n - \Delta t^2 M_h^{-1} A_h^s \underline{u}_h^n - \mathcal{R}_{m+1}(\Delta t^2 \lambda M^{-1} \bar{A}_h^p) \underline{u}_h^n$

Let us add few remarks. As for the (LF), also for (LFC) we must compute the inverse of the mass matrix at each time step $\Delta t$. In elastodynamics, the problems' dimension (in our analysis, we are often in the order of few millions of DoF) requires this operation to be the most efficient as possible, from which the importance of mass-lumping to obtain a diagonal mass matrix, that is trivial to invert. The last remark is related to the application of the scheme to non-linear mechanics. As we have shown in Section 2 and detailed in Appendix A, the analysis of small amplitude elastic wave propagation in a non-linear solid undergoing static large deformation does not affect the applicability of the scheme. The key-point is to rewrite the weak formulation in the so-called *apparent* configuration $\Omega_0$, such that we can isolate the component associated to volumetric deformation and apply the (LFC). In other words, only the shear operator $A^s$ changes with respect to the condition at rest.

## 3.4 Efficiency estimation the LFC scheme

Let us now present a theoretical estimation of the computational gain that we obtain with the (LFC) scheme w.r.t. to the classical Leap-Frog scheme (LF).

We denote $C_S$ the computational time for the application of the shear operator, and $C_P$ the computational time for the application of the pressure operator. The objective is to compare the simulation time of the LF method and of the LFC method, based on $C_S$ and $C_P$ and on the time-steps $\Delta t_{LF}$ and $\Delta t_{LFC}$ derived from the respective stability condition. The quantity that we would like to analyze is the computational gain $G$ defined as follows

$$G = \frac{N_{LF}C_{LF}}{N_{CH}C_{LFC}}$$

where $N_{LF}$ and $N_{LFC}$ are the number of iterations to complete a simulation, while $C_{LF}$ and $C_{LFC}$ are the computational cost of the time step for each method.
We start from the classic LF scheme considering that, with the property that the shear and pressure operators are positive and self-adjoint, we have

$$ \rho(M_h^{-1} A_h^s + \lambda M_h^{-1} \bar{A}_h^p) \le \rho(M_h^{-1} A_h^s) + \rho(\lambda M_h^{-1} \bar{A}_h^p) \quad (37) $$

Then, the CFL condition (19) is respected for the following time-step

$$ \Delta t_{LF}^2 = \frac{4}{\rho(M_h^{-1} A_h^s) + \rho(\lambda M_h^{-1} \bar{A}_h^p)} = \rho(M_h^{-1} A_h^s) \frac{4}{1 + \lambda R} \quad R := \frac{\rho(M_h^{-1} \bar{A}_h^p)}{\rho(M_h^{-1} A_h^s)}. \quad (38) $$

The choice of time step given by (38) may appear suboptimal because (37) is in general a strict inequality. However, for large $\lambda$ the difference between, $\rho(M_h^{-1} A_h^s + \lambda M_h^{-1} \bar{A}_h^p)$ and the sum of the spectral radii is small with respect to $\lambda$. The computational time of a single time step in the LF scheme is trivially given by

$$ C_{LF} = C_S + C_P. \quad (39) $$

For the LFC scheme, we have constructed the stability condition by choosing

$$ \Delta t_{LFC}^2 = \frac{1}{\rho(M_h^{-1} A_h^s)}. \quad (40) $$

Notice that with respect to the time-step defined in (26), we forget the term $(1 - \epsilon)$ for the sake of simplicity. For the evaluation of the computational cost of a time-step with the LFC method, we have

$$ C_{LFC} = C_S + (m + 1)C_P $$

where we take $m + 1$ as the minimal value to guarantee the stability condition (34), and we combine with (40) in order to obtain

$$ C_{LFC} = C_S + \left[ \Delta t_{LFC} \frac{\sqrt{\lambda \rho(M_h^{-1} \bar{A}_h^p)} e^{1/4}}{2} \right] C_P = C_S + \left[ \frac{e^{1/4}}{2} \sqrt{\lambda R} C_P \right]. $$

Now that we have all the elements, we rewrite the computational gain $G$ as

$$ G = \frac{\frac{1}{\Delta t_{LF}} C_{LF}}{\frac{1}{\Delta t_{LFC}} C_{LFC}} = \frac{\Delta t_{LFC} C_{LF}}{\Delta t_{LF} C_{LFC}} $$

We notice that we can write the ratio between $\Delta t_{LFC}$ and $\Delta t_{LF}$ as a function of the quantity $R$ from relations (38) and (40)

$$ \frac{\Delta t_{LFC}}{\Delta t_{LF}} = \frac{\sqrt{1 + \lambda R}}{2} $$

Finally, the computational gain $G$ is given by the following relation

$$ G = \frac{\sqrt{1 + \lambda R}}{2} \frac{C_S + C_P}{C_S + \left[ \frac{e^{1/4}}{2} \sqrt{\lambda R} C_P \right]} \quad (41) $$

Considering the incompressible limit ($\lambda \to \infty$), we have

$$ \lim_{\lambda \to \infty} G = \frac{\sqrt{\lambda R}}{2} \frac{C_S + C_P}{\frac{e^{1/4}}{2} \sqrt{\lambda R} C_P} = e^{-1/4} \left( \frac{C_S}{C_P} + 1 \right) $$

for which we obtain that the LFC scheme is more efficient than the LF scheme (i.e., $G > 1$) if the computational cost of the pressure operator is sufficiently smaller than the shear operator, in particular

$$ C_S/C_P > e^{1/4} - 1 \approx 0.2840. \quad (42) $$
In practice, this condition is always met, for which we show that at the incompressible limit the proposed (LFC) scheme is more efficient than the (LF).

We point out that the strategy developed for the space discretization Section 3.2 has an impact also on the efficiency estimation that we have propose by G. Indeed, the optimization for the computation of the pressure contribution corresponds to a decrease of the computational cost for the pressure operator $C_P$, i.e., a higher efficiency of the scheme through relation (42), given that $C_S$ is the same. In addition, notice that in the case of hyperelastic solids linearized around a loaded condition, the computational cost of the shear operator $C_S$ will probably increase as it includes both material and geometric stiffnesses, while $C_P$ remains constant. Considering our analysis, this increases the efficiency of the (LFC) scheme compared to the (LF) method.

**Remark.** The principal motivation of this work is the treatment of elastodynamic problems with nearly-incompressible solids. Nevertheless, we should notice that the applicability of the time-discretization method based on Chebyshev polynomials can be extended to other double-scale problems in time, where the physical phenomenon at fast scale is relevant only for its impact on the slow-scale one. At the same time, the differential operator associated to the slow-scale phenomenon should be more expensive to compute, or equivalently, the operator at fast-scale is sufficiently cheaper to compute. This evidences our efforts, as mentioned above, in the space discretization strategy to optimize the computation of the volumetric component of the solution, with the definition of the divergence operator (Sec. 3.2.3) and the choice for the FE spaces (Sec. 3.2.4).

# 4 Numerical Results

This section is devoted to the presentation of numerical results in order to validate the properties of the proposed scheme (LFC). Before presenting details of the simulations, we consider the non-dimensionalization of the wave propagation problem. As in [5], we introduce the characteristic length L of the domain $\Omega$, the characteristic observation time T and the shear parameter $\mu$. The non-dimensionalized displacement and pressure field reads respectively:

$$ \underline{u}_a(t_a, \underline{x}_a) := L^{-1}\underline{u}(Tt_a, L\underline{x}_a) \quad \text{and} \quad p_c(t_a, \underline{x}_a) := \mu^{-1}p(Tt_a, L\underline{x}_a) $$

with $t_a \in [0, 1]$. We then derive the non-dimensionalized quantities

$$ \underline{\varrho}_a = \frac{L^2 \underline{\varrho}}{T^2 \mu}, \quad \mathbf{C}_a^d = \mu^{-1} \mathbf{C}^d, \quad \lambda_a = \frac{\lambda}{\mu}, \quad \underline{f}_a = \frac{L}{\mu} \underline{f} $$

and with these, we rewrite our wave propagation problem in a non-dimensional form. For the rest of the paper, we drop the notation $(\cdot)_a$ that represents non-dimensionalized quantities.

In the following results, the source term used to generate the wave propagation has always the same form. It is the product of a spatial and a temporal term. The spatial excitation is a $C^\infty$ function with a compact support, localized in the center of our solid $\underline{x}_C$. For the temporal push, we have a Ricker wavelet with a small offset $t_0$. The source term reads:

$$ f(\underline{x}, t) = \begin{cases} A_m \cdot \exp\left[\frac{1}{2\sigma_s}\left(\frac{r_0^2}{d^2 - r_0^2} + 1\right)\right] \cdot \left[\left(\frac{t-t_0}{\sigma_\tau}\right)^2 - \frac{1}{\sigma_\tau}\right] \exp\left[-\frac{(t-t_0)^2}{\sigma_\tau}\right] & \text{if } d^2 < r_0^2 \\ 0 & \text{otherwise} \end{cases} \quad (43) $$

with $d^2 = \sum_{i=1}^{\text{dim}} (x_i - x_{C,i})^2$ with d the distance from the center, $r_0$ the maximum width of the source, $A_m$ the maximum amplitude, $\sigma_\tau$ and $\sigma_s$ the covariance.

## 4.1 Non-locking mixed formulation

We compare simulation results of the nearly-incompressible wave propagation problem considering the two semi-discretized schemes (NId) and (NIMd) respectively obtained with the displacement-based and mixed method. At this point we do not use the strategy developed with Chebyshev polynomials, but we use the classic LF scheme for the time discretization.
Let us now provide numerical and physical parameters of simulations. We consider a 2D circular domain with non-dimensionalized density $\rho = 1$. We consider isotropic homogeneous behavior, characterized by the constitutive relation (non-dimensionalized)

$$ \underline{\underline{\sigma}} = 2\underline{\underline{\varepsilon}}(\underline{u}) + \lambda \operatorname{tr}(\underline{\underline{\varepsilon}}(\underline{u})), \qquad (44) $$

with the non-dimensionalized second Lamé parameters as $\lambda = 2000$. The function described in (43) is used with $\underline{x}_C = (0, 0)$. A vector field is generated by assuming that the source is polarized in the y-direction: $\underline{f} = f\underline{e}_y$. The other chosen parameters of the source are: $A_m = 10^6$, $r_0 = 0.05$, $\sigma_s = 0.5$, $\sigma_\tau = 2.5 \cdot 10^{-4}$.

For the space discretization, we compare the two schemes by using $Q_4$ elements for the displacement-based method (NId), while for the mixed method we discretize with $Q_4 - Q_2^{disc}$ elements, as presented in (11) and (12). In both cases the mesh is highly refined, with 150k Dofs.

In figure 4 we have the magnitude of the displacement field obtained from the mixed method, in which we can appreciate only the propagation of the shear wave. In an isotropic homogeneous material the velocity of propagation of shear wave ($v_s$) and of pressure wave ($v_p$) are related to the material parameters such that

$$ v_s^2 = \frac{\mu}{\rho}, \qquad v_p^2 = \frac{\lambda + 2\mu}{\rho} $$

With the present choice of material parameters, we have that $v_p \approx 45 v_s$. Then, when looking at the resulting displacement field, we only see a propagating shear wave and not the pressure one (in figure 4), as the pressure is too fast for the present timescale. We show the displacement results only of the mixed method (NIMd), as no qualitative difference is present with respect to the solution of the displacement-based method.

Then, we extract the volumetric component from the solution (which corresponds to $\lambda\nabla\operatorname{div}\underline{u}$) a posteriori, by computing the quantity

$$ -\frac{\underline{u}_h^{n+1} - 2\underline{u}_h^n + \underline{u}_h^{n-1}}{\Delta t^2} - M_h^{-1} A_h^s \underline{u}_h^n $$

for both displacement-based and mixed method. The limitations of the displacement-based method previously mentioned is visible. In figure 5 we observe nonphysical oscillations in the pressure field that propagates with the shear wave, which is a numerical effect due to enforcement of incompressibility. This phenomenon disappears with the implemented mixed method, as we can see from figure 6.

<figure><img src="image_5.png" /><figcaption>Figure 4. Elastic wave propagation in homogeneous isotropic solid. Magnitude of displacement field obtained with the mixed scheme. The 2D circular mesh is discretized with 150k Dofs.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Volumetric component of the elastic wave propagation in homogeneous isotropic solid, with incompressible parameter $\lambda = 2000$. The 2D circular mesh is discretized with 150k DoFs. Solution obtained with displacement based method (NId), characterized by the presence of nonphysical oscillations coupled with the shear wave</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Volumetric component of the elastic wave propagation in homogeneous isotropic solid, with incompressible parameter $\lambda = 2000$. The 2D circular mesh is discretized with 150k Dofs. Solution obtained with mixed method (NIMd). In this case, no spurious pressure wave is observed.</figcaption></figure>

One important remark is that in this simulation we are far from the nearly-incompressible case. In biological tissues, the pressure wave has a velocity that is at least 2 orders of magnitude higher than the shear wave velocity. Yet, in our case we can already see the importance of approximating this type of problems with an appropriate mixed formulation. Nevertheless, these are well-known results when dealing with nearly-incompressible problems [2, 37]. The primary objective of the space discretization strategy was to compute efficiently the pressure contribution, ensuring the stability of the method and the accuracy of the computation for both displacements and pressure.

## 4.2 Leapfrog Chebyshev Methods

### 4.2.1 Efficiency

We perform simulations of elastic wave propagation in isotropic homogeneous solids within different domains: square, circle and cube. For the material, we take $\rho = 1$, $\lambda = 1000$ and the same parameters of the source that we have used in the previous section. In Table 1 we present the associated numerical results. We measure the simulation time $\eta_{LF}$ and $\eta_{LFC}$ with the LF and LFC scheme. The average computational cost of $C_P$ and $C_S$ is measured over 1000 iterations, and the ratio $C_S/C_P$ can be then compared to the measured computational gain $G_{eff} = \eta_{LF}/\eta_{LFC}$. We underline that the effective computational gain $G_{eff}$ is related to simulations with $\lambda = 1000$, while the estimated $G$ corresponds to the limit case of incompressibility $\lambda \to \infty$. In addition, the presented cases of isotropic material linearized at rest condition represent very simplified problems compared to possible realistic application. Yet, we see a significant gain between a factor 2 and 3 for different geometries.
We remark that a 3D simulation with 1,5M of DoFs of elastic wave propagating within a nearly-incompressible isotropic homogeneous cube takes around 6 minutes.

In order to test the scheme with a more complex scenario, we consider a thick hollow sphere, with the shear parameter $\mu$ linearly varying through the thickness, such that $\mu \in [1, 1.5]$. We measure an average ratio $C_S/C_P \simeq 4.24$ and a computational gain $G_{eff} \simeq 3.06$. Looking at results from the circle and the sphere, we observe as expected that non-square domains and a material slightly more complex than the linear isotropic one correspond to an increase of the computational cost for the shear operator $C_S$, and so a higher computational gain $G_{eff}$. In a realistic application such as the biomedical one, biological tissues are often hyperelastic, anisotropic and incompressible, undergoing large load. If we were supposed to model this type of problems, then the shear operator is more complex to compute, such that we expect a higher efficiency of the proposed scheme (LFC) scheme with respect to the classic (LF) scheme.

<table><caption>**Table 1.** Numerical results of time efficiency of LF-Chebyshev Scheme with respect to Leap-Frog scheme, considering different domains: a square, a circle, a cube and a thick hollow sphere. The square, circle and cube are characterized by a linear isotropic behavior, while the thick hollow sphere presents the shear parameter $\mu$ linearly increasing in the radial direction. $\eta_{LF}$ and $\eta_{LFC}$ are the measured simulation time for the computation of the elastic wave propagation in an isotropic homogeneous solid. $G_{eff}$ is the effective computational gain measured as $\eta_{LF}/\eta_{LFC}$. $C_S/C_P$ is the average ratio of the computational cost between the shear and the pressure operator over 1000 of iterations.</caption><thead><tr><th>Shape</th><th>η<sub>LF</sub>[s]</th><th>η<sub>LFC</sub>[s]</th><th>G<sub>eff</sub></th><th>C<sub>S</sub>/C<sub>P</sub></th><th>NDoF</th></tr></thead><tbody><tr><td>Square</td><td>12.24</td><td>5.23</td><td>2.34</td><td>1.99</td><td>51k</td></tr><tr><td>Circle</td><td>41.82</td><td>14.76</td><td>2.83</td><td>2.77</td><td>154k</td></tr><tr><td>Cube</td><td>799.0</td><td>377.18</td><td>2.12</td><td>1.76</td><td>1,594k</td></tr><tr><td>H. Sphere</td><td>427.03</td><td>138.779</td><td>3.06</td><td>4.24</td><td>2,519k</td></tr></tbody></table>

### 4.2.2 Space-time Convergence analysis

As a final validation of the proposed numerical scheme, we present a convergence analysis over space and time, comparing different boundary conditions. We consider a 2D square with $\Omega = [0, 1]^2$, discretized with $N$ elements for each direction, from $N = 8$ to 256. The parameters of the source are the following: $A_m = 10^3$, $r_0 = 0.05$, $\sigma_s = 0.5$, $\sigma_\tau = 0.5 \cdot 10^{-4}$.

We compute the $L^\infty(L^2(\Omega))$ norm and the $L^\infty(H^1(\Omega))$ norm of the displacement error taking the LF scheme as a reference, and we plot it with respect to the space-step $h$. In these simulations we take $\lambda = 1000$. For each image, we consider the convergence of simulations with final times of simulation $T = 0.6$ and $T = 1.0$. In shorter simulations ($T = 0.6$) the shear wave almost reaches the boundary, but is not reflected, while with the longer simulations $T = 1.0$, we have reflection and possible effects of boundary conditions. Figure 7 represents the results obtained with a free boundary ($t = 0$). Results show that second order convergence is maintained for both $L^2(\Omega)$ and $H^1(\Omega)$ norms (in space) for the shortest simulations ($T = 0.6$), (i.e., no reflection). With longer simulations ($T = 1.0$), we see instead a strong deterioration of convergence rate, which is not the same for the two norms: the order is two times smaller for $L^2(\Omega)$ norm and four times smaller for $H^1(\Omega)$. Then, we choose Dirichlet boundary condition $\underline{u} = 0$ on $\Gamma$, and we present the related convergence rate in figure 8. The simulations are obtained with the same material and source parameters as before. The convergence rate is now no longer affected by the boundary condition, as we see that no deterioration appears between simulations with different length in time. In addition, we see that the error is homogeneous between the two space norms $L^2(\Omega)$ and $H^1(\Omega)$, as expected. We point out that with Dirichlet boundary condition, the proposed (LFC) scheme maintains the second order accuracy of the classical (LF) scheme. The loss of accuracy observed when using free surface boundary conditions is due to a boundary layer effect. Such effect is comparable to what is observed in standard penalization strategies used to treat purely incompressible stokes problems (see e.g. [38]). We leave the treatment of such difficulty as a perspective.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Convergence of scheme LFC with respect to the space step $h$ – with $t = 0$ on $\Gamma$ (*free boundary condition*). On the left: we plot the $L^\infty(L^2(\Omega))$-displacement error for a final simulation time $T = 0.6$ (before the shear wave reaches the boundary) and $T = 1.0$ (after reflection on the boundary). On the right we have an equivalent plot but with $L^\infty(H^1(\Omega))$-displacement error.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Convergence of scheme LFC with respect to the space step $h$ – with $u = 0$ on $\Gamma$ (*Dirichlet boundary condition*). On the left, we plot the $L^\infty(L^2(\Omega))$ norm of displacement error. On the right, we plot the $L^\infty(H^1(\Omega))$ norm of the displacement error. In both case we consider a final simulation time $t_f = 0.5$ – before the shear wave reaches the boundary – and $t_f = 1.0$ – after reflection on the boundary.</figcaption></figure>

# 5 Numerical application: elastic wave propagation in an incompressible hyperelastic stretched cube

In the previous section, we presented the numerical scheme for wave propagation problem in nearly-incompressible isotropic solids linearized around at rest condition. We now want to present the applicability of the LFC scheme for the approximation of linear wave propagation problems within solids with large deformations. The formulation obtained in section 2 allows to give a physical interpretation of the wave propagation, depending on material and prestress effect.

In the following applications, the solution of the non-linear nearly-incompressible static problem is approximated by the solution of the pure incompressible problem (45) with a non-homogeneous Dirichlet boundary condition yet unspecified. We solve the following static non-linear problem :
Find $(\underline{u}_0, p_0) \in \mathcal{W} \times \mathcal{L}$ such that

$$ \left\{ \begin{array}{l} \displaystyle \int_{\hat{\Omega}} \underline{\Sigma}_0^d : D_{\underline{e}_0}(\hat{v}) \, d\hat{\Omega} - \int_{\hat{\Omega}} p_0 J_0 \underline{C}_0^{-1} : D_{\underline{e}_0}(\hat{v}) \, d\hat{\Omega} = 0 \quad \forall \hat{v} \in \mathcal{V}(\hat{\Omega}) \\[1.5em] \displaystyle \int_{\hat{\Omega}} \underline{\varrho}(J_0 - 1) \cdot q \, d\hat{\Omega} = 0 \quad \forall q \in \mathcal{L}(\hat{\Omega}) \end{array} \right. \qquad (45) $$

where $p_0$ is now a Lagrangian multiplier.

Let us consider a hyperelastic isotropic material with the Neo-Hookean elastic potential

$$ \mathcal{W}^d = \mu(J_1 - 3) \quad \text{with} \quad J_1 = I_1 I_3^{-\frac{1}{3}} \qquad (46) $$

where $J_1$ is the first reduced invariant presented in (4), and $\mu \in \mathbb{R}^+$ is the material parameter. As presented in the previous section, we consider a non-dimensionalization of the problem with respect to a characteristic length $L$, an observation time $T$ and the material parameter $\mu$.

Solving the non-linear static problem will give us access to the deformed apparent configuration $\Omega_0$ given by $\underline{u}_0$, the related state of stress characterized by $\underline{\sigma}_0$ and the internal pressure $p_0$.

The wave propagation problem is then obtained by linearization around $(\underline{u}_0, p_0)$, with the introduction of a large penalizing constant $\lambda \in \mathbb{R}^+$

For any time $t \in [0, 1]$, find $(\tilde{\underline{u}}(t), \tilde{p}(t)) \in \mathcal{V} \times \mathcal{L}$ such that

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega_0} \varrho_0 \partial_{tt} \tilde{\underline{u}} \cdot \underline{v} \, d\Omega_0 + \int_{\Omega_0} \zeta \underline{\varepsilon}(\tilde{\underline{u}}) : \underline{\varepsilon}(\underline{v}) \, d\Omega_0 + \int_{\Omega_0} \underline{\nabla} \underline{\varepsilon} \underline{\underline{v}} \cdot \tilde{\underline{u}} \underline{\sigma}_0 : \underline{\nabla} \underline{\underline{v}} \underline{v} \\ \qquad + \int_{\Omega_0} \tilde{p} \, \mathrm{div} \, \underline{v} \, d\Omega_0 = \int_{\Omega_0} \tilde{\varrho}_0 \underline{f} \cdot \underline{v} \, d\Omega_0 \quad \forall \underline{v} \in \mathcal{V} \\[1.5em] \displaystyle \int_{\Omega_0} q \, \mathrm{div} \, \tilde{\underline{u}} \, d\Omega_0 = \lambda^{-1} \int_{\Omega_0} \tilde{p} q \, d\Omega_0 \quad \forall q \in \mathcal{L} \end{array} \right. \qquad (47) $$

with

$$ \tilde{\underline{u}}(t=0) = \partial_t \tilde{\underline{u}}(t=0) = 0 \quad \text{in } \Omega_0 $$

where the material tangent stiffness is now given by the scalar $\zeta$, that depends on the loading condition of the non-linear static problem. Notice that the reduction of the material stiffness to a scalar $\zeta$ is related to our choice of the potential $\mathcal{W}^d$ in (46). For the following simulations, we will always take the penalizing parameter $\lambda = 1000$.

For the application of the scheme, we consider a simple example for the non-linear static case, such that the solution of problem (45) can be found analytically. In order to test the scheme and see the effects of the loaded state over the wave propagation, we take a homogeneous cube with a domain $\hat{\Omega} = [0, 1]^3$ and we choose a non-linear solid with isotropic homogeneous property as presented in (46).

The solution is sought as a uniform stretch state given by $\lambda_1, \lambda_2$ and $\lambda_3$, which describe the stretch in the Cartesian coordinates. We have a deformation map $\phi_0(\xi)$, that reads

$$ \phi_0(\xi) = (\lambda_1 \xi_1, \lambda_2 \xi_2, \lambda_3 \xi_3)^T $$

with $\lambda_1 = \bar{\lambda}_1$ (given) and $\lambda_2 = \lambda_3 = 1/\sqrt{\bar{\lambda}_1}$ due to incompressibility.

The state of stress in the reference configuration is given by the following relation

$$ \underline{\Sigma}(\underline{u}, p) = \underline{\Sigma}_d - pJ\underline{C}^{-1} = \frac{\partial \mathcal{W}^d}{\partial \underline{e}} - pJ\underline{C}^{-1} \qquad (48) $$

where the hydrostatic pressure $p$ is now defined as a lagrangian multiplier. We can derive analytically the state of stress in the apparent configuration, taking into account that in the transverse direction there is no constraint. Then, the non-dimensionalized Cauchy tensor is in the form $\underline{\sigma}_0 = \sigma_0 \underline{e}_1 \otimes \underline{e}_1$ with

$$ \begin{aligned} p_0 &= -\frac{2}{3} \left( \bar{\lambda}_1^2 - \frac{1}{\bar{\lambda}_1} \right) \\ \sigma_0 &= \frac{4}{3} \left( \bar{\lambda}_1^2 - \frac{1}{\bar{\lambda}_1} \right) - p_0 = 2 \left( \bar{\lambda}_1^2 - \frac{1}{\bar{\lambda}_1} \right) \end{aligned} $$
<figure><img src="image_10.png" /><figcaption>Figure 9. Examples of loaded configurations analyzed for the application of the numerical scheme</figcaption></figure>

Once we have the solution of the non-linear static problem, we can derive the parameters for the wave propagation problem (47). The geometric stiffness is given by $\underline{\underline{\sigma}}_0$, while the material stiffness is obtained by linearization and it reduces to

$$ \zeta = \frac{4}{3} \text{tr}(\underline{\underline{C}}_0) + 2p_0 = \frac{4}{\bar{\lambda}_1} \qquad (49) $$

From these preliminary results, we can already observe that increasing the imposed stretch $\bar{\lambda}_1$ results in an increase of the geometric stiffness and a reduction of the material stiffness parameter $\zeta$.

The volumetric transient source for the excitation of the wave is defined with the previous formula (43) and polarized in the $e_3$-direction such that $f = f(x_0, t) \cdot e_3$. The spatial excitation is localized in the center of the stretched cube. We take $r_0 = 0.05$, $\sigma_s = 10$ and $A_m = 5$. For the temporal push we have a Ricker wavelet with a small offset $t_0 = 0.04$ and $\sigma_\tau = 200$.

In figure 10, we see an example of propagation of the elastic wave within the incompressible cube. The mesh is discretized with 30 hexahedra in each direction, with $Q_4 - Q_2^{disc}$ elements and 5.3M DoFs. The average time to run these simulations in a normal machine with visual results is of 6 minutes.

In order to see the effect of loaded configuration over the propagation of the wave, we consider three cases of stretch $\bar{\lambda}_1 = [1.0, 1.1, 1.2]$. In figure 11, we present results of wave propagation in the mid $x_1 - x_2$ plane of the cube. In the first case, where the cube has no pre-stretch, the wave propagates with a circular symmetric shape as the one given by the source. By increasing the static stretch $\bar{\lambda}_1$, we observe a deformation of the wave front, for which the wave is faster in the direction $x_1$ of the imposed stretch and slower in the perpendicular direction $x_2$. As we noticed before, the material stiffness has actually decreased since the parameter $\zeta$, related to the Neo-Hookean density function, is inversely proportional to the pre-stretch $\bar{\lambda}_1$. As a consequence, a major role for the variation in the wave propagation is played by the geometric stiffness, which depends on the non-linear variation of strains and the pre-stress which has the only non-zero component in the $\underline{x}$ direction.

To have a final validation of these results, we looked at the elastic wave propagation in the cube artificially turning off and on different terms of the problem, in order to compare the different effects. In figure 12(a) we show the reference case, wave propagation in the cube subjected to a prestretch $\bar{\lambda}_1 = 1.2$, which corresponds to the simulation previously discussed (figure 11(c)). Then, we consider (figure 12(b)) the wave propagation in a domain that corresponds to the one obtained from the non-linear static problem with $\bar{\lambda}_1 = 1.2$, but in the simulation we do not account for any pre-stretch effect on the material and geometric stiffness. The simulation in figure 12(c) show the wave propagation with only the prestress effect due to $\bar{\lambda}_1 = 1.2$, while we "turn off" the material stiffness term taking $\zeta = 4$. Finally, in figure 12(d), we have taken $\sigma_0 = 0$ and $\zeta = 4/\bar{\lambda}_1$ with $\bar{\lambda}_1 = 1.2$, which corresponds to a material
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Simulation of 3D elastic wave propagation within an incompressible isotropic cube. In the images we observe the propagation in time (from left to right) of the magnitude of the displacement field. The transient source is focalized in the center of the cube and polarized in the vertical direction. Magnitude displacement in [0.0, 0.1]</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Increasing effect of pre-stretch over elastic wave propagation in nearly-incompressible cube with non-linear isotropic behavior subjected to uniaxial tension. Displacement magnitude in the middle $x_1 - x_2$ plane. From left $\bar{\lambda}_1 = [1.0, 1.1, 1.2]$. Mesh discretized with 5.3M DoFs. Wave observed at $t = 0.3$</figcaption></figure>

behavior with the linearized material stiffness term and *no prestress* effect is considered.

The resulting displacement fields reflect what we expected by looking at the physical terms of our problem (47). When considering only the linearized material behavior (figure 12(d)), we see that the wave propagates homogeneously slower with respect to the unloaded behavior given in (figure 12(a)). As mentioned before, the material stiffness parameter $\zeta$ is inversely proportional to the pre-stretch (49), which has a visible effect on the propagation of the shear wave. Obviously, this is related to the choice of the non-linear potential. The neo-hookean potential (46) presents a decreasing stiffness (slope) with the increasing loading condition, which corresponds to a softer response of the linearized material.

Still, the overall results that we have observed in figure 11 is an anisotropic response of the solid, with a relevant increase of the wave speed in the direction of the pre-stretch and a decrease in the perpendicular direction. As we can see in figure 12(c), the geometric stiffness (prestress) is the most relevant term in the propagation of the shear wave. Indeed, if we compare this result with the "real" solution of the preloaded wave propagation problem 12(a),
the wave with only the geometric stiffness (c) actually propagates faster. Trivially, this result comes from the sum of the geometric effect with the material effect, where the latter tends to lower the stiffness of the material and so the speed of the wave front.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Comparison of physical effects over wave propagation at time $t = 0.3$. (a) Reference case: linear wave propagation without any preload. (b) Deformed geometry as for $\bar{\lambda}_1 = 1.2$, but no preload effect accounted in the problem. (c) Solution for $\bar{\lambda}_1 = 1.2$ with related material stiffness, while the geometric stiffness term is neglected ($\sigma_0 = 0$). (d) Solution for $\bar{\lambda}_1 = 1.2$ with related geometric stiffness, while the material stiffness term is neglected ($\zeta = 4$).</figcaption></figure>

# 6 Conclusion

In conclusion, we have presented a stable, robust and efficient explicit scheme for the simulation of linear elastic wave propagation in nearly-incompressible solids. The scheme allows the approximation of solids with complex geometries subjected to large static deformations.

The efficiency of the scheme is based on the use of *damped Chebyshev polynomials*, applied here for the first time to treat dynamic problems with nearly-incompressible solids. This method, together with the space discretization strategy, guarantees a reduction of simulation time of at least a factor 2 with respect to the Leap-Frog method, second order convergence and preservation of energy. For slightly more complex problems we observed an improved reduction of the simulation time (around a factor 3). Finally, we have presented the applicability of the scheme for wave propagation in non-linear solids with large deformations. In this latter application, we have considered Dirichlet's boundary conditions on the whole solid, in order to guarantee a second order convergence of the method. Further improvements of the method could focus either on the treatment of the boundary layer effect when considering Neumann conditions (i.e., free boundary surface), in order to guarantee an optimal accuracy also when the shear wave reaches the boundary; or it could focus on extending the proposed strategy to adaptive time steps and non-linear elastodynamics using for instance as a starting point the explicit scheme of [39].

We recall that the method has been developed for the reproduction of dynamic elastographic measurements of biological tissues, a valuable tool to improve early diagnosis and follow-up treatment. In this scenario, materials are often characterized by strong non-linear anisotropic response, complex geometries and viscosity (or in general material with frequency-dependent behavior [40]). While these aspects add further layers of complexity to the problem, they might also improve the efficiency of the proposed method to reduce the computational cost associated to the enforcement of nearly-incompressibility.
References

[1] Theodore Sussman and Klaus-Jürgen Bathe. "A finite element formulation for nonlinear incompressible elastic and inelastic analysis". In: *Computers and Structures* 26.1-2 (1987), pp. 357-409.

[2] Klaus-Jürgen Bathe. *Finite element procedures*. Klaus-Jurgen Bathe, 2006.

[3] J.T. Oden and J.E. Key. "Numerical analysis of finite axisymmetric deformations of incompressible elastic solids of revolution". In: *International Journal of Solids and Structures* 6.5 (1970), pp. 497-518. ISSN: 0020-7683. DOI: https://doi.org/10.1016/0020-7683(70)90027-2. URL: https://www.sciencedirect.com/science/article/pii/0020768370900272.

[4] Robert L. Taylor, Karl S. Pister, and Leonard R. Herrmann. "On a variational theorem for incompressible and nearly-incompressible orthotropic elasticity". In: *International Journal of Solids and Structures* 4.9 (1968), pp. 875-883. ISSN: 0020-7683. DOI: https://doi.org/10.1016/0020-7683(68)90084-X. URL: https://www.sciencedirect.com/science/article/pii/002076836890084X.

[5] Federica Caforio and Sébastien Imperiale. "A conservative penalisation strategy for the semi-implicit time discretization of the incompressible elastodynamics equation". In: *Advanced Modeling and Simulation in Engineering Sciences* 5.1 (2018), p. 30. DOI: 10.1186/s40323-018-0121-8. URL: https://doi.org/10.1186/s40323-018-0121-8.

[6] Patrick Le Tallec. "Numerical methods for nonlinear three-dimensional elasticity". In: *Handbook of numerical analysis* 3 (1994), pp. 465-622.

[7] T. Elguedj et al. "B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements". In: *Computer Methods in Applied Mechanics and Engineering* 197.33 (2008), pp. 2732-2762. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2008.01.012. URL: https://www.sciencedirect.com/science/article/pii/S0045782508000248.

[8] Thomas JR Hughes. "Equivalence of finite elements for nearly incompressible elasticity". In: *J. Appl. Mech*. (1977).

[9] F. Brezzi. "On the Existence, Uniqueness and Approximation of Saddle-Point Problems Arising from Lagrangian Multipliers". In: *Publications mathématiques et informatique de Rennes* S4 (1974), pp. 1-26. URL: http://eudml.org/doc/274605.

[10] Ivo Babuška. "The finite element method with Lagrangian multipliers". In: *Numerische Mathematik* 20.3 (1973), pp. 179-192.

[11] Daniele Boffi and Rolf Stenberg. "A remark on finite element schemes for nearly incompressible elasticity". In: *Computers and Mathematics with Applications* 74.9 (2017). Advances in Mathematics of Finite Elements, honoring 90th birthday of Ivo Babuška, pp. 2047-2055. ISSN: 0898-1221. DOI: https://doi.org/10.1016/j.camwa.2017.06.006. URL: https://www.sciencedirect.com/science/article/pii/S0898122117303541.

[12] Alfio Quarteroni and Alberto Valli. *Numerical approximation of partial differential equations*. Vol. 23. Springer Science and Business Media, 2008.

[13] Franco Brezzi and Juhani Pitkäranta. "On the stabilization of finite element approximations of the Stokes equations". In: *Efficient Solutions of Elliptic Systems: Proceedings of a GAMM-Seminar Kiel, January 27 to 29, 1984*. Springer. 1984, pp. 11-19.

[14] Thomas J.R. Hughes and Leopoldo P. Franca. "A new finite element formulation for computational fluid dynamics: VII. The stokes problem with various well-posed boundary conditions: Symmetric formulations that converge for all velocity/pressure spaces". In: *Computer Methods in Applied Mechanics and Engineering* 65.1 (1987), pp. 85-96. ISSN: 0045-7825. DOI: https://doi.org/10.1016/0045-7825(87)90184-8. URL: https://www.sciencedirect.com/science/article/pii/0045782587901848.
[15] Thomas JR Hughes, Leopoldo P Franca, and Marc Balestra. "A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations". In: *Computer Methods in Applied Mechanics and Engineering* 59.1 (1986), pp. 85–99.

[16] Claudio Baiocchi, Franco Brezzi, and Leopoldo P Franca. "Virtual bubbles and Galerkin-least-squares type methods (Ga. LS)". In: *Computer methods in applied mechanics and engineering* 105.1 (1993), pp. 125–141.

[17] Douglas N Arnold, Franco Brezzi, and Michel Fortin. "A stable finite element for the Stokes equations". In: *Calcolo* 21.4 (1984), pp. 337–344.

[18] Ernst Hairer et al. "Geometric numerical integration". In: *Oberwolfach Reports* 3.1 (2006), pp. 805–882.

[19] Alfio Quarteroni. *Modellistica numerica per problemi differenziali*. Vol. 2. Springer Science Business Media, 2009.

[20] Gary C Cohen. *Higher-order numerical methods for transient wave equations*. Vol. 5. Springer, 2002.

[21] Constantin Carle, Marlis Hochbruck, and Andreas Sturm. "On leapfrog-Chebyshev schemes". In: *SIAM Journal on Numerical Analysis* 58.4 (2020), pp. 2404–2433.

[22] Marcus Grote, Simon Michel, and Stefan Sauter. "Stabilized leapfrog based local time-stepping method for the wave equation". In: *Mathematics of Computation* 90.332 (2021), pp. 2603–2643.

[23] W Ye et al. "A selective mass scaling method for shear wave propagation analyses in nearly incompressible materials". In: *International Journal for Numerical Methods in Engineering* 109.2 (2017), pp. 155–173.

[24] Lars Olovsson, Kjell Simonsson, and Mattias Unosson. "Selective mass scaling for explicit finite element analyses". In: *International Journal for Numerical Methods in Engineering* 63.10 (2005), pp. 1436–1445.

[25] Michel Destrade and Giuseppe Saccomandi. "Waves in Nonlinear Pre-Stressed Materials". In: 2007. URL: https://api.semanticscholar.org/CorpusID:137240942.

[26] Moniba Shams, Michel Destrade, and Ray W Ogden. "Initial stresses in elastic solids: constitutive laws and acoustoelasticity". In: *Wave Motion* 48.7 (2011), pp. 552–567.

[27] Federica Caforio. "Mathematical modelling and numerical simulation of elastic wave propagation in soft tissues with application to cardiac elastography". 2019SACLX001. PhD thesis. 2019. URL: http://www.theses.fr/2019SACLX001/document.

[28] André Dalmora et al. "A time-domain spectral finite element method for acoustoelasticity: modeling the effect of mechanical loading on guided wave propagation". In: *Wave Motion* 129 (2024), p. 103328.

[29] Philippe G Ciarlet. *Three-dimensional elasticity*. Elsevier, 1988.

[30] Philippe G Ciarlet and Giuseppe Geymonat. "Sur les lois de comportement en élasticité non linéaire compressible". In: (1982).

[31] Patrick Joly. "Numerical methods for elastic wave propagation". In: *Waves in nonlinear pre-stressed materials* (2007), pp. 181–281.

[32] Gonçalo Pena. "Spectral element approximation of the incompressible Navier-Stokes equations in a moving domain and applications". PhD thesis. Ecole Polytechnique Fédérale de Lausanne (EPFL), 2009.

[33] Christine Bernardi and Yvon Maday. "Uniform inf-sup conditions for the spectral discretization of the Stokes problem". In: *Mathematical Models and Methods in Applied Sciences* 9.03 (1999), pp. 395–414.
[34] Marc Duruflé, Pascal Grob, and Patrick Joly. "Influence of Gauss and Gauss-Lobatto quadrature rules on the accuracy of a quadrilateral finite element method in the time domain". In: *Numerical Methods for Partial Differential Equations: An International Journal* 25.3 (2009), pp. 526-551.

[35] Franco Brezzi and Richard S Falk. "Stability of higher-order Hood-Taylor methods". In: *SIAM Journal on Numerical Analysis* 28.3 (1991), pp. 581-590.

[36] R. E. Harte. "Spectral Mapping Theorems". In: *Proceedings of the Royal Irish Academy. Section A: Mathematical and Physical Sciences* 72 (1972), pp. 89-107. ISSN: 00358975. URL: http://www.jstor.org/stable/20488717.

[37] Franco Brezzi, Richard S Falk, and L Donatella Marini. "Basic principles of mixed virtual element methods". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 48.4 (2014), pp. 1227-1240.

[38] Daniele Boffi, Franco Brezzi, Michel Fortin, et al. *Mixed finite element methods and applications*. Vol. 44. Springer, 2013.

[39] Frédéric Marazzato et al. "An explicit pseudo-energy conserving time-integration scheme for Hamiltonian dynamics". In: *Computer Methods in Applied Mechanics and Engineering* 347 (2019), pp. 906-927.

[40] Nicolas Favrie and Bruno Lombard. "A hyperbolic generalized Zener model for nonlinear viscoelastic waves". In: *Wave Motion* 116 (2023), p. 103086.

[41] Robert W Penn. "Volume changes accompanying the extension of rubber". In: *Transactions of the Society of Rheology* 14.4 (1970), pp. 509-517.

[42] Jiun-Shyan Chen and Chunhui Pan. "A Pressure Projection Method for Nearly Incompressible Rubber Hyperelasticity, Part I: Theory". In: *Journal of Applied Mechanics* 63 (1996), pp. 862-868. URL: https://apisemanticscholar.org/CorpusID:121749424.

[43] TYP Chang, AF Saleeb, and G Li. "Large strain analysis of rubber-like materials based on a perturbed Lagrangian variational principle". In: *Computational mechanics* 8.4 (1991), pp. 221-233.

# A Derivation of elastic wave propagation in quasi-static loaded nearly-incompressible solid

This section is devoted to the derivation of the weak formulation of the elastic wave propagation problem in the framework of hyperelastic nearly-incompressible solids, given in Eq. (7). For the sake of clarity, we re-present some quantities that were already given in Section 2.

## Fundamental law of dynamics

Let us recall the standard weak formulation associated to the wave propagation problem defined in the current configuration $\Omega$:

For any $t \in [0, T]$, find $u(t) \in W(\Omega)$ such that

$$ \int_{\Omega} \rho (\partial_{tt} u \circ \underline{\phi}^{-1} - \underline{f}) \cdot \underline{v} \, d\Omega + \int_{\Omega} \underline{\sigma} : \underline{\nabla}_{\underline{x}} \underline{v} \, d\Omega = 0 \quad \forall \underline{v} \in \mathcal{V}(\Omega), \quad (50) $$

where $\rho$ is the density in the current configuration and $\underline{\sigma}$ is the Cauchy stress tensor defined in $\Omega(t)$. We remark that the current configuration in which the problem is written is unknown, so that it is classical to express the problem in a total Lagrangian form with respect to the reference configuration. We may neglect some details as it is a classical development, for which we refer again the reader to [29].

As a first step towards the Lagrangian form, we remark that for any quantity $w(\underline{\xi})$ defined in $\hat{\Omega}$, we have that

$$ \int_{\Omega} \underline{w} \circ \phi^{-1}(\underline{x}) \, d\Omega = \int_{\hat{\Omega}} J \underline{w}(\underline{\xi}) \, d\hat{\Omega} \quad (51) $$
and the gradient of a field in the reference configuration can be rewritten in the current configuration using

$$ (\underline{\nabla}_{\underline{\xi}} \underline{w}) \circ \underline{\phi}^{-1} = \underline{\nabla}_{\underline{x}} (\underline{w} \circ \underline{\phi}^{-1}) \cdot (\underline{F} \circ \underline{\phi}^{-1}). \quad (52) $$

We rewrite first the inertial term in $\hat{\Omega}$

$$ \int_{\hat{\Omega}} \rho (\partial_{tt} \underline{u} \circ \underline{\phi}^{-1} - \underline{f}) \cdot \underline{v} \, d\hat{\Omega} = \int_{\hat{\Omega}} \rho \circ \underline{\phi}(\underline{\xi}) (\partial_{tt} \underline{u} - \underline{f} \circ \underline{\phi}(\underline{\xi})) \cdot \underline{v} J d\hat{\Omega} = \int_{\hat{\Omega}} \hat{\rho}(\partial_{tt} \underline{u} - \hat{\underline{f}}) \cdot \hat{\underline{v}} \, d\hat{\Omega}. \quad (53) $$

where we note the mass density in the reference configuration $\hat{\rho} : \hat{\Omega} \to \mathbb{R}$, and is given by $\hat{\rho}(\underline{\xi}) = J \rho \circ \underline{\phi}$. We also define the body force in the reference configuration as $\hat{f} := \underline{f} \circ \underline{\phi}(\underline{\xi})$.

From the relation between the second Piola Kirchhoff tensor $\underline{\Sigma}$ the Cauchy stress tensor given in Eq. (2), we have

$$ \int_{\hat{\Omega}} \underline{\sigma} : \underline{\nabla}_{\underline{x}} \underline{v} \, d\hat{\Omega} = \int_{\hat{\Omega}} J^{-1} \underline{F} \circ \underline{\Sigma} \cdot \underline{F}^T : \underline{\nabla}_{\underline{\xi}} \hat{\underline{v}} \cdot \underline{F}^{-1} \, J \, d\hat{\Omega} = \int_{\hat{\Omega}} \underline{\Sigma} : \underline{F}^T \cdot \underline{\nabla}_{\underline{\xi}} \hat{\underline{v}} \, d\hat{\Omega}. $$

We then introduce the Gateaux derivative of the Green-Lagrange strain tensor $\underline{e}$ at $\underline{u}$ in the direction $\hat{\underline{v}}$ given by the following equation

$$ D_{\underline{e}}(\underline{u}) \cdot \hat{\underline{v}} = \frac{1}{2} \left( \underline{\nabla}_{\underline{\xi}} \hat{\underline{v}}^T \cdot \underline{F}(\underline{u}) + \underline{F}(\underline{u})^T \cdot \underline{\nabla}_{\underline{\xi}} \hat{\underline{v}} \right) \quad (54) $$

for which we recall the definition of the Gateaux derivative

$$ D A(\underline{u}) \cdot \underline{w} = \lim_{\underline{h} \to 0} \frac{A(\underline{u} + \underline{h}\underline{w}) - A(\underline{u})}{\underline{h}}. $$

As $\underline{\Sigma}$ is a symmetric tensor, we can consider the contraction with the symmetrized version of $\underline{F} \cdot \underline{\nabla}_{\underline{\xi}}$, that corresponds to

$$ \int_{\hat{\Omega}} \underline{\Sigma} : \underline{F}^T \cdot \underline{\nabla}_{\underline{\xi}} \hat{\underline{v}} \, d\hat{\Omega} = \int_{\hat{\Omega}} \underline{\Sigma} : \frac{1}{2} \left( \underline{F}^T \cdot \underline{\nabla}_{\underline{\xi}} \hat{\underline{v}} + \underline{\nabla}_{\underline{\xi}} \hat{\underline{v}}^T \cdot \underline{F} \right) d\hat{\Omega} = \int_{\hat{\Omega}} \underline{\Sigma} : D_{\underline{e}}(\underline{u}) \cdot \hat{\underline{v}} \, d\hat{\Omega}. \quad (55) $$

Combining (53) and (55) we obtain the weak formulation in the reference configuration, that reads:

For any $t \in [0; T]$, find $\underline{u} \in W(\hat{\Omega})$ such that

$$ \int_{\hat{\Omega}} \hat{\rho}(\partial_{tt} \underline{u} - \hat{\underline{f}}) \cdot \hat{\underline{v}} \, d\hat{\Omega} + \int_{\hat{\Omega}} \underline{\Sigma} : D_{\underline{e}}(\underline{u}) \cdot \hat{\underline{v}} \, d\hat{\Omega} = \int_{\Gamma_N} \hat{\underline{t}} \cdot \hat{\underline{v}} \, d\hat{\Omega} \quad \forall \hat{\underline{v}} \in V(\hat{\Omega}) \quad (56) $$

In the framework of hyperelastic materials, the state of stress within the solid depends on the hyperelastic potential $W^e = W^e(\underline{C})$ and, in the case of isotropic material, we have

$$ W^e = W^e(I_1, I_2, I_3) \quad \text{s.t.} \quad \underline{\Sigma} = \frac{\partial W^e}{\partial \underline{e}}(I_1, I_2, I_3) = 2 \sum_{i=1}^{3} \frac{\partial W^e}{\partial I_i} \frac{\partial I_i}{\partial \underline{C}} $$

where $I_1, I_2, I_3$ are the eigenvalues of $\underline{C}(\underline{u})$ known as invariants

$$ I_1 = \text{tr} \underline{C}, \quad I_2 = \text{tr}(\text{cot} \underline{C}) = \frac{1}{2} ((\text{tr} \underline{C})^2 - \text{tr}(\underline{C}^2)), \quad I_3 = \det \underline{C} = J^2. $$

The application of this work concerns hyperelastic nearly-incompressible materials, like soft tissues characterized by high water content. A classical decomposition of the Second Piola Kirchhoff tensor $\underline{\Sigma}$ is given by the definition of the hydrostatic pressure $p$ and the deviatoric component $\underline{\Sigma}^d$

$$ p = -\frac{1}{3J} \underline{\Sigma} : \underline{C} \quad \text{and} \quad \underline{\Sigma}^d = \underline{\Sigma} + pJ\underline{C}^{-1} \quad \text{s.t.} \quad \underline{\Sigma} = \underline{\Sigma}^d - pJ\underline{C}^{-1}. \quad (57) $$
The stress tensor is thereby decomposed in its deviatoric and volumetric term, and in particular we have $\Sigma^d : C = 0$ (volumetric deformation do not generate any work from the deviatoric component).

The nearly-incompressibility is taken into account with an additional strain energy density that constrains volume changes. This term can be interpreted as a dilatation strain energy function that depends only on $I_3$, for which a small variation in volume is penalized by a large bulk modulus. However, to achieve a separation in the energy density function between the component related to distortion and dilatation deformation [41, 42], we cannot use $I_1$ and $I_2$ as they are not fully isochoric. We introduce the reduced invariants

$$J_1 = I_1 I_3^{-1/3}, \quad J_2 = I_2 I_3^{-2/3}, \quad J_3 = I_3^{1/2} = J \quad (58)$$

and we do the classical splitting

$$\mathcal{W}^e(I_1, I_2, I_3) = \mathcal{W}^d(J_1, J_2) + \mathcal{W}^v(J),$$

where $\mathcal{W}^d(J_1, J_2)$ and $\mathcal{W}^v(J_3)$ correspond respectively to the elastic potential referred to deviatoric and volumetric deformations. We suggest the reader to look in [43] for further details. From definitions (57), we derive that

$$p(J) = -\frac{d\mathcal{W}^v(J)}{dJ} \quad \text{and} \quad \Sigma^d = 2 \sum_{i=1}^{2} \frac{\partial \mathcal{W}^d}{\partial J_i} \frac{\partial J_i}{\partial C}.$$

Consequently, we can specify an equivalent formulation of (56), with the additional term related to volumetric deformation. It reads:

For any $t \in [0; T]$ find $\underline{u}(t) \in W(\hat{\Omega})$ such that

$$\int_{\hat{\Omega}} \hat{\varrho} \partial_{tt} \underline{u} \cdot \hat{v} d\hat{\Omega} + \int_{\hat{\Omega}} \left( \Sigma^d(\underline{e}(\underline{u})) - p J C^{-1} \right) : D\underline{e}(\underline{u}) \cdot \hat{v} d\hat{\Omega} \\ = \int_{\hat{\Omega}} \hat{\varrho} \hat{f} \cdot \hat{v} d\hat{\Omega} + \int_{\Gamma_N} \hat{t} \cdot \hat{v} d\hat{s} \quad \forall \hat{v} \in V(\hat{\Omega}). \quad (59)$$

We recall that in (59) the terms $p, J$ and $C$ are function of the unknown solution $\underline{u}$, while $\hat{\varrho}, \hat{f}$ and $\hat{t}$ are given parameters.

## A.1 Elastic wave propagation in quasi-static loaded solid

Let us now consider that the load applied to the solid is defined as follows

$$\hat{f}(\xi, t) = \hat{f}_0(\xi) + \delta\hat{f}(\xi, t)$$

The solid undergoes large deformation due to the application of the quasi-static load $\hat{f}_0(\xi)$, and it is subjected to a transient small-amplitude perturbation $\delta\hat{f}(\xi, t)$ with $\delta$ as a small parameter, that generates the wave propagation. Solutions can be written as

$$\underline{u}(\xi, t) = \underline{u}_0(\xi) + \delta\underline{u}(\xi, t) + \mathcal{O}(\delta^2),$$

where $\underline{u}_0$ is the governing solution of the quasi-static non-linear mechanics and $\delta\underline{u}$ is the small-amplitude wave propagating in the loaded configuration.

### A.1.1 Wave propagation problem in the reference configuration

With the external load defined as in (5) we derive the governing equation for the small perturbation $\tilde{u}$, i.e., the wave propagation, by linearizing (56) around the static loaded configuration given by the displacement $\underline{u}_0$.

We start the process of linearization with the component related to the deviatoric component of the internal work

$$\int_{\hat{\Omega}} \Sigma^d(\underline{e}(\underline{u})) : D\underline{e}(\underline{u}) \cdot \hat{v} d\hat{\Omega} \quad (60)$$
First, we linearize the derivative of the Green-Lagrange strain tensor

$$D\underline{e}(\underline{u}) \cdot \hat{v} = D\underline{e}(\underline{u}_0) \cdot \hat{v} + \delta D^2 \underline{e}(\tilde{u}, \hat{v}) + \mathcal{O}(\delta^2) \quad (61)$$

with

$$D^2 \underline{e}(\tilde{u}, \hat{v}) := \frac{1}{2} \left( \nabla_{\underline{\xi}} \tilde{u}^T \nabla_{\underline{\xi}} \hat{v} + \nabla_{\underline{\xi}} \hat{v}^T \nabla_{\underline{\xi}} \tilde{u} \right)$$

For the linearization of the deviatoric component of the stress tensor $\underline{\Sigma}^d$, let us first introduce some useful notation that we will use for the rest of the paper

$$D\underline{e}_0(\underline{w}) := D\underline{e}_0(\underline{w}) \circ \underline{\phi}(\underline{\xi}) := D\underline{e}_0(\underline{u}_0) \cdot \underline{w} = \frac{1}{2} \left( \nabla_{\underline{\xi}} \underline{w}^T \underline{F}_0 + \underline{F}_0^T \nabla_{\underline{\xi}} \underline{w} \right) \quad (62)$$

$$C_0^d : \underline{s} := D\underline{\Sigma}^d(\underline{e}(\underline{u}_0)) : \underline{s} \quad (\text{Gateaux derivative}) \quad (63)$$

where $\underline{F}_0 = \underline{F}(\underline{u}_0)$ and the first order derivative of $\underline{e}$ evaluated in $\underline{u}_0$ in the direction $\underline{w}$ is now denoted with $D\underline{e}_0(\underline{w})$. The 4th-order tensor $C_0^d$ is characterized by minor and major symmetries, and it accounts for the deviatoric behavior of the solid at the quasi-static configuration defined by $\underline{u}_0$. The linearization of the deviatoric stress reads

$$\underline{\Sigma}^d(\underline{e}(\underline{u})) = \underline{\Sigma}^d_0 + \delta\underline{\Sigma}^d + \mathcal{O}(\delta^2) \quad \text{with} \quad \begin{cases} \underline{\Sigma}^d_0 := \underline{\Sigma}^d(\underline{e}(\underline{u}_0)) \\ \underline{\Sigma}^d := C_0^d : D\underline{e}_0(\tilde{u}) \end{cases} \quad (64)$$

Finally, due to the symmetry of $\underline{\Sigma}^d$ and $C_0^d$, we write

$$\int_{\hat{\Omega}} \underline{\Sigma}^d(\underline{e}(\underline{u})) : D\underline{e}(\underline{u}) \cdot \hat{v} d\hat{\Omega} = \\ = \int_{\hat{\Omega}} \underline{\Sigma}^d_0 : D\underline{e}_0(\hat{v}) d\hat{\Omega} + \delta \left[ \int_{\hat{\Omega}} D\underline{e}_0(\tilde{u}) : C_0^d : D\underline{e}_0(\hat{v}) d\hat{\Omega} + \int_{\hat{\Omega}} \underline{\Sigma}^d_0 : \nabla_{\underline{\xi}} \tilde{u}^T \nabla_{\underline{\xi}} \hat{v} d\hat{\Omega} \right] + \mathcal{O}(\delta^2). \quad (65)$$

Let us now present the linearization of the volumetric stress term

$$\int_{\hat{\Omega}} p J \underline{C}^{-1} : D\underline{e}(\underline{u}) \cdot \hat{v} d\hat{\Omega}.$$

We linearize all terms dependent on $\underline{u}$ as follows

$$\underline{C}^{-1} = \underline{C}^{-1}_0 + \delta\underline{\tilde{G}} \quad \text{with} \quad \begin{cases} \underline{C}^{-1}_0 = \underline{C}(\underline{u}_0)^{-1} \\ \underline{\tilde{G}} = -2 \underline{C}^{-1}_0 \cdot D\underline{e}_0(\tilde{u}) \cdot \underline{C}^{-1}_0 \end{cases} \\ \\ J = J_0 + \delta\tilde{J} \quad \text{with} \quad \begin{cases} J_0 = [\det(\underline{C}_0)]^{1/2} \\ \tilde{J} = J_0 \underline{C}^{-1}_0 : D\underline{e}_0(\tilde{u}) \end{cases} \\ \\ p = p_0 + \delta\tilde{p} \quad \text{with} \quad \begin{cases} p_0 = p(J_0) = -\frac{d\mathcal{W}^v}{dJ}(J_0) \\ \tilde{p} = \frac{dp}{dJ}(J_0)\tilde{J} = -\frac{d^2\mathcal{W}^v}{dJ^2}(J_0) J_0 \underline{C}^{-1}_0 : D\underline{e}_0(\tilde{u}). \end{cases}$$
From these results and the previous linearization of the deformation tensor, the volumetric term reads

$$ \int_{\hat{\Omega}} pJ \underline{C}^{-1} : D_{\underline{e}}(\underline{u}) \cdot \hat{v} d\hat{\Omega} = \int_{\hat{\Omega}} p_0 J_0 \underline{C}_0^{-1} : D_{\underline{e}_0}(\hat{v}) d\hat{\Omega} + \\ + \delta \left[ \int_{\hat{\Omega}} (\tilde{p}J_0 + p_0 \tilde{J}) \underline{C}_0^{-1} : D_{\underline{e}_0}(\hat{v}) d\hat{\Omega} + \int_{\hat{\Omega}} p_0 J_0 \tilde{\underline{G}} : D_{\underline{e}_0}(\hat{v}) d\hat{\Omega} \right. \\ \left. + \int_{\hat{\Omega}} p_0 J_0 \underline{C}_0^{-1} : \underline{\nabla}_{\underline{\xi}} \tilde{u}^T \underline{\nabla}_{\underline{\xi}} \hat{v} d\hat{\Omega} \right] + \mathcal{O}(\delta^2). \quad (66) $$

The linearization of the inertial term takes into account the static condition of the non-linear solution, hence we have:

$$ \int_{\hat{\Omega}} \hat{\varrho} \partial_{tt} \underline{u} \cdot \hat{v} d\hat{\Omega} = \delta \int_{\hat{\Omega}} \hat{\varrho} \partial_{tt} \tilde{u} \cdot \tilde{v} d\hat{\Omega}. $$

Let us now separate the non-linear static part from the linear part of the problem, in order to distinguish the different physics of the problem for which we will need different methods of approximation. We also assume no surface traction for the sake of simplicity ($t = 0$).

**Quasi-static problem** The non-linear problem in the reference configuration reads: Find $\underline{u}_0 \in W(\hat{\Omega})$ such that

$$ \left\{ \begin{array}{l} \displaystyle \int_{\hat{\Omega}} \left( \underline{\Sigma}_0^d - p_0 J_0 \underline{C}_0^{-1} \right) : D_{\underline{e}_0}(\hat{v}) d\hat{\Omega} = \int_{\hat{\Omega}} \hat{\varrho} \underline{f}_0 \cdot \hat{v} d\hat{\Omega} \quad \forall \hat{v} \in \mathcal{V}(\hat{\Omega}) \\ \\ p_0 = - \frac{d\mathcal{W}^v}{dJ}(J_0) \end{array} \right. \quad (67) $$

We recall that the Second Piola Kirchhoff evaluated at $\underline{u}_0$, that will be referred to as the prestress is

$$ \underline{\Sigma}_0 := \underline{\Sigma}(\underline{u}_0) = \underline{\Sigma}_0^d - p_0 J_0 \underline{C}_0^{-1}. \quad (68) $$

Once the non-linear solution is given, we obtain the wave propagation in the quasi-static loaded solid by solving the linearized problem, following

**Wave equation in reference configuration** For any $t \in [0, T]$, find $\tilde{u}(t) \in \mathcal{V}(\hat{\Omega})$ such that

$$ \left\{ \begin{array}{l} \displaystyle \int_{\hat{\Omega}} \hat{\varrho} \partial_{tt} \tilde{u} \cdot \hat{v} d\hat{\Omega} + \int_{\hat{\Omega}} D_{\underline{e}_0}(\tilde{u}) : \underline{C}_0^d : D_{\underline{e}_0}(\hat{v}) d\hat{\Omega} \\ \\ \displaystyle \qquad + 2 \int_{\hat{\Omega}} p_0 J_0 \underline{C}_0^{-1} \cdot D_{\underline{e}_0}(\tilde{u}) : \underline{C}_0^{-1} \cdot D_{\underline{e}_0}(\hat{v}) d\hat{\Omega} \\ \\ \displaystyle \qquad + \int_{\hat{\Omega}} \lambda^* \underline{C}_0^{-1} : D_{\underline{e}_0}(\tilde{u}) \underline{C}_0^{-1} : D_{\underline{e}_0}(\hat{v}) d\hat{\Omega} \\ \\ \displaystyle \qquad + \int_{\hat{\Omega}} \underline{\Sigma}_0 : \underline{\nabla}_{\underline{\xi}} \tilde{u}^T \underline{\nabla}_{\underline{\xi}} \hat{v} d\hat{\Omega} = \int_{\hat{\Omega}} \hat{\varrho} \underline{\tilde{f}} \cdot \hat{v} d\hat{\Omega} \end{array} \right. \quad (69) $$

where

$$ \lambda^* := -J_0^2 \frac{dp}{dJ}(J_0) - J_0 p_0 = J_0 \left( J_0 \frac{d^2 \mathcal{W}^v}{dJ^2}(J_0) + \frac{d \mathcal{W}^v}{dJ}(J_0) \right) $$

with initial condition

$$ \tilde{u}(0) = 0, \quad \partial_t \tilde{u}(0) = 0 $$

We remark that the parameter $\lambda^*$ will have a high value as the bulk modulus of the potential $\mathcal{W}^v$ is taken very large to enforce the incompressibility.
### A.1.2 Wave propagation problem in the apparent configuration

In an application framework, the wave propagation is measured in the loaded configuration described by $\underline{u}_0$. It is useful, indeed, to reformulate the linearized problem in the quasi-static loaded configuration denoted with $\Omega_0$, whose set of points $\underline{x}_0$ can be retrieved with the bijective deformation map $\underline{\phi}_0(\underline{\xi})$:

$$ \underline{\phi}_0 : \underline{\xi} \to \underline{x}_0(\underline{\xi}) = \underline{\phi}_0(\underline{\xi}) $$

$$ \underline{\phi}_0(\underline{\xi}) = \underline{u}_0(\underline{\xi}) + \underline{\xi} $$

We also introduce the functional space $\mathcal{V}(\Omega_0) := \{\underline{v} \in H^1(\Omega_0)^d\}$ where we seek the solution of the linear solution. Given a function $\underline{w}$ defined in $\Omega_0$, we can define the linear strain tensor in the apparent configuration as

$$ \underline{\varepsilon}_{\underline{x}_0}(\underline{w}) = \frac{1}{2}(\underline{\nabla}_{\underline{x}_0}\underline{w} + \underline{\nabla}_{\underline{x}_0}\underline{w}^T). \qquad (70) $$

We now want to refer the mechanical quantities of the wave equation (69) with respect to $\Omega_0$. We combine (52) and (62) to write the first order derivative of Green-Lagrange strain tensor evaluated in $\underline{u}_0$ w.r.t. the apparent configuration $\Omega_0$

$$ (\underline{D}_{\underline{\underline{e}}_0}(\underline{w} \circ \underline{\phi}_0) \circ \underline{\phi}_0^{-1})(\underline{x}_0) = \frac{1}{2} \left( \underline{F}_{\underline{\underline{e}}_0}^T \cdot \underline{\nabla}_{\underline{x}_0} \underline{w}^T \cdot \underline{F}_{\underline{\underline{e}}_0} + \underline{F}_{\underline{\underline{e}}_0}^T \cdot \underline{\nabla}_{\underline{x}_0} \underline{w} \cdot \underline{F}_{\underline{\underline{e}}_0} \right) = \underline{F}_{\underline{\underline{e}}_0}^T \cdot \underline{\varepsilon}_{\underline{x}_0}(\underline{w}) \cdot \underline{F}_{\underline{\underline{e}}_0} $$

From this relation one can derive the two equations below that are used to treat the terms related to the linearization of the volumetric work, in the apparent configuration we have

$$ (\underline{C}_{\underline{\underline{e}}_0}^{-1} : \underline{D}_{\underline{\underline{e}}_0}(\underline{w} \circ \underline{\phi}_0)) \circ \underline{\phi}_0^{-1} = \underline{1} : \underline{\varepsilon}_{\underline{x}_0}(\underline{w}) = \mathrm{div}_{\underline{x}_0} \underline{w}, $$

$$ (\underline{C}_{\underline{\underline{e}}_0}^{-1} \cdot \underline{D}_{\underline{\underline{e}}_0}(\tilde{\underline{u}}) : \underline{C}_{\underline{\underline{e}}_0}^{-1} \cdot \underline{D}_{\underline{\underline{e}}_0}(\underline{w} \circ \underline{\phi}_0)) \circ \underline{\phi}_0^{-1} = \underline{\varepsilon}_{\underline{x}_0}(\tilde{\underline{u}} \circ \underline{\phi}_0^{-1}) : \underline{\varepsilon}_{\underline{x}_0}(\underline{w}). $$

With the relation between the Cauchy tensor and the second Piola Kirchhoff (2), we bring in the apparent configuration the term depending on the pre-stress $\underline{\Sigma}_{\underline{\underline{e}}_0}$

$$ (\underline{\Sigma}_{\underline{\underline{e}}_0} : \underline{\nabla}_{\underline{\underline{\xi}}} \tilde{\underline{u}}^T \underline{\nabla}_{\underline{\underline{\xi}}} (\underline{w} \circ \underline{\phi}_0)) \circ \underline{\phi}_0^{-1} = J_0 \underline{\nabla}_{\underline{x}_0} (\tilde{\underline{u}} \circ \underline{\phi}_0^{-1}) \cdot \underline{\sigma}_{\underline{\underline{e}}_0} : \underline{\nabla}_{\underline{x}_0} \underline{w} $$

$$ \text{with } \underline{\sigma}_{\underline{\underline{e}}_0} := \underline{\sigma}(\underline{u}_0) \circ \underline{\phi}_0^{-1} = (J_0^{-1} \underline{F}_{\underline{\underline{e}}_0} \cdot \underline{\Sigma}_{\underline{\underline{e}}_0} \cdot \underline{F}_{\underline{\underline{e}}_0}^T) \circ \underline{\phi}_0^{-1} $$

For the sake of lightening the reading, we will now drop the notation $\cdot \circ \underline{\phi}_0^{-1}$.

We now use (51) considering the deformation map $\underline{\phi}_0$ between $\hat{\Omega}$ and $\Omega_0$, in order to obtain the formulation of the elastic wave propagation problem in the apparent configuration, which reads:
Wave equation in apparent configuration $\Omega_0$: For any $t \in [0, T]$, find $\tilde{u}(t) \in \mathcal{V}(\Omega_0)$ such that

$$ \left\{ \begin{array}{ll} \displaystyle \int_{\Omega_0} \varrho_0 \partial_{tt} \tilde{u} \cdot \underline{v} \, d\Omega_0 + \int_{\Omega_0} \underline{\varepsilon}_{\underline{x}_0} (\tilde{u}) : \tilde{\mathbf{C}}_0 : \underline{\varepsilon}_{\underline{x}_0} (\underline{v}) \, d\Omega_0 & \\ \\ \displaystyle + \int_{\Omega_0} \underline{\nabla}_{\underline{x}_0} \tilde{u} \underline{\sigma}_0 : \underline{\nabla}_{\underline{x}_0} \underline{v} \, d\Omega_0 & \forall \underline{v} \in \mathcal{V}(\Omega_0) \quad (71) \\ \\ \displaystyle + \int_{\Omega_0} \lambda_0^* \mathrm{div}_{\underline{x}_0} (\tilde{u}) \mathrm{div}_{\underline{x}_0} (\underline{v}) \, d\Omega_0 = \int_{\Omega_0} \varrho_0 \tilde{f} \cdot \underline{v} \, d\Omega_0 & \end{array} \right. $$

where

$$ \begin{align} \tilde{\mathbf{C}}_0 : \underline{s} &:= J_0^{-1} \underline{F}_0 \left( \mathbf{C}_0^d : \underline{F}_0^T \underline{s} : \underline{F}_0 \right) \underline{F}_0^T + 2 p_0 \underline{s} \tag{72} \\ \lambda_0^* &:= J_0^{-1} \lambda^* = J_0 \frac{d^2 \mathcal{W}^v}{dJ^2}(J_0) + \frac{d \mathcal{W}^v}{dJ}(J_0) \nonumber \\ \varrho_0 &:= J_0^{-1} \tilde{\varrho} \nonumber \end{align} $$

with initial condition

$$ \tilde{u}(0) = 0 \qquad \partial_t \tilde{u}(0) = 0 $$