F Treyssède
email: fabien.treyssede@ifsttar.fr

K L Nguyen

A.-S Bonnet-Bendhia

C Hazard

Finite element computation of trapped and leaky elastic waves in open stratified waveguides

Keywords: Leaky modes, Open waveguides, Perfectly matched layers, Finite element

Elastic guided waves are of interest for inspecting structures due to their ability to propagate over long distances. In numerous applications, the guiding structure is surrounded by a solid matrix that can be considered as unbounded in the transverse directions. The physics of waves in such an open waveguide significantly differs from a closed waveguide, i.e. for a bounded cross-section. Except for trapped modes, part of the energy is radiated in the surrounding medium, yielding attenuated modes along the axis called leaky modes. These leaky modes have often been considered in non destructive testing applications, which require waves of low attenuation in order to maximize the inspection distance. The main difficulty with numerical modeling of open waveguides lies in the unbounded nature of the geometry in the transverse direction. This difficulty is particularly severe due to the unusual behavior of leaky modes: while attenuating along the axis, such modes exponentially grow along the transverse direction. A simple numerical procedure consists in using absorbing layers of artificially growing viscoelasticity, but large layers may be required. The goal of this paper is to explore another approach for the computation of trapped and leaky modes in open waveguides. The approach combines the so-called semi-analytical finite element method and a perfectly matched layer technique. Such an approach has already been successfully applied in scalar acoustics and electromagnetism. It is extended here to open elastic waveguides, which raises specific difficulties. In this paper, two-dimensional stratified waveguides are considered. As it reveals a rich structure, the numerical eigenvalue spectrum is analyzed in a first step. This allows to clarify the spectral objects calculated with the method, including radiation modes, and their dependency on the perfectly matched layer parameters. In a second step, numerical dispersion curves of trapped and leaky modes are compared to analytical results.

Introduction

Context and objectives

Elastic guided waves are of interest for inspecting elongated structures due to their ability to propagate over long distances. When the guiding structure is embedded into a solid matrix, the waveguide can usually be considered as open (unbounded in its transverse direction). The physics of waves in open waveguides significantly differs from closed waveguides. The modes that are generally considered in such structures are either trapped or of radiation type. Another category of slightly more artificial objects, called leaky modes, can also be considered.

Trapped modes can be seen as perfectly guided waves in the sense that on the one hand, they propagate without attenuation along the waveguide axis (their axial wavenumber is real) and on the other hand, they decay exponentially in the transverse directions (their transverse wavenumber is pure imaginary). Such waves are confined near the core or the surface of the structure and do not loose energy in the infinite surrounding medium. They constitute a discrete set, possibly empty. For scalar waveguides (i.e. for acoustic, SH or electromagnetic waves), trapped modes only exist when the velocity of waves in the surrounding medium is greater than in the core [START_REF] Collin | Field Theory of Guided Waves[END_REF]. For elastic waveguides, this existence criterion is more complex because two kinds of waves, longitudinal and shear waves, are present [START_REF] Thurston | Elastic waves in rods and clad rods[END_REF]. In particular, the coupling between both can lead to the presence of surface waves, called Stoneley waves [START_REF] Stoneley | Elastic waves at the surface of separation between two solids[END_REF], which are trapped at the interface between different materials. Besides, a phenomenon analogous to the scalar case occurs for body waves: trapped waves are present when the shear wave velocity is greater in the surrounding medium, but no longer when this velocity is lower. The latter configuration is widely encountered in civil engineering applications (tendons, rockbolts, buried pipes, railways), the outer medium usually consisting of soft materials (cement, concrete, soil).

Unlike trapped modes, radiation modes have an oscillating behavior in the transverse directions (their transverse wavenumber is real) and can be either oscillating or evanescent in the longitudinal direction (real or pure imaginary axial wavenumber). They constitute a continuous spectrum of standing 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]. This continuum is intrinsic to the unbounded nature of the cross-section of the waveguide. It allows to represent all radiation phenomena which may occur in the surrounding medium. In the context of non-destructive evaluation (NDE), radiation modes are of less interest because the continuum hides the key property which is required: the ability to propagate over long distances. Actually, such an information is condensed in the notion of leaky modes.

As opposed to trapped modes, leaky modes can be interpreted as imperfectly guided waves in the sense that they are waves propagating with attenuation along the waveguide axis. This attenuation is due to the leakage of energy into the surrounding medium. Owing to these radiation losses, the axial wavenumber becomes complex. In practice, such losses can strongly limit the application of guided wave techniques. Hence, an accurate determination of leaky modes appears to be a necessary step for the NDE of embedded civil structures based on guided waves.

Leaky modes are true solutions of the equations. However, their behavior is somewhat unusual at infinity: while exponentially decreasing along the axis, leaky modes exponentially grow along the transverse direction. This behavior can be readily deduced from the dispersion relation,

k 2 + k 2 x = ω 2 /c 2 ∞ ,
where k denotes the axial wavenumber, k x is the transverse wavenumber, ω is the angular frequency and c ∞ denotes the speed of waves inside the infinite surrounding medium (one assumes a scalar waveguide for simplicity). The imaginary part of the dispersion relation then yields:

k ′ x k ′′ x = -k ′ k ′′ ,
where ′ and ′′ denote real and imaginary parts respectively. For instance, considering a leaky wave propagating and attenuating in the positive direction (k ′ > 0, k ′′ > 0) yields the inequality k ′ x k ′′ x < 0, which means that the wave is necessarily exponentially growing (k ′′ x < 0) as it radiates toward a given transverse direction (k ′ x > 0). The exponential growth of leaky modes is rather well-known in electromagnetism (see [START_REF] Collin | Field Theory of Guided Waves[END_REF] for instance). This particularity has sometimes been mentioned in elastodynamics [START_REF] Alsop | Reflection and transmission of inhomogeneous waves with particular application to Rayleigh waves[END_REF][START_REF] Simmons | Leaky axisymmetrical modes in infinite clad rods .1[END_REF][START_REF] Hladky-Hennion | Conical radiating waves from immersed wedges[END_REF][START_REF] Schroder | On the complex conjugate roots of the Rayleigh equation: The leaky surface wave[END_REF][START_REF] Vogt | The scattering of guided waves in partly embedded cylindrical structures[END_REF]. From a mathematical point of view, leaky modes are complex poles that do not belong to the proper eigenvalue spectrum because such poles are located on a forbidden Riemann sheet [START_REF] Chapman | Lamb's problem and comments on the paper "on leaking modes" by Usha Gupta[END_REF][START_REF] Collin | Field Theory of Guided Waves[END_REF]. These modes constitute a discrete set which is actually not part of the complete modal set constituted by the discrete set of trapped modes and the continuous spectrum of radiation modes representing the contribution of body 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] Harris | Comment on "on the complex conjugate roots of the Rayleigh equation: The leaky surface wave[END_REF][START_REF] Margerin | Generalized eigenfunctions of layered elastic media and application to diffuse fields[END_REF][START_REF] Hu | Understanding leaky modes: slab waveguide revisited[END_REF]. However, leaky modes have the potential to give a good representation of the continuum over an area close to the core or the surface of the waveguide. From a physical point of view, leaky modes are then able to resemble the transverse gradual increase of a wave beam, which occurs over a limited transverse distance inside the surrounding region, as the beam propagates far from the excitation zone [START_REF] Alsop | Reflection and transmission of inhomogeneous waves with particular application to Rayleigh waves[END_REF][START_REF] Simmons | Leaky axisymmetrical modes in infinite clad rods .1[END_REF][START_REF] Hu | Understanding leaky modes: slab waveguide revisited[END_REF].

The concept of leaky modes have turned out to be meaningful for the NDE of solid waveguides [START_REF] Lowe | Comparison of the modal porperties of a stiff layer embedded in a solid medium with the minima of the plane-wave reflection coefficient[END_REF][START_REF] Lowe | The influence of the modal properties of a stiff layer embedded in a solid medium on the field generated in the layer by a finite-sized transducer[END_REF][START_REF] Beard | Non-destructive testing of rock bolts using guided ultrasonic waves[END_REF][START_REF] Beard | Ultrasonic guided waves for inspection of grouted tendons and bolts[END_REF][START_REF] Ryden | Guided wave propagation in three-layer pavement structures[END_REF][START_REF] Castaings | Finite element model for waves guided along solid systems of arbitrary section coupled to infinite solid media[END_REF], where the problem consists in looking for modes having low attenuation in order to maximize the inspection distance. Analytical tools have been developed for modeling open elastic waveguides, based on the Thomson-Haskell method [START_REF] Thomson | Transmission of elastic waves through a stratified solid medium[END_REF][START_REF] Haskell | Dispersion of surface waves on multilayered media[END_REF] or the global matrix method [START_REF] Lowe | Plate waves for the NDT of diffusion bonded titanium[END_REF][START_REF] Pavlakovic | Leaky guided ultrasonic waves in NDT[END_REF]. Methods based on Debye series [START_REF] Laguerre | Ultrasonic transient bounded-beam propagation in a solid cylinder waveguide embedded in a solid medium[END_REF] have also been proposed. Analytical techniques are yet limited to simple geometries (plates, cylinders).

For complex geometries, a classical approach relies on the finite element (FE) discretization of the eigenproblem in the transverse direction, which is often referred to as the semi-analytical finite element (SAFE) method (see for instance [START_REF] Gavric | Computation of propagative waves in free rail using a finite element technique[END_REF][START_REF] Liu | Strip element method to analyze wave scattering by cracks in anisotropic laminated plates[END_REF][START_REF] Hayashi | Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example[END_REF][START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF]). Yet for open waveguides, a difficulty arises because the geometry is unbounded in the transverse direction. This difficulty is particularly severe due to the exponential growth of leaky modes. A simple procedure that can be used with existing codes has been proposed in [START_REF] Castaings | Finite element model for waves guided along solid systems of arbitrary section coupled to infinite solid media[END_REF][START_REF] Fan | Torsional waves propagation along a waveguide of arbitrary cross section immersed in a perfect fluid[END_REF], which consists in using absorbing layers (AL) of artificially growing viscoelasticity. In practice, large layers may be required in order to reduce artificial reflections by the absorbing layer. Other methods have been specifically developed in order to deal with the FE modeling of open waveguides. Hladky-Hennion et al. [START_REF] Hladky-Hennion | Finite element modelling of radiating waves in immersed wedges[END_REF] have presented a procedure based on a non-reflecting condition together with an original change of variable, but the approach has been limited to structures immersed in fluids. Recently, Mazzoti et al. [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] have coupled the SAFE method with a boundary element technique, which yields nonlinear eigenproblems and requires higher computation times.

Instead of AL, an alternative technique consists in using perfectly matched layers (PML). Contrary to the approach developed in [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], the PML approach has the advantage to lead to a linear eigenvalue problem (same for AL). This technique has already been applied to the scalar wave equation for acoustic waveguides [START_REF] Bonnet-Bendhia | A multimodal method for non-uniform open waveguides[END_REF][START_REF] Pelat | A coupled modal-finite element method for the wave propagation modeling in irregular open waveguides[END_REF] as well as electromagnetic waveguides [START_REF] Huang | The perfectly matched layer boundary condition for modal analysis of optical waveguides: Leaky mode calculations[END_REF][START_REF] Ould Agha | On the use of PML for the computation of leaky modes. An application to microstructured optical fibres[END_REF]. Closely linked to such problems, a mathematical analysis has been reported in [START_REF] Kim | The computation of resonances in open systems using a perfectly matched layer[END_REF] for the FE computation of acoustic resonances in open cavities with PMLs.

The main goal of this work is to compute trapped and leaky modes in open waveguides by applying a SAFE-PML approach to the equilibrium equations of elastodynamics (non-scalar). The numerical eigenspectrum, including radiation modes, is analyzed in order to clarify the influence of the various PML parameters. Theoretically, the PML technique allows to strongly attenuate any type of wave without reflection, thanks to an analytical continuation of equations into complex spatial coordinates. Compared to AL, one expects that the perfectly matched property will allow reduction of the artificial layer size. Another difference between the two approaches is seldom mentioned: computing leaky modes with PMLs is mathematically relevant, since both leaky modes and PMLs are defined through analytic extensions. Conversely, the ability of AL to approximate leaky modes has, up to our knowledge, no theoretical explanation.

In this paper, one-dimensional modal problems are considered, i.e., two-dimensional elastic waveguides corresponding to stratified planes. For clarity, Sec. 1.2 recalls the mathematical definition of the various kind of modes in the continuous problem. In Sec. 2, the SAFE-PML approach is presented. The eigenvalue spectrum of the method is analyzed in Sec. 3. In Sec. 4, numerical dispersion curves are computed and validated thanks to analytical results of the literature.

Mathematical definition of the various kinds of modes

One considers a two-dimensional waveguide S × R in the ( x, z) plane, where z is the waveguide axis, x is the transverse direction and S ⊂ R represents the cross-section of the waveguide. The tilde notation will be explained later in section 2.3 with the introduction of the PML method. The waveguide is filled by a linearly elastic material and we assume small plane strains and displacements in the ( x, z) plane. The time-harmonic dependence is chosen as e -iωt . As this paper is restricted to the study of propagation modes, acoustic sources and external forces are dropped. Consider the two-dimensional elastodynamic equations satisfied by the displacement vector U( x, z) :

∇ • σ(U) -ω 2 ρ U = 0 in S × R,
where ρ is the mass density and σ(U) is the stress tensor, which is related to the strain tensor ǫ(U) = (∇U+(∇U) T )/2 by the linear relation σ(U) = C : ǫ(U) where C is the stiffness tensor. The superscript T denotes the matrix transpose. When S R, a free or fixed boundary condition is added to the above equations.

We assume that both C and ρ depend only on the transverse coordinate x, which means that we consider uniform rectilinear waveguides. The modes we are interested in are particular solutions of the above equations which have the form:

U( x, z) = ũ( x) e ikz for k ∈ C. (1)
The equations which are satisfied by ũ are simply deduced from those of U by replacing derivatives ∂/∂z with products by ik. We are thus led to a one-dimensional problem for ũ which can be written formally as a linear equation:

L(k) ũ = 0, (2)
where L(k) has a quadratic dependence on k.

In the case of a closed waveguide, that is, a bounded cross-section S , it is well known that this problem has nonzero solutions ũ only for a discrete infinite set of values of k ∈ C. For real values of k, the associated solutions ũ correspond to guided waves which propagate without attenuation in the z-direction.

On the other hand, the case of an open waveguide, that is, an unbounded cross-section S , is more involved. The different kinds of modes that can be considered depend on the behavior of ũ( x) at infinity (as | x| → ∞). This behavior can be easily expressed if we assume that outside a bounded part of the cross-section, our waveguide is made of an isotropic homogeneous medium. In this case, following the Helmholtz decomposition method (see [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF] for instance), the displacement field appears as a superposition of longitudinal and shear waves:

ũ = ∂/∂ x ik φl + -ik ∂/∂ x φs (3)
where φl and φs are scalar potentials which satisfy:

∂ 2 φl ∂ x2 +       ω 2 c 2 l -k 2       φl = 0 and ∂ 2 φs ∂ x2 + ω 2 c 2 s -k 2 φs = 0,
where c l and c s are the longitudinal wave and shear wave velocities respectively. Hence φl is a linear combination of two exponential functions exp(±ik l x) where k 2 l = ω 2 /c 2 lk 2 , and the same for φs with

k 2 s = ω 2 /c 2 s -k 2 : φl/s ( x) = A + l/s exp(+ik l/s x) + A - l/s exp(-ik l/s x),
where A ± l/s ∈ C. If no condition on φl and φs is imposed, then for all k ∈ C, one can always find a non-zero solution ũ to [START_REF] Thurston | Elastic waves in rods and clad rods[END_REF]. For non-real values of k l and k s , this solution is in general exponentially increasing unless the associated coefficients A ± l/s vanish. It is generally admitted that an exponentially increasing behavior of the field is not physically relevant, so that the potentials φl and φs are naturally chosen either oscillating or exponentially decreasing, depending on whether k l and k s take real values or not. The solutions which are exponentially decreasing are referred to as trapped modes since they remain confined near the core of the waveguide. All other bounded solutions are called radiation modes for they represent waves which propagate outside the core. The values of k associated with trapped modes constitute a finite set which represents the point spectrum of the transverse elastodynamic operator L(k), whereas those associated with radiation modes compose its continuous spectrum.

But for non-real values of k l and k s , there is another category of modes which are solutions to Eq. ( 2) and grow exponentially at infinity. These solutions, called leaky modes, are of physical interest. They correspond to the case where each potential φl/s only involves one exponential function (i.e., either A + l = 0 or A - l = 0, and the same for A ± s ) and at least one of these potentials is exponentially increasing. The values of k associated with such leaky modes make up a discrete set but are not, strictly speaking, spectral objects. Their mathematical characterization derives from a complicated analytic continuation process which is out of the goal of the present paper. We refer for instance to the analytic dilation technique (also called complex scaling) described in [START_REF] Hislop | Introduction to Spectral Theory[END_REF]. The latter technique is actually very close to the PML method used here for the numerical computation of modes.

Numerical method

Variational formulation for elastodynamics

The starting point of the SAFE method consists in the variational formulation of the one-dimensional elastodynamic problem [START_REF] Thurston | Elastic waves in rods and clad rods[END_REF]. For the sake of simplicity, we use stress and strain vectors instead of tensors. It can be shown that the variational formulation of (2) can be written as:

S δǫ T σd x -ω 2 S ρδ ũT ũd x = 0, (4)
with the stress-strain relationship σ = Cǫ and the strain-displacement relations:

ǫ = L x ∂ ∂ x + ik L z ũ, δǫ = L x ∂ ∂ x -ik L z δ ũ. (5)
L x =           1 0 0 0 0 1           , L z =           0 0 0 1 1 0           .

SAFE modeling of closed waveguides

In this subsection, the SAFE formulation is briefly recalled for closed waveguides (i.e. waveguides in vacuum with a bounded cross-section S ). In this case, the tilde notation introduced in the previous subsection is useless and can be omitted. For details on the SAFE method, the reader might refer to [START_REF] Gavric | Computation of propagative waves in free rail using a finite element technique[END_REF][START_REF] Liu | Strip element method to analyze wave scattering by cracks in anisotropic laminated plates[END_REF][START_REF] Hayashi | Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example[END_REF][START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF][START_REF] Treyssède | Elastic waves in helical waveguides[END_REF] for instance.

The SAFE method consists in applying a 1D finite element method along the transverse direction x to the variational problem (4). The displacement is then expressed on each element e as:

u(x) = N e (x)U e , (6)
where N e (x) is a matrix of nodal interpolating functions and U e contains nodal displacements on the element e.

The FE discretization of the variational formulation (4) finally yields:

{K 1 -ω 2 M + ik(K 2 -K T 2 ) + k 2 K 3 }U = 0 (7)
with the following elementary matrices:

K e 1 = e (L x N e ,x ) T CL x N e ,x dx, K e 2 = e (L x N e ,x ) T CL z N e dx, K e 3 = e N eT L T z CL z N e dx, M e = e ρN eT N e dx. ( 8
)
where the notation N e ,x is used for ∂N e /∂x. For each frequency ω, Eq. ( 7) corresponds to an eigenvalue problem whose eigensolutions yields the propagation modes (k j , U j ).

SAFE-PML method

Open waveguides, that is, unbounded cross-sections S , are now considered. We assume that outside a possibly inhomogeneous region representing the core of the waveguide, the medium is homogeneous. In order to account for the presence of unbounded layers, the basic idea consists in closing the waveguide by replacing the unbounded homogeneous region with a PML, truncated at some finite distance. For clarity, Fig. 1 gives the example of a bilayer open waveguide. Throughout this paper, h will denote the PML thickness and d will denote the distance of the PML interface from the core. a is the core thickness.

The principle of PMLs is easily understood in our one-dimensional situation. As shown in the introduction, the modes we are interested in can be expressed outside the core by means of scalar potentials φl and φs which are linear combinations of exponential functions exp(±ik l/s x), where k l/s denotes the transverse wavenumber, longitudinal or shear. Here x represents the original transverse coordinate, which is real. Consider for instance the case of longitudinal waves for x > 0. As mentioned in the introduction, the imaginary part of k l is negative for leaky modes propagating and attenuating as z → +∞, so that exp(ik l x) increases exponentially as x tends to +∞. The key point is to notice that this function extends to an entire function for complex values of x. Hence, instead of considering real x, one can choose in the complex plane a particular path x(x) parameterized by a real variable x ∈ (0, +∞) and such that exp(ik l x(x)) decays exponentially as x tends to +∞, which amounts to saying that arg x(x) >arg k l for large enough x (in this paper, arg denotes the principal argument defined to lie in the interval ] -π, +π]). The principle of PMLs consists in replacing the initial equations satisfied by ũ( x) for real x with those satisfied by ũ( x(x)), which can be seen as a simple (but unusual) change of variable x → x. This amounts to considering an infinite PML. As explained in Ref. [START_REF] Hislop | Introduction to Spectral Theory[END_REF] for a scalar equation, the effect of such a PML is on the one hand, to rotate in the complex plane the two parts of the continuous spectrum, associated respectively with longitudinal and shear radiation modes, and on the other hand, to reveal the leaky modes located in the region covered by the displacement of the continuous spectrum.

From a numerical point of view, the PML domain is truncated at a finite distance (here x = a + d + h) where an arbitrary boundary condition, usually of Dirichlet type, is imposed. As leaky modes decay exponentially in the PML, this truncation should not perturb strongly their computation provided the PML is thick enough. The truncation of the PML yields a discretization of the rotated continuous spectrum, and the associated radiation modes are oscillating in the PML. These facts will be made clear in Sec. 3.

In the present paper, we choose:

x(x) = x 0 γ(ξ)dξ, (9)
where γ is a complex-valued function of x satisfying:

• γ(x)=1 for x a + d; • Im{γ(x)} > 0 for x > a + d.
Hence, from Eq. ( 9), the change of variables x → x yields for any function f :

∂ f ∂ x = 1 γ ∂ f ∂x , d x = γdx, (10)
where f ( x) = f ( x(x)) = f (x). For instance, the displacement and the mass density are rewritten as ũ( x) = u(x) and ρ( x) = ρ(x). Note that ρ as well as the other physical features involved in matrix C are constant in the PML region.

Hence their respective analytic continuations are simply these constants. From Eqs. ( 10) and ( 5), the strain-displacement relations are expressed as:

ǫ = 1 γ L x ∂ ∂x + ikL z u, δǫ = 1 γ L x ∂ ∂x -ikL z δu.
With Eq. ( 6), the one-dimensional FE discretization of the variation formulation (4) along x gives exactly the same form of eigenproblem as Eq. ( 7), but with the following elementary matrices:

K e 1 = e 1 γ (L x N e ,x ) T CL x N e ,x dx, K e 2 = e (L x N e ,x ) T CL z N e dx, K e 3 = e γN eT L T z CL z N e dx, M e = e ργN eT N e dx.
The only difference with the elementary matrices (8) lies in the γ terms, appearing in the expression of K e 1 , K e 3 and M e .

Solving the eigenvalue problem

The eigenproblem ( 7) is quadratic for k, which cannot be directly held by standard eigensolvers. This eigensystem can be transformed into the more suited generalized linear eigensystem [START_REF] Baronian | Couplage des méthodes modale et éléments finis pour la diffraction des ondes élastiques guidées : Application au contrôle non destructif (Coupling of modal and finite elements methods for the diffraction of guided elastic waves -application to non destructive testing[END_REF]:

(A -kB) Û = 0 with: A = 0 I -(K 1 -ω 2 M) -i(K 2 -K T 2 ) , B = I 0 0 K 3 , Û = U kU .
Note that in the presence of PMLs, K 1 , K 3 and M are complex symmetric but not hermitian matrices so that neither A nor B are Hermitian, which makes the numerical treatment of the above eigensystem complicated.

Other linearization forms of the quadratic eigenproblem are possible, as reviewed in Ref. [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. Our choice is justified by the fact that B is non-singular (B is also hermitian and positive-definite without PML). Although it could be replaced with any invertible matrix (e.g. K 3 or M), the identity matrix is chosen in the first row of A and B in order to reduce memory allocation.

In this paper, the above eigensystem is solved using the ARPACK library [START_REF] Lehoucq | ARPACK User's Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods[END_REF], based on the implicitly restarted Arnoldi method, which allows to compute a specified number of eigenvalues that are nearest in absolute value to a user-defined shift. The non-hermitian character of matrices with PML may raise problems with the norm of the resolvent, which can be quite large far away from poles and yield spurious eigenvalues [START_REF] Kim | The computation of resonances in open systems using a perfectly matched layer[END_REF]. Yet this drawback can be reduced by setting the PML not too far from the core, as will be shown in Sec. 3.2. The question of eigensolver accuracy due to non-hermitian matrices is left beyond the scope of this paper. In this paper, the possible loss of accuracy of the eigensolver appears to be of less importance than that of FE discretization or PML parametrization, since it will be shown in Secs. 3 and 4 that leaky and trapped modes can be accurately computed with the SAFE-PML method.

It is worth noting that, due to the symmetry of K 1 , K 3 and M, if k is an eigenvalue, then -k is also an eigenvalue. Hence, the eigenproblem has two sets of eigensolutions (k + j , U + j ) and (-k + j , U - j ) ( j = 1, ..., n), representing n positivegoing and n negative-going wave modes.

PML parameters

From the definition of the complex coordinate (9), the PML depends on three user-defined parameters: the thickness h, the distance d and the non-dimensional attenuating function γ(x).

A proper choice of the attenuating function is crucial for an efficient wave absorption inside the layer. In timeharmonic source problems (excluding modal problems), this function has often been chosen as 1 + iσ(x)/ω, where σ(x) is a continuous real function increasing inside the PML [START_REF] Harari | Studies of FE/PML for exterior problems of time-harmonic elastic waves[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Bermudez | An optimal perfectly matched layer with unbounded absorbing function for timeharmonic acoustic scattering problems[END_REF]. This particular choice derives naturally from the use of time-dependent PMLs. However, the problem slightly differs for the computation of modes in open waveguides.

As explained in section 2.3, if k l/s denotes the transverse wavenumber (shear or longitudinal) of a given mode, the effect of an infinite PML is to transform its actual behavior exp(ik l/s x) for real x into an exponentially decaying behavior exp(ik l/s x(x)) as x tends to +∞, where x(x) denotes the complex path defined in [START_REF] Schroder | On the complex conjugate roots of the Rayleigh equation: The leaky surface wave[END_REF]. Hence a truncated PML will act properly if exp(ik l/s x(x)) is sufficiently small at the end of the PML so that only small reflections occur. From Eq. ( 9), we have x(a

+ d + h) = a + d + γh where: γ = 1 h a+d+h a+d γ(s)ds
denotes the average value of γ(x) inside the PML. As a consequence, the total attenuation from the interface to the end of the PML is given by: exp -Im(k l/s γh) = exp -|k l/s || γ|h sin(arg k l/s + arg γ) .

This formula tells us on the one hand that the PML can play its attenuating role for leaky modes such that arg k l/s > arg γ. The attenuation vanishes if arg k l/s becomes close toarg γ. Of course, increasing arg γ increases the region of the complex plane were leaky modes can be computed. Besides, the PML attenuation increases with | γ|h. But for a given h, the value of | γ| cannot be chosen arbitrarily large. Otherwise the attenuation would be too fast, as well as the wavelength in the PML region too small (increasing | γ|h also increases Re(k l/s γh)), to be approximated by the FE discretization. High values of | γ|h indeed require refined meshes. This drawback is well-known in PML problems (see [START_REF] Bermudez | An optimal perfectly matched layer with unbounded absorbing function for timeharmonic acoustic scattering problems[END_REF] for instance).

As far as trapped modes are concerned (arg k l/s = π/2), waves remain exponentially decaying provided that arg γ < π/2. However, the PML can enhance this decay if γ is such that Im(k l/s γ) > Im(k l/s ), or equivalently

| γ| cos arg γ = Reγ > 1.
When considering open waveguide problems, a constant complex function has generally been adopted [START_REF] Ould Agha | On the use of PML for the computation of leaky modes. An application to microstructured optical fibres[END_REF][START_REF] Bonnet-Bendhia | A multimodal method for non-uniform open waveguides[END_REF][START_REF] Pelat | A coupled modal-finite element method for the wave propagation modeling in irregular open waveguides[END_REF]:

γ(x) = 1 if x ≤ a + d, γ if x > a + d,
where γ is a user-defined parameter. Such a choice somewhat simplifies the mathematical analysis of the problem. In this paper, the following continuous parabolic function for both the real and imaginary parts of γ(x) is also considered and appears to provide better results as will be shown in Secs. 3.2.2 and 4.2:

γ(x) =            1 if x ≤ a + d, 1 + 3( γ -1) x -a -d h 2 if x > a + d.
Note that any frequency dependence is dropped in order to avoid the computation of SAFE matrices at each frequency.

Analysis of PML effects on the eigenvalue spectrum

As it reveals a rich structure, a better understanding of the numerical spectrum is required before computing the dispersion curves of trapped and leaky modes. This section is hence devoted to the study of the eigenspectrum obtained by the SAFE-PML method. First, a simplified analytical model is proposed to outline the main effects of PMLs on the elastodynamic spectrum. Second, numerical experiments are performed with the SAFE-PML method. The influence of PML parameters on leaky and radiation modes is highlighted.

The case of trapped mode is not discussed in this section. As trapped modes naturally decrease in the transverse direction, their computation is less critical. Trapped modes will be considered in Sec. 4.2.

Analytical solution

Let us assume a semi-infinite isotropic homogeneous medium in the transverse direction x ∈ [0, +∞[. This problem can be seen as an open waveguide problem (with a = 0). Let us close the problem with a PML of thickness h and position d. L = d + h denotes the total thickness of the waveguide. The equations to be solved are the equilibrium equations of elastodynamics together with the boundary conditions at x = 0 and L, or equivalently (see [START_REF] Schroder | On the complex conjugate roots of the Rayleigh equation: The leaky surface wave[END_REF]) at x = 0 and L given by:

L = x(L) = L 0 γ(ξ)dξ = d + h γ. ( 11
)
which can be referred to as the complex thickness of the medium. As in the introduction, the Helmholtz decomposition tells us that the displacement field ũ can be expressed as the superposition (3) of longitudinal and shear waves, where the potentials φl and φs satisfy respectively:

∂ 2 φl ∂ x2 +       ω 2 c 2 l + λ       φl = 0 and ∂ 2 φs ∂ x2 + ω 2 c 2 s + λ φs = 0, (12)
where λ = -k 2 , c l and c s are the longitudinal wave and shear wave velocities respectively.

For free or fixed boundary conditions, the potentials φl and φs are coupled. However, mixed boundary conditions yield uncoupled problems leading to analytical solutions. Let us apply the following mixed conditions: ũx | 0, L = 0, σxz | 0, L = 0. These conditions yield the uncoupled boundary conditions for potentials:

∂ φl ∂ x 0, L = 0, φs | 0, L = 0, (13)
and the eigenvalues of Eqs. ( 12) and ( 13) can be readily obtained as:

λ (n) l = - ω 2 c 2 l + n 2 π 2 L2 , λ (m) s = - ω 2 c 2 s + m 2 π 2 L2 ,
together with the associated eigenfunctions:

φ(n) l ( x) = A (n) l cos nπ x L , φ(m) s ( x) = A (m) s sin mπ x L , (14)
where n ≥ 0 and m > 0 denote integers. A (n) l and A (m) s are normalization factors. The above derivation of the modes may not seem rigorous since the Helmholtz decomposition (3) and the differential equations ( 13) concern a complex domain { x(x) ∈ C; 0 < x < L}. If we rewrite the equations in terms of the real variable x ∈ (0, L), it is clear that the above formal derivation yields particular solutions of the initial elastodynamic problem. And these are the only solutions, which is easily justified here, since the problem is one-dimensional, hence we know in advance the dimension of the space of solutions of the various differential equations.

From this result, it can be concluded that:

• two spectra are present for elastodynamics instead of one for acoustics: they are associated respectively to compressional waves and shear waves;

• in the complex plane (Reλ, Imλ), the poles associated with the compressional wave spectrum and the shear wave spectrum are located on a half line of origin -ω 2 /c 2 l and -ω 2 /c 2 s respectively, both lines having a rotation angle equal to -2 arg L with respect to the axis Im λ = 0;

• the spectra are getting denser by increasing | L|, and become continuous in the limit | L| → ∞.

For clarity, Fig. 2a plots an example of spectrum calculated for d = 0.1, h = 0.9 and γ = 1 + 2i.

Note that the initial open problem is without PML and can be recovered in the limit L = L → ∞, yielding two continua of modes both located on horizontal half-lines (the rotation angle being zero). It can be deduced that the introduction of an infinite PML rotates each continuum in the complex plane, and that the truncation of the layer discretizes each continuum. In the truncated PML problem, each of the two spectra can then be viewed as a discretized continuum of radiation modes, representing standing waves in the transverse direction.

It can be noticed that the eigenvalues and the rotation angle only depend on the global PML parameter L, function of d, h and γ. From Eq. ( 11), the angle can be expressed as:

θ = -2 arg (d + h γ). ( 15
)
This result coincides with the mathematical derivation obtained in Ref. [START_REF] Goursaud | étude mathématique et numérique de guides d'ondes ouverts non uniformes, par approche modale (mathematical and numerical study of non uniform open waveguides, modal approach[END_REF] for a bilayer scalar waveguide with a constant γ function, where d should be understood as the distance between the PML and the core of the waveguide. For a given value γ, the eigenvalues and the rotation angle are not affected by the absorbing function profile γ(x). However, the eigenfunctions ( 14) are strongly dependent on that profile. Figure 2b depicts the real part of the eigenfunction ϕ (5) l (x) = φ(5) l ( x(x)), where x has been calculated from Eq. ( 9) both for a constant and a parabolic profile (d = 0.1, h = 0.9 and γ = 1 + 2i). With a constant γ function, it can be observed that a slope discontinuity occurs at the PML interface. With a continuously increasing function, no slope discontinuity occurs but the oscillations are significantly stretched near the interface and narrowed near the PML end. This will have important consequence in the SAFE-PML spectrum, as discussed in Sec. 3.2.2.

Note that the present analytical solution does not imply any trapped or leaky modes. In particular, no Rayleigh wave occurs because the boundary conditions are mixed at the surface x = 0 (Rayleigh waves only occur with a free boundary condition [START_REF] Rayleigh | On waves propagated along the plane surface of an elastic solid[END_REF]). PML effects on leaky modes are numerically studied in the next subsection.

SAFE-PML numerical experiments 3.2.1. Test case description

In the following, numerical tests are performed for a three layer waveguide taken from the work of Lowe [START_REF] Lowe | Plate waves for the NDT of diffusion bonded titanium[END_REF]. Figure 3a depicts the structure, a thin core of Alpha case of 0.1mm embedded on both sides by Titanium half-spaces. Material characteristics are given in Table 1. In this test case, no trapped modes occur. The eigenvalue spectrum is computed for a single frequency Ω = 17.68, where Ω = ωa/c s0 is the dimensionless angular frequency and c s0 denotes the shear wave velocity of the core (Alpha case).

In SAFE-PML simulations, finite elements are quadratic (three-node elements). Both constant and parabolic γ functions are considered. Choosing the optimized value of γ is a difficulty inherent to PML problems: its modulus should not be set too low for a proper absorption of waves, and not too high otherwise numerical reflection will occur. In this section, the averaged value of absorbing functions is set to γ = 1 + 2i.

Figures 4a and4b compare the dimensionless spectrum for λa 2 obtained for γ = 1 + 0.2i and γ = 1 + 0.5i. A constant γ function is used. The PML thickness is set to h = 0.9a. The FE size, denoted ∆x, is equal to 0.05a. From these figures, two kinds of modes can be observed. The first kind of modes forms two branches starting from (-ω 2 a 2 /c 2 s , 0) and (-ω 2 a 2 /c 2 l , 0) respectively, where c s and c l denote the bulk wave velocities of Titanium. Except far from the real axis, each branch resembles a straight line whose rotation angle is approximately equal to -2 arg(d+h γ) (about -120 • in Fig. 4a and -45 • in Fig. 4b). By analogy with Sec. 3.1, these branches can be interpreted as a discretized continuum of longitudinal waves and a discretized continuum of shear waves. Modes belonging to these branches are hence radiation modes, mainly oscillating in the surrounding medium, which are strongly dependent on the PML parameters.

c l (m/s) c s (m/s) ρ (kg/m 3 )
The second kind of modes corresponds to leaky modes. These modes, which are intrinsic to the physics, can be compared with analytical solutions obtained from a Thomson-Haskell matrix approach [START_REF] Thomson | Transmission of elastic waves through a stratified solid medium[END_REF][START_REF] Haskell | Dispersion of surface waves on multilayered media[END_REF]. As shown for γ = 1 + 2i (Fig. 4a), numerical results for leaky modes are in good agreement with Thomson-Haskell results. Yet for γ = 1 + 0.5i (Fig. 4b), some leaky modes cannot be computed by the SAFE-PML method: as already explained in Sec. 2.5, the argument of γ is not large enough to attenuate the transversal growth of leaky modes. Indeed, the rotation of each continuum allows to discover leaky modes (hidden in the original problem without PML), the number of discovered modes increasing with the rotation angle.

Comparison between constant and parabolic PML functions

Figures 5a and5b show the spectrum obtained with a constant PML function and a parabolic profile respectively. In the latter case, the accuracy of leaky modes is improved. This result is due to the fact that the perfectly matched property satisfied for the continuous problem no longer holds after FE discretization, which yields artificial reflections. Moreover, the discontinuity of the constant PML function at the interface between the physical and the PML regions yields more artificial reflection than continuously increasing absorbing functions. As expected from Eq. ( 15), note that the rotation angles remain identical to the constant γ case, because h, d and γ have been kept constant for this comparison.

Concerning higher order radiation modes (further from the real axis), it can be observed that each continuum deviates from straight lines. This phenomenon is due to the FE approximation, as will be shown in Sec. 3.2.4. With γ constant, the deviation of the discretized continuum of shear waves occurs a little further than for compressional waves and is not shown in Fig. 5a.

Note that with a parabolic profile, the deviation of radiation modes from straight lines occurs at a significantly closer distance from the real axis. As shown in Sec. 3.1, choosing an increasing absorbing profile affects the mode- shapes of radiation modes. In particular, the transverse oscillations are narrowed near the PML end, which the FE discretization may no longer approximate.

PML thickness effects

Figure 6 gives the spectrum computed by increasing the PML thickness to h = 1.8a, both for the constant and the parabolic profiles.

Concerning radiation modes, it can be observed that increasing the PML thickness yields a denser discretization of continua, as expected from Sec. 3.1 (increasing h increases L). Note that the deviation limit from straight lines is left unchanged by the parameter h.

As a side remark, increasing the thickness h while leaving γ unchanged may not improve the accuracy of leaky modes when a constant profile is chosen (compare Figs. 5a and6a). While the parabolic function gets smoother as h is increased (in the sense that the discontinuity of the second derivative of γ at x = a +d decreases), the constant function keeps the same discontinuity at the PML interface, yielding artificial reflections that are not reduced by increasing h.

Let us emphasize that the thickness of the PML should not be taken too small, otherwise spurious reflections at the end of the layer will pollute the results. For a given PML (a given h and a given absorbing profile), the reflection coefficients associated to these spurious reflections (which can be easily computed analytically) only depend on λ. Their magnitude are roughly speaking exponentially decreasing when λ moves away from the continua. On the contrary, computing a leaky mode which is located very close to the continua requires a priori a larger h (together with an appropriate discretization).

Mesh size effects

Figure 7 gives the spectrum obtained by refining the FE mesh to ∆x = 0.02a. Two effects can be observed. First, refining the FE mesh tends to restore the perfectly matched property: the accuracy of leaky modes increases (compare Figs. 6 with 7). Second, the deviation occurs at greater distance from the real axis (the continua are also better approximated). The deviation phenomenon is indeed essentially due to the FE discretization, unable to approximate higher order modes as in conventional eigenvalue FE problems.

PML position effects

Figure 8 shows the influence of increasing the distance d between the core and the PML. As d increases, it can be observed that the deviation from straight lines approaches the real axis and that the computation of leaky modes is deteriorated. For instance for γ constant and d = 3a, only two leaky modes are obtained, the other ones being polluted by sub-branches emanating from continua.

Regarding source PML problems, this phenomenon is unusual. Yet for scalar modal problems [START_REF] Goursaud | étude mathématique et numérique de guides d'ondes ouverts non uniformes, par approche modale (mathematical and numerical study of non uniform open waveguides, modal approach[END_REF][START_REF] Kim | The computation of resonances in open systems using a perfectly matched layer[END_REF], it has been mathematically explained by the increase with d of the norm of the resolvent of the eigenproblem, due to the non-hermitian character of matrices involved in the discretized problem. From a physical point of view, setting the PML too far from the core allows leaky modes to exponentially grow in the transverse direction before entering the PML, which can deteriorate their calculation. Therefore, a general requirement for the computation of leaky modes is that the PML interface should be as close as possible to the core.

As observed in Fig. 8b, results obtained from a parabolic function are less affected by the parameter d than with a constant profile. But for d = 3a, Fig. 8d shows again the pollution by sub-branches emanating from continua. It should also be noticed that the rotation angle of continua (measured near the real axis) is significantly decreased compared to previous figures. This is also well explained by Eq. ( 15): increasing d decreases |2 arg(d + h γ)|.

Computation of dispersion curves

This section focuses on the computation of dispersion curves of leaky and trapped modes. In order to assess the efficiency of the SAFE-PML method, three test cases taken from the literature are considered. Results are compared with analytical solutions obtained from a Thomson-Haskell matrix approach.

Modal filtering of radiation modes

As shown in Sec. 3, the SAFE-PML method does not only provide trapped and leaky modes, but also many radiation modes related to the rotated continuous spectrum of the initial physical problem (see Sec. 3.1). These numerous modes are non intrinsic to the physics because they mainly resonate inside the PML. It is desirable to get rid of them in the visualization of dispersion curves. A post-processing step must hence be applied, consisting in identifying and removing modes resonating mainly in the artificial layers.

Our filtering process is based on an adapted notion of kinetic energy. By analytical continuation into the complex coordinate x, the averaged kinetic energy in the region [x 1 ; x 2 ] can be defined as E = 1 2 x2 x1 ρṽ T ṽd x, where ṽ = d ũ/dt is the velocity vector, bars denote time averaging over one period.

From the change of variable [START_REF] Vogt | The scattering of guided waves in partly embedded cylindrical structures[END_REF], the relations v = -iωu and v T v = ω 2 2 u * u, the kinetic energy expression becomes:

E = ω 2 4 x 2 x 1 ρu * uγdx. (16)
The superscript * refers to the complex conjugate transpose.

Let us define the total kinetic energy E tot from Eq. ( 16) with bounds x 1 and x 2 corresponding to the whole crosssection (including PMLs). From the interpolation of the displacement vector [START_REF] Alsop | Reflection and transmission of inhomogeneous waves with particular application to Rayleigh waves[END_REF], the FE expression of E tot is:

E tot = ω 2 4 U * MU. (17)
Without PML (γ=1), this expression coincides with the results standardly found for closed waveguide problems [START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF][START_REF] Treyssède | Mode propagation in curved waveguides and scattering by inhomogeneities: Application to the elastodynamics of helical structures[END_REF]. With PMLs, M is complex and E tot is no longer real. We also define the PML kinetic energy E PML , obtained from Eq. ( 16) with bounds x 1 and x 2 restricting the integration to the PML region only. The filtering criterion used for our tests is based on the ratio of kinetic energy in the PML region over the kinetic energy in the whole domain. Physical modes are then identified if this criterion is smaller than a user-defined value ρ max :

|E PML | |E tot | < ρ max . (18)
Note that since the kinetic energy is complex, the modulus of kinetic energy is used.

Results

Numerical tests are realized for three multilayer waveguides, shown in Fig. 3 and taken from the work of Lowe [START_REF] Lowe | Plate waves for the NDT of diffusion bonded titanium[END_REF][START_REF] Lowe | Comparison of the modal porperties of a stiff layer embedded in a solid medium with the minima of the plane-wave reflection coefficient[END_REF][START_REF] Lowe | The influence of the modal properties of a stiff layer embedded in a solid medium on the field generated in the layer by a finite-sized transducer[END_REF]. The first one has already been considered in Sec. 3 and only involves leaky modes. The second test case is an Epoxy layer of 100 µm on an Aluminum half-space, for which the first two modes are trapped. The third one corresponds to a thin Alpha case layer of 50 µm on a Titanium half-space. This last case is of particular interest because the first mode is trapped in a low frequency range and becomes leaky mode at higher frequencies. Material characteristics are given in Table 1.

The thickness of PML layers is set to h = 0.9a. Following the suggestion of Sec. 3.2.5, the PML interface is set close to the core (d = 0.1a) in order to reduce the effects of the exponential growth of leaky modes. Note that it could be possible to use d = 0 instead. However a small but non zero value is preferred as it yields a more general configuration that will typically occur when extending the method to three-dimensional waveguides of arbitrary crosssection. A Dirichlet condition is chosen at the PML end (zero displacement in both directions). Numerical tests, not presented here for conciseness, have shown that Neumann or mixed boundary conditions yield negligible differences on leaky and trapped modes (such modes mainly resonate outside the PML). Finite elements are quadratic (threenode elements) and sized in order to satisfy the criterion ∆x λ s /5, where λ s = c s / f max denotes the minimum shear wavelength. SAFE-PML dispersion curves are presented after modal filtering and compared with the reference solutions obtained from the Thomson-Haskell matrix method. The user-defined filtering threshold ρ max is specified for each test. Note that a proper choice of ρ max depends on the absorbing profile, since the PML kinetic energy depends on γ(x).

Let us first consider the first test case (Ti-Alpha-Ti). Figure 9 compares the phase velocity and attenuation curves computed with γ constant and γ parabolic for γ = 1 + 4i. In order to reduce the number of eigenvalues to solve, the computation has been centered around the wavenumber ω/c l , where c l denotes the longitudinal velocity in the surrounding medium (this partly avoids the computation of opposite-going modes). 50 modes have been solved.

Figure 9 shows that the parabolic function gives better results than the constant function, as can be observed in the high frequency range. This confirms the trends already found in Sec. 3.2. Therefore in the remainder, all PML results will be obtained with γ parabolic. It can be observed that numerical results obtained with γ constant yield additional modes in a higher frequency range. Though unfiltered by the criterion (18), these modes are not of leaky type. In order to illustrate the efficiency of the filtering method based on Eq. ( 18), Fig. 10 depicts the numerical results without filtering (ρ max =1). The radiation modes yield many dispersion curves and prevent a proper visualization of leaky mode curves. Another benefit of the filtering process is that leaky modes can be filtered out in the low-frequency range, indicating waves insufficiently absorbed by the PML.

Figure 11 represents the phase velocity as a function of frequency for the Epoxy-Al waveguide. 30 modes have been solved, centered around ω/c s . One focuses on the first two trapped modes (k ′ l/s = 0, k ′′ l/s > 0, k ′′ = 0), represented by continuous lines in figures. Numerical results also show other modes, of less interest here (these modes are leaky modes of high axial attenuation or higher order trapped modes).

Note that in the complex plane (Reλ, Imλ), trapped modes in the original problem (without PML) would be located on the real axis. Hence, trapped modes do not need to be discovered by spectrum PML rotation (increasing the imaginary part of γ is useless). For a proper computation of trapped modes, a first idea is to increase the value of d in order to let the trapped modes naturally decrease in the transverse direction. Yet increasing d could lead to numerical problems as previously mentioned in Sec. 3. An interesting alternative technique consist in increasing Reγ instead (see Sec. 2.5), which is equivalent to stretch the transverse direction and therefore to increase the PML thickness. Comparing the results obtained with γ = 4 + i and 20 + i confirms that increasing Reγ give more accurate Figure 12a represents the attenuation as a function of frequency for the Alpha-Ti bilayer with γ = 10 + 10i. 20 modes have been computed around ω/c l . In this test case, one focuses on the computation of the first mode (the less attenuated one). As shown in Fig. 12b, this mode is trapped up to a dimensionless frequency close to 4 and then becomes leaky beyond. Therefore, the real part of γ is chosen as high as the imaginary part in this example. Good agreement with the analytical solution is found.

From Fig. 12a, it is observed that the second mode is filtered out at low frequencies, due to insufficient attenuation by the PML. Besides, the highest mode is not accurately computed. Numerical tests (not shown for conciseness) have shown that setting γ to a lower value, γ = 4 + 4i, gives more satisfying results for higher order modes. Such a value is yet less suitable for the first mode.

Conclusion

A SAFE-PML technique has been proposed to compute the elastic propagation modes of open waveguides. The method yields three kinds of modes: leaky and trapped modes, which are intrinsic to the physics, and radiation modes, which are located along two branches corresponding to a continuum of compressional waves and a continuum of shear waves respectively. The PML rotates both branches from the real axis, with a rotation angle approximately equal to -2 arg L where L is the so-called complex thickness of the semi-infinite layer. The truncation of the PML discretizes both continua and also controls the modal density along each branch (increasing | L| yields a denser discretization). Furthermore, the rotation angle must be large enough in order to discover the leaky modes of interest.

In this study, parabolic absorbing functions have shown to yield more accurate results than constant functions. Besides, the PML interface should be set close to the core in order to avoid the deteriorating effects of the exponential growth of leaky modes. A comparison with analytical dispersion curves obtained from a Thomson-Haskell approach has been conducted. The SAFE-PML method gives satisfying results both for leaky and trapped modes. An efficient filtering criterion based on the ratio of kinetic energy inside the PML over the whole domain has been proposed in order to eliminate radiation modes in the visualization of dispersion curves.

The main drawback of the approach is the optimization of PML parameters (inherent to any PML method). For a given PML thickness and a given mesh size, the PML absorbing function cannot be arbitrarily large due to the FE approximation. However, relatively thin PML layers allow to obtain accurate results.

Figure 1 :

 1 Figure 1: (a) Example of open waveguide, (b) closure by a truncated PML.

Figure 2 :

 2 Figure 2: (a) Analytical spectrum (ω = 17.68, c l = 1.706, c s = 0.909), (b) real part of ϕ(5) l (x) calculated with γ constant (continuous line) and γ parabolic (dashed line). PML parameters are: d = 0.1, h = 0.9 and γ = 1 + 2i.

Figure 3 :

 3 Figure 3: Description of test cases, (a) Titanium-Alpha-Titanium, (b) Epoxy-Aluminium and (c) Alpha-Titanium waveguides.

Figure 4 :

 4 Figure 4: Spectrum of Ti-alpha-Ti waveguide computed with γ constant, d = 0.1a, h = 0.9a, ∆x = 0.05a (Ω = 17.68) for: (a) γ = 1 + 2i, (b) γ = 1 + 0.5i. Crosses: SAFE-PML results, circles: analytical solutions for leaky modes.

Figure 5 :

 5 Figure 5: Spectrum of Ti-alpha-Ti waveguide computed with d = 0.1a, h = 0.9a, ∆x = 0.05a ( γ = 1 + 2i, Ω = 17.68) for: (a) γ constant, (b) γ parabolic. Crosses: SAFE-PML results, circles: analytical results for leaky modes.

Figure 6 :

 6 Figure 6: Same caption as Fig. 5 but with h = 1.8a.

Figure 7 :

 7 Figure 7: Same caption as Fig. 5 but with h = 1.8a and ∆x = 0.02a.

Figure 8 :

 8 Figure 8: Same caption as Fig. 5 but with: (a) d = a and γ constant, (b) d = a and γ parabolic, (c) d = 3a and γ constant, (d) d = 3a and γ parabolic.

Figure 9 :

 9 Figure 9: Comparison between SAFE-PML (circles: γ constant, crosses: γ parabolic) and analytical (continuous lines) solutions for Ti-Alpha-Ti waveguide, (a) phase velocity and (b) attenuation. PML parameters are: d = 0.1a, h = 0.9a, γ = 1 + 4i, ρ max = 0.8 for γ constant and ρ max = 0.9 for γ parabolic (∆x = 0.05a).

Figure 10 :

 10 Figure 10: Same caption as Fig. 9, but without filter (ρ max = 1, γ parabolic only).

Figure 11 :

 11 Figure 11: Comparison between SAFE-PML (crosses) and analytical (continuous lines) phase velocities for Epoxy-Al waveguide for (a) γ = 4 + i (ρ max = 0.6), (b) γ = 20 + i (ρ max = 0.9). PML parameters are: d = 0.1a, h = 0.9a, parabolic function (∆x = 0.2a).

Figure 12 :

 12 Figure 12: (a) Comparison between SAFE-PML (crosses) and analytical (continuous lines) attenuation for Alpha-Ti waveguide, (b) zoom on the first mode. PML parameters are: d = 0.1a, h = 0.9a, γ = 10 + 10i, ρ max = 0.8, parabolic function (∆x = 0.1a).

Table 1 :

 1 Material characteristics

	Titanium	6060	3230	4460
	Alpha	6666	3553	4460
	Aluminium	6370	3170	2700
	Epoxy	2610	1100	1170

Acknowledgements

This work has been partially funded by the French Agence Nationale de la Recherche in frame of the Programme Inter Carnot Fraunhofer PICF of the Franco-German FilameNDT project.