Numerical modeling of wave propagation in anisotropic viscoelastic laminated materials in transient regime: Application to modeling ultrasonic testing of composite structures

Composite laminate structures remain an important family of materials used in cutting‐edge industrial areas. Building efficient numerical modeling tools for high‐frequency wave propagation in order to represent ultrasonic testing experiments of these materials remains a major challenge. In particular, incorporating attenuation phenomena within anisotropic plies, and thin intermediate isotropic layers between the plies often represent significant obstacles for standard numerical approaches. In our work, we address both issues by proposing a systematic study of the fully discrete propagators associated to the Kelvin‐Voigt, Maxwell, and Zener models, and by incorporating effective transmission conditions between plies using the mortar element method. We illustrate the soundness of our approach by proposing intermediate one‐dimensional and two‐dimensional numerical evidence, and we apply it to a more realistic configuration of a curved laminate composite structure in a three‐dimensional setting.

quantitative models of UT inspections. Essentially, it entails that direct solvers should address many specific physical behaviors, while providing essential results in a reasonable time.
The development of models fitting these requirements is the subject of numerous contributions. Although it is not the purpose to review them herein, one can at least mention the class of high-frequency asymptotic methods, [8][9][10] finite difference methods, 11 boundary element methods, 12 finite element methods in transient regime, 13,14 or more recent developments, see, for example, References 15,16. In this communication, the primary focus of our attention is the class of time-domain finite element methods. In order to achieve sufficient performances, a popular approach using finite element methods is to consider dedicated solvers built for specific hardware architectures. [17][18][19] Alternatively, authors in Reference 20 have proposed a specific block-structured strategy based on spectral finite elements and a domain decomposition approach, yielding enhanced performances on standard computers. Initially used locally around a defect to compute the diffracted field and echo response 21 within the CIVA platform, 22 this approach has been used in the context of the UT inspection of an entire laminate composite specimen. 23,24 In the context of time-domain wave propagation modeling within composite materials, one of the main challenges is to incorporate efficiently the attenuation effects in the numerical scheme. These effects are mainly due to the viscosity of the epoxy matrix and the multiple scattering of the incident wave by the reinforcing constituents embedded in this matrix. In addition, another challenge associated to the UT of laminates, is the modeling of the structural noise inherent to the stratified structure of the material. Usually, this aspect is taken into account by considering thin intermediate layers (with thickness of the order of several microns), much smaller than the characteristic wavelength of the source wave. A typical numerical approach, applied, for example, in Reference 20, is to address this issue by using finite elements with a lower order of approximation in the stratification direction. However, as long as the intermediate layers are meshed, the stability condition on the time step of any explicit time scheme necessarily deteriorates as their thickness decreases. This aspect may potentially leads to substantial computational costs in practical configurations.
Our work is dedicated to proposing relevant numerical strategies in order to address both challenges. In particular, the contributions of our work are twofold: • We compare fully discrete propagators of standard viscoelastic models of Kelvin-Voigt, Maxwell, and Zener 25 obtained after finite element space discretization and explicit time discretization. To do so, we propose a common calibration strategy for these models, applicable in the general case of anisotropic constitutive laws. For each model, we use energy related arguments to derive the stability condition of the fully discrete scheme. In the specific case of the Zener propagator, this approach is identical to the one proposed in Reference 26, however, we use a slightly different formulation yielding a more concise proof. From the expression of their stability conditions, we are able to analyze and compare the computational and memory costs of the three discrete propagators. This study enables us to address the question of choosing the relevant viscoelastic model, depending on the specifics of the configuration at hand while considering performance issues.
• In order to efficiently account for the intermediate layers, we adapt the mortar element method 27,28 to incorporate effective transmission conditions (ETCs) between viscoelastic models. We prove that the stability conditions obtained for each discrete propagators are not modified by the effective parameters of the interphase.
The outline of the article is as follows. In Section 2, we derive and compare the fully discrete propagators of the viscoelastic models. In Section 3, we study the numerical scheme obtained when using the mortar element method to incorporate ETCs between viscoelastic propagators. In Section 4 we illustrate our strategy in the context of the UT of curved carbon fiber reinforced polymer composite structures in a three-dimensional (3D) setting. Finally, we give some conclusions and perspectives in Section 5.

COMPARING FULLY DISCRETE VISCOELASTIC PROPAGATORS
In this section, we compare fully discrete propagators, that is, after time and space discretization, of the Kelvin-Voigt, Maxwell, and Zener viscoelastic models. 25 First, we start by proposing, for each model, a simple calibration procedure. In essence, this calibration aims at identifying complex stiffness modulii at a frequency of interest, and can be applied to general anisotropic materials. Second, we prove the stability conditions of the explicit time schemes associated to each model using energy arguments. Finally, we derive qualitative comparaison of the computational cost associated to the discrete propagators.

Identified viscoelastic properties at a central frequency of interest
We are interested in the propagation of high-frequency elastic waves within a given solid material Ω, an open subset of R d with d ∈ {1, 2, 3}. The waves are typically generated from an input excitation source, for example, a piezoelectric actuator, usually represented by a boundary load. In most of practical configurations, it is reasonable to assume that this excitation follows a relatively narrow-band signal, with a central angular frequency denoted by ⋆ . From a prior theoretical analysis and/or experimental campaign, we assume that it is possible to determine the mechanical material properties of the solid. 29 This characterization may lead to the estimation of the phase velocity and attenuation of incident plane waves, around the frequency of interest and for varying incident angles and polarizations, see, for instance, . This information enables to identify [34][35][36] the elastic and viscous behavior of the material. In the frequency domain, it amounts to obtaining the real and imaginary parts of the stiffness tensor appearing in the constitutive law, locally in the frequency spectrum, that is, around ⋆ . In addition, one may possess an estimation of the dependency of the attenuation w.r.t. the angular frequency, the variation of the phase velocity being sufficiently limited to be neglected in most cases. This aspect is usually leading the choice of the relevant viscoelastic relaxation tensor. 25,37 Throughout this article, we will consider exclusively the standard Kelvin-Voigt, Maxwell, and Zener models, leading to quadratic, constant, and locally linear attenuation laws, respectively, see Reference 25 or Appendix A. The mechanical properties of the material identified at ⋆ are: ⋆ the mass density,  ⋆ and  ⋆ the real and imaginary parts of the stiffness tensor. They are symmetric and positive definite fourth-order tensors, and we assume that they satisfy a "low-loss" assumption where  is the fourth-order identity tensor, and  ⋆ is referred to as the quality factor of the material. 25,38

Generalities on space discretization using lumped finite elements
In this section, we consider viscoelastic formulations written in terms of displacement vector field, or displacement and symmetric stress tensor. In the scope of our work, we restrict ourselves to primal-dual formulations. 39 The displacement variables are sought in V = [H 1 (Ω)] d , and the symmetric stresses tensor variable in W = [L 2 (Ω)] d ′ , with d ′ = 1 2 (d + 1)d. While the construction of a discrete space W h ⊂ W is naturally obtained using discontinuous elements, we assume that standard continuous finite elements 40 are used in order to build a H 1 -conform discrete space V h . Assuming that a conform tessellation  h of Ω is given, V h is defined in the following manner where F K is the geometric transformation of a reference finite element into the element K ∈  h , andV h is a finite dimensional functional space defined on the reference element. We focus our attention on traditional finite element procedures that considerV h as the space of (potentially high-order) Lagrange polynomials defined on simplexes, for example, triangles or tetrahedra, or on parallelotopes, for example, quadrangles or hexahedra. In the following, we will denote by N h and N ′ h the dimensions of V h and W h , respectively. For any function u * h ∈ V h , we denote by ⃗ U * ∈ R N h the vector regrouping the coefficients of the decomposition of u * h onto the basis of V h , idem for any functions in W h . In the perspective of applying such Lagrange finite element procedure to a time-domain wave propagation problem, we make the additional assumption that we can find a suitable quadrature formulae that enables consistent mass lumping. 14,[41][42][43] On this assumption, denoting by a strictly positive function and by  a symmetric positive fourth-order tensor field, the bilinear forms of the type are represented in the corresponding discrete spaces V h × V h and W h × W h by strictly positive diagonal matrices, trivially invertible.

Explicit energy related time discretization
In our work, we focus our attention on proposing explicit time schemes, meaning that, at every time steps, obtaining the next iterate only requires the inversion of one or several diagonal matrices. Generally, this class of time schemes are at best conditionally stable. The stability condition is often referred to as the CFL condition, and in practice it corresponds to an upper bound on the time step. Deriving the stability condition for this class of time discretizations may be performed in several ways. In the scope of the present communication, we will systematically favor the approach based on discrete energy arguments. This approach, extensively detailed in References 14,43, can be decomposed in the following steps: (1) identify a discrete functional satisfying a conservation or decay property; (2) prove that the functional is an equivalent norm on the discrete solution on a specific condition on the time step.
As an illustrative example, let us consider the case of inviscid elastodynamics, that is,  ⋆ = 0. Assuming free surface boundary conditions, the field equations read where ⋅ is the vector divergence operator, and (⋅) is the linearized Green-Lagrange tensor. The previous equations are completed with nonzero initial conditions. The weak form of this problem consists in finding, for any time t > 0, the solution u ∈ V such that for any u * ∈ V we have The mass bilinear form is defined as in (3) with = ⋆ , and the stiffness bilinear form reads In this context, a popular example of explicit time scheme is the so-called leap-frog scheme where the mass matrix M V ⋆ and stiffness matrix K  ⋆ are the discrete representations of m V ⋆ (⋅, ⋅) and k  ⋆ (⋅, ⋅) on V h × V h . Note that, using the mass lumping technique, the mass matrix is diagonal, leading to a fully explicit time scheme. As a first step for deriving the stability condition, we recall the following standard conservation property, obtained by multiplying (5), then the discrete energy functional satisfies the following conservation relation The second step is to find the condition ensuring that the functional (6) is an equivalent norm for the discrete solution of (5). To this goal, we recall the popular result: 14,39 Theorem 1. Denoting by r(⋅) the spectral radius of a given matrix, the scheme (5) is stable on the following condition Proof. Remarking that, we can recast the expression of the discrete energy functional (6) in the following manner .
Hence, on the condition (8), this functional is positive which necessarily implies the stability of the numerical scheme due to the conservation property (7). ▪ In order to better interpret this stability condition, and in the perspective of the comparative study between the various time schemes detailed in this section, we recall an estimation related to spectral radii as appearing in (8). Using the decomposition of the bilinear forms in (4) on the elements of the tessellation  h , and invoking arguments of change of variables from the reference finite element it is possible to prove, see, for example, Reference 39, that In the previous inequalityĈ( ⋆ , ⋆ ) is a constant depending exclusively on the physical parameters, the choice of the reference finite element and the order of the local Lagrange basis functions. Note that, in this context, h is a global upper bound on the size of the elements in  h . Using (9), a less sharp stability condition for the numerical scheme (5) reads and we recover the well-known result that the stability condition of explicit time schemes applied to hyperbolic systems depends linearly on the space step. For the numerical experiments conducted in this communication we use the sharper condition (8), the computations of the spectral radii being performed using a standard power iteration algorithm.

Kelvin-Voigt model
The Kelvin-Voigt model 25 consists in finding, for any time t > 0, the displacement field u ∈ V satisfying model, it will require some additional efforts to be extended to the Maxwell and Zener models. From the constitutive law in (11), we obtain̂( is the time-domain Fourier transform. Hence, using this relation at ⋆ we obtain relations (12). ▪ For any u * ∈ V, the weak formulation associated to (11) reads (14). The energy functional defined as satisfies the following decay property Proof. Multiplying (14) by 1 Introducing an order one estimate of the discrete velocity we can recast the last term of (17) in the following fashion Remarking that we can transform (17) in order to obtain the expression of the discrete energy defined in (15), yielding the sought decay property. ▪ From the decay property (16) satisfied by the discrete solution, we can derive the following stability condition on the time step.

Theorem 2. Defining r
, the time scheme (14) is stable on the following condition which implies that the time step should satisfy Proof. We aim at finding the sufficient condition for functional (15) to be an equivalent discrete norm. Using the following relation where || ⋅ ||  kv is the (semi)-norm associated to K  kv , the energy functional reads proving that (18) is a necessary condition for its positivity. Computing this only positive root of the second-order inequality, we obtain (19). ▪ Using the estimate (9) on the spectral radii and (18), a less sharp, albeit easier to interpret, upper bound on the time step is whereĈ  kv andĈ  kv depends only on the material constitutive law, the mass density and the choice of the reference finite element. Hence, the time step required for the numerical scheme (14) to be stable is asymptotically, that is, for small values of h, a quadratic function of the space step. This type of stability condition is similar to the one obtained when using explicit time scheme for parabolic systems. It is only natural to find a stability condition of this nature, since one can derive an explicit parabolic scheme satisfied by the velocity unknown when omitting the elastic part of (14). While the sensitivity to the space step is important, another key point of this stability condition is that it appears to depend linearly w.r.t. the quality factor. Indeed, using the definition (12), we have where C  ⋆ is a constant depending only on the quality factor defined in (1). Note that, the current analysis has been performed with a first-order uncentered scheme for the viscous terms, however, one can expect to obtain similar behaviors for higher order uncentered schemes. In practical configurations, these aspect represents a major difficulty since, for a fixed simulation time window, it can lead to significantly more time steps. However, one asset of (14) is that its implementation is very similar to the inviscid case (5) and only requires a second application of the "stiffness-like" matrix K  kv .

Maxwell model
The Maxwell system 25 are field equations written in terms of displacement and stress. Completed with nonzero initial conditions, it reads In (21),  mx and  mx are homogeneous to a compliance tensor and the inverse of a viscosity tensor, respectively. In other words, we define  mx et  mx two symmetric and positive tensors such that The weak primal-dual formulation associated to the Maxwell model consists in finding, for any time t > 0, the solution (u, ) ∈ V × W such that ∀(u * , * ) ∈ V × W we have The corresponding formulation, semidiscretized in space, reads Inspired from the previous work of Reference 26 carried out for the Zener model, we consider the following scheme , (26). The energy functional defined as satisfies the following decay property Proof. Multiplying the first relation in (26) by 1 2Δt ( ⃗ U n+1 − ⃗ U n−1 ) ⊺ , we obtain a first conservation relation In order to obtain an equivalent expression of the second term, we multiply the second relation in (26) by Computing the average of this last relation at times t n+1/2 and t n−1/2 yields 1 4Δt Following the recombination given in Reference 44, we can extract the following expression which leads to ) .
Finally, we conclude by combining the previous relation with (29). ▪ As for the previous case of the Kelvin-Voigt model, a stability condition for the numerical scheme (26) is obtained by deriving the condition for the energy functional (27) to be a positive functional.

Theorem 3. Defining the following stiffness matrix
the time discrete scheme (26) is stable on the following condition Proof. To start with, we use the following relation in order to separate the energy functional (27) in the following manner Note that the first term is a positive term and that the second term can be written in the following form where || ⋅ ||  mx is the discrete (semi)-norm associated to the matrix K  mx . We use the definition (31) of K  mx and the square of binomials in order to rewrite the previous expression as Hence, the energy function is a positive function under the sought condition on the time step. ▪ Using the estimate (9), we can verify that the stability condition obtained for the Maxwell model varies linearly with the space step, as for the inviscid case (10). In addition, using the calibration (23), one can remark that Hence, from the estimate on the spectral radii (9), the time step related to the Maxwell discrete propagator behaves as so that the difference compared with the inviscid case is limited in a low-loss context. In terms of computational cost, the major computational load in (26) consists in applying B and B ⊺ at every time steps, which is roughly equivalent to the application of a single stiffness matrix. In this regards, the computational cost required in order to obtain the next iterate in the Maxwell scheme is similar to the inviscid case. It requires, however, the storage of the stress unknowns in addition to the displacement unknowns, thus increasing its memory load.

Zener model
The Zener model 25,26 is defined in terms of displacement and stress unknowns, satisfying the following field equations completed with nonzero initial conditions. In (34),  zn et  zn are symmetric positive fourth-order tensors representing elasticity and viscosity, respectively. The parameter zn > 0 corresponds to a relaxation time of the Zener model. 25 Lemma 6. Let ( ⋆ ,  ⋆ , ⋆ ) be the mechanical properties obtained at a frequency of interest ⋆ , we define a suitable calibration of the Zener model (34) as Proof. Applying the Fourier transform in time domain to the second equation in (34) we obtain or, equivalently, Thus, choosing zn = −1 ⋆ yieldŝ( Identifying the real and imaginary parts with the target elasticity and viscosity tensors, we obtain the following pair of relations leading to the calibration result (35). ▪ Note that, on the low-loss assumption (1), one can verify that the tensor  zn expressed in (35) is a positive tensor. Furthermore, the tensor  zn −  zn is also positive. This condition has been identified in Reference 26 as a necessary condition in order to ensure the well posedness of the Zener model. As in Reference 26 we introduce an internal variable s ∈ W defined as the purely viscous part of the stress tensor, namely, The Zener model written in terms of displacement and internal variable reads where  zn and  zn are two symmetric positive tensors such that Hence, the weak primal-dual formulation associated to the Zener model consists in finding, for any time t > 0, the , and m W  zn (⋅, ⋅) are defined as in (3), k  zn (⋅, ⋅) is defined as in (4), and b(⋅, ⋅) is defined in (25). As for the Maxwell model, and also inspired from Reference 26, we consider the following time scheme It is a second-order scheme, the first and second equation being centered in t n and t n+ 1 2 , respectively. As previously, we start by giving a first result of discrete energy decay in order to derive in a second step the stability condition associated to (39).

Theorem 4.
Denoting by || ⋅ ||  zn and || ⋅ ||  zn the discrete norms associated to the matrices M W  zn and M W  zn , respectively, we define the following functionals , (39). The energy functional defined as satisfies the following identity Proof. Multiplying the first relation in (39) After multiplication of the second relation of (39) by 1 2 ( ⃗ S n+1 + ⃗ S n ) ⊺ , we compute the mean values obtained at times t n+1/2 and t n−1/2 , then, using the identity (30) (with replacement of the stress unknown by the internal variable unknown), we obtain which, combined with (43), leads to the relation (42). ▪ Note that, even though its proof differs, the previous result is identical to Theorem 4.1 of Reference 26. Using this energy identity, we can deduce the following stability result. Theorem 5. Defining the following stiffness matrix the time scheme (39) is stable on the following condition Proof. Replacing stress unknowns by internal variables unknowns in relation (33) in order to reformulate P such thatẼ We use the square of binomials in order to recast the second term in the previous expression intõ Finally, using the relation (20) applied to P n+ 1 2  zn , we obtaiñ where || ⋅ ||  zn stands for the (semi)-norm associated to K  zn . Hence, the discrete energy is positive under the condition (45). ▪ Using the estimate (9), we can verify that the stability condition obtained for the Zener model varies linearly with the space step, as for the Maxwell case and inviscid case. In addition, using the calibration (35), one can remark that so that, from (9), we have While the effect of the viscosity on the time step is expected to be limited in a low-loss context, it is interesting to remark that the Zener discrete propagator appears to be more sensitive to the viscosity compared with the Maxwell case. In addition, note that the main part of the computation load in (39) comes from the application of K  zn , B, and B ⊺ at every time steps. This corresponds roughly to the application of two stiffness-like matrices. Hence, the computation of a next iterate in the Zener scheme is similar to the Kelvin-Voigt case, and is two times more expensive than the Maxwell and the inviscid case. However, the memory load required for (39) is similar to the Maxwell discrete scheme.
comparison of fully discrete explicit viscoelastic propagators In Table 1, we propose to summarize various characteristics of each fully discrete propagators. As a direct consequence, for a given set of target mechanical properties, one can expect to have the following relation between the stability conditions As expected, these stability conditions are identical when  ⋆ = 0. Note that, in the specific case of  ⋆ = , we should expect to have Δt mx = Δt zn .

Numerical illustration on a one-dimensional configuration
The goal of this section is twofold. First, we provide numerical evidence of the specific behaviors, detailed in Table 1, of the stability conditions of the numerical schemes. Second, we illustrate how the numerical schemes can adequately represent the various characteristics of the Kelvin-Voigt, Maxwell, and Zener models in terms of attenuation laws and phase velocities. For the latter item, we wish to rely on comparisons with analytical results, and to do so we restrict ourselves to the one-dimensional (1D) case. Let us consider the following domain Ω =]0 ; L[, with L = 30 mm. We set the frequency of interest to f ⋆ = 4 MHz, with ⋆ = 2 f ⋆ , the mass density to ⋆ = 8 g ⋅ cm −3 , and the target elasticity modulus to C ⋆ = 288 GPa. The velocity of waves in the inviscid case is set to v P = 6 mm ⋅ s −1 , which roughly corresponds to the velocity of longitudinal bulk waves in steel materials.

Sensitivity of the stability conditions w.r.t. space discretization and quality factor
In Figure 1A, the mesh step is set to h = 0.1 mm, and we use polynomials of order 5 within each elements. We allow the target imaginary part of the stiffness coefficient to vary in the interval D ⋆ ∈ [10 GPa; 288 GPa]. For each values, we plot the associated CFL conditions in log scale. We clearly observe that the hierarchy (46) between the time steps is respected. The CFL conditions in the inviscid, Maxwell, and Zener cases are close to each other, while the ones associated to the Kelvin-Voigt propagator are significantly lower and following a linear tendency as expected. Note that, in the particular case of  ⋆ = , one can verify that Δt mx = Δt zn . In Figure 1B, the target imaginary part is set to D ⋆ = 10 GPa, the mesh step varies in the interval [0.015; 0.1] mm, and for each values we plot in log scale the stability conditions of each propagators. We can clearly observe the linear behaviors in the inviscid, Maxwell, and Zener cases, while the Kelvin-Voigt CFL conditions follow a quadratic dependency.

Attenuation laws and phase velocities
In the context of 1D viscoelastic propagation models, we denote by k ∈ C a complex wave number related to a monochromatic 1D plane wave. We decompose k = − i , so that is directing the attenuation of the plane waves. The dependency of w.r.t. to the angular frequency ∈ R + is referred to as the attenuation law of the viscoelastic model. In the following, we propose to evaluate numerically these attenuation laws for each discrete propagators, along with the phase velocities. To do so, we set the target viscosity to D ⋆ = 10 GPa. We impose the solution at x = 0 to be a wide-band signal, typically a Ricker's wavelet, centered in the frequency domain at f ⋆ . At x = L, we use absorbing boundary conditions to avoid undesired reflections.
By performing a Fourier analysis of the numerical solution at an observation point, set to x = L∕2, we are able to evaluate the "numerical" attenuation law obtained from the various fully discrete propagators. Using the expression of the exact attenuation law provided in (A2), we compare in Figure 2A the expected attenuation (in solid, dashed, and dotted lines) with the numerical ones (in square markers). A first remark is that the curves for the three-viscoelastic models intersect exactly at the frequency of interest, proving the validity of our calibration procedure. Second, we have an excellent agreement between the expected and the numerical attenuation laws, validating the consistency of our complete numerical procedure with respect to the continuous equations. Indeed, we are able to reproduce the quadratic, constant, and locally linear attenuation law corresponding to the Kelvin-Voigt, Maxwell, and Zener, as recalled in Appendix A.
In the same way, we recover the phase velocity of the discrete propagators. In Figure 2B, we compare the dispersion curves obtained numerically with the exact dispersion curves using the expression of the exact phase velocities for viscoelastic models provided in (A2). We remark that the limited dispersion of the discrete propagators is in very good agreement with the theoretical ones and coincide at the frequency of interest. However, we should point out the specific case of the Kelvin-Voigt propagator that appears to incorporate spurious numerical dispersion at high frequencies, most likely due to the fact that it is a first-order time scheme.

Choosing the relevant viscoelastic propagator
To the best of the authors' knowledge, a popular choice of viscoelastic model is the Kelvin-Voigt model (11). A possible explanation for this tendency is the rather simple form of the constitutive law, yielding a direct link with the complex stiffness coefficients obtained from the material characterization. Note that this direct link clearly appears in the straightforward calibration given in (12). Another potential argument that may explain the popularity of the Kelvin-Voigt model is its simplicity of implementation. In practice, it simply amounts to applying a second stiffness operator, usually already implemented for performing iterations of the inviscid model (5). Nonetheless, the detailed analysis proposed in this section has shed light to significant issues related to the computational performances of the Kelvin-Voigt fully discrete propagator. Hence, as a conclusion of this section, we discuss the rationale for choosing the relevant viscoelastic model for a given configuration. An important discriminant characteristic of the viscoelastic models is their related attenuation law, as depicted in Figure 2A. For the Zener and Kelvin-Voigt model, the attenuation grows, respectively, as a (locally) linear and a quadratic function of the frequency. As a consequence, this particular aspect may entail a shift of the ultrasonic wave field toward lower frequencies. We show a concrete example of this effect in the Figure 10A of the final result Section 4. These shifts may become more visible as the distance of propagation increases, or equivalently as the quality factor decreases. Hence, if the main expected results of the numerical modeling is related to the spectral content of the outputs, then the choice of the viscoelastic model should be driven by its attenuation law.
Otherwise, a more practical argument may be used. At this point, it is important to emphasize that, by construction, the calibration procedure proposed in this section ensures that the three models recover an identical attenuation at the frequency of interest. Hence, assuming that it is sufficient, one can choose the viscoelastic model depending on the performances associated to the fully discrete scheme. In this regard, we refer to Table 1 to infer a strict hierarchy between propagators. With its low dependency of the CFL condition w.r.t. the quality factor, and with only one stiffness operation per time step, the Maxwell propagator (26) is nearly as effective as the inviscid propagator (5). The Zener propagator also exhibits a robust CFL condition w.r.t. the quality factor. However, it costs another stiffness operation, yielding roughly a doubled iteration cost compared with the Maxwell propagator. Similarly, the Kelvin-Voigt model presents two stiffness operations per time step. More importantly, Figure 1A,B illustrate the strong sensitivity of its CFL condition w.r.t. the quality factor and the mesh step. Hence, the timestep of the Kelvin-Voigt model can be in practice significantly smaller, leading to an increase number of iterations to reach a fixed final time. Due to the explicit uncentered treatment of the viscous term in (14), we have obtained a first-order time scheme. Thus, in addition to its cost, the Kelvin-Voigt propagator may lead to less accurate results in some cases. This hierarchy, illustrated in a practical case of 3D computations in the Tables 3 to 5, clearly suggests that, if the choice of the viscoelastic propagator is driven by computational performances, the best choice is the Maxwell model.

FULLY DISCRETE PROPAGATORS OF ANISOTROPIC VISCOELASTIC LAMINATED MATERIALS
The goal of this section is to detail a fully discrete propagator for a class of laminated materials. In each layer (or ply) of the laminate we assume that the material follows an anisotropic viscoelastic constitutive law. Between each pair of plies, we consider thin intermediate layers of viscoelastic materials. In practice, these intermediate layers play an important roles in the modeling of laminated structures, since they can be used to represent structural noise or specific bonding conditions between plies.
From a computational standpoint, addressing this type of configurations is challenging, especially in a 3D context. Apart from the construction of the viscoelastic discrete propagators for the plies, subject treated in the previous section, one of the main difficulty is the efficient incorporation of the thin intermediate layers within the numerical solver. A reasonable approach, proposed, for instance, in Reference 20, is to use finite elements with anisotropic orders to limit the cost overhead.
In the scope of our work, we focus our attention on an alternate approach based on ETCs. [45][46][47] In this context, the effect of a thin intermediate layer (or interphase) is taken into account through so-called spring-mass conditions. In practice, this approach circumvent the difficulty of meshing the interphase medium, potentially yielding significant improvements in computational costs. The key points of this section are: • After recalling the standard elements of ETCs and their calibration w.r.t. the target complex stiffness modulii, we show how they can be incorporated to viscoelastic discrete propagators using the mortar element method. 27,28 • We prove that the incorporation of these ETCs does not affect the stability condition computed (independently) for each ply. In particular, this condition on the time step does not depend on the material properties of the interphase.
To simplify the presentation of the coupling scheme obtained from the mortar element method, we consider the specific case of two plies sharing one interphase, and we assume that each plies support a Zener formulation (34). The extensions to Kelvin-Voigt or Maxwell models with multiple intermediate layers are straightforward.

Plane and curved geometry of a laminate material
and each subdomain geometries are defined by In practical cases, the material at hand has a curved geometry. This geometry is usually obtained during manufacturing processes by deformation of a plane geometry. A convenient way to model these configurations is to consider a given mapping Φ, yielding the resulting geometries In order to ensure a geometry suitable for standard numerical procedures, we make some regularity assumption on this transformation. In particular, we assume that Φ is a diffeomorphism, that is, it is a differentiable bijection and its inverse is also differentiable, and that Φ is a conformal map, that is, it locally preserves orientations. On these assumptions, we can define in Ω the local orthonormal basis {e k } d k=1 satisfying where • is the gradient operator w.r.t. the coordinate system in Ω • applied to vector fields, and || ⋅ || is the Euclidian norm. With these notations, let  • i be a fourth-order tensor defined on a ply Ω • i , the resulting tensor  i defined in the curved geometry is given by

Upper surface, lower surface, and mid-surface of intermediate layers
We denote the boundaries of the plies in contact with the intermediate layer by ETCs aim at linking together the mechanical fields, that is, normal stresses and traces, at the mid-surface of the intermediate layer. In the following we denote by Γ this mid-surface, and we define ± ∶ Γ → Γ ± , the geometrical transformations from the mid-surface to the upper and lower surfaces. Note that, from the regularity assumptions on Φ, these transformations are diffeomorphisms, and in the particular case where Φ is the identity transformation, they boil down to simple translations of ± 1 2 •, in the stratification direction.

Material properties of the plies and intermediate layers
In the following, the mechanical properties of the plies and the intermediate layers identified at ⋆ are denoted by where ⋆,i (or ⋆ ) are the mass densities, the real part  • ⋆,i (or  •, ⋆ ) and imaginary part  • ⋆,i (or  •, ⋆ ) of the stiffness tensors are symmetric and positive definite fourth-order tensors, satisfying the low-loss assumption (1) and defined on the prior geometry Ω • . We consider general anisotropic stiffness tensors for the plies, and we denote by  ⋆,i and  ⋆,i the resulting tensors in (47). Concerning the intermediate layers we make the assumption of isotropic materials. Their corresponding mechanical properties in the curved configuration are simply defined by the Lamé coefficients C, ⋆ , C, ⋆ for the real part, and D, ⋆ , D, ⋆ for the imaginary part.

Definition of standard spring-mass ETCs
For every subdomains Ω i , we denote by i the corresponding stress tensor field, and we associate to each interphase the following mechanical fluxes where n is the normal vector field of the mid-surface Γ oriented from Γ − to Γ + . The mean value and jump of the mechanical fluxes at the mid-surface are denoted by Introducing the functional space X = L 2 (Γ), we assume enough regularity on the stress tensor fields so that ⟨ ⟩ ∈ X, and [ ] ∈ X.
In the same manner, considering the displacement field u i associated to each subdomains Ω i , we associate to the interphase between the plies the mean value and jump of the trace For the inertial terms authors in References 49,52 consider the following definition Proposing an expression of these coefficients, typically using an asymptotic expansion w.r.t. the ratio of the thickness over the wavelength, is still an active field of research. For instance, spring-mass models can be used to represent complex configurations such as damaged interfaces 53 or distributions of interfacial flaws. 54,55 It should be noted, however, that the ETCs given in (49) are an approximation of the rigorous asymptotic model, derived in Reference 56. More precisely, they correspond to neglecting terms related to the surface gradient of the mean and jump of the trace, which we can expect to be a reasonable approximation in the case of small contrasts.

ETCs through mortar elements
The mortar element method 27,28,57,58 is a domain decomposition approach aiming at coupling independent space discretizations. It has been successfully applied to elastodynamics and wave propagation problems. [59][60][61][62][63] More recently, 20 this method has been used to couple different formulations, for example, elastodynamics, acoustics, and corresponding perfectly matched layers [64][65][66] formulations. In the present work, we extend this approach to viscoelastic formulations coupled with spring-mass ETCs. To start with, let us introduce the following Lagrange multipliers We define the ETCs (49) written in terms of Lagrange multipliers as where I is the identity second-order tensor. Interestingly enough, we remark that this system is identical, modulo the dimension of the tensors, to the system (24) (36). Prior to propose the weak formulation of the mortar-coupling scheme, we start with the following lemma expressing the boundary terms appearing in the Green identity. Lemma 8. Let us denote by n ± the normal vector fields of the boundaries Γ ± , outgoing from Ω 1 and Ω 2 , and by dΓ ± the surface measures. Assuming that the intermediate layer is sufficiently thin so that the following developments can be reasonably truncated in order to keep their leading-order terms exclusively, then we have for any u * 1 ∈ V 1 and u * 2 ∈ V 2 .
Proof. For any u * 1 ∈ V 1 , discarding terms in O( ) in (55) yields Remarking that the definition (52) of the Lagrange multiplier entails we obtain the first relation in (56). For the second relation, by discarding the terms in O( ) in (55), we obtain and we conclude using again (57). ▪ After multiplication of the Equation (37) by test functions, integration, and application of Green identities with the boundary terms expressed using (56), we obtain the weak formulation of the coupled system. For the volume fields, it consists in finding (u 1 , s

for every t > 0 and for any test functions
in corresponding functional spaces, we have In addition, multiplying (53) by test functions and integrating, we derive the weak formulation satisfied by the Lagrange multipliers ( , ) ∈ X × X for any ( * , * ) ∈ X × X, We distinguish several types of bilinear forms emanating from the imperfect transmission conditions. We define the coupling bilinear forms between the Lagrange multiplier and the volume unknowns as In addition, we consider the following bilinear forms operating solely on the Lagrange multipliers Let us introduce X h ⊂ X a discrete approximation space for the Lagrange multiplier, and we denote by M h its dimension. The previous bilinear forms are represented in their corresponding discrete approximation spaces by the following finite element matrices the set of vectors regrouping the components of discrete Lagrange multipliers in the basis of X h . Using these notations, we consider the following fully discrete scheme for (58) which is centered in t n and t n+ 1 2 for the displacements and internal variables equations, respectively. In addition, for (58), we propose the following scheme, centered in t n , and we consider the following energy functional The solution of (62) and (63) satisfies the following energy identity Proof. To start with, we follow the demonstration leading to the identity energy (42) for each subdomains. Thus, multiplying the first and second relation in (62) by 1 2 ) ⊺ , respectively, and multiplying the second relation in (63) by 1 First, let us remark that Second, using the first relation in (63), we can give the following expression for the last terms in (67) 1 4 leading to (66) using the square of binomials. ▪ Note that the energy functional (65) representing the kinetic and potential energy associated to the interphase is unconditionally positive w.r.t. the time step. Hence, deriving a stability condition for the discrete coupling scheme can be simply done by insuring that energy functionals associated to each subdomains are positive, which is exactly the result provided in Theorem 5. More precisely we can give the following stability result. Theorem 6. For a subdomain index i ∈ {1, 2}, we denote by K  zn i the matrices associated to the discrete Zener propagators in Ω i , defined as in (44). The coupling scheme (62) and (63) is stable on the following condition In other words, the introduction of the ETCs using the mortar element method does not influence the stability condition of each subdomains, for any values of the effective mechanical parameters of the interphase.

Resolution algorithm using the Schur complement method
In order to give practical solving steps for the numerical scheme (62) and (63), we start by separating the computation of the internal variables from the displacement unknowns. This is possible since the internal variables appear in the displacement equations in an explicit manner. The coupled system satisfied by the displacement unknowns and the Lagrange multipliers can be expressed as where the various right-hand sides of this system are defined as We define the following matrices representing the block matrices of the system (69). With these notations, the Schur complement matrix is defined by Introducing two intermediate variables ⃗ U ⋆ 1 and ⃗ U ⋆ 2 , we can obtain, from applying a Schur complement method, the following solving steps: .
Note that, using mass lumping technique the matrix A is diagonal thus trivially invertible. However, for this strategy to be relevant, we need to verify that the Schur matrix S is invertible.

Lemma 10. The Schur matrix defined in (70) is symmetric and positive definite.
Proof. Remarking that B = − Δt 2 2 C, we have proving that the Schur matrix is symmetric. The expression of this matrix is where ⋅ represents the block matrix deduced by symmetry. Let ⃗ Λ and ⃗ M be two nonzero vectors of R M h and define we can verify that

Lumped Schur complement matrix in the conform case
In practice, we define the discrete space X h as one of the volume trace space, that is, This approximation is similar to (55), and corresponds to neglecting the effect of the deformation − on the surface mesh. Introducing the following bilinear form and its corresponding finite element representation M Γ ∈  M h ,M h (R), the definition (72) of X h entails some simplifications of the expressions of the bilinear forms (60). In particular, denoting by T 1 ∈  M h ,N 1,h (R) the matrix selecting the values on Γ − of a finite element vector in V h,1 , we have Furthermore, on the assumption of conform interfaces, that is, we can extend the previous expression to bilinear forms applied on elements of V h,2 . More precisely, if we define T 2 ∈  M h ,N 2,h (R) the matrix selecting the values on Γ + of a finite element vector in V h,2 , we have With this specific decomposition of the coupling bilinear forms, and assuming that a mass lumping technique is used when assembling M Γ , we can verify, from the expression (71), that the Schur matrix is a block diagonal matrix, each block being of size 2d × 2d. Therefore, it can be efficiently inverted, each block problems in parallel. In other words, in the conform case, the computational overhead of the introduction of imperfect transmission conditions using the mortar element method is very limited.

3.7
Numerical illustrations on a two-dimensional configuration where ≤ 1 is a contrast coefficient. Denoting by v ⋆,L and v ⋆,T the velocity of longitudinal and transverse waves in the material, and by v ⋆,L and v ⋆,T the corresponding velocities in the interphase, the previous expression (73) of the Lamé coefficients were derived so that On the upper surface of the domain, we impose a normal surface load with a centered support of 2 mm with the following time dependency with f ⋆ = 6 MHz, t 0 = 0.65μs and f = 0.2 s −1 . In each subdomain, we set the identical mesh steps of h x = h y = 0.125 mm and use finite elements of order 4. For clarity purposes, we illustrate this simple configuration in Figure 3A.
To start with, we consider the case of a relatively high contrast of = 1 5 and an interphase thickness of = 10 μm. We show in Figure 3B a snapshot of the solution obtained at the arrival time of the first longitudinal wave packet. In this snapshot, the intermediate layer is meshed with one element of order 2 in the thickness direction. At the abscissa x = 8 mm we consider an observation line. In Figure 4A,B, we compare the first and second components of the displacement field obtained either by meshing the intermediate layer or by applying ETCs. We observe a very good agreement between the two solutions Figure 5.
In a second numerical experiment, we impose the same contrast coefficient but we increase the thickness of the intermediate layers to = 100 μm. We show in Figure 3C the snapshot of the solution obtained with a meshed intermediate layer using finite elements of order 3 in the thickness direction. Compare to the previous case we observe a significant longitudinal reflection and some contributions trapped within the interphase. We also extract the displacement fields at an observation line at x = 8 mm and compare the numerical solution obtained from meshed layers and ETCs. While we keep a reasonable agreement for the second component, the first component, gathering most of the transverse contributions of higher frequency, are not adequately captured using ETCs. In particular, the solution obtained from ETCs does not manage to represent the spatial variation around the thin layer location, around y = 5 mm.
To better illustrate the evolution of the error introduced by the effective transmission conditions, we propose to consider the following quantity depending on the contrast coefficient and the thickness of the intermediate layer where u n h is the numerical solution obtained with meshed interphase and u ,n h is the one obtained using ETCs. In Figure 6, we plot this relative error for ∈ [10 −3 ; 10 −1 ] and ∈ { 1 5 , 1 3 , 1}. We clearly observe that, for a fixed thickness, as the contrast coefficient decreases, the relative error increases. Interestingly enough, for every tested values of the contract coefficient, we observe a linear asymptotic behavior w.r.t. . Note that authors in Reference 56 derive an asymptotic . Nonetheless it requires surface gradient terms, which are more involved computationally and, in particular, it would necessarily leads to sparse Schur matrices, even with conform interfaces.

NUMERICAL ILLUSTRATION: APPLICATION TO UT OF A LAMINATED COMPOSITE
In this section, we illustrate our modeling strategy by considering a more realistic 3D configuration. In this configuration, we assume that the specimen under study is composed of 20 anisotropic viscoelastic plies and 19 intermediate layers. The geometry of the domain, prior to its deformation, is characterized by the following lengths L 1 = L 2 = 10 mm, L 3 = 4.995 mm, •, = 5 ⋅ 10 −3 mm, so that the plies have a thickness of 245 ⋅ 10 −3 mm in the "plane" geometry. The mass density of the plies is set to ⋆,i = 1.6 g ⋅ cm −3 . For the subdomain with odd or even indexes, the real part of the target constitutive law reads, using the Voigt notation, The global transformation Φ of the geometry is a simple cylindrical transformation of radius 10 mm and angle 4 rad. This transformation is affecting the definition of the mechanical properties of the plies in virtue of (47). On the upper surface of the domain, we impose a normal surface load with a centered circular support of radius 2 mm. The time dependency of this surface load is as in (74), with F = 5 MHz, t 0 = 0.65 μs, and f = 0.25 s −1 . We consider a maximal time of 2 μs, which roughly corresponds to the time required for the principal wave packet to reach the lower surface of the domain. Concerning the space discretization of the domain, in the (x, y)-plane we set the mesh step to 285 ⋅ 10 −3 mm and we use finite elements of order 4. Within each ply's thickness we use two elements of order 4, and in the case of meshed intermediate layers we use one element of order 2.
As a first step, we neglect the viscous part of the target constitutive law and we compare the inviscid model with either meshed intermediate layers or ETCs. We show in Figure 7 a snapshot of the numerical solution obtained in the former case, and in Figure 8A we plot the third component of both solutions over the line x = y = 5 mm. We observe a very good agreement between the two numerical solutions. In Table 2, we compare various characteristics of the simulations, carried out on a laptop PC Intel (R) Core TM i7 − 8850H at 2.6GHz. We observe that we have close to a factor three between the timesteps, leading to a factor three in terms of computational time. In terms of memory load we observe, however, a factor two when using ETCs, which is due to the storage of the inverse of the Schur matrices at each interface between plis.
In a second step, we compute the solution obtained with a Kelvin-Voigt, a Maxwell, and a Zener model. In the three cases we use ETCs to model the intermediate layers between plies. In Tables 3 to 5, we gather the various computation costs associated to each viscoelastic propagator. We observe, on this practical case, the hierarchy described in Section 2.8: the Maxwell propagator is as efficient as the inviscid case; the Zener model exhibits a doubled response time, mainly due to the second application of a stiffness operator; the Kelvin-Voigt is roughly four times more expensive than the Maxwell propagator due to two stiffness operations, and a smaller timestep.
We show in Figure 7 the snapshots obtained from the Zener model at specific times, to be compared with the inviscid case. In Figure 8B  In addition, we consider an observation point of coordinates (x, y, z) = (5, 5, 0) situated at the bottom of the observation line. After increasing the maximal time to 4 μs, we plot over time in Figure 9A the solutions of the inviscid and the three-viscoelastic models. For visibility purposes, we plot in Figure 9B only the signals obtained from the viscous propagators. These plots demonstrate the validity of our calibration procedure since the viscoelastic solutions give similar results.
To better understand their differences, we analyze the frequency content of the main wave packet. To do so, we select the numerical signals up to time t = 3 μs and apply a simple zero padding. On these windowed signals, we apply a   Figure 10A. We clearly see the global effect of the attenuation, which is equivalently recovered by the three-viscoelastic models. At the frequency of interest-5 MHz for this specific case-the corresponding three curves match. Away from the frequency of interest, we observe some frequency shifts. As detailed in Section 2.8, these shifts are due to the frequency dependency of the attenuation law associated to each model, namely: constant, (locally) linear and quadratic functions for the Maxwell, Zener, and the Kelvin-Voigt case, respectively, see Figure 2A. From the frequency content of the signals, we are able to compute by a simple spectral division the "numerical" attenuation laws associated to each model for the main wave packet. We present them in Figure 10B. The three curves match at the frequency of interest, and we recover the theoretical dependencies of the attenuation laws of the three models.

CONCLUSION AND PERSPECTIVES
In this work, we have proposed a systematic approach to compare explicit fully discrete wave propagators of standard viscoelastic models. By proposing a common and simple calibration procedure for these models, applicable in the general context of anisotropic constitutive laws, we were able to provide qualitative analysis of the corresponding numerical schemes. In particular, we have studied how the stability conditions may evolve w.r.t the quality factor or the spatial discretization. This general analysis has been illustrated with various numerical solutions.
In the second part of this communication, we have shown how one can integrate ETCs to represent isotropic intermediate layers between curved anisotropic viscoelastic domains using the mortar element method. Based on the stability analysis of the viscoelastic models carried out in the first part, we have proven that the introduction of the mortar elements does not modify the stability condition of the coupled scheme, this for any effective parameters of the transmission conditions. In addition, even though this approach is valid for nonconform space discretization, in the specific case of conform interfaces we have obtained an explicit algorithm. We have proposed 2D illustrations providing numerical evidence of the validity of this strategy, and a numerical study on the convergence of the solution with respect to thickness of the intermediate layer. In a 3D context, we have considered the case of wave propagation within a curved viscoelastic laminated material. In this case, the use of ETCs to represent the epoxy intermediate layers has yielded a significant improvement, in terms of computational costs, compared with the standard numerical solution obtained by meshing these layers.
Several perspectives related to the presented work may be envisioned. Concerning the comparison of the discrete propagators, we may consider the generalized Zener model, which in essence is an extension of the model (34) with multiple relaxation mechanisms. In this context, the main difficulty would be to propose an explicit calibration of this model applicable for the anisotropic case. In relation to this issue, we may refer to the recent advances of References 68,69. Concerning the incorporation of ETCs, it seems possible, albeit more involved, to incorporate more elaborate asymptotic models such as the one proposed in Reference 56. Using this type of transmission conditions, yielding better convergence w.r.t. the intermediate layer thickness, may be interesting in order to address a larger range of UT configurations. A rather more direct perspective is to consider the coupling between a fluid domain (acoustic formulation) and a solid domain (elastodynamic formulation) with standard ETCs. A typical example of application is the simulation of UT of immersed and coated materials. Note that, in this perspective, being able to use nonconforming space discretization thanks to the mortar element method is a major asset. Another natural perspective is the use of these ETCs to model delamination flaws, typically appearing after impacts on laminated composite structures. The major goal of this approach would be to incorporate flaws of arbitrary shapes without penalizing the computational cost of the numerical solver. As a more general perspective, it seems promising to rely on the fact that the coupling scheme does not depend on the effective parameters of the interphase. This robustness may be put to good use within an inverse or imaging process. For example, using specific ETCs to represent delamination flaws, we may rely on variational data-assimilation approaches [70][71][72] or topological gradient imaging [73][74][75] in order to reconstruct specific parameters of the flaws.
Using (A5), we obtain  1 v c ( ) and we conclude by remarking that |M( )| = M 1 ( ) From (A2), we can verify that the ratio between the attenuation and the (real) wavenumber is so that we recover the result given in the Section 2.3 of Reference 25. Using the expression (A2) of the attenuation law, we can propose the following characterization of the standard viscoelastic models.
Proof. From (A2), the "low-loss" limit of the attenuation law is The complex modulus associated to the Kelvin-Voigt model reads where  kv > 0 and  kv > 0 are scalar values in a 1D context. Using (A7) we have proving the quadratic dependence of the attenuation w.r.t. the frequency. For the Maxwell case, the complex modulus is given by yielding the following real and imaginary parts and quality factor M mx 1 ( ) = 2 mx  mx 1 + ( mx ) 2 , M mx 2 ( ) =  mx 1 + ( mx ) 2 , Q mx ( ) = mx .
The expression of the limit (A7) reads Furthermore, using a high-frequency asymptotic we observe that the Maxwell model follows a constant attenuation law Finally, concerning the Zener model, the complex modulus is M zn ( ) =  zn + i zn  zn 1 + i zn , with real and imaginary parts defined as follow The corresponding quality factor is given by The high-frequency asymptotic of the real part of the complex modulus is constant w.r.t. the frequency, Hence, assuming low-loss and high-frequency asymptotics, the attenuation law of the Zener model satisfies The derivative of the quality factor w.r.t. the angular frequency is dQ zn d ( ) = zn ( ( zn ) 2  zn −  zn )