Matthieu Gallezot
email: matthieu.gallezot@ifsttar.fr

Fabien Treyssede
email: fabien.treyssede@ifsttar.fr

Laurent Laguerre
email: laurent.laguerre@ifsttar.fr

F Treyssède

A modal approach based on perfectly matched layers for the forced response of elastic open waveguides

Keywords: open waveguide, forced response, perfectly matched layer, orthogonality, PML modes, leaky modes

come

A modal approach based on perfectly matched layers for the forced response of elastic open waveguides

Introduction

Context and state-of-the-art

Elastic guided waves are interesting for many applications involving elongated structures (e.g. nondestructive evaluation (NDE), structural health monitoring (SHM), exploration geophysics...), because of their ability to propagate over large distances. When the structure (the core) is embedded into a large solid matrix, it can be considered as an open waveguide (unbounded in the transverse direction). Such a configuration typically occurs in civil engineering and in geophysics for example. Contrary to waveguides in vacuum (closed waveguides), most of the waves in open waveguides are attenuated by leakage in the surrounding medium as they propagate. Moreover, the underlying physics is deeply transformed.

Indeed, in open waveguides, three main kinds of modes are distinguished: trapped modes, radiation modes, and leaky modes. Trapped modes propagate without leakage attenuation along the waveguide axis and decay exponentially in the transverse direction. These waves are confined in the core of the waveguide or at the interface. Their existence depends on material contrasts between the core and the surrounding medium [START_REF] Thurston | Elastic waves in rods and clad rods[END_REF]. Radiation modes are standing waves in the transverse direction, which are propagative or evanescent in the axis direction [START_REF] Marcuse | Theory of Dielectric Optical Waveguides[END_REF][START_REF] Malischewsky | Surface Waves and Discontinuities[END_REF]. Finally, leaky modes propagate with leakage attenuation along the waveguide axis. These modes dramatically grow exponentially in the transverse direction [START_REF] Marcuse | Theory of Dielectric Optical Waveguides[END_REF][START_REF] Hladky-Hennion | Conical radiating waves from immersed wedges[END_REF][START_REF] Vogt | The scattering of guided waves in partly embedded cylindrical structures[END_REF].

Modes are computed by considering the source-free problem. This can be done by analytical methods (e.g. the Thomson-Haskell [START_REF] Thomson | Transmission of elastic waves through a stratified solid medium[END_REF][START_REF] Haskell | The dispersion of surface waves on multilayered media[END_REF], the stiffness matrix [START_REF] Kausel | Stiffness matrices for layered soils[END_REF] or the global matrix [START_REF] Pavlakovic | Disperse: A general purpose program for creating dispersion curves[END_REF] methods), which are yet limited to canonical waveguide geometries. Numerical methods are well-suited for more complex geometries. The idea is then to discretize only the cross-section of the waveguide while describing analytically the direction of the wave propagation. In closed waveguides, this approach yields a linear eigenvalue problem. It has been applied under various names in the literature, such as the extended Ritz technique [START_REF] Nelson | Vibrations and waves in laminated orthotropic circular cylinders[END_REF][START_REF] Aalami | Waves in prismatic guides of arbitrary cross section[END_REF] ; the thin layer method (TLM) for stratified waveguides in geophysics [START_REF] Kausel | An explicit solution for the Green functions for dynamic loads in layered media[END_REF][START_REF] Park | Wave motion in finite and infinite media using the thin-layer method[END_REF] ; the strip-element method [START_REF] Liu | Strip element method to analyze wave scattering by cracks in anisotropic laminated plates[END_REF], the Semi-Analytical Finite Element method (SAFE) [START_REF] Marzani | Time-transient analysis response for ultrasonic guided waves propagating in damped cylinders[END_REF][START_REF] Damljanovic | Forced response of a cylindrical waveguide with simulation of the wavenumber extraction problem[END_REF] or more recently the Scaled Boundary Finite Element Method (SBFEM) [START_REF] Gravenkamp | A numerical approach for the computation of dispersion relations for plate structures using the Scaled Boundary Finite Element Method[END_REF][START_REF] Gravenkamp | The computation of dispersion relations for three-dimensional elastic waveguides using the scaled boundary finite element method[END_REF] in ultrasonics. In this paper, this approach will be referred to as waveguide formulation to avoid the use of acronyms.

Extending the waveguide formulation to open waveguides is not straightforward because of the unbounded nature of the problem in the transverse direction. This difficulty is enhanced by the transverse growth of leaky modes. Therefore, the waveguide formulation must be coupled to other techniques to numerically compute the modes of open waveguides.

The first class of methods avoids the discretization of the embedding medium by using appropriate boundary conditions. The waveguide formulation has been combined with the boundary element method to model three-dimensional waveguides immersed in fluids [START_REF] Mazzotti | Dispersion analysis of leaky guided waves in fluid-loaded waveguides of generic shape[END_REF][START_REF] Mazzotti | Ultrasonic leaky guided waves in fluid-coupled generic waveguides: hybrid finite-boundary element dispersion analysis and experimental validation[END_REF] or embedded in solids [START_REF] Mazzotti | A coupled SAFE-2.5D BEM approach for the dispersion analysis of damped leaky guided waves in embedded waveguides of arbitrary cross-section[END_REF]. Similarly, exact boundary conditions have been proposed for two-dimensional waveguides (plates and cylinders) immersed in fluids [START_REF] Gravenkamp | Numerical modeling of elastic waveguides coupled to infinite fluid media using exact boundary conditions[END_REF]. All of these boundary conditions lead to a highly non-linear eigenproblem that is difficult to solve. The latter can be linearized in the case of two-dimensional plates immersed in perfect fluids [START_REF] Hayashi | Calculation of leaky Lamb waves with a semi-analytical finite element method[END_REF][START_REF] Inoue | Transient analysis of leaky Lamb waves with a semi-analytical finite element method[END_REF]. In the case of high-contrast solid waveguides, the waveguide formulation can also be coupled to an approximate condition (the so-called dashpot boundary condition) [START_REF] Gravenkamp | Computation of dispersion curves for embedded waveguides using a dashpot boundary condition[END_REF]. With this approximation, the eigenproblem remains linear.

The second class of methods requires a discretization of the surrounding medium, which must be truncated. The eigenproblem remains linear. To avoid spurious reflections due to truncation, the waveguide formulation have been combined with non-reflecting [START_REF] Hladky-Hennion | Finite element modelling of radiating waves in immersed wedges[END_REF] and continued-fraction absorbing [START_REF] Astaneh | Dispersion analysis of composite acousto-elastic waveguides[END_REF] boundary conditions in fluids, or paraxial approximation in solids [START_REF] Park | Wave motion in finite and infinite media using the thin-layer method[END_REF][START_REF] Kausel | Physical interpretation and stability of paraxial boundary conditions[END_REF].Absorbing layers of artificially growing damping can also be used to simulate fluid [START_REF] Fan | Torsional waves propagation along a waveguide of arbitrary cross section immersed in a perfect fluid[END_REF] or solid [START_REF] Castaings | Finite element model for waves guided along solid systems of arbitrary section coupled to infinite solid media[END_REF] infinite media. Another technique consists in using a Perfectly Matched Layer (PML) to model the infinite surrounding medium (solid or fluid) [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF][START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF][START_REF] Zuo | Numerical modeling of embedded solid waveguides using SAFE-PML approach using a commercially available finite element package[END_REF][START_REF] Zuo | SAFE-PML approach for modal study of waveguides with arbitrary cross sections immersed in inviscid fluid[END_REF][START_REF] De Oliveira Barbosa | Perfectly matched layers in the thin layer method[END_REF].

Contrary to absorbing layers, the PML avoids most spurious reflections from the layer, which allows its thickness to be greatly reduced. Moreover, it have been shown that the computation of leaky modes with a PML is mathematically relevant (see Ref. [START_REF] Kim | The computation of resonances in open systems using a perfectly matched layer[END_REF] for scalar wave problems).

As far as the forced response problem is concerned, modal expansion methods have been widely applied in closed waveguides [START_REF] Auld | Acoustic fields and waves in solids[END_REF]. However, their application to open elastic waveguides is more intricate and has been barely considered in the literature. With a numerical approach, the particular case of two-dimensional plates immersed in fluids has recently been handled by a waveguide formulation with exact boundary conditions [START_REF] Inoue | Transient analysis of leaky Lamb waves with a semi-analytical finite element method[END_REF]. The case of a stratified plate over or between half-spaces has been treated using a PML in Ref. [START_REF] De Oliveira Barbosa | Perfectly matched layers in the thin layer method[END_REF].

Theoretically, the forced response of an open waveguide can be expanded on trapped modes and radiation modes [START_REF] Marcuse | Theory of Dielectric Optical Waveguides[END_REF][START_REF] Malischewsky | Surface Waves and Discontinuities[END_REF][START_REF] Collin | Field theory of guided waves[END_REF], such that the displacement field can be symbolically written as:

u(r, ω) = trapped + radiation modes (1)
where, in addition to trapped modes, complex poles of backward type are also likely to occur depending on the problem type [START_REF] Tamir | Guided complex waves. Part 1: Fields at an interface[END_REF][START_REF] Nedospasov | Unusual energy properties of leaky backward Lamb waves in a submerged plate[END_REF]. Let us briefly recall the origin of Eq. ( 1). In the wavenumber domain, the solution of the problem is a multivalued function owing its dependence on the transverse wavenumber of the unbounded medium. The transverse wavenumber is indeed the square root of a complex number on a twosheeted Riemann surface. To evaluate analytically the inverse spatial Fourier transform of the solution, a branch cut is defined separating the proper Riemann sheet (where trapped modes occur) from the improper Riemann sheet (where leaky modes occur). Hence, the inverse transform integration is performed only on the proper sheet and gives rise [see Eq. ( 1)] to the discrete sum of trapped modes and to the continuum of radiation modes, which represents the branch cut contribution. This continuum is characteristic of the unboundedness of the modal problem and is difficult to manipulate from a mathematical point of view. In elastodynamics, there are two continua instead of one because two transverse wavenumbers occur (longitudinal and shear waves) [START_REF] Maupin | The radiation modes of a vertically varying half-space: a new representation of the complete Green's function in terms of modes[END_REF][START_REF] Margerin | Generalized eigenfunctions of layered elastic media and application to diffuse fields[END_REF]. The continua can actually be approximated with a convenient discrete set of leaky modes, e.g. using the steepest descent method [START_REF] Collin | Field theory of guided waves[END_REF][START_REF] Tamir | Guided complex waves. Part 1: Fields at an interface[END_REF]. This approximation is valid in a zone restricted near the core, in which leaky modes can provide useful practical information such as axial attenuation and travelling velocity of waves packets [START_REF] Pavlakovic | High-frequency low-loss ultrasonic modes in imbedded bars[END_REF].

In a recent work [START_REF] Gallezot | Contribution of leaky modes in the modal analysis of unbounded problems with perfectly matched layers[END_REF], the authors have shown that when the surrounding medium is modified by an infinite PML, the forced response can theoretically be obtained with a modal expansion on trapped modes, revealed leaky modes and two continua of radiation modes related to the PML branch cuts. With numerical methods, the PML must yet be truncated to a finite thickness. The continua are then transformed into discrete sets of PML modes. It remains necessary to investigate the use of the so-obtained numerical modes to compute the forced response, based on an appropriate orthogonality relationship to guarantee the uniqueness of the solution. It is of particular interest to derive a general orthogonality relationship that remains applicable for a wide class of problems, including viscoelastic materials, full anisotropy and 3D waveguides of arbitrary cross-section. Moreover, it is necessary to clarify the contribution and the physical relevance of leaky and PML modes in the solution of the problem.

The aim of this article is to compute the forced response of elastic open waveguides with a modal approach based on a waveguide formulation combined with a PML. This approach will be termed PML waveguide formulation in this paper. It is briefly recalled in Sec. 1.2. Then, the modal approach used to compute the forced response is presented in Sec. the ability of the approach to handle three-dimensional waveguide problems.

PML waveguide formulation for mode computation

This section briefly recalls the PML waveguide formulation to compute the eigenmodes of an open waveguide [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]. One considers an elastic open waveguide, where z is the propagation axis of the waveguide.

The core of the waveguide has an arbitrary cross-section. It is embedded in an infinite medium. Small strains and displacements are assumed. A time-harmonic dependence e -jωt is chosen. The elastodynamic equilibrium equations are written as:

∇.σ(ũ) + ρω 2 ũ = f , ( 2
)
where ũ is the displacement vector field. The tilde notation is explained in the following with the PML introduction. σ(ũ) is the stress tensor. The stress-strain relation is σ(ũ) = C : ǫ(ũ), where C is the stiffness tensor (complex for viscoelastic materials). The strain-displacement relation is ǫ(ũ) = (∇ũ + (∇ũ) T )/2.

The superscript T denotes the transpose. f is the vector of volumic forces.

For numerical purpose, the infinite surrounding medium is truncated by a finite PML of thickness h. A radial PML is introduced by analytic continuation [START_REF] Chew | A 3D perfectly matched medium from modified maxwells equations with stretched coordinates[END_REF] of the elastodynamic equations (2) into the complex transverse coordinate r:

r(r) = r 0 γ(ξ)dξ. ( 3
)
Note that a Cartesian PML could be used instead of a radial PML [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF], as done in Sec. 4.5. The function γ(r) is a user-defined complex-valued function for absorbing outgoing waves in the surrounding medium, such that:

• γ(r) = 1 outside the PML region (r < d),

• Im(γ(r)) > 0 inside the PML region (d < r < d + h).

d is the position of the PML interface. A boundary condition is arbitrarily applied at the end of the PML. In this paper, a Dirichlet condition is used. Eq. ( 2) must be transformed to go back to the real radial direction r. The change of variable r → r for any function g(r) gives:

g(r) = g(r), dr = γ(r)dr, ∂g ∂r = ∂g ∂r 1 γ(r) . ( 4
)
Following the waveguide formulation in closed waveguides [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF], the strain-displacement relation can be written as:

ǫ = (L S + L z ∂ ∂z )u, ( 5
)
where ǫ is the strain vector and u the displacement vector. L S is the operator containing all terms but derivatives with respect to the z-axis, and L z is the operator of z-derivatives. The operator L S includes derivatives with respect to r and hence it depends on γ(r) and r [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]. A FE discretization is then applied to the cross-section (spectral elements can be equally used [START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF]). The displacement on each element can be written u e (x, y, z, ω) = N e (x, y)U e (z, ω), where N e (x, y) is the matrix of two-dimensional interpolation functions and U e (z, ω) is the vector of nodal displacements. The derivation of the weak formulation of Eq. ( 2) yields the global matrix system:

(K 1 -ω 2 M)U + (K 2 -K T 2 )U ,z -K 3 U ,zz = F. (6)
with the element matrices:

K e 1 = L T S CL z rγ r dS, K e 2 = L T S CL z rγ r dS, K e 3 = L T z CL z rγ r dS, M e = ρN eT N e rγ r dS, F e = N eT f rγ r dS. ( 7
)
where C is the stiffness matrix defined from Voigt notation of the tensor C.

A spatial Fourier transform in the z-direction is applied to Eq. ( 6), which leads to the PML waveguide formulation [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]:

K 1 -ω 2 M + jk(K 2 -K T 2 ) + k 2 K 3 Û = F, (8)
where Û(k, ω) is the vector of nodal displacements, F is the vector of nodal forces and k is the axial wavenumber with the following notation for any function g(z):

ĝ(k) = +∞ -∞ g(z)e -jkz dz. ( 9
)
Further details of the derivation can be found in Refs. [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]. For paper self-containedness, additional steps are given in Appendix A for an axisymmetric open waveguide, which will be used for most of the numerical results presented in Secs. 3 and 4.

The source term in Eq. ( 8) is dropped (i.e. F = 0). The source-free problem gives a quadratic eigenvalue problem in terms of k. The eigenproblem is linearized (see Appendix B) and solved for each angular frequency ω. The matrices K 1 , K 3 and M are symmetric. From this property, it can be shown that the eigenspectrum is made of symmetric pairs of eigenvalues and eigenvectors, (k m , U m ) and (-k m , U - m ), (m = 1, . . . , M ), representing M positive-going modes and M negative-going modes respectively. In the following, the negative-going modes will be denoted (k -m , U -m ).

As recalled in Sec. 1.1, the PML waveguide formulation gives three kinds of eigenmodes: trapped modes (if any), leaky modes and the so-called PML modes. PML parameters have a strong influence on the accuracy of trapped and leaky modes. Theoretically (see e.g. Refs. [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]), leaky modes are attenuated by the PML if:

arg(k l/s ) > -arg(γ), (10)
where k l/s is the longitudinal or shear transverse wavenumber and γ is the average value of γ(r) inside the PML region. When the condition ( 10) is fulfilled, only small reflections occur at the end of a thick enough PML. As usual with PML, the attenuation rate should not be too high to be properly approximated by the FE discretization. Concerning trapped modes, the PML can accentuate their natural transverse decay if Re(γ) > 1 [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]. As for PML modes, these modes mainly resonate inside the PML and fully depend on PML parameters.

Forced response of an elastic open waveguide

Any arbitrary field in a waveguide can be expanded on a unique set of modes, provided that these modes are orthogonal [START_REF] Auld | Acoustic fields and waves in solids[END_REF]. In elastic open waveguides, orthogonality relationships theoretically hold between trapped modes and radiation modes, see e.g. Ref. [START_REF] Malischewsky | Surface Waves and Discontinuities[END_REF][START_REF] Maupin | The radiation modes of a vertically varying half-space: a new representation of the complete Green's function in terms of modes[END_REF]. However, the continua of radiation modes are difficult to manipulate. Conversely, leaky modes form a discrete set, which may approximate the radiation field [START_REF] Sammut | Leaky modes on a dielectric waveguide: orthogonality and excitation[END_REF][START_REF] Lee | On leaky mode approximations for modal expansion in multilayer open waveguides[END_REF][START_REF] Hu | Understanding leaky modes: slab waveguide revisited[END_REF] and which is easier to handle. However, they cannot satisfy orthogonality relationships [START_REF] Vogt | The scattering of guided waves in partly embedded cylindrical structures[END_REF][START_REF] Sammut | Leaky modes on a dielectric waveguide: orthogonality and excitation[END_REF], as they theoretically grow to infinity in the transverse direction (as recalled in Eq. ( 1), leaky modes do not belong to the theoretical modal basis). In this section, a general orthogonality relationship is derived when a finite PML is introduced to close the open problem.

Effects of a PML on the modal basis

In a recent work, the authors have studied the modal basis of an elastic open waveguide when the surrounding medium is replaced by an infinite PML [START_REF] Gallezot | Contribution of leaky modes in the modal analysis of unbounded problems with perfectly matched layers[END_REF]. It has been shown that the introduction of an infinite PML of constant attenuation profile γ modifies the transverse wavenumbers and yields a rotation of branch cuts by an angle of -arg(γ) [START_REF] Gallezot | Contribution of leaky modes in the modal analysis of unbounded problems with perfectly matched layers[END_REF]. The rotation of the branch cuts induces a partial access to the initial improper sheets where leaky modes occur. Therefore with a PML, a discrete set of revealed leaky modes must be included in the modal basis. The set of trapped modes is left unchanged and the solution field with an infinite PML can hence be written:

u(r, ω) = trapped modes + revealed leaky modes + PML branch cuts. (11)
The solution field described by Eq. ( 11) remains physically valid only in the region out of the PML (r < d).

With numerical methods, the PML must be truncated to a finite thickness. The continuous sums related to the PML branch cuts are then discretized [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF][START_REF] Olyslager | Discretization of continuous spectra based on perfectly matched layers[END_REF]. These discretized sets are PML modes, non-intrinsic to the physics since they fully depend on the parameters of the truncated PML and on the boundary condition applied at the end of the PML. Therefore, the role of these sets in the convergence of the modal expansion might be questionable at first sight. This topic will be addressed in Sec. 3 through numerical results.

Derivation of a biorthogonality relationship

Let us start from the linear eigenvalue problem (see Appendix B), written for a mode m:

(A -k m B)x m = 0, (12)
with x T m = [U T m k m U T m ] T .
A is not symmetric, and then left eigenvectors must be considered [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. The left eigenvector y n satisfies the eigensystem:

y T n (A -k n B) = 0. (13)
The left and right eigenvectors are related by the following general biorthogonality relationships:

y T n Bx m = b m δ mn , y T n Ax m = k m b m δ mn , ( 14
)
where b m is a normalization factor and δ mn is the Kronecker' symbol. B is symmetric, and from Eq. ( 13), the left eigenvector y m satisfies:

(A T -k m B)y m = 0. (15)
It can be checked that the above eigensystem is satisfied with:

y T m = [U T -m k m U T -m ] T . ( 16
)
The left eigenvectors can then be directly related to the right eigenvectors. This avoids to solve the left eigenproblem. Hence, from Eq. (B.2) the biorthogonality relationship y T n Bx m = b m δ mn can be expanded as:

U T -n (K 1 -ω 2 M)U m + k n k m U T -n K 3 U m = b m δ mn . ( 17
)
From Eq. (B.3) written for a given mode m, one gets

(K 1 -ω 2 M)U m = -jk m (K 2 -K T 2 )U m -k 2 m K 3 U m .
Eq. ( 17) can then be rewritten as:

jk m U T -n (K T 2 + jk m K 3 )U m -jk m U T -n (K 2 -jk n K 3 )U m = b m δ mn . ( 18
)
In Eq. ( 18), one can identify the so-called modal forces [START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF]:

F m = (K T 2 + jk m K 3 )U m , (19)
such that after rearrarengement, Eq. ( 18) yields:

jω 4 (U T m F -n -U T -n F m ) = Q m,-m δ mn , ( 20
)
with Q m,-m = -ωb m /4k m . Eq. ( 20) provides a general orthogonality relationship involving the whole cross-section of the waveguide (i.e including the finite PML). It is applicable to any kind of modes without distinction (trapped, leaky and PML modes).

This orthogonality relationship is identical to the orthogonality relationship derived for lossy closed waveguides in [START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF]. It remains valid for 3D waveguides of arbitrary cross-section, as well as for fully anisotropic and viscoelastic materials (it is noteworthy that Eq. ( 20) is also applicable with absorbing viscoelastic layers, see e.g. results in Ref. [START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF]). As shown in this reference, Eq. ( 20) can be readily expressed in an analytical form:

Q m,-n = jω 4 S (u m .t -n -u -n .t m )dS = 0 if m = n, (21)
where for each mode, t m is the traction vector applied on the section S. Actually, Eq. ( 21) corresponds to Auld's real orthogonality relationship [START_REF] Auld | Acoustic fields and waves in solids[END_REF], which could therefore be used for the analytical modeling of open waveguides truncated by a finite PML. Note that except for real wavenumber modes (perfectly guided modes or trapped modes), Q m,-m differs from the averaged power flow P m [START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF]. The averaged power flow can yet be computed from the Poynting vector of each mode from Eq. ( 29).

As a side remark, a mathematically exact orthogonality relationship has been proposed in Ref. [START_REF] Sammut | Leaky modes on a dielectric waveguide: orthogonality and excitation[END_REF] for leaky modes in dielectric waveguides by deforming the radial coordinate into the complex plane. The approach of the present paper, applied to elastodynamics, is quite similar.

Modal solution of the problem

Let us consider the excitation of the waveguide by a force. The linear formulation [see Eq. (B.1)] is premultiplied by y T n and the vector x is decomposed on 2M modes as follows:

x = M -M α m x m (m = 0).
The system is then expanded and, thanks to the biorthogonality relation [START_REF] Mazzotti | Ultrasonic leaky guided waves in fluid-coupled generic waveguides: hybrid finite-boundary element dispersion analysis and experimental validation[END_REF], it can be checked that the modal coefficients are given by

α m = -ωU T -m F/[4Q m,-m (k m -k)].
Hence, the displacement solution to the forced response problem is:

Û = - +M m=-M m =0 ωU T -m F 4Q m,-m (k m -k) U m . ( 22
)
Then, the inverse Fourier transform U(z, ω) = 1/(2π)

+∞ ∞ Û(k, ω)e jkz dk is used to go back to the real space (z, ω). With the help of Cauchy residue theorem, it can be shown that the field outside the source region going towards the positive z direction is [START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF]:

U(z, ω) = M m=1 E m F(k m )e jkmz (23)
where E m is the generalized modal excitability matrix:

E m = jω 4Q m,-m U m U T -m . ( 24
)
(E m ) ij is the displacement amplitude of mode m at degree of freedom (dof) i when excited by a unit force at dof j. If one is interested in the field going towards the negative z direction, then the summation in Eq. ( 23) must be performed on negative-going modes. The +j factor in Eq. ( 24) is replaced by -j and the subscripts m and -m are switched.

Finally, the time-transient displacement is obtained by the inverse Fourier transform 1/(2π) +∞ ∞ U(z, ω)e -jωt dω.

Determination of the travelling direction

The solution given by Eq. ( 23) requires to determine the travelling direction of each mode. Modal characteristics can be computed straightforwardly from the eigensolutions (k m , U m ) and the finite element matrices [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]. The phase velocity of a mode m is:

v pm = ω Re(k m ) . ( 25
)
Its attenuation in dB/m is:

η m = 8, 686Im(k m ). ( 26
)
The group velocity v gm = ∂ω/∂k m is computed as follows:

v gm = Re 2ωU T -m MU m U T -m (j(K 2 -K T 2 ) + 2k m K 3 )U m -1 . ( 27
)
The phase velocity does not necessarily give the travelling direction owing to the existence of backward modes, which have phase and group velocites of opposite sign [START_REF] Auld | Acoustic fields and waves in solids[END_REF]. The direction of propagation of the modes can be properly determined according to the following criterion:

• the sign of Im(k m ) if Im(k m ) = 0, such that if Im(k m ) > 0,
the mode is a positive-going mode, and if Im(k m ) < 0 the mode is negative-going,

• the sign of v gm if k m ∈ R (that is for trapped modes, when they exist).

The case of trapped modes is somehow particular. The case of purely trapped modes never occurs in viscoelastic media since the wavenumbers always have a non zero imaginary part. In purely elastic waveguides, trapped waves exist only if the shear velocity is greater in the core, unless longitudinal and shear waves couple into interface trapped waves, called Stoneley waves [START_REF] Thurston | Elastic waves in rods and clad rods[END_REF].

Remarks on the energy velocity

In lossy waveguides, it has been shown that the group velocity can yield non-physical values [START_REF] Bernard | Guided waves energy velocity in absorbing and non-absorbing plates[END_REF], so that the energy velocity is usually preferred in closed waveguides. However, the energy velocity must be averaged on the whole cross-section, which is theoretically unbounded considering open waveguides.

For trapped modes, which are non-attenuated waves that exponentially decay in the transverse direction, the group velocity is equal to the energy velocity. Hence, the criterion proposed in Section 2.4 to determine their travelling direction is consistent (it does not matter whether the group velocity or the energy velocity is used).

For leaky modes, the energy velocity cannot be theoretically defined because these modes grow to infinity in the transverse directions. In practice, the energy velocity integration of leaky modes is usually restricted to the core, which can give a reasonable approximation of the velocity of wave packets [START_REF] Pavlakovic | High-frequency low-loss ultrasonic modes in imbedded bars[END_REF][START_REF] Bernard | Guided waves energy velocity in absorbing and non-absorbing plates[END_REF]. The energy velocity restricted to the core cross-section can be calculated as follows:

v em = Re(P m ) Re(T m ) + Re(V m ) , ( 28
)
where for each mode m, P m is the normal component of the Poynting vector integrated on the core crosssection, T m is the core cross-section and time-averaged kinetic energy, and V m is the core cross-section and time-averaged potential energy. These quantities can be post-processed from [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]:

P m = -jω 2 U * m (K T 2 + jk m K 3 )U m , (29)
T m = ω 2 4 U * m MU m , (30)
V m = 1 4 U * m (K 1 + jk m K 2 -jk * m K T 2 + k * m k m K 3 )U m . ( 31
)
where the * superscript denotes conjugate transpose and the overbar stands for the dof restriction to the core cross-section.

It is noteworthy that the integration restriction to the core is quite arbitrary. Therefore, the so-defined energy velocity cannot be considered as a general indicator to determine the travelling direction of any type of modes (of PML modes, in particular). In the following of this paper, the above-defined energy velocity will be only used for the post-processing of dispersion curves of leaky modes and their comparison with literature results.

The homogeneous test case

In this section, the forced response of a homogenous medium excited by a point load is computed by the numerical method and compared to an analytical solution. A homogeneous medium can be viewed as an open waveguide without a core. This test case is actually the worst scenario for the numerical modal approach since there is no contrast of impedance: only the discrete set of PML modes occurs (there is neither trapped nor leaky modes). As outlined in Section 2.1, the contribution of PML modes to the solution has to be clarified because these modes are not intrinsic to the physics (they mainly oscillate in the PML region, which mainly depends on user-defined parameters).

Description

Let us assume a homogeneous, isotropic and unbounded medium, made of cement grout material (see Fig. 1). Material properties are given in Table 1. In this section, the medium is pure elastic (η l = η s = 0).

The medium is excited by a point force in the z-direction at the origin: The transient axial displacement along the z-axis (i.e. at r = 0) at a distance z to the source is given by the following analytical solution (see e.g. Ref. [54]):

f = F (t)δ(x)e z (f = F (t) 2πr δ(r)δ(z)
u ref z (z, t) = 1 2πρz 3 z/cs z/c l τ F (t -τ )dτ + 1 4πρc 2 l z F (t - z c l ). (32)
The first term of Eq. ( 32) is a near-field contribution. The second term is the far-field contribution of the P-wave (along the z-axis, there is no far-field contribution of the S-wave [54]).

The problem is solved with an axisymmetrical PML waveguide formulation (see Appendix A). The PML starts at r = d = 10mm. Following Refs. [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF], the attenuation function of the PML is chosen as parabolic and is given by:

γ(r) =      1 if r ≤ d 1 + 3(γ -1) r -d h 2 if r > d, (33)
where γ = 1

e(z) = |u ref z (z, ω) -u num z (z, ω)| 2 dω |u ref z (z, ω)| 2 dω . ( 34
)
where u ref z (z, ω) is obtained from the analytic Fourier transform of Eq. ( 32), given in Appendix C (from Parseval's equality, the relative error has the same value both in frequency and time domains).

Results

Figure 2 represents the discrete modal spectrum in the complex wavenumber kd-plane (Re(kd) > 0, Im(kd) > 0) at the centre frequency, computed for γ = 4 + 4j and the following three cases: (h = 4d, M = 50), (h = 4d, M = 30) and (h = 3d, M = 50). This enables to study the influence of two key features on the relative error of the forced response:

• the complex thickness d + γh,

• the number of modes M retained in the modal expansion.

The number of dofs of the quadratic eigenvalue problem is equal to 68 dofs for h = 3d and 84 dofs for h = 4d. These values must be multiplied by two for the linearized eigenvalue problem (see Appendix B).

As observed in Fig. 2, the poles close to the real axis (Im(kd) = 0) lie along two hyperbolas, corresponding to the branch cuts rotated by the PML. As briefly recalled in Sec. 2.1, these branch cuts are rotated by an angle equal to -arg(γ) in the case of an infinite PML (for further details, see Ref. [START_REF] Gallezot | Contribution of leaky modes in the modal analysis of unbounded problems with perfectly matched layers[END_REF]). For a PML of finite thickness, this angle becomes equal to -arg(d + hγ) as shown in Ref. [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF]. Note that the deviation from the theoretical hyperbolas observed here for high order poles is explained by the FE approximation: this deviation can be moved further away from the real axis by reducing the size of elements (see Ref. [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF]).

The relative error of the forced response is plotted in Fig. 3 as a function of the distance z. It can be observed that a good accuracy is reached for the three cases in the range z ∈ [20d, 90d], with a relative error lower than 0.5%. This shows the ability of PML modes, when the PML is properly parametrized, to reconstruct the analytical solution. In particular, the geometrical decay of the field is well recovered by the superposition of PML modes, which are by nature exponentially decaying.

Further away from the source (z > 90d), it turns out that the error gets lower by increasing the complex thickness (compare red and blue results with the black ones). As shown in previous studies [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF], the discretization of the continua gets denser as the complex thickness increases. This effect can be observed in Fig. 2 (see e.g. the distance labelled A and B between two poles for two different complex thicknesses).

Indeed, only modes with low |Im(kd)| can significantly contribute to the forced response at far distances. It is then necessary to have a denser spectrum to reassemble a geometrically decaying solution.

Yet in the near field (z < 20d), even highly attenuated PML modes (with high |Im(kd)|) can significantly contribute to the solution. For a given modal density (i.e. a given complex thickness), the error can then be reduced by increasing the maximum imaginary part of the computed spectrum, that is to say, the number of modes M retained in the modal expansion (compare blue with red results). Conversely, for a fixed number of modes M , the error in the near field could also be reduced by decreasing the complex thickness.

It has to be mentioned that the numerical results of Fig. 3, obtained in this paper for the elastodynamic problem, are consistent with the theoretical proof of convergence of the PML mode expansion derived in Ref. [START_REF] Olyslager | Discretization of continuous spectra based on perfectly matched layers[END_REF] for the Helmholtz equation. As an example, Fig. 4 compares the transient numerical and analytical solutions calculated at the distance z = 175d. Good agreement is found, which confirms the ability of the numerical approach to accurately compute the forced response even in the fully homogeneous configuration (worst scenario). Note that the distance z = 175d will be of interest when considering a core waveguide example as discussed later in Sec. 4.

The argument of the complex thickness, denoted as θ (θ ∈]0, π/2[), also plays a role on the accuracy of the numerical results.

Application

A steel core is now included in the homogeneous medium previously studied (see Fig. 6). First, the core has a circular cross-section. This configuration enables to highlight the contribution of leaky modes, compared to that of PML modes, as a function of the distance to the source. As before, this open waveguide is modelled with the axisymmetrical PML waveguide formulation. Then, a square cross-section is considered as an example of three-dimensional waveguide.

Circular bar test case description

The material properties are given in Table 1. The shear wave velocity in the core (steel) is greater than in the embedding medium (cement grout). For this configuration, there is no trapped modes [START_REF] Thurston | Elastic waves in rods and clad rods[END_REF][START_REF] Pavlakovic | High-frequency low-loss ultrasonic modes in imbedded bars[END_REF], so that only leaky modes and PML modes occur. The core has a radius a = 10 mm. Both the elastic and viscoelastic cases will be considered. The embedding medium is truncated by a PML of thickness h starting at a distance r = d with the same attenuation profile as Eq. [START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF]. In order to avoid spurious eigenvalues and limit the transverse growth of leaky modes, it is preferable to set the PML close to the core [START_REF] Kim | The computation of resonances in open systems using a perfectly matched layer[END_REF] (d -a should be small). The PML interface has been set to d = a.

Dispersion curves

For the visualisation of dispersion curves, only the modes intrinsic to the physics are of interest (leaky and trapped modes). It is then necessary to filter out PML modes. A filtering criterion has been proposed in previous studies [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF] to identify and remove PML modes. A slightly modified version is presented here. The proposed filtering criterion is based on the ratio of the imaginary part over the modulus of the kinetic energy, defined from Eq. ( 30) but integrating the whole cross-section (i.e. including the PML, so that the overbar has to be discarded in Eq. ( 30)). The modes are retained in the visualisation if they fulfil the following criterion:

1 - Im(T m ) |T m | > η min (35)
where 0 < η min < 1 is a user-defined parameter.

The filtered dispersion curves obtained for the viscoelastic case are shown in Fig. 7. The figure shows 12 modes, of leaky type, corresponding to longitudinal modes (these modes are usually labelled L(0,n) in the literature). Also shown are the results of Ref. [START_REF] Pavlakovic | High-frequency low-loss ultrasonic modes in imbedded bars[END_REF], which confirms the ability of the PML waveguide formulation to accurately compute leaky modes [START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF][START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF][START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF].

Excitation parameters

The forced response of the open waveguide is computed with the same type of force as in Sec. 3, with an amplitude of 1N, considering the following two signals:

• f c = 60 kHz with a 5 cycles Hanning window. At this frequency, only one leaky mode propagates (see point A in Fig. 7). This low-frequency excitation allows a direct comparison with the results of Sec. 3 and aims to clarify the influence of leaky modes on the solution compared to that of PML modes.

The problem is solved on the frequency range [f min , f max ], with f min = 12 kHz and f max = 120 kHz, divided into 323 frequency steps. The size of finite elements is set to ∆r = 0.25a, yielding 84 dofs for the quadratic eigenvalue problem (84x2 for the linearized one).

• f c = 2.284 MHz with a 10 cycles Hanning window. This high-frequency excitation involves several leaky modes and is of interest for NDE because the mode of lowest attenation is excited (see point B in Fig. 7). For this high-frequency excitation, the problem is solved for [f min , f max ], with f min = 1.14

MHz and f max = 3.42 MHz, divided into 1823 steps. The size of finite elements is set to ∆r = 0.019a, yielding 240 dofs for the quadratic eigenvalue problem (240x2 for the linearized one).

Low frequency response

For the low-frequency excitation (f c = 60 kHz), the PML parameters are set according to the results of Sec. 3 in order to ensure a good accuracy of the PML mode contribution (h = 4a, γ = 4 + 4j). As in Sec. 3,

0.2 0.4 0.6 0.8 1 1.2 t (ms) -3 -2 -1 0 1 2 3 U z (mm)

×10 -12

Leaky mode 1 PML modes contribution the axial component u z of the forced response is computed along the z-axis. 50 PML modes and 1 leaky mode are now included in the modal superposition.

Figure 8 shows the transient response at z/a = 175 both for the elastic and the viscoelastic cases. As expected, the amplitude of the response is reduced in the viscoelastic case (but the times of flight are nearly unchanged). Three wave packets can be observed. The first one is the contribution of the leaky mode.

The second and third packets are the contribution of PML modes (it has been checked that these packets disappear if PML modes are not retained in the modal expansion). Interestingly, the second and third packets can be compared to the results of Fig. 4 obtained for the homogeneous test case (i.e. without steel core): it can be noticed that the times of flight coincide each other. Therefore, the contribution of PML modes observed in Fig. 8, inside the steel core, can be interpreted as the trace of the field radiated into the surrounding medium.

To confirm this interpretation, Fig. 9 plots the quantity proportional to the energy, |u z (r = 0, z)| 2 , as a function of the distance z computed at the centre frequency (elastic case). Up to a certain distance (z/a = 170), the leaky mode contribution prevails: the energy decreases exponentially. But further away from the source, the leaky mode contribution becomes negligible and the PML mode contribution prevails:

the energy decreases with a geometrical decay. In between, there is a transition zone (around z/a = 175)

where both types of modes contribute nearly equally as can be observed in Fig. 8.

Similar results have been obtained analytically in the case of optical scalar waveguides [START_REF] Hu | Understanding leaky modes: slab waveguide revisited[END_REF], showing that radiation modes first reassemble the leaky mode contribution (exponentially decaying) and then yields the long-term diffraction (geometrically decaying). Provided that the PML parameters are set properly, Fig. 8 and 9 show that the PML technique enables to compute the long-term diffraction phenomenon. This is not the case with absorbing viscoelastic layers [START_REF] Hu | Understanding leaky modes: slab waveguide revisited[END_REF].

High frequency response

The NDE of open waveguides requires that enough energy remains inside the core. At high frequencies, the energy can be concentrated into the core and the leakage attenuation reduced [START_REF] Pavlakovic | High-frequency low-loss ultrasonic modes in imbedded bars[END_REF]. In Fig. 7, the lowest attenuation is reached for the 12th leaky mode, with a characteristic series of lobes above 20 MHz.mm (the minimum attenuation is reached for the second lobe, at point B, and equals 159 dB.mm/m). At higher frequencies, the attenuation of higher order modes increases because of greater losses due to viscoelasticity than to leakage [START_REF] Pavlakovic | High-frequency low-loss ultrasonic modes in imbedded bars[END_REF]. Since it can propagate over greater distances than the others, the 12th mode is hence quite attractive.

In the context of NDE, only leaky modes are of interest. The contribution of PML modes (long-term diffraction, geometrically decaying) does not need to be considered because it occurs when a lot of energy has been already lost from the core (that is to say, when the problem can somehow no longer be considered as a waveguide problem). The computational cost can then be reduced here because, for the high-frequency excitation (f c = 2.284 MHz), an acceptable convergence of leaky modes can actually be reached with a thin PML. After numerical tests, the PML thickness has been reduced to h = 0.1a with γ = 1 + 2j.

Figure 10b shows the transient solution at z/a = 175. M = 20 leaky modes have been used to compute the forced response. Thanks to the modal approach used in this paper, the contribution of each mode m to the total forced response at a given distance of propagation can be computed. This modal contribution is given by the vector E m F(k m )e jkmz (see Sec. 2.3). The so-obtained values are plotted in Fig. 10a for the dof corresponding to the axial displacement at r = 0, over the energy velocity curves. It shows that the first wave packet is mainly due to the 12th mode with a smaller contribution of the 11th mode. The wave packet of the second arrival is composed of weaker modal contributions, corresponding to parts of dispersion curves having higher attenuation and slower velocity. As expected, the energy is essentially contained in the first arrival because the centre frequency of the excitation coincides with the maximum energy velocity and the minimum attenuation of the 12th mode (point B in Fig. 7).

Excitability of leaky modes

In the analysis of elastic waveguides, the excitability of a given mode is an important feature, defined by the displacement-force ratio, which typically enables to optimise the type and location of sensors to be used [START_REF] Wilcox | Mode and transducer selection for long range lamb wave inspection[END_REF]. The expression of excitability arises from the orthogonality of modes. Yet from a theoretical point of view and as mentioned earlier, leaky modes cannot satisfy an orthogonality relationship because the wave fields of such modes do not vanish at infinity in the transverse direction. To circumvent this problem, some authors have proposed to restrict the integration over the core [START_REF] Vogt | The scattering of guided waves in partly embedded cylindrical structures[END_REF].

However, the PML technique used in this paper enables to calculate the excitability, as defined from Eq. ( 24), by integrating over the whole cross-section including the PML. Like the biorthogonality relation [START_REF] Mazzotti | Ultrasonic leaky guided waves in fluid-coupled generic waveguides: hybrid finite-boundary element dispersion analysis and experimental validation[END_REF], the so-defined excitability is applicable to any type of modes, and in particular to leaky modes. The excitability is found to be strong for odd modes (modes 1, 3, 5, 7, 9, 11) and weak for even modes up to the 12th mode. As previously mentioned, the 12th mode is somehow specific as it is composed of a series of lobes. These lobes also exhibit a strong variation of excitability (contrary to the energy velocity and the attenuation). For instance, the excitability of the second lobe (point B) is much higher than the excitability of the first lobe (point C), by an amount of 20 log 10 |E B /E C | = 42 dB.

These differences of excitability can be explained by inspecting the mode shapes. As an example, Fig. 12 shows the mode shapes of the first two lobes of the 12th mode. The modal displacement of the first lobe is significantly lower at the centre of the core (i.e. where the point force is applied), hence reducing the excitability at this point.

A square bar example

A three-dimensional waveguide with a core cross-section of square shape is now considered. As previously mentioned, the theory of Sec. 2 remains fully applicable. However, the computational cost significantly increases because the cross-section is now two-dimensional (instead of one-dimensional in the previous test cases). Furthermore, the eigenspectrum contains flexural and torsional leaky modes (in addition to longitudinal modes) and much more PML modes as well. In the following, the forced response of the square bar is computed in a high-frequency range and compared to the circular case.

Square bar test case description

The dispersion curves of a square bar open waveguide have already been obtained with a Cartesian PML waveguide formulation discretized with spectral elements in Ref. [START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF]. The square bar is made of viscoelastic steel and is embedded into viscoelastic cement grout. Materials are given in Table 1. The half-width of the square is given by a = 10mm. The Cartesian PML is introduced by an analytic continuation of the elastodynamic equations (2) into the complex plane (x, ỹ). As in the circular bar case, its attenuation is chosen parabolic. The PML parameters are then given in both directions x and y by the mean attenuations γx and γy , the thickness h x and h y and the distances d x and d y . In the following, d x = d y = a. Details about the implementation of a Cartesian PML can be found in Ref. [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF] and are not recalled in this paper for conciseness. As for spectral elements discretization, the specific element matrices and details about the convergence of the method can be found in Ref. [START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF]. To reduce the number of computed modes, the eigensolver shift is set to k 0 = ω/c l0 , where c l0 is the longitudinal bulk velocity in the core. This enables to compute mainly longitudinal modes, while avoiding the computation of PML modes or highly attenuated leaky modes of minor interest [START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF]. Note that opposite-going modes also have to be computed (using a -k 0 shift) to obtain the forced response.

The waveguide is excited by a unit force of 1 N applied at the centre of the square shape in z direction.

The excitation signal is a 10 Hanning window with f c = 2.247 MHz. At this frequency, the mode with the lowest attenuation is excited. The problem is then solved for [f min , f max ], with f min = 1.797 MHz and f max = 2.697 MHz, divided into 542 steps. The half-width of the square is meshed with 4 spectral elements of order 8 (see Fig. 13), yielding 20 451 dofs for the quadratic eigenvalue problem (20 451x2 for the linearized one).

High-frequency response

As shown in Ref. [START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF], low-leakage longitudinal modes also occurs in square bars at high frequencies.

Herein, the lowest attenuation is reached for the second lobe of the 12th longitudinal mode at f = 2.247

MHz and equals 147 dB.mm/m. These values are in good agreement with the results of Ref. [START_REF] Treyssède | Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides[END_REF]. They have been obtained with the following PML parameters: h x = h y = 0.25a and γx = γy = 1 + 2j. Note that the wave packets are of slightly smaller amplitudes than in the circular bar. Although the attenuation of the modes is a bit smaller in the square bar, the energy spreads over a wider cross-section area. Moreover, the energy repartition among the modes is also modified. Contrary to the circular case, each lobe of the 12th mode significantly contributes to the response (see Fig. 14a). These differences can be explained by the mode shapes of the square bar. The latter are displayed in Fig. 15a and Fig. 15b for the first two lobes at their lowest attenuation. At the centre of the circular bar, the modal displacement is roughly ten times smaller for the first lobe than for the second lobe (Fig. 12). In the square bar, the modal displacement at the centre (x = y = 0) is only twice smaller for the first lobe than for the second lobe, hence increasing the modal contribution of the first lobe in the response.

Conclusion

A modal approach based on PML has been proposed to compute the forced response of elastic open waveguides. Solving the eigenproblem of the PML waveguide formulation leads to three discrete sets of modes: trapped modes, revealed leaky modes and PML modes. The latter comes from a continuous set of modes, discretized by the truncation of the PML, and fully depends on the PML parameters. A biorthogonality relationship has been derived, applicable to any types of modes, which enables to obtain the forced response with a unique modal superposition, as well as to define the excitability of leaky modes.

A numerical study in a homogeneous elastic medium has shown that, despite their dependence on PML parameters and their exponential decay, the contribution of the discrete set of PML modes enables to reassemble the exact solution (geometrically decaying). In a waveguide, the PML modes enable to obtain the long-term diffraction phenomenon, which occurs at great distances to the source when most energy has leaked into the surrounding medium. However, leaky modes dominate over a wide range of distance inside the core. In this range, leaky modes are sufficient to obtain a good approximation of the solution. Furthermore, the leaky mode contribution remains accurate even with a small PML thickness, hence reducing the computational cost.

Acknowledgements

This work was supported by the Région Pays de la Loire.

Appendix A. Axisymmetrical PML waveguide formulation

Let us go back to Eq. ( 2) written in the cylindrical coordinate system (r, θ, z). The waveguide is axisymmetric, and then the cross-section is reduced to the radial direction r. The core of the waveguide has a radius r = a, and is made of one or several layers. Only axisymmetric solutions are considered, so that ∂(.)/∂θ = 0 and ũθ = 0. Hence the displacement field in Eq. ( 2) is expressed as ũ = [ũ r ũz ] T , and only depends on the variables r and z.

Based on Eq. ( 3), the change of variable ( 4) is applied to Eq. ( 2). In the strain-displacement relation [START_REF] Vogt | The scattering of guided waves in partly embedded cylindrical structures[END_REF], the operators L S and L z are given by:

L S =     1 γ     1 0 0 0 0 0 0 1     ∂ ∂r +     0 0 1 0 0 0 0 0     1 r     , L z =     0 0 0 0 0 1 1 0     . (A.1)
Then, a one-dimensional finite element discretization is applied along the radial direction, yielding the following interpolation on each element: u e (r, z, ω) = N e (r)U e (z, ω). (A.2)

N e (r) is the matrix of interpolation one-dimensional functions. U e (z, ω) is the vector of nodal displacements.

The weak formulation is established by integrating by parts the second order derivatives with respect to r in Eq. ( 2). Regardless of the presence of a PML, detailed steps of the integration can be found in the literature (see [START_REF] Marzani | Time-transient analysis response for ultrasonic guided waves propagating in damped cylinders[END_REF] for instance) and are not repeated here. It yields Eq. ( 6), where the element matrices are given by Eq.( 7) with dS = 2πrdr. Under the assumption used in this appendix, C is reduced to a four-by-four matrix. For instance, the matrix C for an isotropic material is:

C =     λ + 2µ λ λ 0 λ λ + 2µ λ 0 λ λ λ + 2µ 0 0 0 0 µ     .
(A.3)

  2. A general biorthogonality relationship is derived for the modes of the PML waveguide formulation. A validation test case is presented by considering a fully homogenous medium in Sec. 3. Results are compared with an analytic solution, which allows understanding how the PML parameters act on PML modes together with the accuracy of the solution. Section 4 shows numerical results obtained for a circular bar waveguide (widely encountered in civil engineering). Time-domain results are obtained, which enables to clearly identify the distinctive role of PML modes and leaky modes and their ability to approximate the total field. An open square bar example is finally considered, in order to show

  e z in cylindrical coordinates). The signal F (t) is a sinus of centre frequency f c = 60 kHz with an amplitude of 1 N, modulated over 5 cycles by a Hanning window. Material ρ (kg/m 3 ) c l (m/s) c s (m/s) η l (Np/wavelength) η s (Np/wavelength)

Figure 1 :

 1 Figure 1: (Colour online) Sketch of the excited homogeneous medium truncated by a finite PML.

  )dξ is the average value of γ(r) inside the PML region. γ and h are user-defined parameters. The solution is computed in the frequency range [0, f max ], with f max = 120 kHz, divided into 176 frequency steps. The cross-section is discretized with three-nodes line elements of length ∆r = 0.25d ≈ λ min /5, where λ min = c s /f max . The eigenproblem (B.4) is solved with the ARPACK library[START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. Based on an implicit restarted Arnoldi method, this library enables to compute a user-defined number of eigenvalues (denoted as 2M ) that are the closest in absolute value to a user-defined shift, set to 0 in this paper (except in Sec. 4.5). The displacement at the FE node located at r = 0 is then computed from Eq. (23) and denoted as u num z (z, ω). The modal expansion[START_REF] Hayashi | Calculation of leaky Lamb waves with a semi-analytical finite element method[END_REF] contains PML modes only. The accuracy of the numerical solution is assessed as a function of the distance z with the following relative error in the L2-norm:

Figure 2 :

 2 Figure 2: (Colour online) Wavenumber spectrum of an elastic homogeneous medium computed with the PML waveguide formulation at f = fc. PML parameters: γ = 4 + 4j; blue triangle: h = 4d, M = 50 ; red circle: h = 4d, M = 30 ; black cross: h = 3d, M = 50. Dashed black line: theoretical branch cut rotated by -arg(d + hγ) with h = 4d.

Figure 3 :Figure 4 :

 34 Figure 3: (Colour online) Relative error as a function of the distance z between PML waveguide formulation and analytical solutions. Same parameters and color legend as in Fig. 2.

Figure 5

 5 shows the relative error for different values of θ by keeping constant the modulus of the complex thickness, set to |d + hγ| = 25 (h = 4d, M = 50). This means that γ now varies as a function of θ, such that Im(γ) increases as θ increases and conversely. Not too far from the source, the error remains low regardless of the value of θ (less than 1.2% for z/d ∈ [0, 150d]). For greater distances, the error increases with a rate depending on θ. The error remains small (below 2%) in a wide range of values (θ ∈ [30, 60] deg), with a minimum reached around θ = 45 deg. Outside this range, the accuracy of the solution decreases because the PML attenuation is not strong enough.More precisely, when the argument is small (e.g. θ = 20 deg), waves are reflected at the end of the PML (Im(γ) is too low), which deteriorates the solution. Equivalently, θ must not be too close to π/2 either (see e.g. θ = 70 deg), otherwise the rotated branch cuts become too close to the real axis. In this case, the problem becomes equivalent to the computation of trapped modes (these modes occur along the real axis also), and as shown in Ref.[START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF], a proper transverse attenuation of such modes requires to increase Re(γ) instead of Im(γ).

Figure 5 :Figure 6 :

 56 Figure 5: (Colour online) Relative error as a function of distance between PML waveguide formulation and analytical solutions for a complex thickness 25e jθ . Dotted orange line: θ = 20 deg ; crossed blue line: θ = 30 deg ; black line: θ = 45 deg ; dashed-dotted red line: θ = 60 deg ; dashed green line: θ = 70 deg.

Figure 7 :

 7 Figure 7: (Colour online) PML waveguide formulation results (black dots) and results of Ref. [44] (red line) for (a) the energy velocity and (b) the attenuation (steel circular bar embedded into cement grout, viscoelastic case). The PML parameters are: h = a, γ = 1 + 2j. Filtering criterion: η min = 0.6. Point A: f = 0.60 MHz.mm ; Point B: f = 22.84 MHz.mm; Point C: f = 20.92 MHz.mm.

Figure 8 :

 8 Figure 8: (Colour online) Transient axial displacement uz(r = 0, z = 175a, t) in a steel-cement grout circular bar open waveguide, in the pure elastic case (blue line) and in the viscoelastic case (red line), for a low-frequency excitation (fc = 60 kHz). PML parameters: γ = 4 + 4j, h = 4a, M = 51.

Figure 9 :

 9 Figure 9: (Colour online) Energy as a function of distance in a steel-cement grout circular bar open waveguide (elastic case). PML parameters: γ = 4 + 4j, h = 4a. Blue line: modal superposition over 50 PML modes and 1 leaky mode (M = 51) ; dashed red line: contribution of the leaky mode only (exponential decay).

Figure 10 :

 10 Figure 10: (Colour online) (a) Energy velocity curves with the modal contribution of axial displacement at r=0 and z=175a at each frequency (colormap) and (b) Transient axial displacement uz(r = 0, z = 175a, t) in a steel-cement grout circular bar open waveguide (viscoelastic case) for a high-frequency excitation (fc = 2.284 MHz). PML parameters: γ = 1 + 2j, h = 0.1a, M = 20.

Figure 11 :Figure 12 :

 1112 Figure 11: Excitability in a steel-cement grout circular bar open waveguide (viscoelastic case) for a point-force excitation at the core centre oriented in the axial direction. PML parameters: h = a, γ = 1 + 2j.

Figure 11 represents

 11 Figure 11 represents the modulus of the excitability over a wide frequency range for an axial point force applied at the centre of the core. For clarity of the figure, the excitability of each mode is shown only in a frequency range where their maximum energy velocity and their minimum attenuation occur (these regions are of low modal dispersion).

Figure 13 :

 13 Figure13: Spectral element mesh (8th order elements) of a square bar (viscoelastic steel) embedded into viscoelastic cement grout. Black: square bar ; gray: embedding medium (PML).

Figure 14 :

 14 Figure 14: (Colour online) (a) Energy velocity curves with the modal contribution of axial displacement at r=0 and z=175a at each frequency (colormap) and (b) Transient axial displacement uz(r = 0, z = 175a, t) in a steel-cement grout square open waveguide (viscoelastic case) for a high-frequency excitation (fc = 2.247 MHz). PML parameters: γx = γy = 1 + 2j, hx = hy = 0.25a, M = 50.

Figure 14b shows the

  Figure 14b shows the transient axial displacement at z/a = 175 and the modal contribution are displayed in Fig. 14a. M = 50 modes have been used to compute the forced response. Because of the symmetry of the source, only longitudinal modes are excited. As in the circular case (Fig. 10b), two wave packets propagate in the square bar. The first wave packet is due to the arrivals of the 11th and 12th leaky modes, which are the fastest and the less attenuated ones. The second wave packet contains weaker modal contributions of slower velocities.

Figure 15 :

 15 Figure 15: (Colour online) Mode shape Um in the z direction of the 12th longitudinal leaky mode for (a) the first lobe (f = 2.06 MHz) (b) the second lobe (f = 2.247 MHz) in a steel-cement grout square open waveguide.

Table 1 :

 1 Material properties

Appendix B. Linearization of the quadratic eigenvalue problem

The system [START_REF] Kausel | Stiffness matrices for layered soils[END_REF] can be linearized [START_REF] Marzani | Time-transient analysis response for ultrasonic guided waves propagating in damped cylinders[END_REF][START_REF] Tisseur | The quadratic eigenvalue problem[END_REF][START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF], and yields for instance:

with:

The free response problem is considered by dropping the source term:

and yields the generalized eigenvalue problem:

Because of the PML, the matrices are complex. A and B are not Hermitian.

Appendix C. Fourier transform of Eq. (32)

This appendix gives the analytic Fourier transform

)e +jωt dt of Eq. [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF].

First, the near-field term is written as a convolution product:

with:

The Fourier transform of a(t) * F (t) is equal to A(ω)F (ω), where F (ω) is the spectrum of the excitation force, and A(ω) is obtained from an integration by parts, yielding:

Finally, the Fourier transform of Eq. ( 32) is: