# A time-domain spectral finite element method for acoustoelasticity: modeling the effect of mechanical loading on guided wave propagation

Andre Luiz Dalmora, Alexandre Imperiale, Sébastien Imperiale, Philippe Moireau

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

## To cite this version:

Andre Luiz Dalmora, Alexandre Imperiale, Sébastien Imperiale, Philippe Moireau. A time-domain spectral finite element method for acoustoelasticity: modeling the effect of mechanical loading on guided wave propagation. Wave Motion, 2024, 129, 103328 (23 p.). 10.1016/j.wavemoti.2024.103328. hal-04417255

HAL Id: hal-04417255

https://inria.hal.science/hal-04417255v1

Submitted on 25 Jan 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# A time-domain spectral finite element method for acoustoelasticity: modeling the effect of mechanical loading on guided wave propagation

André Dalmora¹,²,³, Alexandre Imperiale¹✉, Sebastien Imperiale²,³, Philippe Moireau⁵,²,³

¹Université Paris-Saclay, CEA, List, F-91120, Palaiseau, France; ²Inria, Team-MΞDISIM, Inria-Saclay Ile de France 91120 Palaiseau, France; ³LMS, Ecole Polytechnique, CNRS – Institut Polytechnique de Paris

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

For correspondence: alexandre.imperiale@cea.fr

**Funding:** This research was funded by the following project: "GW4SHM"(gw4shm.eu) project from the European Union's Horizon 2020 Research and Innovation program under the Marie Skłodowska-Curie, grant number 860104.

**Present address:** CEA Saclay, LIST - DIN, Bât. 565 - PC 120, 91191 Gif-sur-Yvette Cedex - France

**Competing interests:** The author declare no competing interests.

## Abstract

Ultrasonic testing techniques such as guided wave-based structural health monitoring aim to evaluate the integrity of a material with sensors and actuators that operate in *situ*, i.e. while the material is in use. Since ultrasonic wave propagation is sensitive to environmental conditions such as pre-deformation of the structure, the design and performance evaluation of monitoring systems in this context is a complicated task that requires quantitative data and the associated modeling effort. In our work, we propose a set of numerical tools to solve the problem of mechanical wave propagation in materials subjected to pre-deformation. This type of configuration is usually treated in the domain of acoustoelasticity. A relevant modeling approach is to consider two different problems: a quasi-static nonlinear problem for the large displacement field of the structure and a linearized time-domain wave propagation problem. After carefully reviewing the modeling ingredients to represent the configurations of interest, we propose an original combination of numerical tools that leads to a computationally efficient algorithm. More specifically, we use 3D shell elements for the quasi-static nonlinear problem and the time-domain spectral finite element method to numerically solve the wave propagation problem. Our approach can represent any type of material constitutive law, geometry or mechanical solicitation. We present realistic numerical results on 3D cases related to the monitoring of both isotropic and anisotropic materials, illustrating the genericity and efficiency of our method. We also validate our approach by comparing it to experimental data from the literature.

## 1 Introduction

In many high-end industries, safety regulations include assessments of the condition of structures and materials used for operations. Typical examples of application areas are nuclear energy, petrochemical industry, transportation, or aeronautics. In fact, any industry where the safety of critical components is of paramount importance. Over time, the need to ensure or assess the health of these components has led to the development of advanced Nondestructive Testing (NDT) techniques such as Structural Health Monitoring (SHM). SHM systems aim to evaluate the structural integrity of the material in question using sensors and processing units *in situ*. In
other words, the monitoring systems are attached to the structure once and for all, so that they can continuously generate relevant field data. One means of implementing such systems is to rely on ultrasonic Guided Waves (GW) [50, 65]. These waves have the advantage of propagating over longer distances – compared to bulk waves, which are arguably more common in the NDT community – and thus probe a larger volume of material. Nonetheless, GW propagation is potentially more complex. Namely, it involves dispersion effects that vary from one propagating mode to another, and in the particular case of SHM systems with *in situ* control processes, it is known to be affected by Environmental and Operational Conditions (EOCs) [61]. One of the most general types of EOC that we consider in our work is mechanical solicitations sustained by the structure during its use. These mechanical conditions lead to potential deformations and internal stresses that affect GW propagation. This is the so-called acoustoelastic effect [39, 40]. To provide meaningful assistance in the development and evaluation of SHM system performance, our work addresses the challenge of incorporating acoustoelastic effects into a wave propagation modeling tool that is both generic – *i.e.* applicable to any geometry, constitutive law or mechanical solicitation – and efficient.

There are numerous previous works in the literature that address this issue but, to our knowledge, they are mainly based on semi-analytical approaches. For instance, in [41, 48, 51, 52, 53, 56, 55], analytical formulations are proposed and can be solved at low cost. Nevertheless, due to assumptions regarding the geometry or the deformation, such formulations can be limited when precisely modeling complex and realistic cases. For complex shapes and loading configurations, in [32, 42, 57, 60, 62] the authors propose semi-analytical methods that encompass stress inhomogeneity in one or two dimensions. These methods treat the problem in the frequency domain and efficiently compute dispersion curves of guided waves in specified wave guides, given the stress profile through the thickness. In our work, we aim at proposing a time-domain approach for this problem, while considering arbitrary geometries, hyperelastic laws, and inhomogeneous stress profiles caused by mechanical loading.

Traditionally, due to the different characteristic times and amplitudes of the structural strain, one can identify two (unidirectionally coupled) problems: (1) a quasi-static nonlinear mechanical problem for the structural deformation; (2) a time-dependent wave problem resulting from a linearization procedure around the quasi-static deformation. In this work, we aim to extend the reach of acoustoelastic GW propagation modeling by combining a set of generic numerical tools to solve both problems. For the first problem, the main challenge (apart from mastering the inherent nonlinearities) is that classical finite element methods (FEMs) give inadequate results for thin, elongated geometries due to numerical locking [30, 38]. To overcome this difficulty we use 3D shell elements [38]. They are based on an asymptotic expansion of the displacement field in the thickness direction and, combined with the so-called Mixed Interpolation of Tensorial Components (MITC) method [14, 20, 24] they are locking-free. Moreover, they are readily available for any kind of constitutive law. The output of this first computational stage is the (large and quasi-static) displacement field sustained by the structure. This output is then used as an input for the wave propagation problem. For this second problem, we use a different type of discretization, namely the Spectral Finite Element Method (SFEM) in the time domain [23, 25, 33]. Since the unknown of each problem are fundamentally different in nature, it is natural to resort to different discretization schemes, in particular different meshes. This entails the use of an interpolation operator from the first problem's mesh onto the second one. The SFEM is based on high-order Lagrangian elements defined on Gauss-Lobatto points. This avoids the Runge phenomenon and achieves spectral-like convergence.
Using a consistent mass lumping procedure, they lead to a very efficient explicit fully discrete scheme – with a stability condition that depends on the mechanical deformation of the structure. To illustrate the validity and efficiency of our approach, we provide meaningful 3D simulation results associated with SHM configurations, as well as experimental validation results with data from the literature [41].

This article is divided into two main sections. In section 2 we explain the methodology we use in our work. In particular, we first recall the general formalism of nonlinear elastodynamics, followed by the linearization procedure around a quasi-static displacement field. In this way, we can accurately define the wave propagation problem involving the acoustoelastic effects. After comparing our modeling approach with the relevant previous work on acoustoelasticity, we conclude this section with the dedicated time and space numerical schemes applied in order to solve the one-way coupled problems, namely the quasi-static nonlinear problem and the linearized wave propagation problem. The section 3 is devoted to the detailed 3D numerical illustrations of our approach. First, we provide illuminating examples addressing the potential issue of stability of the fully discrete wave propagation model as a function of the type of deformation and the constitutive law. Then we provide illustrations of SHM configurations, namely a tube of isotropic material subjected to a four-point bending test and a stratified anisotropic plate subjected to torsion. These examples illustrate the performances and versatility of our approach, while in the last part of the section we present experimental validations showing the reliability of our modeling toolchain.

# 2 Methodology for modeling the effect of mechanical loading on elastic wave propagation

## 2.1 Nonlinear elastodynamics with quasi-static surface traction & time dependent ultrasonic actuator

This section is dedicated to deriving a model representing the propagation of “low amplitude – high frequency” waves, within a material subject to a “high amplitude – (very) low frequency” mechanical loading. This model enters the traditional framework of finite deformations of hyperelastic solids. Readers may refer to the reference textbooks [12, 17, 21] and references therein for an exhaustive presentation of this generic formalism. We do not seek here to deviate from the very well-understood formalism of hyperelasticity. Rather, we aim at specifying the corresponding modeling elements, in particular the various forces that act on the material of interest. This aspect is of paramount importance through our work, since it leads to a linearization procedure and its subsequent wave propagation model.

Let us denote by $t \ge 0$ the time variable with $T > 0$ the maximal time of interest, and by $\Omega(t) \subset \mathbb{R}^3$ the volume occupied by the material during its deformation, with boundary $\Gamma(t)$. We assume that there exists $\Gamma_D \subset \Gamma(t)$ a clamped part of the boundary. In the following, we consider the functional space $V(\Omega(t))$ of admissible virtual displacements vanishing on $\Gamma_D$.

### Quasi-static surface traction.

We consider a sub-part of the boundary $\Gamma_0(t) \subset \Gamma$ that is the area of action of a high amplitude and quasi-static force, and such that $\Gamma_0(t) \cap \Gamma_D = \emptyset$. More precisely, let us denote by $f_0$ a prescribed scalar-valued function representing the amplitude of a surface traction. The associated virtual work amounts to

$$ \int_{\Gamma_0(t)} f_0(x) n_0(x, t) \cdot w(x) \, d\Gamma, \quad \forall w \in V(\Omega(t)), $$
where $n_0$ is the outgoing unit normal vector field of $\Gamma_0$. The quasi-static hypothesis represents the fact that the time variation of this surface traction is assumed to be significantly lower than the ultrasonic source term described subsequently.

### Time dependent ultrasonic actuator.

We also consider a third part of the boundary $\Gamma_a(t) \subset \Gamma(t)$ where a low amplitude and time-dependent ultrasonic actuator is vibrating, thus generating propagating wave packets. We assume that $\Gamma_a(t) \cap \Gamma_0(t) = \emptyset$ and $\Gamma_a(t) \cap \Gamma_D = \emptyset$. We denote by $f_a^\delta$ a given source term representing this ultrasonic actuator. The superscript $\delta$ represents the amplitude ratio of the ultrasonic source term compared to $f_0$. In the following, we extensively use the assumption that this ratio is small, i.e.

$$ \delta = \frac{\sup_{(x,t) \in \Omega \times [0;T]} \|f_a^\delta(x,t)\|_3}{\sup_{x \in \Omega} |f_0(x)|} \ll 1. \quad (1) $$

In the following, we assume – without loss of generality – that there exists a vector field $f_a$ that corresponds to a rescaled source term, namely $f_a^\delta = \delta f_a$. The associated virtual work becomes

$$ \delta \int_{\Gamma_a(t)} f_a(x, t) \cdot w(x) \, d\Gamma, \quad \forall w \in \mathcal{V}(\Omega(t)), \quad \forall t \in [0; T]. $$

### Weak formulations of the nonlinear elastodynamic problem.

Let us introduce $\rho$ the mass density of the material. With these notations and the previous definitions of the various boundary conditions, the principle of virtual works states that the velocity field $v(t)$ and the Cauchy stress tensor $\sigma(t)$ satisfy, for any $w \in \mathcal{V}(\Omega(t))$

$$ \int_{\Omega(t)} \rho (\partial_t v + (\nabla_x v) \cdot v) \cdot w \, d\Omega + \int_{\Omega(t)} \sigma : \nabla_x w \, d\Omega \\ = \int_{\Gamma_0(t)} f_0 n_0 \cdot w \, d\Gamma + \delta \int_{\Gamma_a(t)} f_a \cdot w \, d\Gamma, \quad (2) $$

with given initial conditions $v|_{t=0}$ and $\sigma|_{t=0}$. Albeit being the standard weak form of the equilibrium relation, relation (2) deserves some comments. It is not an explicit form because the deformed domain $\Omega$ is unknown, nor is it a closed form since it lacks the constitutive law linking the Cauchy stress tensor with the kinematic unknown. It should be noted that this constitutive law is rarely written on a deformed configuration, but rather on a so-called reference configuration. Also, from a purely practical viewpoint, the definition of the source term $f_a$ representing the ultrasonic actuator is more likely to be available on this reference configuration.

From this line of arguments, it is natural to derive a formulation written on a reference coordinate system. Let us denote by $\hat{\Omega}$ the volume occupied by the material in this reference setting, and by $\xi \in \hat{\Omega}$ the associated coordinates. For any time $t > 0$, the deformation is represented by the bijective mapping,

$$ \phi(t) : \hat{\Omega} \longrightarrow \Omega(t) \\ \xi \longrightarrow x(\xi, t) = \phi(\xi, t), $$

and we define the displacement field $u(t)$ as

$$ u(\xi, t) = x(\xi, t) - \xi, \quad \forall \xi \in \hat{\Omega}. $$
<figure><img src="image_4.png" /><figcaption>Figure 1. Reference and deformed configurations and their associated notations.</figcaption></figure>

We denote by $\hat{\Gamma}$ the boundary of $\hat{\Omega}$, and $\hat{\Gamma}_D$, $\hat{\Gamma}_0$, and $\hat{\Gamma}_a$ the sub-parts of the boundary corresponding to $\Gamma_D$, $\Gamma_0(t)$, and $\Gamma_a(t)$ respectively. See Figure 1 for an illustration of both coordinate systems and domains.

The deformation gradient and its determinant are denoted by

$$F(\xi, t) = \nabla_{\xi}\phi(\xi, t) = I + \nabla_{\xi}u(\xi, t) \quad \text{and} \quad J(\xi, t) = \det F(\xi, t), \quad \forall \xi \in \hat{\Omega},$$

respectively. With these notations, we derive the equivalent of (2) set on the reference configuration, thus explicitly showing the dependency of the displacement field – which is the chosen kinematic unknown in the following. To do so, we apply a standard change of variables from $\Omega$ to $\hat{\Omega}$ – see [21] and references therein for more details. Upon the straightforward relation $v \circ \phi = \partial_t u$, we start with the inertial terms, namely, for all $w \in \mathcal{V}(\hat{\Omega})$,

$$\int_{\Omega(t)} \varrho (\partial_t v + (\nabla_x v) \cdot v) \cdot (w \circ \phi^{-1}) \, d\Omega = \int_{\hat{\Omega}} J (\varrho \circ \phi) \partial_{tt}^2 u \cdot w \, d\hat{\Omega} = \int_{\hat{\Omega}} \hat{\varrho} \partial_{tt}^2 u \cdot w \, d\hat{\Omega},$$

where $\hat{\varrho} = (\varrho \circ \phi)J$ is the reference mass density. Concerning the stiffness term, we have, for all $w \in \mathcal{V}(\hat{\Omega})$,

$$\int_{\Omega(t)} \sigma : \nabla_x (w \circ \phi^{-1}) \, d\Omega = \int_{\hat{\Omega}} (\sigma \circ \phi) : (\nabla_{\xi} w \cdot F^{-1}) \, J d\hat{\Omega} = \int_{\hat{\Omega}} T : \nabla_{\xi} w \, d\hat{\Omega},$$

where $T = (\sigma \circ \phi) \cdot F^{-1}J$ is the first Piola-Kirchhoff tensor. Upon the assumption of hyperelasticity [17, 21] there exists a potential function $\mathcal{W}(F)$ that depends exclusively on the deformation gradient and such that

$$T(F) = D_F \mathcal{W}(F).$$

Note that prescribing a closed form of this energy potential function actually defines the constitutive relation of the material – a matter addressed in the subsequent section. The first Piola-Kirchhoff tensor is non-symmetric, and in the following, we resort to its symmetric counterpart, which can be introduced as follows,

$$\int_{\hat{\Omega}} T : \nabla_{\xi} w \, d\hat{\Omega} = \int_{\hat{\Omega}} (F^{-1} \cdot T) : (F^{\top} \cdot \nabla_{\xi} w) \, d\hat{\Omega} = \int_{\hat{\Omega}} \Sigma : D_u e(u) w \, d\hat{\Omega},$$

with $\Sigma = F^{-1} \cdot T$ the (symmetric) second Piola-Kirchhoff tensor, and

$$e(u) = \frac{1}{2} (\nabla_{\xi} u + \nabla_{\xi} u^{\top} + \nabla_{\xi} u^{\top} \cdot \nabla_{\xi} u), \quad D_u e(u) w = \frac{1}{2} (F^{\top} \cdot \nabla_{\xi} w + \nabla_{\xi} w^{\top} \cdot F),$$
the Green-Lagrange tensor and its differential w.r.t. the displacement field. For hyperelastic materials, the second Piola-Kirchhoff tensor is usually defined exclusively as a function of the Green-Lagrange tensor, namely $\Sigma = \Sigma(e(\mathbf{u}))$ in the previous expression, and

$$ \Sigma(e) = D_e \mathcal{W}(e). $$

For the quasi-static surface traction, one resorts to the formula for transport of normal vector field to derive the following relations

$$ \forall \mathbf{w} \in \mathcal{V}(\hat{\Omega}), \quad \int_{\Gamma_0(t)} f_0 \mathbf{n}_0 \cdot (\mathbf{w} \circ \phi^{-1}) \, d\Gamma = \int_{\hat{\Gamma}_0} \hat{f}_0 (\co F \cdot \hat{\mathbf{n}}_0) \cdot \mathbf{w} \, d\hat{\Gamma}, $$

where $\hat{f}_0 = f_0 \circ \phi$, $\hat{\mathbf{n}}_0$ is the outgoing unit normal vector field of $\hat{\Gamma}_0$, and $\co F = JF^{-\top}$ is the co-factor of the deformation gradient. Finally, concerning the ultrasonic source term, we have

$$ \forall \mathbf{w} \in \mathcal{V}(\hat{\Omega}), \quad \int_{\Gamma_a(t)} \mathbf{f}_a \cdot (\mathbf{w} \circ \phi^{-1}) \, d\Gamma = \int_{\hat{\Gamma}_a} \hat{\mathbf{f}}_a \cdot \mathbf{w} \, \hat{\gamma} \, d\hat{\Gamma}, $$

where $\hat{\gamma} = \|\co F \cdot \hat{\mathbf{n}}_a\|_3$ with $\hat{\mathbf{n}}_a$ the unit normal vector field of $\hat{\Gamma}_a$. This term represents the effect of the deformation on the surface measure. In the previous relation, $\hat{\mathbf{f}}_a$ is the time and space dependency of the actuator on the reference configuration. Its expression depends on the chosen modeling of the actuator. Without being completely exhaustive, we can at least provide three main cases of interest:

* The ultrasonic source is independent of the deformation of the material. For instance, one could imagine configurations where a hydrophone is emitting an ultrasonic beam toward the immersed specimen of interest, subject to the deformation. In that case, we simply have $\hat{\mathbf{f}}_a = \mathbf{f}_a \circ \phi$.

* The actuator is considered to be generating a time-dependent source term with polarization in the normal direction of the material. This corresponds to the case where there exists a scalar-valued function $f_a$ such that $\mathbf{f}_a = f_a \mathbf{n}_a$, where $\mathbf{n}_a$ is the outgoing unit normal vector field of $\Gamma_a$. In that case, we define

$$ \hat{\mathbf{f}}_a = (f_a \circ \phi) \frac{\co F \cdot \hat{\mathbf{n}}_a}{\hat{\gamma}}. $$

* The polarization lies within the local tangent plane of the surface $\Gamma_a$. More precisely if we denote by $(\tau_{1,a}, \tau_{2,a})$ two vectors field generating this tangent plane, then there exists two scalar-valued functions $f_{1,a}$ and $f_{2,a}$ such that $\mathbf{f}_a = f_{1,a} \tau_{1,a} + f_{2,a} \tau_{2,a}$, entailing

$$ \hat{\mathbf{f}}_a = \sum_{i=1,2} (f_{i,a} \circ \phi) \frac{\mathbf{F} \cdot \hat{\tau}_{i,a}}{\|\mathbf{F} \cdot \hat{\tau}_{i,a}\|_3}, \quad \tau_{i,a} \circ \phi = \mathbf{F} \cdot \hat{\tau}_{i,a}, \quad \forall i \in \{1,2\}. $$

We can now give the weak formulation of the problem of interest written in the reference configuration: for any time $t \in [0; T]$, we aim at finding the displacement field $\mathbf{u}(t)$ such that, for any test function $\mathbf{w} \in \mathcal{V}(\hat{\Omega})$, we have

$$ \int_{\hat{\Omega}} \hat{\rho} \partial_{tt}^2 \mathbf{u} \cdot \mathbf{w} \, d\hat{\Omega} + \int_{\hat{\Omega}} \Sigma(e(\mathbf{u})) : D_u e(\mathbf{u}) \mathbf{w} \, d\hat{\Omega} \\ = \int_{\hat{\Gamma}_0} \hat{f}_0 (\co F \cdot \hat{\mathbf{n}}_0) \cdot \mathbf{w} \, d\hat{\Gamma} + \delta \int_{\hat{\Gamma}_a} \hat{\mathbf{f}}_a \cdot \mathbf{w} \, \hat{\gamma} \, d\hat{\Gamma}, \quad (3) $$

<footnote n="173">173 with given initial conditions $\mathbf{u}|_{t=0}$ and $\partial_t \mathbf{u}|_{t=0}$.</footnote>
## 2.2 Linearization around a quasi-static displacement field

Upon the assumption of a quasi-static mechanical loading and a low amplitude time-dependent ultrasonic source term – as depicted in (1) – we propose to decompose the solution of (3) following the Ansatz

$$u(\xi, t) = u_0(\xi) + \delta u_1(\xi, t) + O(\delta^2), \quad \forall \xi \in \hat{\Omega}, \quad \forall t \in [0; T]. \quad (4)$$

This (natural) decomposition amounts to seeking a large quasi-static deformation $u_0$ generated by the surface traction term, and a low amplitude time-dependent contribution $u_1$ emanating from the ultrasonic source term. Intuitively, one would expect $u_0$ to satisfy a nonlinear static problem, and $u_1$ a wave equation. In this section, we follow standard linearization arguments to derive these two problems. Moreover, following the existence and uniqueness results given in [17], we define the space

$$\mathcal{V}(\hat{\Omega}) = \{w \in [H^1(\hat{\Omega})]^3 \mid w|_{\Gamma_D} = 0\},$$

and the space $\mathcal{W}(\hat{\Omega})$ of admissible displacements in the reference configuration

$$\mathcal{W}(\hat{\Omega}) = W^{2,p}(\hat{\Omega})^3 \cap \mathcal{V}(\hat{\Omega}), \text{ for some } p > 3. \quad (5)$$

We now comment briefly this choice. Imbeddings results (see [27, 17]) show that, for any $w^* \in \mathcal{W}(\hat{\Omega})$, first, $\nabla_{\xi} w^*$ belongs to $W^{1,p}(\hat{\Omega})^{3 \times 3}$ that is an algebra (in particular products of functions in this space do belong to the same space), second, $\nabla_{\xi} w^*$ is continuous and bounded up to the boundary. We deduce that

$$e(w^*) \in W^{1,p}(\hat{\Omega})^{3 \times 3} \quad \text{and} \quad \forall w \in \mathcal{V}(\hat{\Omega}), \quad D_u e(w^*) w \in L^2(\hat{\Omega})^{3 \times 3}.$$

Finally, assuming sufficient smoothness of the potential function $\mathcal{W}(F)$ it can be inferred that $\Sigma(e(w^*)) \in W^{1,p}(\hat{\Omega})^{3 \times 3}$ hence, it is continuous and bounded. A similar smoothness property holds for the fourth order tensor $D_e \Sigma(e(w^*))$. The properties given above are implicitly used to give meanings to the integrals introduced in the section below.

### Stiffness and surface traction operators and associated differentials.
Let us denote by

$$A : \mathcal{W}(\hat{\Omega}) \to \mathcal{V}(\hat{\Omega})' \\ w^* \to A(w^*),$$

the nonlinear stiffness operator such that, for any $w \in \mathcal{V}(\hat{\Omega})$,

$$\langle A(w^*), w \rangle = \int_{\hat{\Omega}} \Sigma(e(w^*)) : D_u e(w^*) w \, d\hat{\Omega},$$

where $\langle \cdot, \cdot \rangle$ denotes the duality product in $\mathcal{V}(\hat{\Omega})$. The Gâteaux differential of A at $w^*$ is denoted by $D_u A(w^*) \in \mathcal{L}(\mathcal{V}(\hat{\Omega}), \mathcal{V}(\hat{\Omega})')$, and reads, for any $v, w \in \mathcal{V}(\hat{\Omega})$

$$\langle D_u A(w^*) v, w \rangle = \int_{\hat{\Omega}} D_u e(w^*) v : D_e \Sigma(e(w^*)) : D_u e(w^*) w \, d\hat{\Omega} \\ + \int_{\hat{\Omega}} \Sigma(e(w^*)) : D_u^2 e(v, w) \, d\hat{\Omega}, \quad (6)$$

where

$$D_u^2 e(v, w) = \frac{1}{2} (\nabla_{\xi} v^T \cdot \nabla_{\xi} w + \nabla_{\xi} w^T \cdot \nabla_{\xi} v).$$
Note that the differential (6) is obtained from the differential of a product combined with a chain rule, and the definition of the Green-Lagrange tensor. It actually defines a bilinear form of particular interest in the following. More precisely, for a given $\mathbf{w}^* \in W(\hat{\Omega})$ defining

$$ a(\mathbf{w}^*; \cdot, \cdot) : V(\hat{\Omega}) \times V(\hat{\Omega}) \mapsto \mathbb{R} \\ (\mathbf{v}, \mathbf{w}) \mapsto a(\mathbf{w}^*; \mathbf{v}, \mathbf{w}) = \langle D_u A(\mathbf{w}^*) \mathbf{v}, \mathbf{w} \rangle, \quad (7) $$

then one can verify that $a(\mathbf{w}^*; \cdot, \cdot)$ is a symmetric bilinear form. Note that $D_u A$ is a contribution of two terms: the first one, involving $D_u \Sigma$, is related to constitutive nonlinearities since it involves second order derivatives of the potential $W$ with respect to the Green-Lagrange tensor $e$; the second one, involving $\Sigma$, is related to geometrical nonlinearities since it comes exclusively from the linearization of the Green-Lagrange tensor. Another more compact form – but less prone to interpretation – of the Gâteaux differential of A can be deduced from the expression

$$ \langle A(\mathbf{w}^*), \mathbf{w} \rangle = \int_{\hat{\Omega}} \mathbf{T}(F(\mathbf{w}^*)) : \nabla_{\xi} \mathbf{w} \, d\hat{\Omega}, $$

where $F(\mathbf{w}^*) = \mathbf{I} + \nabla_{\xi}\mathbf{w}^*$. Computing the Gâteaux differential of A at $\mathbf{w}^*$ gives the alternative expression

$$ \begin{aligned} a(\mathbf{w}^*; \mathbf{v}, \mathbf{w}) = \langle D_u A(\mathbf{w}^*) \mathbf{v}, \mathbf{w} \rangle &= \int_{\hat{\Omega}} \nabla_{\xi} \mathbf{v} : D_F \mathbf{T}(F(\mathbf{w}^*)) : \nabla_{\xi} \mathbf{w} \, d\hat{\Omega} \\ &= \int_{\hat{\Omega}} \nabla_{\xi} \mathbf{v} : D_F^2 W(F(\mathbf{w}^*)) : \nabla_{\xi} \mathbf{w} \, d\hat{\Omega}. \end{aligned} \quad (8) $$

In the following, we also denote by

$$ \begin{aligned} \mathbf{B} : W(\hat{\Omega}) &\mapsto V(\hat{\Omega})' \\ \mathbf{w}^* &\mapsto \mathbf{B}(\mathbf{w}^*), \end{aligned} $$

the nonlinear surface traction operator such that, for any $\mathbf{w} \in V(\hat{\Omega})$,

$$ \langle \mathbf{B}(\mathbf{w}^*), \mathbf{w} \rangle = \int_{\hat{\Gamma}_0} \hat{f}_0 (\co F(\mathbf{w}^*) \cdot \hat{n}_0) \cdot \mathbf{w} \, d\hat{\Gamma}. $$

The Gâteaux differential of B at $\mathbf{w}^*$ is denoted by $D_u B(\mathbf{w}^*) \in \mathcal{L}(V(\hat{\Omega}), V(\hat{\Omega})')$, and reads, for any $\mathbf{v}, \mathbf{w} \in V(\hat{\Omega})$

$$ \langle D_u B(\mathbf{w}^*) \mathbf{v}, \mathbf{w} \rangle = \int_{\hat{\Gamma}_0} \hat{f}_0 ((D_u \co F(\mathbf{w}^*)) \mathbf{v}) \cdot \hat{n}_0) \cdot \mathbf{w} \, d\hat{\Gamma}, \quad (9) $$

where

$$ D_u \co F(\mathbf{w}^*) \mathbf{v} = ((F^{-\top}(\mathbf{w}^*) : \nabla_{\xi} \mathbf{v}) \mathbf{I} - F^{-\top}(\mathbf{w}^*) \cdot \nabla_{\xi} \mathbf{v}) \cdot \co F(\mathbf{w}^*). $$

Note that the expression of the differential of the cofactor matrix derives from the following relations

$$ D_u J(\mathbf{w}^*) \mathbf{v} = \co F(\mathbf{w}^*) : \nabla_{\xi} \mathbf{v}, \quad D_u (F^{-\top}(\mathbf{w}^*)) \mathbf{v} = -F^{-\top}(\mathbf{w}^*) \cdot \nabla_{\xi} \mathbf{v} \cdot F^{-\top}(\mathbf{w}^*). $$

Compared to (6), no particular property – such as symmetry – can be derived from the bilinear form associated to the differential (9). In order to link with the definition of the weak formulation (3), we introduce the difference of these two nonlinear operators, namely

$$ \begin{aligned} A_b : W(\hat{\Omega}) &\mapsto V(\hat{\Omega})' \\ \mathbf{w}^* &\mapsto A_b(\mathbf{w}^*) = A(\mathbf{w}^*) - B(\mathbf{w}^*), \end{aligned} \quad (10) $$

such that its associated Gâteaux differential, for any $\mathbf{w}^* \in W(\hat{\Omega})$, simply reads

$$ D_u A_b(\mathbf{w}^*) = D_u A(\mathbf{w}^*) - D_u B(\mathbf{w}^*) \in \mathcal{L}(V(\hat{\Omega}), V(\hat{\Omega})'). $$
212 Quasi-static and linearized wave problems.

Let us denote by $\mathcal{H}(\hat{\Omega}) = [L^2(\hat{\Omega})]^3$ and the associated mass-weighted inner product by

$$ (v, w)_{\mathcal{H}} = \int_{\hat{\Omega}} \hat{\rho} v \cdot w \, d\hat{\Omega}, \quad \forall v, w \in \mathcal{H}. $$

We introduce the nonlinear operator

$$ L(t, \cdot) : \begin{aligned} & \mathcal{W}(\hat{\Omega}) \mapsto \mathcal{V}(\hat{\Omega})' \\ & w^* \mapsto L(t, w^*), \end{aligned} $$

213 such that, for any $w \in \mathcal{V}(\hat{\Omega})$,

$$ \langle L(t, w^*), w \rangle = \int_{\hat{\Gamma}_a} \hat{f}_a(w^*, t) \cdot w \hat{\gamma}(w^*) \, d\hat{\Gamma}. \quad (11) $$

From the definition of the operators (10) and (11), we remark that the weak formulation (3) can be written in the following form

$$ \frac{d^2}{dt^2}(u, w)_{\mathcal{H}} + \langle A_b(u), w \rangle = \delta\langle L(t, u), w \rangle, \quad \forall w \in \mathcal{V}(\hat{\Omega}). $$

Using the Ansatz (4), and the time invariance of $u_0$ leads to

$$ \langle A_b(u_0), w \rangle + \delta \left( \frac{d^2}{dt^2}(u_1, w)_{\mathcal{H}} + \langle D_u A_b(u_0) u_1, w \rangle - \langle L(t, u_0), w \rangle \right) = O(\delta^2). $$

214 Hence, satisfying the relation (3) up to second order terms w.r.t $\delta$ amounts to satisfying the following two problems:

215
216
i) Find $u_0 \in W(\hat{\Omega})$ the solution, for any $w \in \mathcal{V}(\hat{\Omega})$, of
$$ \langle A_b(u_0), w \rangle = 0. $$

218
219
ii) Find $(0, T) \ni t \mapsto u_1(t) \in \mathcal{V}(\hat{\Omega})$ satisfying in $\mathcal{D}'(0, T)$ - a.k.a in the sense of distribution:

$$ \frac{d^2}{dt^2}(u_1(\cdot), w)_{\mathcal{H}} + \langle D_u A_b(u_0)u_1(\cdot), w \rangle = \langle L(\cdot, u_0), w \rangle, \quad w \in \mathcal{V}(\hat{\Omega}) \quad (12) $$

220 with given initial conditions $u_1(0) \in \mathcal{H}((\hat{\Omega}))$ and $\frac{du_1}{dt}(0) \in \mathcal{V}(\hat{\Omega})'$.

221 Approximated linearized wave problem.

222 A major drawback of the formulation (12) of the wave problem is that it contains the tangent of the surface traction operator, namely $D_u B(u_0)$, defined in (9). It is associated with a non-symmetric bilinear form, with an undefined sign and implying the surface gradient of its first argument. In addition to these unfortunate mathematical properties, one can interpret this term as the effect of the surface traction onto the wave field. In the applications of interest, where the loading area is restricted compared to the complete surface domain, we expect this term to bear little effect. Combining these arguments, we neglect the contribution of the tangent of the surface traction operator in the wave propagation problem, and we consider the following pair of problems:

231 i) Find $u_0 \in \mathcal{V}(\hat{\Omega})$ the solution, for any $w \in \mathcal{V}(\hat{\Omega})$, of

$$ \langle A_b(u_0), w \rangle = 0. \quad (13) $$
$$ \text{ii) Find } (0, T) \ni t \mapsto \tilde{u}(t) \in \mathcal{V}(\hat{\Omega}) \text{ satisfying in } \mathcal{D}'(0, T): \frac{d^2}{dt^2}(\tilde{u}(\cdot), w)_{\mathcal{H}} + a(\mathbf{u}_0; \tilde{u}(\cdot), w) = \langle L(\cdot, \mathbf{u}_0), w \rangle, \quad \forall w \in \mathcal{V}(\hat{\Omega}) \quad (14) $$

with given initial conditions $\tilde{u}(0) \in \mathcal{H}((\hat{\Omega}))$ and $\frac{d\tilde{u}}{dt}(0) \in \mathcal{V}(\hat{\Omega})'$, and $a(\mathbf{u}_0; \cdot, \cdot)$ defined in (7).

One important remark that can be made at this point is that the two problems previously defined can be addressed successively. From a practical viewpoint this means that one can first solve (13), extract the solution $\mathbf{u}_0$ to compute the tangent operator and the source term in (14) and solve the wave equation to obtain $\tilde{u}$.

Well-posedness of the wave propagation problem and energy relation. Since the focus of our work is on wave propagation modeling, we do not discuss the questions of the existence and uniqueness of the nonlinear static problem (13). Theoretical elements on this matter can be found for instance in [12, 17], they usually entail a so-called “poly-convexity” property of $\mathcal{W}$. We assume in the following that there exists at least one solution of (13), and rather discuss the well-posedness of the wave propagation problem of interest (14). We rely on the standard abstract results that can be found e.g. in the reference textbooks [9, 11, 19]. In particular, since the symmetry of the bilinear form is easily verified from its definition given in (7), we only need the following assumptions to apply Theorem 4.1 in [11]:

* The right-hand side of (14) satisfies $L(\cdot, \mathbf{u}_0) \in \mathcal{C}^1([0, T]; \mathcal{V}(\hat{\Omega})')$,

* The bilinear form $a(\mathbf{u}_0; \cdot, \cdot)$ defined in (7) is continuous and coercive, i.e. there exists two strictly positive constants $C_0$ and $\alpha_0$ such that

$$ |a(\mathbf{u}_0; v, w)| \le C_0 \|v\|_{\mathcal{V}} \|w\|_{\mathcal{V}}, \quad a(\mathbf{u}_0; v, v) \ge \alpha_0 \|v\|_{\mathcal{V}}^2, \quad \forall v, w \in \mathcal{V}(\hat{\Omega}). \quad (15) $$

Upon these assumptions the problem (14) admits a unique solution such that

$$ \tilde{u} \in L^2((0, T); \mathcal{V}(\hat{\Omega})), \quad \frac{d\tilde{u}}{dt} \in L^2((0, T), \mathcal{H}), \quad \frac{d^2\tilde{u}}{dt^2} \in L^2((0, T), \mathcal{V}(\hat{\Omega})') $$

with additionally, from continuous injection results,

$$ \tilde{u} \in \mathcal{C}^0([0, T]; \mathcal{H}(\hat{\Omega})) \cap L^\infty([0, T]; \mathcal{V}(\hat{\Omega})), \quad \frac{d\tilde{u}}{dt} \in \mathcal{C}^0([0, T]; \mathcal{V}(\hat{\Omega})') \cap L^\infty([0, T]; \mathcal{H}(\hat{\Omega})). $$

Note that the coercivity assumption could be replaced by a weaker hypothesis, assuming the coercivity up to a term proportional to $\|v\|_{\mathcal{H}}^2$. It is clear from the expression of the bilinear form that the coercivity constant depends intricately on the choice of the potential function $\mathcal{W}$ and on the linearization point $\mathbf{u}_0$. In particular, in section 3 we provide illuminating examples of numerical solutions where, depending on the solution of the static problem, we observe a loss of coercivity leading to an ill-posed wave propagation problem. In addition to the previous well-posedness result, we recall the standard energy relation which is an important framework to devise a stable time discretization procedure for the problem (14). Let us denote the kinematic, potential and total energies associated with the unique solution of (14) by

$$ \mathcal{K}(t) = \frac{1}{2} \| \frac{d\tilde{u}}{dt} \|_{\mathcal{H}}^2, \quad \mathcal{P}(t) = \frac{1}{2} a(\mathbf{u}_0; \tilde{u}, \tilde{u}), \quad \mathcal{E} = \mathcal{K} + \mathcal{P}. $$

Assuming a null source term in (14), namely $L(t; \mathbf{u}_0) = 0$, then, by using formally the time derivative of the solution as a test function in (14), we recover the following energy conservation

$$ \frac{d}{dt} \mathcal{E}(t) = 0. \quad (16) $$
258 Energy potentials.

259 Our numerical investigations mainly concern frame-indifferent hyperelastic material with no residual stress. The material being frame-indifferent the hyperelastic potential 260 – assumed smooth – is only dependent on the Green-Lagrange tensor e 261

$$
\mathcal{W}(\mathbf{F}) \equiv \mathcal{W}(e).
$$

262 First, we consider isotropic materials. Since the material has no residual stress, we 263 necessarily have (see [17]),

$$
\mathcal{W}(e) = \frac{\lambda}{2}(\operatorname{tr} e)^2 + \mu \operatorname{tr} e^2 + \alpha_1(\operatorname{tr} e)^3 + \alpha_2(\operatorname{tr} e)\operatorname{tr} e^2 + \alpha_3\operatorname{tr} e^3 + o(|e|^3). \quad (17)
$$

264 Note that in the expression above, if the material is heterogeneous, the parameter $\lambda, \mu$ and the $\{\alpha_i\}_{i=1}^3$ may depend on $\xi$. The coefficients associated with third-order terms, 265 namely $\{\alpha_i\}_{i=1}^3$, are usually referred to as the Third Order Elastic Constants (TOECs). 266 The two first terms of this expansion correspond to the so-called Saint-Venant-Kirchhoff 267 (SVK) potential 268

$$
\mathcal{W}^{\text{SVK}}(e) = \frac{\lambda}{2}(\operatorname{tr} e)^2 + \mu \operatorname{tr} e^2,
$$

269 while, neglecting higher-order terms, yields the potential

$$
\mathcal{W}^{3\mathrm{rd}}(e) = \mathcal{W}^{\mathrm{SVK}}(e) + \alpha_1(\mathrm{tr} \, e)^3 + \alpha_2(\mathrm{tr} \, e)\mathrm{tr} \, e^2 + \alpha_3\mathrm{tr} \, e^3, \quad (18)
$$

such potential can be written using Landau’s coefficients [15]

$$
A = 3\alpha_3, \quad B = \alpha_2 \quad \text{and} \quad C = 3\alpha_1.
$$

270 Note that these two potentials are not poly-convex in general, hence stability of the 271 nonlinear mechanical problem – and therefore of its linearized version – may be an 272 issue. The most simple elastic potential mitigating this problem is the Compressible 273 Neo-Hookean (CNH) potential,

$$
\mathcal{W}^{\text{CNH}}(e) = \frac{\lambda}{2} \left( \sqrt{\det(I + 2e)} - 1 \right)^2 + \mu \operatorname{tr} e - \frac{\mu}{2} \log \left( \det(I + 2e) \right).
$$

274 The potential $\mathcal{W}^{\text{CNH}}$ satisfies the following expansion,

$$
\mathcal{W}^{\text{CNH}}(e) = \frac{\lambda}{2}(\operatorname{tr} e)^2 + \mu \operatorname{tr} e^2 + \frac{\lambda}{2}(\operatorname{tr} e)^3 - \lambda(\operatorname{tr} e)\operatorname{tr} e^2 - \frac{4\mu}{3}\operatorname{tr} e^3 + o(|e|^3),
$$

275 showing that, up to third order term in $|e|^3$ the CNH constitutive law corresponds to 276 (18) with $\alpha_1 = \lambda/2$, $\alpha_2 = -\lambda$ and $\alpha_3 = -4\mu/3$. The hyperelastic potential of frame- 277 indifferent isotropic material can be re-written in terms of the first three invariants of 278 the Cauchy-Green deformation tensor $C = F^T F$,

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

279 It can be verified that

$$
\mathcal{W}^{\text{SVK}} = \frac{(\lambda + 2\mu)}{8} I_1^2 - \frac{(3\lambda + 2\mu)}{4} I_1 - \frac{\mu}{2} I_2 + \frac{(9\lambda + 6\mu)}{8} \quad (19)
$$

280 and

$$
\mathcal{W}^{\text{CNH}} = \frac{\lambda}{2}(\sqrt{I_3} - 1)^2 + \frac{\mu}{2}(I_1 - 3 - \log(I_3)). \quad (20)
$$
Murnaghan [2] proposed different TOECs to characterize the second-order elastic deformations, namely $l, m$ and $n$,

$$ \mathcal{W}^{3\text{rd}} = \mathcal{W}^{\text{SVK}} + \frac{l}{24}(I_1 - 3)^3 + \frac{m}{12}(I_1 - 3)(I_1^2 - 3I_2) + \frac{n}{8}(I_1 - I_2 + I_3 - 1). \quad (21) $$

In the context of acoustoelasticity, Hughes and Kelly related these TOECs with changes in bulk wave velocities for prestressed isotropic materials and experimentally assessed them in [3]. The relations between $(l, m, n)$ and $(\alpha_1, \alpha_2, \alpha_3)$ appearing in the expression (18) are given by

$$ \alpha_1 = l - \frac{m}{3} + \frac{n}{6}, \quad \alpha_2 = m - \frac{n}{2}, \quad \alpha_3 = \frac{n}{3}. $$

We also investigate an orthotropic configuration where one plane of symmetry is lost, orthogonal to a principal axis $a$, in that case, the material is said to be transversely isotropic. To fully describe its constitutive behavior we need to add two more invariants

$$ I_4 = a \cdot C \cdot a, \quad I_5 = a \cdot C^2 \cdot a. $$

One example of the hyperelastic potential for transversely isotropic materials is the one from [22],

$$ \mathcal{W}^{\text{CNH-TI}} = \mathcal{W}^{\text{CNH}} + [\alpha + \beta \ln J + \gamma(I_4 - 1)](I_4 - 1) - \frac{\alpha}{2}(I_5 - 1), \quad (22) $$

where $\alpha, \beta$ and $\gamma$ are additional parameters that can be calibrated to match observed transversely isotropic behavior.

### Acoustoelastic effects.

The change in the behavior of elastic waves when propagating in a deformed or prestressed medium is known as acoustoelastic effect. This phenomenon is extensively studied, see e.g. [1, 2, 3]. In [13], the author analyses the structure and properties of incremental elastic deformations models, including the unicity and stability aspects with respect to the constitutive law. As far as isotropic materials are considered, the behavior of elastic waves – bulk wave velocities, in particular – has been mathematically studied by the authors in [39]. We can relate their considerations to stability properties and positivity of the bilinear form $a(u_0; \cdot, \cdot)$. Indeed, from the expression (8), assuming that the solution behaves locally as a plane wave, with polarization $m$ propagating in the direction $n$, and assuming smoothness of the parameters and of $u_0$, we arrive at the necessary condition for the positivity of $a(u_0; \cdot, \cdot)$,

$$ \forall (m, n) \in \mathbb{R}^3 \times \mathbb{R}^3, \quad m \otimes n : D_F^2 \mathcal{W}(F(u_0)) : m \otimes n \ge 0. \quad (23) $$

Note that the inequality above should be verified for each $\xi \in \tilde{\Omega}$ since $u_0$ is a function of $\xi$ even when the material is homogeneous. Now, following [39], this amounts to checking the positivity of the so-called (symmetric) Christoffel tensor $Q(u_0, n) \in \mathbb{R}^{3 \times 3}$, indeed,

$$ \forall (m, n) \in \mathbb{R}^3 \times \mathbb{R}^3, \quad m \otimes n : D_F^2 \mathcal{W}(F(u_0)) : m \otimes n = m \cdot Q(u_0, n) \cdot m. $$

To further analyse it, we use the Christoffel tensor $Q(u_0, n)$ given in [39] in the isotropic case and for the reference configuration. It depends on $\mathcal{W}_i$ and $\mathcal{W}_{ij}$, the first and second-order derivatives of the potential with respect to the invariant $I_i$ and $I_j$ respectively. Note that the dependence in $u_0$ of $\mathcal{W}_i$ and $\mathcal{W}_{ij}$ is temporarily omitted for the sake of clarity. In the case where

$$ \mathcal{W}_{12} = \mathcal{W}_{22} = \mathcal{W}_{13} = \mathcal{W}_{23} = 0, $$
which is verified for the Saint-Venant-Kirchhoff and Compressible Neo-Hookean constitutive laws, the Christoffel tensor is given by

$$ \begin{aligned} Q(\mathbf{u}_0, \mathbf{n}) = & (4I_3^2\mathcal{W}_{33} + 2I_3\mathcal{W}_3)\hat{\mathbf{n}} \otimes \hat{\mathbf{n}} + 2(\mathcal{W}_1 + I_1\mathcal{W}_2)\mathbf{I} \\ & + (4\mathcal{W}_{11} + 2\mathcal{W}_2)(F(\mathbf{u}_0) \cdot \mathbf{n}) \otimes (F(\mathbf{u}_0) \cdot \mathbf{n}) - 2\mathcal{W}_2(\|F(\mathbf{u}_0) \cdot \mathbf{n}\|_3^2\mathbf{I} + \mathbf{B}(\mathbf{u}_0)), \end{aligned} \quad (24) $$

where $\mathbf{B}(\mathbf{u}_0) = F(\mathbf{u}_0)F^\top(\mathbf{u}_0)$ and $\hat{\mathbf{n}} = F^{-\top}(\mathbf{u}_0) \cdot \mathbf{n}$. The expression above takes a simple form when considering the CNH constitutive law

$$ Q^{\text{CNH}}(\mathbf{u}_0, \mathbf{n}) = (\lambda I_3(\mathbf{u}_0) + \mu)\hat{\mathbf{n}} \otimes \hat{\mathbf{n}} + \mu\mathbf{I}. $$

When considering the SVK constitutive law, we get

$$ \begin{aligned} Q^{\text{SVK}}(\mathbf{u}_0, \mathbf{n}) = & \frac{1}{2}(I_1(\mathbf{u}_0)\lambda - 3\lambda - 2\mu + 2\mu\|F(\mathbf{u}_0) \cdot \mathbf{n}\|^2)\mathbf{I} \\ & + (\lambda + \mu)(F(\mathbf{u}_0) \cdot \mathbf{n}) \otimes (F(\mathbf{u}_0) \cdot \mathbf{n}) + \mu\mathbf{B}(\mathbf{u}_0). \end{aligned} \quad (25) $$

Considering pure dilatation deformation, i.e., $\mathbf{u}_0 = \xi(d_0 - 1)$ and $F(\mathbf{u}_0) = d_0\mathbf{I}$ with $d_0 \ge 1$ we have,

$$ Q^{\text{SVK}}(\xi(d_0 - 1), \mathbf{n}) = \frac{(d_0^2 - 1)}{2}(3\lambda + 2\mu)\mathbf{I} + d_0^2[(\lambda + \mu)\mathbf{n} \otimes \mathbf{n} + \mu\mathbf{I}] $$

and

$$ Q^{\text{CNH}}(\xi(d_0 - 1), \mathbf{n}) = (\lambda d_0^4 + \mu d_0^{-2})\mathbf{n} \otimes \mathbf{n} + \mu\mathbf{I}, $$

One can see that the eigenvalues of $Q^{\text{SVK}}(\xi(d_0 - 1), \mathbf{n})$, are increasing functions of $d_0$, meaning that waves tend to propagate faster – at least locally. The CNH constitutive law behaves similarly when $2\lambda \ge \mu$. This behavior contradicts observed experiments [3] where velocities of the shear and pressure waves decrease for some materials – aluminum typically – that undergoes an increase in volume. In this case, proper behavior can be modeled using high-order elasticity models such as the Murnaghan’s, despite some expected and observed stability issues (see the numerical experiments after). In an hypothetical case of pure stretching in one direction, we have locally

$$ F(\mathbf{u}_0) = \begin{pmatrix} d_0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, $$

that is to say $\mathbf{u}_0 = e_1(e_1 \cdot \xi)(d_0 - 1)$, with $d_0 \ge 1$, we have

$$ Q^{\text{SVK}}(\mathbf{u}_0, e_1) = \left( \frac{d_0^2 - 1}{2}(\lambda + 2\mu) + \mu \right) \mathbf{I} + (d_0^2(\lambda + 2\mu) - \mu)e_1 \otimes e_1, $$

and

$$ Q^{\text{CNH}}(\mathbf{u}_0, e_1) = (\lambda + d_0^{-2}\mu)e_1 \otimes e_1 + \mu\mathbf{I}. $$

In this particular case, the behavior of SVK and CNH constitutive laws differs. Indeed one can observe that for increasing values of $d_0$ the eigenvalues of $Q^{\text{SVK}}$ are all increasing while the eigenvalues of $Q^{\text{CNH}}$ associated to the eigenvector $e_1$ is given by $\lambda + (d_0^{-2} + 1)\mu$, hence it decreases with $d_0$. This behavior is observed in our numerical results.
In the acoustoelasticity literature, a more general formulation comprising anisotropic behavior is usually presented. The expansion of the potential up to third-order terms are written in the indicial notation as

$$ \mathcal{W}(e) \approx C_{ijkl}e_{ij}e_{kl} + C_{ijklmn}e_{ij}e_{kl}e_{mn}, \qquad (26) $$

which is a generalization of the expansion (17). The fourth-order tensor $C_{ijkl}$ is the traditional elasticity tensor, with components related to the second order elastic constants from Hooke's law, and the sixth-order tensor $C_{ijklmn}$ has its components related to the third-order elastic constants. One may note that this description of the constitutive behavior complexifies significantly the study of poly-convexity.

Obtaining material parameters for suitable constitutive models is an essential and difficult task. Many works concern the experimental assessment of these parameters. For example, the seminal work presented in [3] gives the TOECs for polystyrene, iron, and Pyrex glass. In [4, 5, 6, 8, 7, 10, 26, 36, 41, 44, 45], the TOECs are obtained for aluminum and other metals. Transversely isotropic materials are specifically treated in [18, 46, 49, 53]. Most of the presented methods rely on pitch-catch time of flight measurements, which is potentially subject to poor robustness [36]. Alternatively, more sophisticated methods based on laser vibrometer aided measurements and space-time Fourier Transform can be applied to enhance the quality of the parameteres reconstruction [59]. Note that, in practice, it appears that experimental set up uncertainties can have a large influence on the estimation of the parameters, as studied in [36]. Additionally, fabrication processes may significantly modify the acoustoelastic response of the material [10, 16].

## 2.3 Dedicated time & space numerical schemes

We now move on to proposing a discretization strategy that addresses the problem of solving the one-way coupled problems “i) ⇒ ii)” defined in (13) and (14). Each problem, namely the nonlinear static problem and the wave propagation problem, presents specific numerical challenges – that have been tackled independently in the literature – and we propose here to aggregate the suitable numerical methods to derive a complete and efficient fully discrete algorithm.

### Solving the quasi-static problem using 3D shell elements.

The nonlinearities associated to problem (13) are handled using Newton-Raphson iterations: for a given $u_0^k$, find $u_0^{k+1} \in W(\hat{\Omega})$ such that

$$ \langle D_u A_b(u_0^k)(u_0^{k+1} - u_0^k), w \rangle = \langle A_b(u_0^k), w \rangle, \quad \forall w \in V(\hat{\Omega}), $$

until convergence. In some cases, typically when the maximal amplitude of the surface traction is large, this iterative process may not converge. To increase robustness in these cases, we apply a continuation algorithm – i.e. successive reductions of the loading amplitude, see [47] and references therein for more details. Writing the increment $u = u_0^{k+1} - u_0^k$, we need to solve at each iteration a linearized elasticity problem in the form of

Find $u \in W(\hat{\Omega})$ the solution, for any $w \in V(\hat{\Omega})$ of

$$ a_k(u, w) - b_k(u, w) = \ell_k(w), \qquad (27) $$

where $a_k(\cdot, \cdot) = a(u_0^k; \cdot, \cdot)$ defined in (7), and

$$ b_k(u, w) = \langle D_u B(u_0^k)u, w \rangle, \quad \ell_k(w) = \langle A_b(u_0^k), w \rangle. $$
It is well-known that for configurations where $\hat{\Omega}$ is a thin elongated geometry, standard finite element procedures have very poor performances, or can even fail to provide meaningful results due to locking phenomena [30, 38]. To circumvent this difficulty we use 3D shell elements [38] that enter the framework of general shell elements [14, 20]. These elements possess desirable properties in the scope of our work: (1) they are locking-free; (2) they rely on the full 3D variational formulations of the type of (27), hence, they can be applied to any type of constitutive law; (3) they rely on a quadratic expansion of the displacement unknown, involving nodes on the outer surface. This makes them readily suitable for coupling with other elements or with themselves, i.e. stacking of shell elements. This aspect is particularly useful to represent stratified materials.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Definition of the MITC9 tying points on the reference hexahedron – extracted from [24] – with $\alpha = 1/\sqrt{3}$ and $\beta = \sqrt{3}/\sqrt{5}$. The tying points for the remaining tensorial components are deduced by symmetry.</figcaption></figure>

Let us assume for illustration purposes that the material at hand is homogeneous and that the domain of interest can be expressed as $\hat{\Omega} = \hat{\omega} \times [-\eta; \eta]$, where $\hat{\omega} \subset \mathbb{R}^2$ is the mid-surface, and $\eta > 0$ is the half-thickness of the domain. Since 3D shell elements rely on a quadratic expansion of the displacement unknown, after defining quadrangular quadratic finite elements on the mid-surface, they amount to considering standard $Q^2$ Lagrange elements on the full 3D domain [38]. More precisely, let us define $\mathcal{T}_H$, a tesselation of $\hat{\Omega}$ with mesh step $H$ and composed of hexahedra with one element in the thickness. Denoting by $\hat{K} = [-1; 1]^3$ the reference element, any hexahedron $K \in \mathcal{T}_H$ is the deformation of the reference element through a $C^1$-diffeomorphism denoted by $\mathcal{F}_K$. We consider the following approximation space

$$V_H(\hat{\Omega}) = \{w_H \in C^0(\hat{\Omega}) \mid \forall K \in \mathcal{T}_H, \exists \hat{w} \in Q^2(\hat{K}), w_H|_K = \hat{w} \circ \mathcal{F}_K^{-1}\},$$

and its vectorial counterpart including the relevant boundary conditions

$$\mathcal{V}_H(\hat{\Omega}) = \{w_H \in [V_H(\hat{\Omega})]^3, w_H|_{\Gamma_D} = 0\} \subset \mathcal{V}(\hat{\Omega}).$$

Considering a direct application of an approximation of problem (27) using $\mathcal{V}_H(\hat{\Omega})$ would not alleviate the locking phenomena. To do so, one can resort to the so-called Mixed Interpolation of Tensorial Components (MITC) method [14, 20, 24]. In essence, this approach consists of interpolating the (discrete) Green-Lagrange tensor on so-called
tying points. The location of these points is different from the Lagrange interpolation points defining the discrete displacement unknown and is also different from one component of the interpolated tensor to another. In our work, we use the MITC9 point distribution since we consider quadratic elements – see Figure 2 for the positioning of these points. Let $i, j \in \{1, 2, 3\}$ be indexes of the components of the Green-Lagrange tensor, and $n_{\text{TP}}^{ij}$ the total number of the tying points denoted by

$$ \{\hat{\xi}_{\text{TP},k}^{ij}\}_{k=1}^{n_{\text{TP}}^{ij}} \subset \hat{K}, $$

and defined on the reference hexahedron. The MITC process applied to a component $e_{ij}$ of the Green-Lagrange tensor is denoted by $\mathcal{I}_{\text{TP}}(e_{ij})$ and reads

$$ \forall K \in \mathcal{I}_H, \quad \forall \xi \in K, \quad \mathcal{I}_{\text{TP}}(e_{ij})|_K(\xi) = \sum_{k=1}^{n_{\text{TP}}^{ij}} \hat{\lambda}_{\text{TP},k}(\mathcal{F}_K^{-1}(\xi)) e_{ij}(\mathcal{F}_K(\hat{\xi}_{\text{TP},k}^{ij})), $$

where $\{\hat{\lambda}_{\text{TP},i}\}_{i=1}^{n_{\text{TP}}}$ are the Lagrange shape functions associated with the tying points. In the following, we denote by $\bar{e} = \mathcal{I}_{\text{TP}}(e)$ the MITC process applied to tensors, in the sense that each component is interpolated on the tying points, namely

$$ \forall i, j \in \{1, 2, 3\}, \quad \bar{e}_{ij} = \mathcal{I}_{\text{TP}}(e)_{ij} = \mathcal{I}_{\text{TP}}(e_{ij}). $$

In the same fashion, considering a discrete linearization point $w_H^* \in \mathcal{V}_H(\hat{\Omega})$, we introduce, for any $v_H, w_H \in \mathcal{V}_H(\hat{\Omega})$,

$$ D_u \bar{e}(w_H^*) \cdot w_H = \mathcal{I}_{\text{TP}}(D_u e(w_H^*) \cdot w_H), \quad D_u^2 \bar{e}(v_H, w_H) = \mathcal{I}_{\text{TP}}(D_u^2 e(v_H, w_H)), $$

the interpolation on the tying points of the first and second differential of the Green-Lagrange tensor. With these notations, the discrete problem associated with (27) reads:

$$ \begin{array}{l} \text{Find } \mathbf{u}_H \in \mathcal{V}_H(\hat{\Omega}) \text{ the solution, for any } \mathbf{w}_H \in \mathcal{V}_H(\hat{\Omega}), \text{ of} \\ \qquad \bar{a}_k(\mathbf{u}_H, \mathbf{w}_H) - b_k(\mathbf{u}_H, \mathbf{w}_H) = \ell_k(\mathbf{w}_H), \end{array} $$

with

$$ \bar{a}_k(\mathbf{u}_H, \mathbf{w}_H) = \int_{\hat{\Omega}} D_u \bar{e}(\mathbf{u}_{0,H}^k) \mathbf{u}_H : D_e \mathbf{\Sigma}(\bar{e}(\mathbf{u}_{0,H}^k)) : D_u \bar{e}(\mathbf{u}_{0,H}^k) \mathbf{w}_H \mathrm{d}\hat{\Omega} \\ + \int_{\hat{\Omega}} \mathbf{\Sigma}(\bar{e}(\mathbf{u}_{0,H}^k)) : D_u^2 \bar{e}(\mathbf{u}_H, \mathbf{w}_H) \mathrm{d}\hat{\Omega}. $$

In the previous expression of the discrete problem, $\mathbf{u}_{0,H}^k \in \mathcal{V}_H(\hat{\Omega})$ represents the discrete displacement field obtained from the prior Newton-Raphson iteration, and is also used as a linearization point for the bilinear and linear forms $b_k(\cdot, \cdot)$ and $\ell_k(\cdot)$.

Spectral elements and explicit leap-frog scheme for the wave propagation problem. We now provide key elements for efficient space and time discretization of the wave problem (14) relying on the combination of spectral finite elements, the mass-lumping technique and the explicit leapfrog scheme. Readers may refer to [23, 25, 33] for a more detailed description of these subjects.

To start with, let us denote by $\mathbf{p} = (p_1, p_2, p_3)^\top \in \mathbb{N}^{*3}$ an integer triplet representing orders of polynomial on every direction of the reference hexahedron. We define the space of polynomial with “anisotropic orders” as

$$ Q^p(\hat{K}) = \operatorname{span} \{ \hat{s}_1^i \hat{s}_2^j \hat{s}_3^k, \quad \forall (i, j, k) \in \mathbb{N} : i \le p_1, j \le p_2, k \le p_3 \}, $$
where $\hat{s} = (\hat{s}_1, \hat{s}_2, \hat{s}_3)^\top$ are the local coordinates on the reference element. From this local space, we can build the global conforming approximation space that supports the semi-discretization procedure of the wave problem (14). To do so, we consider $\mathcal{T}_h$ a hexahedral mesh of $\hat{\Omega}$ – different from $\mathcal{T}_H$ – with mesh step $h$ and such that any hexahedron $K \in \mathcal{T}_h$ is the deformation of the reference element through a $\mathcal{C}^1$- diffeomorphism denoted by $\mathcal{G}_K$. The scalar-valued approximation space reads

$$V_h(\hat{\Omega}) = \{w_h \in \mathcal{C}^0(\bar{\Omega}) \mid \forall K \in \mathcal{T}_h, \exists \hat{w} \in \mathcal{Q}^p(\bar{K}), w_h|_K = \hat{w} \circ \mathcal{G}_K^{-1}\},$$

and, as previously, its vectorial counterpart including relevant boundary conditions is

$$\mathcal{V}_h(\hat{\Omega}) = \{w_h \in [V_h(\hat{\Omega})]^3, w_h|_{\Gamma_D} = 0\} \subset \mathcal{V}(\hat{\Omega}).$$

The local Lagrange (polynomial) basis functions are defined on the reference element from Gauss-Lobatto points. This is a key point in order to avoid Runge's phenomenon, and to render spectral-like convergence as the order of the local polynomials increases. In practice, these points are set on the reference edge $\hat{e} = [-1; 1]$ for any order of approximation (see e.g. [33] for an explicit definition and related examples of these points), and the location of points are deduced in the reference cube by a simple extrusion process.

In addition to this specific choice of approximation space, we use a quadrature formula to numerically compute the integrals involved in the bilinear forms $(\cdot, \cdot)_H$ and $a(\mathbf{u}_0; \cdot, \cdot)$ appearing in (14). The formula is defined on the reference hexahedron. We denote the $\hat{n}_q$ associated weights and points respectively by

$$\{\hat{\omega}_i\}_{i=1}^{\hat{n}_q}, \{\hat{s}_i\}_{i=1}^{\hat{n}_q}.$$

For any function $f$ at least continuous over each element in $\mathcal{T}_h$, the numerical integration process reads

$$\int_{\Omega} f \, d\Omega \approx \sum_{K \in \mathcal{T}_h} \sum_{k=1}^{\hat{n}_q} \hat{\omega}_k f(\mathcal{G}_K(\hat{s}_k)) |\det \nabla \hat{s} \mathcal{G}_K(\hat{s}_k)|.$$

The quadrature points are identical to the Gauss-Lobatto points, and the quadrature weights are the strictly positive Gauss-Lobatto weights. Matching the quadrature points with the points associated with the local Lagrange basis functions is referred to as “mass lumping” [23, 25, 33]. It leads to a diagonal approximated mass matrix. It has been the subject of particular interest in the context of finite element methods applied to time domain wave propagation problems since it enables fully explicit time discretization schemes. To conserve stability and accuracy, specific conditions must be enforced on the quadrature points and weights. For quadrangles and hexahedra, these conditions are met at any order of approximations by the Gauss-Lobatto rules [33, 35]. Following this numerical integration procedure, we denote by

$$(\cdot, \cdot)_H \approx (\cdot, \cdot)_h, \quad a(\mathbf{u}_0; \cdot, \cdot) \approx a_h(\mathbf{u}_0; \cdot, \cdot), \qquad (28)$$

the approximated bilinear forms so that the semi-discrete formulation of the wave problem (14) reads

For any time $t \in [0; T]$, find $\tilde{u}_h(t) \in \mathcal{V}_h(\hat{\Omega})$ the solution, for any $w_h \in \mathcal{V}_h(\hat{\Omega})$, of

$$\frac{d^2}{dt^2}(\tilde{u}_h, w_h)_h + a_h(\mathbf{u}_{0,H}; \tilde{u}_h, w_h) = \langle L(t, \mathbf{u}_{0,H}), w_h \rangle. \qquad (29)$$

with given initial conditions $\tilde{u}_h|_{t=0}$ and $\partial_t \tilde{u}_h|_{t=0}$.
Note that in (29) $u_{0,H}$ is the discrete solution – after convergence of the Newton-Raphson iterations – of problem (13). Moving on to the time discretization of the wave propagation problem, we define $\Delta t > 0$ as the time step. We denote by $t^n = n\Delta t$, and we assume that there exists $N \in \mathbb{N}^*$ such that $T = N\Delta t$. To obtain the fully discrete scheme associated with (29), we apply an explicit second-order time scheme centered at $t^n$, a.k.a. the leapfrog scheme [25, 28]. Namely, denoting by $\tilde{u}_h^n$ the approximation of $\tilde{u}_h(t^n)$ solution of (29), we define the following fully discrete scheme:

$$ \text{For any } n \in \llbracket 1; N \rrbracket, \text{ find } \tilde{u}_h^{n+1} \in \mathcal{V}_h(\hat{\Omega}) \text{ the solution, for any } w_h \in \mathcal{V}_h(\hat{\Omega}), \text{ of} \\ \left( \frac{\tilde{u}_h^{n+1} - 2\tilde{u}_h^n + \tilde{u}_h^{n-1}}{\Delta t^2}, w_h \right)_h + a_h(\mathbf{u}_{0,H}; \tilde{u}_h^n, w_h) = \langle L(t^n, \mathbf{u}_{0,H}), w_h \rangle. \quad (30) $$

with given initial steps $\tilde{u}_h^0$ and $\tilde{u}_h^1$.

From this expression of the fully discrete problem, it is possible to recover a discrete energy conservation property, which is the exact translation of (16). Denoting by $\|\cdot\|_h$ the norm associated with the discrete inner-product $(\cdot, \cdot)_h$, we define the kinematic and potential discrete energies associated with the unique solution of (30) as

$$ \mathcal{H}_h^{n+1/2} = \frac{1}{2} \left\| \frac{\tilde{u}_h^{n+1} - \tilde{u}_h^n}{\Delta t} \right\|_h^2, \quad \mathcal{P}_h^{n+1/2} = \frac{1}{2} a_h(\mathbf{u}_{0,H}; \tilde{u}_h^{n+1}, \tilde{u}_h^n), $$

and the total energy as $\mathcal{E}_h^{n+1/2} = \mathcal{H}_h^{n+1/2} + \mathcal{P}_h^{n+1/2}$. Assuming a null source term in (30), namely $L(t^n, \mathbf{u}_{0,H}) = 0$, then, taking $w_h = \frac{1}{2\Delta t}(\tilde{u}_h^{n+1} - \tilde{u}_h^{n-1})$ as test function and using the symmetry of the stiffness bilinear form, we have the following energy conservation

$$ \frac{\mathcal{E}_h^{n+1/2} - \mathcal{E}_h^{n-1/2}}{\Delta t} = 0. \quad (31) $$

This conservation property is critical since it leads to the conditions – on the stiffness bilinear form and on the time step – for the fully discrete (30) to be stable. Indeed, note that the total energy can be re-written in the following form

$$ \mathcal{E}_h^{n+1/2} = \frac{1}{2} \left\| \frac{\tilde{u}_h^{n+1} - \tilde{u}_h^n}{\Delta t} \right\|_h^2 - \frac{\Delta t^2}{8} a_h \left( \mathbf{u}_{0,H}; \frac{\tilde{u}_h^{n+1} - \tilde{u}_h^n}{\Delta t}, \frac{\tilde{u}_h^{n+1} - \tilde{u}_h^n}{\Delta t} \right) \\ + a_h \left( \mathbf{u}_{0,H}; \frac{\tilde{u}_h^{n+1} + \tilde{u}_h^n}{2}, \frac{\tilde{u}_h^{n+1} + \tilde{u}_h^n}{2} \right). $$

Hence, upon the assumption that

* The stiffness bilinear form in (30) satisfies the coercivity property (15),

* The time step satisfies

$$ \Delta t \le 2 \left( \sup_{w_h \in \mathcal{V}_h(\hat{\Omega})} \frac{a_h(\mathbf{u}_{0,H}; w_h, w_h)}{\|w_h\|_h^2} \right)^{-\frac{1}{2}}, \quad (32) $$

we can ensure that the total energy is a positive functional and that it actually corresponds to an equivalent norm on the discrete solution. Following arguments provided e.g. in [28, 43, 63] one can obtain the following stability result: there exists a constant $C > 0$ that depends on the initial steps $\tilde{u}_h^0$ and $\tilde{u}_h^1$, the source term $L(t^n; \mathbf{u}_{0,H})$ and the maximal time $T$, such that the solution of (30) satisfies

$$ \|\tilde{u}_h^{n+1}\|_h \le C(T, L, \tilde{u}_h^0, \tilde{u}_h^1), \quad \forall n \in \llbracket 1; N \rrbracket. $$
Remarks on practical considerations.

Prior to giving substantiated numerical results proving the soundness of the presented modeling strategy, let us gather some important remarks.

First, in the definition of the discrete space $V_h(\hat{\Omega})$ we consider anisotropic orders of approximation. This aspect appears to be particularly convenient for thin and potentially stratified materials. Indeed, it allows for high-order elements in the mid-surface of the domain – with low numerical dispersion or numerical anisotropy [25, 31, 34] – while maintaining a reasonable number of Degrees of Freedom (DoF) by using lower order in the thickness of the material.

Second, the use of approximated numerical integration as denoted in (28) implies in particular a lumped, i.e. diagonal, mass matrix. This is an important aspect in the context of high-frequency wave propagation since the leapfrog scheme (30) thus becomes fully explicit. At each time step, only a diagonal matrix needs to be inverted. Additionally, using numerical integration with points matching the ones defining the local Lagrange polynomial leads to very sparse discrete gradient computations. Hence, the approximation $a(\mathbf{u}_0; \cdot, \cdot) \approx a_h(\mathbf{u}_0; \cdot, \cdot)$ also implies significant performance gains – see e.g. [58] for more details on these practical aspects.

Moreover, in practice, the mesh step for the locking-free shell elements depends mostly on the geometrical characteristics of the domain and the material variations – typically one element per material layer in the case of stratified materials. For the wave problem, to these constraints is added the constraint of having enough discretization points per wavelength. With these qualitative meshing directives in mind, in the context of high-frequency wave propagation, one can expect that the $\mathcal{T}_h$ is significantly more refined than $\mathcal{T}_H$, i.e. $h \ll H$.

In light of the previous remark, in order to solve the fully discrete one-way coupled problems, one needs to resort to a mean of interpolating $\mathbf{u}_{0,H} \in V_H(\hat{\Omega})$ onto $V_h(\hat{\Omega})$. Let us denote by $\mathcal{I}_h(\cdot)$ the nodal interpolator from $V_H(\hat{\Omega})$ to $V_h(\hat{\Omega})$, without entering into the details of the exact expression of this interpolator, let us simply remark that the stiffness bilinear form actually reads

$$a_h(\mathbf{u}_{0,H}; \cdot, \cdot) = a_h(\mathcal{I}_h(\mathbf{u}_{0,H}); \cdot, \cdot).$$

Finally, the condition (32) on the time step is often referred to as the Courant – Friedrichs – Lewy (CFL) condition of the fully discrete scheme. Note that in the particular cases of interest in our work, this CFL condition naturally depends on the linearization point, i.e. the displacement field $\mathbf{u}_{0,H}$. Furthermore, in terms of stability, in addition to this CFL condition, the coercivity assumption is capital – a condition that is also required for the well-posedness of the continuous problem.

## 3 Numerical results and experimental validation

In this section, we apply our model and numerical methods for different configurations. First, we assess the stability issue caused by a potential failure to ensure coercivity. Second, we apply our methods to realistic cases in the SHM context. At last, we present how we validate our model with experimental data available in the literature. The computations were done using a laptop workstation equipped with a Intel i9-9880H CPU and 32Gb of RAM. The wave propagation simulations are performed within the SFEM kernel of the CIVA platform<fnref n="1" /> developed at CEA.

<footnote n="1">¹https://www.extende.com</footnote>
<table><caption>**Table 1.** Material parameters for the metals used in this work.</caption><thead><tr><th>Material</th><th>Density(g/cm³)</th><th>λ(GPa)</th><th>μ(GPa)</th><th>l(GPa)</th><th>m(GPa)</th><th>n(GPa)</th></tr></thead><tbody><tr><td>Aluminum</td><td>2.704</td><td>54.308</td><td>27.174</td><td>-181.0</td><td>-289.0</td><td>-336.0</td></tr><tr><td>Steel</td><td>7.850</td><td>111</td><td>82.1</td><td>-459.0</td><td>-461.0</td><td>-358.0</td></tr></tbody></table>

## 3.1 Illustration of potential stability issues of the fully discrete scheme

As discussed previously, the tangent stiffness operator in (14) must satisfy the coercivity assumption (15), otherwise the well-posedness of the problem cannot be guaranteed. This assumption cannot be satisfied for any deformation or constitutive law. In this section, we aim to illustrate this problem in a test model. At a discrete level, coercivity is equivalent to positiveness of the operator, which can be assessed by computing its eigenvalues. In the zero-deformation case, the presented model equates to linear elasticity (Hooke's law), and $a(0; \cdot, \cdot)$ is known to be coercive [17]. Here, we assess stability for deformed cases by analyzing the evolution of the system energy and the eigenvalues of the local stiffness matrices.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Configuration and deformations for illustrating the potential stability issues.</figcaption></figure>

If instability would appear by the deformation, it could be interpreted as a shift of tangent stiffness operator eigenvalues towards negative values. In this case, the negative eigenvalues that are related to the instability of the scheme are expected to be among those of the smallest amplitude, i.e. low spatial frequency. Computing the lowest eigenvalues for the stiffness operator may have prohibitive computational cost due to the size of the problems at hand. To circumvent this computational difficulty, we analyze instead the eigenvalues of the local stiffness matrix related to each finite element. Although negative eigenvalues for the local stiffness matrices do not necessarily mean negative eigenvalues for the global matrix, it does give an insight into the problematic regions and related severity. In addition to the study of the lowest eigenvalue, we analyze the evolution in time of the total system's kinetic and potential energy for different initial conditions, showing the nature of such stability issues.

The model configuration used is an aluminum cuboid of 25mm × 25mm × 100mm depicted in Figure 3a. Dirichlet boundary condition is considered at $x=0$. The material parameters are present in Table 1. We analyze different constitutive laws and types of mechanical solicitations. The analysis is done using the aforementioned SVK, CNH, and Murnaghan's constitutive laws. For each case, three different types of mechanical solicitations as body forces will be considered:

* **Stretching**: 300N/mm³ (Y direction),

* **Bending**: 7.5N/mm³ (-Z direction),
• **Compreming:** 135N/mm³ (-Y direction).

511

512

513

514

515

The associated deformations are illustrated in Figure 3b for the CNH case. The SVK law has slightly different deformation with the same nature and amplitudes. Computing the quasi-static deformation using Murnaghan's constitutive law did not converge due to its particularly ill-posed tangent operator.

516

517

518

519

520

521

522

523

524

After computing the deformations caused by these forces for each constitutive law, the displacement field is given as input to our wave propagation solver. The same constitutive law is used for the quasi-static problem and the wave propagation problem, except in Murnaghan case. For the Murnaghan case, we compute the deformation with CNH instead. The deformation field is interpolated from the quasi-static problem mesh to the wave propagation problem mesh. After discretizing the functional space using spectral finite elements, we assemble the local tangent stiffness matrices and compute their eigenvalues. The minimum eigenvalues, for each element, are shown in Figure 4 for different types of deformations and constitutive laws.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Minimum (local) tangent stiffness matrix eigenvalues for different types of deformation and constitutive laws.</figcaption></figure>

525

526

527

528

529

530

531

The CNH constitutive law presents less severity, i.e. lowest amplitude negative eigenvalues. The SVK and CNH cases are more unstable for compressed cases – cases associated with buckling phenomena. The Murnaghan’s law presents the highest negative eigenvalues with critical behavior in the stretched case. The stability issue can be also assessed by observing the evolution of the system’s potential and kinetic energy. We run the wave propagation problem, without source-term, from initial conditions with specific spatial frequency signature, namely for $n = 1, 2, 3$

$$ \tilde{u}_h^0 = - \left[ \sin \left( \frac{25-x}{25} 2\pi n \right) \sin \left( \frac{25-y}{25} 2\pi n \right) \sin \left( \frac{25-z}{25} 2\pi n \right) \right] (e_x + e_y + e_z), \\ \tilde{u}_h^1 = -\tilde{u}_h^0. $$

532

533

534

It allows us to solicit the low-frequency modes and assess the instability, relating the frequency content to the issue. Running simulations with these initial conditions, we monitor the total, kinetic and potential energies as shown in Figure 5.
<figure><img src="image_8.png" /><figcaption>(a) Stable behavior for a stable case (Stretched; CNH).</figcaption></figure>

<caption>(b) Unstable case (Bended; CNH) for different initial conditions.</caption>

<caption>**Figure 5.** Evolution of total energies (a.u.) for the wave propagation problem.</caption>

As the total energy is conserved, if both kinetic and potential energies are positive, the system's energy is bounded. However, if any of the energies become negative, the solution may diverge. For the stable cases, the energy transfer between kinetic and potential is kept throughout the simulation. For the unstable cases, as soon as the potential energy presents negative values, the system starts to diverge. The excitation of lower frequencies presented instability in fewer time steps as they excite lower frequency modes, i.e. modes related to lower eigenvalues. For the presented forces, the CNH case is unstable only under compression, SVK is unstable for bending and compression and Murnaghan's is unstable for all presented deformations. In general, the stability issues are not only dependent of buckling phenomena but also on the choice of constitutive law.

## 3.2 Illustrations related to Structural Health Monitoring

We now illustrate the generic aspect of our work by applying it to different cases. The first case is a steel pipe subjected to the 4-point bending test, inspired by an experiment done in [54]. We simulate ultrasonic propagation with different deformation intensities and compare the acoustoelastic effects in the signal with ones caused by a crack through the thickness of the pipe. In the second case, we model ultrasonic wave propagation in a 16-layer stratified composite plate under torsion forces. For such modeling, we use the transversely isotropic law (22) for each of the 16 layers associated with a principal fiber direction. These examples illustrate the application of our method to realistic and complex loading configurations, and the simulation results could be used to assess the efficiency of a SHM system.

### Illustration on a pipe subject to a 4-point bending test.

We present here a case of wave propagation in a steel pipe subjected to a 4-point bending test. This configuration is usually used to perform fatigue tests in welded pipes. The experimental setup is inspired by the experiment done in [54]. We model and simulate the wave propagation in the pipe subjected to such bending forces with different intensity levels. We use Murnaghan's constitutive law with material parameters for steel (Table 1). For these levels of deformation and material parameters the wave propagation problem does not present stability issues. The surface force is
applied at two regions close to the pipe extremities. The configuration is depicted in Figure 6, with a 2.94m long, 8mm thick pipe and an inner radius of 102.65mm. The resulting deformation for a total force of 220kN is shown in Figure 7.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Numerical setup and mesh for the quasi-static problem.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Deformation of the pipe for a 220kN loading. (Deformation scaled 20x for visualization).</figcaption></figure>

Given the different deformations as input, we interpolate them in the wave propagation problem mesh and compute the wavefield generated by a unique source term. The wave problem mesh, source location and measurement point are depicted in Figure 8, as well as the position of a potential crack. The ultrasonic waves are excited close to one edge of the pipe and measured at the other edge. For a case with no cracks, we simulate the wavefield and extract the signal from the modeled transducer at the measurement point and plot in Figure 9 for different levels of mechanical loadings. Snapshots of the wavefield are shown in Figure 10.

In Figure 9, we can see that the effect of bending forces results in a change in time of flight, decreasing the wave velocity, and also a change in amplitude of the signals with higher forces.

To evaluate how the effect of mechanical loading on the signal compares with the one caused by a crack, we set up and simulate three numerical configurations:

* S: non deformed pipe (considered as baseline; no defects)

Simulation is done considering no deformation with a pristine pipe,

* $S_p$: deformed pipe (without defect)

Simulation is done considering the pipe under bending forces of 220kN,

* $S_d$: pipe with defect (without deformation)

Simulation is done considering the undeformed pipe with a crack of 19mm wide through half the thickness.

From each case, we extract the signal of a modeled transducer at the measurement point. The defect through half the thickness was modeled by disconnecting elements. By
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Transducer, potential crack and measurement point positions for the wave propagation problem. On the right-hand side, the mesh for the wave propagation problem characterized by the Gauss-Lobatto points.</figcaption></figure>

<table><caption>**Table 2.** Material parameters for CFRP. [29]</caption><thead><tr><th>Material</th><th>Density(g/cm<sup>3</sup>)</th><th>C<sub>11</sub></th><th>C<sub>22</sub></th><th>C<sub>12</sub></th><th>C<sub>23</sub></th><th>C<sub>66</sub></th></tr></thead><tbody><tr><td>CFRP</td><td>1.6</td><td>155.17</td><td>13.22</td><td>5.34</td><td>6.98</td><td>5.77</td></tr></tbody></table>

subtracting the signals from the baseline, we can isolate the individual contribution to the signal of the mechanical loading and the crack. The baseline signal $S$ is subtracted from $S_p$ and $S_d$ and the results are plotted in Figure 11.

Regarding performances, the quasi-static problem is 148k DoFs and takes 3 minutes to run 5 Newton-Raphson iterations while using 2.2Gb of RAM. The wave propagation problem requires around 3.86M mesh points to ensure at least two elements of fourth-order per wavelength, resulting in 11.58M DoFs. The timestep is $\Delta t = 0.11\mu s$ to comply with the CFL condition (32). Using an efficient implementation of the spectral finite elements, the numerical solver uses 2.3Gb of RAM, 73 seconds for initialization (starting the model, initial outputs and interpolating the structure displacement field) and 40 minutes to compute approximately 4336 timesteps (including wavefield outputs).

The effect on the signal caused by the deformation alone is comparable to the effect caused by a crack for such configuration and type of measurement. Thus, not considering these effects when analysing wave propagation data for this level of mechanical loading results in less reliable detection and evaluation of such type of damages.

### Illustration on a stratified anisotropic plate.

Carbon fiber-reinforced polymers are largely employed in aerospace industries, where ensuring safety is essential. They are made of a polymer matrix reinforced with stacked layers of fibers. Usually, each fiber layer has a principal direction, presenting a transversely isotropic response to mechanical excitations. Here, we model the wave propagation in a 16-layer Carbon Fiber Reinforced Polymer (CFRP) plate with dimensions of 900 × 300 × 2 (16 × 0.125) mm³ under torsion as illustrated in Figure 6.

Each layer is modeled using the hyperelastic potential (22) where the direction related to the invariants $I_4$ and $I_5$ were set accordingly to the fiber directions in Figure 12. The material parameters are shown in Table 2. The deformation was caused by setting the Dirichlet boundary condition at $x = 0$ and applying a body force with the following form

$$f = \frac{y - 150}{300} 10^{-3} \text{N/mm}^3, \quad \forall y \in [0, 300].$$
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Transducer signal at the measurement point for different bending forces.</figcaption></figure>

The resulting deformation by solving the quasi-static problem is shown in *Figure 13*.
Snapshots of the wavefield computed solving the wave propagation problem can be
seen in *Figure 14*.

The quasi-static problem has 270k DoFs and takes 50 minutes to run 10 Newton-Raphson iterations while using 12Gb of RAM. The wave propagation problem requires around 8.3M mesh points, having at least two elements of fourth-order per wavelength, resulting in 25M DoFs. Satisfying the CFL condition (32), the timestep is $\Delta t = 0.0159\mu s$. The wave propagation solver requires 4Gb of RAM, 9 minutes for initialization (starting the model, initial outputs and interpolating the structure displacement field) and 2 hours and 54 minutes to compute approximately 18867 timesteps (including wavefield outputs). This computation is particularly expensive as we need to satisfy the CFL condition for the mesh step related to the thickness discretization (16 elements for 2mm). This problem could be mitigated by using implicit-explicit schemes as in [64].

This illustration shows the capability of our simulation tools to model thin structures
with the possibility of choosing a hyperelastic law that best suits the modeling objectives.
As discussed in the paragraph about acoustoelasticity, the acoustoelastic effect is
controlled by the TOECs. Here, we choose to use a hyperelastic law that matches only
second-order elastic constants although the implementation of laws that include the
TOECs can also be considered.

**3.3 Experimental validation on an isotropic plate**

Finally, in this section, we model and simulate an experimental setup to validate our
model in the isotropic case. The experiment is presented in [41], where the authors
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Snapshots of the wavefield computed in the steel pipe subjected to the 4-point bending test at t = 70, 140, 230, 300, 390 and 480µs (left to right, top to bottom).</figcaption></figure>

<table><caption>**Table 3.** Material parameters for aluminum, fitted to experimental data.</caption><thead><tr><th>Material</th><th>Density (g/cm³)</th><th>λ (GPa)</th><th>μ</th><th>l</th><th>m</th><th>n</th></tr></thead><tbody><tr><td>Aluminum (fitted)</td><td>2.704</td><td>54.308</td><td>27.174</td><td>-181.0</td><td>-289.0</td><td>-336.0</td></tr></tbody></table>

apply axial forces to an aluminum plate and measure the change in wave speed caused by the load-induced deformation, for different directions and guided wave modes. Using Murnaghan's constitutive law, the authors propose an eigenproblem to obtain guided wave modes and their wave speeds in biaxially stressed isotropic materials. By applying canonical deformations and measuring the change in wave speed for different angles, the authors calibrated the TOECs of Murnaghan's constitutive law using experimental data. The obtained parameters are presented in Table 3. Here, we aim at modeling and simulating this experiment using the hyperelastic potential (21) with these parameters. We compare the simulated results with their experimental data to validate our model and numerical tools. The experiment is depicted in Figure 15.

In the experimental setup, eleven levels of axial loading are applied from 0MPa to 57.5MPa and the ultrasonic wave speeds are measured for 9 different angles of propagation with transducers positioned at 109mm from the center in the undeformed configuration [37]. For the validation procedures, we use their measured changes in wave speed for the S0 mode excited with a 250kHz 5-cycle Hanning windowed cosine. The wave speed is obtained by retrieving the ToF between transducer pairs for the S0 mode at a given load. To avoid effects of dispersion on the ToF measurement, we identify the 4th zero-crossing of the signal as the least affected by dispersion. As our solution is computed in the reference configuration, there is no need for correcting the distance change between transducers caused by the deformation. The authors in [41] fitted their parameters while neglecting this distance difference when computing the wave speed using the ToF. The computed displacement between the axial pair (#1,#6)
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Transducer signal at measure point subtracted from the baseline for the case with defect-only and load-only.</figcaption></figure>

has a maximum of 0.173mm at 57.5MPa. It represents a change in the total pair distance of around 0.08%, or, in terms of wave speed, around 4.2m/s when considering the mean wave speed of 5300m/s.

We simulate and obtain, for the different load levels, the changes in time of flight (ToF) for the angles 90°, 45°, 0°, -22.5°, -67.5°, corresponding to the transducer pairs (#1,#6), (#5,#10), (#4,#9), (#2,#7) and (#3,#8). The Y axis is considered to be at an angle of 90°. The changes in wave speed in the loading axis with respect to the load intensity are plotted in *Figure 16a*. Here, simulated results with different constitutive laws (SVK, CNH and Murnaghan's) are compared to experimental data and the analytical reference solution proposed in [41]. The behavior of SVK and CNH constitutive laws reflect the analysis done in the paragraph about acoustoelastic effects. With the SVK law, increasing axial load leads to increasing incremental wave speeds, what is opposite to the observed effect. Although the CNH models the right tendency of decreasing the wave speed, it underestimate its effects. Making use of the TOECs, the Murnaghan's law is able to properly model the experimental data. For the maximum load of 57.5MPa, the changes in wave speed with respect to the angle of propagation are plotted in *Figure 16b*. Here, our method implemented with the Murnaghan's law is also compared with the experimental data and the analytical proposition when obtaining the angular-dependency of the changes in wave speed. With our simulation, we retrieve the expected bell-shaped curve, effect of the load-induced anisotropy also obtained from the ToF measurements from experimental and simulated data.

Summarizing our validation achievement, the experiment is modeled and simulated starting from the mechanical loading that results in the deformation of the aluminum plate. This deformation is used as input to our linearized wave propagation model, and, using the proposed constitutive law and material parameters we are able to extract the ultrasonic signals of the same nature as in the experiment. Our simulation reliably matches the experimental results while computing the whole wave field. We also implement different constitutive laws, obtaining their expected behavior. This further validates each modeling component and solvers used.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Configuration for the quasi-static problem (upper left). Zoom showing the stratified aspect (upper right). Ply configuration (bottom).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Deformation of the CFRP plate under torsion. (visualization scaled 2x).</figcaption></figure>

# 4 Conclusion & perspectives

In our work, we have proposed a set of numerical tools addressing the problem of ultrasonic guided-wave propagation in materials subject to pre-deformation. Modeling of such configurations includes two different problems: a quasi-static (time-independent) nonlinear problem satisfied by the large displacement field, and a linearized time-domain wave propagation problem. Our approach is generic in the sense that it can represent any type of constitutive law and geometries. For each problem, we have carefully chosen a dedicated numerical method – namely the 3D shell elements for the first problem and the time-domain SFEM for the second one – leading to computationally efficient algorithms. We have presented numerical results on 3D cases related to SHM configurations, thus illustrating the performances of our method for both isotropic and anisotropic materials. We also have compared our results against experimental data extracted from the literature, thus validating the overall modeling and numerical toolchain.

Thanks to the presented approach, we were able to illustrate how the effects of a pre-deformation can alter the ultrasonic signals used to monitor a specimen under inspection. It makes this work in a good position to get closer to real applications in industrial contexts, and to provide meaningful modeling tools for potential industrial partners aiming at designing SHM systems. In order to enhance the attractiveness of these modeling tools, it would be interesting to include a mean of automatically removing the bias induced by the pre-deformation or the EOCs in general. This
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Wavefield computed in a 16 layers CFRP plate under torsion at $t = 45, 90, 180, 240, 270$ and $300\mu s$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Experimental setup for calibrating the TOECs using the acoustoelastic effect.</figcaption></figure>

708 could typically be done by considering an inverse problem where we would aim at 709 reconstructing the pre-deformation sustained by the structure from the guided-wave 710 measurements.

## References

711 [1] M. A. Biot. "The Influence of Initial Stress on Elastic Waves". In: *Journal of Applied Physics* 11.8 (1940), pp. 522-530.

712 [2] F. D. Murnaghan. *Finite Deformation of an Elastic Solid*. John Wiley & Sons, 1951.

713 [3] D. S. Hughes and J. L. Kelly. "Second-Order Elastic Deformation of Solids". In: *Physical Review* 92.5 (1953), pp. 1145-1149.

714 [4] A. Hikata, R. Truell, A. Granato, B. Chick, and K. Lücke. "Sensitivity of Ultrasonic Attenuation and Velocity Changes to Plastic Deformation and Recovery in Aluminum". In: *Journal of Applied Physics* 27.4 (1956), pp. 396-404.

715 [5] R. H. Bergman and R. A. Shahbender. "Effect of Statically Applied Stresses on the Velocity of Propagation of Ultrasonic Waves". In: *Journal of Applied Physics* 29.12 (1958), pp. 1736-1738.

716 [6] R. T. Smith, R. Stern, and R. W. B. Stephens. "Third-Order Elastic Moduli of Polycrystalline Metals from Ultrasonic Velocity Measurements". In: *The Journal of the Acoustical Society of America* 40.5 (1966), pp. 1002-1008.

717 [7] J. R. Asay and A. H. Guenther. "Ultrasonic Studies of 1060 and 6061-T6 Aluminum". In: *Journal of Applied Physics* 38.10 (1967), pp. 4086-4088.

718 [8] D. I. Crecraft. "The measurement of applied and residual stresses in metals using ultrasonic waves". In: *Journal of Sound and Vibration* 5.1 (1967), pp. 173-192.
<figure><img src="image_19.png" /><figcaption>(a) Load-induced wavespeed changes for the S0 guided wave mode at 250kHz.</figcaption></figure>

<caption>(b) Load-induced wavespeed changes for the S0 mode at 250kHz. (Murnaghan's law)</caption>

**Figure 16.** Comparison between simulated data (for different constitutive laws) and experimental data, validating our proposed model.

[9] J. L. Lions. *Optimal control of systems governed by partial differential equations*. Vol. 170. Springer, 1971.

[10] N. N. Hsu. "Acoustical birefringence and the use of ultrasonic waves for experimental stress analysis". In: *Experimental Mechanics* 14.5 (1974), pp. 169–176.

[11] G. Duvaut and J. L. Lions. *Inequalities in mechanics and physics*. Springer-Verlag, 1976.

[12] J. E. Marsden and T. J. R. Hughes. "Topics in the mathematical foundations of elasticity". In: (1978).

[13] R. W. Ogden. "Non-linear elastic deformations". In: *Engineering Analysis* 1.2 (1984), p. 119.

[14] K.-J. Bathe and E. N. Dvorkin. "A formulation of general shell elements—the use of mixed interpolation of tensorial components". In: *International journal for numerical methods in engineering* 22.3 (1986), pp. 697–722.

[15] L. D. Landau, E. M. Lifšic, and L. D. Landau. *Theory of elasticity*. 2. ed. Course of theoretical physics Vol. 7. Oxford: Pergamon Press, 1986.

[16] R. B. Thompson, S. S. Lee, and J. F. Smith. "Angular dependence of ultrasonic wave propagation in a stressed, orthorhombic continuum: Theory and application to the measurement of stress and texture". In: *The Journal of the Acoustical Society of America* 80.3 (1986), pp. 921–931.

[17] P. G. Ciarlet. *Mathematical elasticity. Volume I, Three-dimensional elasticity*. Amsterdam, New York: North-Holland; Sole distributors for the U.S.A. and Canada, Elsevier Science Pub. Co., 1988.

[18] W. H. Prosser and R. E. Green. "Characterization of the Nonlinear Elastic Properties of Graphite/Epoxy Composites Using Ultrasound". In: *Journal of Reinforced Plastics and Composites* 9.2 (1990), pp. 162–173.

[19] R. Dautray and J.-L. Lions. *Evolution problems I, volume 5 of Mathematical analysis and numerical methods for science and technology*. 1992.
[20] M. L. Bucalem and K.-J. Bathe. "Higher-order MITC general shell elements". In: *International Journal for Numerical Methods in Engineering* 36.21 (1993), pp. 3729-3754.

[21] P. Le Tallec. "Numerical methods for nonlinear three-dimensional elasticity". In: *Handbook of Numerical Analysis*. Vol. 3. Elsevier, 1994, pp. 465-622.

[22] J Bonet and A. Burton. "A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations". In: *Computer methods in applied mechanics and engineering* 162.1-4 (1998), pp. 151-164.

[23] D. Komatitsch, J.-P. Vilotte, R. Vai, J. M. Castillo-Covarrubias, and F. J. Sanchez-Sesma. "The spectral element method for elastic wave equations-application to 2-D and 3-D seismic problems". In: *International Journal for numerical methods in engineering* 45.9 (1999), pp. 1139-1164.

[24] K.-J. Bathe, A. Iosilevich, and D. Chapelle. "An evaluation of the MITC shell elements". In: *Computers & Structures* 75.1 (2000), pp. 1-30.

[25] G. C. Cohen. *Higher-order numerical methods for transient wave equations*. 1st edition. Scientific computation. Berlin Heidelberg: Springer-Verlag, 2002.

[26] A. Santos Jr and D. Bray. "Comparison of Acoustoelastic Methods to Evaluate Stresses in Steel Plates and Bars". In: *Journal of Pressure Vessel Technology-transactions of The Asme - J PRESSURE VESSEL TECHNOL* 124 (2002).

[27] R. A. Adams and J. J. Fournier. *Sobolev spaces*. Elsevier, 2003.

[28] P. Joly. "Variational Methods for Time-Dependent Wave Propagation Problems". In: *Topics in Computational Wave Propagation: Direct and Inverse Problems*. Lecture Notes in Computational Science and Engineering. Berlin, Heidelberg: Springer, 2003, pp. 201-264.

[29] S. Lonné. *Modélisation de la propagation ultrasonore dans les matériaux composites obtenus par le procédé de fabrication RTM (Resin Transfer Molding)*. 2003.

[30] K.-J. Bathe. *Finite Element Procedures*. Klaus-Jurgen Bathe, 2006.

[31] J. D. Basabe. "Grid dispersion and stability criteria of some common finite-difference and finite-element methods for acoustic and elastic wave propagation". In: *Seg Technical Program Expanded Abstracts* (2007).

[32] F. Chen and P. D. Wilcox. "The effect of load on guided wave propagation". In: *Ultrasonics* 47.1 (2007), pp. 111-122.

[33] P. Joly. "Numerical Methods for Elastic Wave Propagation". In: *Waves in Nonlinear Pre-Stressed Materials*. CISM Courses and Lectures 495. Vienna: Springer, 2007, pp. 181-281.

[34] G. Seriani and S. P. Oliveira. "Dft modal analysis of spectral element methods for acoustic wave propagation". In: *Journal of Computational Acoustics* 16.04 (2008), pp. 531-561.

[35] M. Duruflé, P. Grob, and P. 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* 25.3 (2009), pp. 526-551.

[36] D. D. Muir, T. E. Michaels, and J. E. Michaels. "One-sided liquid coupled ultrasonic method for recovery of third order elastic constants". In: *AIP Conference Proceedings* 1096.1 (2009), pp. 1275-1282.
805

806

[37] Gandhi, M. *Determination of dispersion curves for acoustoelastic Lamb Wave propagation*. 2010.

807

808

809

[38] D. Chapelle and K.-J. Bathe. *The Finite Element Analysis of Shells - Fundamentals*. 2nd ed. Computational Fluid and Solid Mechanics. Berlin Heidelberg: Springer-Verlag, 2011.

810

811

812

[39] M. Shams, M. Destrade, and R. Ogden. "Initial stresses in elastic solids: Constitutive laws and acoustoelasticity". In: *Wave Motion* 48.7 (2011), pp. 552–567.

813

814

[40] Z. Abiza, M. Destrade, and R. W. Ogden. "Large acoustoelastic effect". In: *Wave Motion* 49.2 (2012), pp. 364–374.

815

816

817

[41] N. Gandhi, J. E. Michaels, and S. J. Lee. "Acoustoelastic Lamb wave propagation in biaxially stressed plates". In: *The Journal of the Acoustical Society of America* 132.3 (2012), pp. 1284–1293.

818

819

[42] P. Loveday, C. Long, and P. Wilcox. "Semi-analytical finite element analysis of the influence of axial loads on elastic waveguides". In: 2012.

820

821

822

[43] J. Chabassier and S. Imperiale. "Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string". In: *Wave Motion* 50.3 (2013), pp. 456–480.

823

824

825

[44] F. Shi, J. E. Michaels, and S. J. Lee. "In situ estimation of applied biaxial loads with Lamb waves". In: *The Journal of the Acoustical Society of America* 133.2 (2013), pp. 677–687.

826

827

828

[45] F. Song, G. L. Huang, and G. K. Hu. "Coupled piezo-elastodynamic modeling of guided wave excitation and propagation in plates with applied prestresses". In: *Journal of Intelligent Material Systems and Structures* 24.5 (2013), pp. 598–611.

829

830

831

[46] A. Bulletti and L. Capineri. "Interdigital Piezopolymer Transducers for Time of Flight Measurements with Ultrasonic Lamb Waves on Carbon-Epoxy Composites under Pure Bending Stress". In: *Journal of Sensors* 2015 (2015), pp. 1–11.

832

833

834

[47] S Léger, J Deteix, and A Fortin. "A Moore-Penrose continuation method based on a Schur complement approach for nonlinear finite element bifurcation problems". In: *Computers & Structures* 152 (2015), pp. 173–184.

835

836

837

[48] A. Pau and F. Lanza di Scalea. "Nonlinear guided wave propagation in prestressed plates". In: *The Journal of the Acoustical Society of America* 137.3 (2015), pp. 1529–1540.

838

839

840

841

[49] A. C. Kubrusly, N. Pérez, T. F. de Oliveira, J. C. Adamowski, A. M. B. Braga, and J. P. von der Weid. "Mechanical Strain Sensing by Broadband Time Reversal in Plates". In: *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control* 63.5 (2016), pp. 746–756.

842

843

[50] M. Mitra and S Gopalakrishnan. "Guided wave based structural health monitoring: A review". In: *Smart Materials and Structures* 25.5 (2016), p. 053001.

844

845

846

[51] N. Pei and L. J. Bond. "Higher order acoustoelastic Lamb wave propagation in stressed plates". In: *The Journal of the Acoustical Society of America* 140.5 (2016), pp. 3834–3843.

847

848

849

[52] B. Dubuc, A. Ebrahimkhanlou, and S. Salamone. "Effect of pressurization on helical guided wave energy velocity in fluid-filled pipes". In: *Ultrasonics* 75 (2017), pp. 145–154.
850

851

852

[53] B. Dubuc, A. Ebrahimkhanlou, and S. Salamone. "The effect of applied stress on the phase and group velocity of guided waves in anisotropic plates". In: *The Journal of the Acoustical Society of America* 142.6 (2017), pp. 3553-3563.

853

854

855

[54] K. Tschöke, B. Weihnacht, E. Schulze, T. Gaul, L. Schubert, and R. Neubeck. "Determination of Defect Sizes with the help of Structural-Health-Monitoring Methods based on Guided Waves". In: *2017-12*. 2017.

856

857

858

859

[55] B. Dubuc, A. Ebrahimkhanlou, and S. Salamone. "Computation of propagating and non-propagating guided modes in nonuniformly stressed plates using spectral methods". In: *The Journal of the Acoustical Society of America* 143.6 (2018), pp. 3220-3230.

860

861

862

[56] M. Mohabuth, A. Kotousov, C.-T. Ng, and L. R. F. Rose. "Implication of changing loading conditions on structural health monitoring utilising guided waves". In: *Smart Materials and Structures* 27.2 (2018), p. 025003.

863

864

865

[57] K. Peddeti and S. Santhanam. "Dispersion curves for Lamb wave propagation in prestressed plates using a semi-analytical finite element analysis". In: *The Journal of the Acoustical Society of America* 143.2 (2018), pp. 829-840.

866

867

868

869

870

[58] Alexandre Imperiale and E. Demaldent. "A macro-element strategy based upon spectral finite elements and mortar elements for transient wave propagation modeling. Application to ultrasonic testing of laminate composite materials". In: *International Journal for Numerical Methods in Engineering* 119.10 (2019), pp. 964-990.

871

872

873

[59] M. Ponschab, D. A. Kiefer, and S. J. Rupitsch. "Towards an Inverse Characterization of Third Order Elastic Constants Using Guided Waves". In: *2019 IEEE International Ultrasonics Symposium (IUS)*. 2019, pp. 1264-1268.

874

875

876

877

[60] C. Cheng, H. Mei, V. Giurgiutiu, S. Yuan, F. Fang, and R. James. "Simulation of guided wave under varying temperature and load conditions". In: *Health Monitoring of Structural and Biological Systems IX*. Online Only, United States: SPIE, 2020, p. 19.

878

879

880

[61] R. Gorgin, Y. Luo, and Z. Wu. "Environmental and operational conditions effects on Lamb wave based structural health monitoring systems: A review". In: *Ultrasonics* (2020), p. 106114.

881

882

883

[62] A. Abderahmane, A. Lhémery, and L. Daniel. "Effects of multiaxial pre-stress on Lamb and shear horizontal guided waves". In: *The Journal of the Acoustical Society of America* 149.3 (2021), pp. 1724-1736.

884

885

886

887

[63] J. Chabassier and S. Imperiale. "Construction and convergence analysis of conservative second order local time discretisation for linear wave equations". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 55.4 (2021), pp. 1507-1543.

888

889

890

[64] H. Methenni. "Modélisation mathématique et méthode numérique pour la simulation du contrôle santé intégré par ultrasons de plaques composites stratifiées". PhD thesis. 2021.

891

892

893

[65] F. Ricci, E. Monaco, N. Boffa, L. Maio, and V. Memmolo. "Guided waves for structural health monitoring in composites: A review and implementation strategies". In: *Progress in Aerospace Sciences* 129 (2022), p. 100790.