Vahan Baronian

Anne-Sophie Bonnet-Ben Dhia

Sonia Fliss
email: sonia.fliss@ensta-paristech.fr

Antoine Tonnoir

Iterative methods for scattering problems in isotropic or anisotropic elastic waveguides

Keywords: elastic waveguide, diffraction, modal expansion, domain decomposition method, iterative methods

We consider the time-harmonic problem of the diffraction of an incident propagative mode by a localized defect, in an infinite elastic waveguide. We propose several iterative algorithms to compute an approximate solution of the problem, using a classical finite element discretization in a small area around the perturbation, and a modal expansion in the unbounded straight parts of the guide. Each algorithm can be related to a socalled domain decomposition method, with an overlap between the domains. Specific transmission conditions are used, so that at each step of the algorithm only the sparse finite element matrix has to be inverted, the modal expansion being obtained by a simple projection, using a bi-orthogonality relation. The benefit of using an overlap between the finite element domain and the modal domain is emphasized. An original choice of transmission conditions is proposed which enhances the effect of the overlap and allows us to handle arbitrary anisotropic materials. As a by-product, we derive transparent boundary conditions for an arbitrary anisotropic waveguide. The transparency of these new boundary conditions is checked for twoand three-dimensional anisotropic waveguides. Finally, in the isotropic case, numerical validation for twoand three-dimensional waveguides illustrates the efficiency of the new approach, compared to other existing methods, in terms of number of iterations and CPU time.

Introduction

The development of non destructive testing techniques using ultrasonic guided waves (see [START_REF] Lowe | Long Range Guided Wave Inspection Usage Current Commercial Capabilities and Research Directions[END_REF] and the references herein) motivates the improvement of existing numerical methods of simulation. In particular, efficient methods are required to compute the scattering of guided waves by arbitrary defects in elastic waveguides. Classically, the waveguide is supposed to be infinite and perfectly uniform, except in a bounded area containing the defect. A natural objective is to reduce the finite element computations to a region as close as possible to this perturbed area. The difficulty is then to handle the artificial boundaries of the finite element domain in order to avoid spurious reflections. This is an old problem [START_REF] Givoli | Numerical Methods for Problems in Infinite Domains[END_REF] which has been satisfactorily solved in case of scalar equations, but still raises open questions for vectorial equations arising in electromagnetic or elastic waveguides. For scalar problems, two classes of methods can be used. Let us explain for each of them what are the specific issues when considering elastic waveguides.

A first class of methods consists in putting on each side of the perturbed area a perfectly matched absorbing layer (PML), so that the computed diffracted field almost vanishes at the end of the layer. This technique is easy to integrate in any Finite Element codes (no specific implementation is needed) and leads us to solve a classical sparse linear system. Unfortunately, it is well-known that PMLs do not work in elastic waveguides [START_REF] Skelton | Guided elastic waves and perfectly matched layers[END_REF] because of the existence in some range of frequencies of backward modes, whose group and phase velocities are of opposite signs. In such configuration, the PMLs do not select the correct outgoing solution. A remedy has been proposed and analyzed in [START_REF] Bonnet-Ben Dhia | On the use of perfectly matched layers in the presence of long or backward guided elastic waves[END_REF] where the physical solution is reconstructed a posteriori by combining several wrong fields computed with PMLs. An alternative consists in using adiabatic viscoelastic absorbing layers [START_REF] Drozdz | Efficient Numerical Modelling of Absorbing Regions for Boundaries Of Guided Waves Problems[END_REF] which are not perfectly matched and need to be sufficiently large to avoid spurious reflections. The main drawback of this approach is then its computational cost. Also let us point out that absorbing layer techniques (perfectly matched or not) require a fine adjustment of some parameters, which may limit their systematic use in an industrial context. Let us finally mention a more recent method for elastic waveguides based on Hardy space infinite elements [START_REF] Halla | Hardy space infinite elements for time-harmonic two-dimensional elastic waveguide problems[END_REF] whose development is still in progress.

A second possibility consists in using the modal decomposition of the field outside the perturbed area to derive transparent boundary conditions on the artificial boundaries of the finite element domain. The advantage is that such conditions are exact (and with an exponentially small error at the discrete level if enough modes are kept in the modal expansion). Different ways to implement such conditions have been proposed in the literature. In [START_REF] Benmeddour | Numerical modeling of guided wave interaction with non-axisymmetric cracks in elastic cylinders[END_REF][START_REF] Treyssède | Mode propagation in curved waveguides and scattering by inhomogeneities: application to the elastodynamics of helical structures[END_REF], a formulation involving both finite elements and modal unknowns is derived, while only one type of unknowns is generally kept, modal unknowns in [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF][START_REF] Gravenkamp | The simulation of Lamb waves in a cracked plate using the scaled boundary finite element method[END_REF] and finite elements unknowns in [START_REF] Baronian | Transparent boundary conditions for the harmonic diffraction problem in an elastic waveguide[END_REF]. For the latter, the difficulty is due to a lack of orthogonality of the displacement fields associated with elastic modes. As a consequence, it is not possible to obtain a diagonal expression of the natural Dirichlet to Neumann operator (relating the normal stress to the displacement) in elastic waveguides. In the isotropic case, an alternative has been proposed in [START_REF] Baronian | Transparent boundary conditions for the harmonic diffraction problem in an elastic waveguide[END_REF] (see [START_REF] Baronian | Couplage des méthodes modale et éléments finis pour la diffraction des ondes élastiques guidées[END_REF] for more details and [START_REF] Baronian | Simulation of non-destructive inspections and acoustic emission measurements involving guided waves[END_REF][START_REF] Baronian | Hybrid modal/FE simulation of guided waves inspections[END_REF] for applications), where the authors derived transparent boundary conditions relating hybrid displacement/stress vectors. This work is based on a bi-orthogonality relation, mixing displacement and stress components (which has been derived first by Fraser in 2D [START_REF] Fraser | Orthogonality relation for the Rayleigh-Lamb modes of vibration of a plate[END_REF] and then extended to the 3D case [START_REF] Gregory | A note on bi-orthogonality relations for elastic cylinders of general cross section[END_REF]). Let us point out that a scalar Lagrange multiplier has to be introduced on the artificial boundaries, because these transparent conditions are not naturally compatible with the variational formulation in the perturbed area. The method is very accurate but requires a specific implementation and leads to a partially dense linear system. Such system can be difficult to invert, in particular for elastic 3D configurations. Moreover, this approach cannot be used in the general anisotropic case.

In our work, we intend to gather advantages of both classes of methods. In other words, we would like to design a method using transparent boundary conditions based on modal expansions and such that only simple and sparse systems have to be inverted. A natural idea is to use iterative algorithms instead of direct ones to solve a system involving transparent boundary conditions. This framework has been already applied to several problems set in unbounded domains [START_REF] Belgacem | On the Schwarz algorithms for elliptic exterior boundary value problems[END_REF][START_REF] Gmati | Numerical study of some iterative solvers for acoustics in unbounded domains[END_REF][START_REF] Belgacem | Convergence bounds of GMRES with Schwarz preconditioner for the scattering problem[END_REF]. One of the main features is that the system to invert (or equivalently the preconditioner) is chosen as a sparse part of the complete system. As a consequence, the dense part of the matrix coming from the non-local transparent condition is involved only in the matrix-vector product step. It is also instructive to relate these iterative algorithms to domain decomposition methods. The specificity here is that only two subdomains are introduced, a bounded one and an unbounded one where the equation can be explicitly solved, using an analytical representation (a modal expansion or an integral representation for instance). Different algorithms are derived by modifying the transmission conditions between the two subdomains, which can overlap or not. It is now well known that algorithms of this type do not converge in general for time harmonic wave equations [START_REF] Belgacem | Convergence bounds of GMRES with Schwarz preconditioner for the scattering problem[END_REF] but they can be used to design preconditioned Krylov methods like GMRES [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving non-symmetric linear systems[END_REF]. The main criterion to discriminate between different algorithms is the rate of convergence of the associated GMRES algorithms.

In the present paper, we want to adapt these ideas to the case of elastic waveguides, which has not been considered in the literature. Combining such point of view with the ideas of [START_REF] Baronian | Transparent boundary conditions for the harmonic diffraction problem in an elastic waveguide[END_REF], we derive new families of modal transparent boundary conditions. The benefit of adding an overlap is emphasized. Besides classical effects (for instance improvement of the convergence rate of associated iterative algorithms), it allows the construction of a particularly efficient transparent boundary condition described in Section 4. Last but not least, this condition can be used for anisotropic waveguides. This is the main contribution of the paper.

The paper is organized as follows. In Section 2, the main notions concerning elastic modes in the isotropic case are summarized. Then, we introduce a well-posed boundary value problem in a semi-infinite straight waveguide, with a modal expansion of its solution. Finally, the scattering problem we are interested in is defined. In Section 3, domain decomposition algorithms are introduced. The subdomains are two semi-infinite straight waveguides and a bounded domain containing the defect. We explain how to choose transmission conditions between the subproblems such that the iterative algorithm is consistent and computable (which means that the solution at each step can be computed with a standard finite element code). However, this first approach is not completely satisfactory: it is restricted to the case of isotropic materials, and the existence of spurious resonances of the overlapping area may affect the numerical computation. To end this section, we explain the link between these domain decomposition methods and modal transparent boundary conditions. This provides a new point of view for the choice of transmission conditions. Up to now, only local conditions have been used, and the new idea, which is investigated is Section 4, is to consider non-local modal transmission conditions. We propose a particular choice which completely removes the resonances of the overlapping area. Moreover, this method can be extended to anisotropic waveguides. In Section 5, the convergence of the Schwarz algorithm is analyzed and we compare the different formulations we proposed, for a simple model case. It is shown that the behaviour of the algorithm for high order modes is highly improved by means of the overlap between the subdomains. However it appears that the convergence cannot be ensured in general, so we are lead to consider the previous domain decomposition method as a preconditioner for a GMRES algorithm. Finally, two-and three-dimensional numerical results are presented in Section 6. The superiority of the algorithm of Section 4, in terms of convergence rate, is pointed out.

Modal analysis and definition of the scattering problem

Modes classification and Fraser bi-orthogonality relation

Let us consider a homogeneous isotropic elastic waveguide of section S in x s plane -S is a bounded domain of R (x s = x) in the 2D case or R 2 (x s = (x, y)) in the 3D case-and of axis z, which occupies the domain S × R of R 2 or R 3 . The density of the material denoted by ρ and the Lamé's coefficients denoted by λ and µ are supposed to be independent of z, but they may depend on x s . Finally, we suppose that the waveguide has a stress-free boundary ∂S × R. In time harmonic regime (of pulsation ω > 0), the propagation in the waveguide is modeled by the following equations

-div σ(u) -ω 2 ρ u = 0 in S × R, σ(u) • ν = 0 on ∂S × R, (1)
where u represents the displacement field (u = (u x , u z ) in the 2D case and u = (u x , u y , u z ) in the 3D case) and σ(u) the stress tensor which is related to the strain tensor ε(u) = 1/2 (∇u + ∇ T u) by Hooke's law

σ(u) = λ div (u) I + 2 µ ε(u). (2)
The vector ν denotes the outward unitary normal to ∂S × R.

We recall in this section the convenient formalism introduced in [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF][START_REF] Pagneux | Scattering matrix properties with evanescent modes for waveguides in fluids and solids[END_REF][START_REF] Pagneux | Lamb wave propagation in elastic waveguides with variable thickness[END_REF] for the 2D case and extended in [START_REF] Baronian | Couplage des méthodes modale et éléments finis pour la diffraction des ondes élastiques guidées[END_REF][START_REF] Baronian | Transparent boundary conditions for the harmonic diffraction problem in an elastic waveguide[END_REF] to the 3D case. First we denote by u s the transverse part of the displacement field u, t s the transverse part of the axial stress σ(u) • e z , σ s the transverse part of the stress tensor and ε s the transverse part of the strain tensor: in 2D:

u s = u x , t s = σ zx , σ s = σ xx , ε s = ε xx , in 3D: u s = u x u y , t s = σ zx σ zy , σ s = σ xx σ xy σ yx σ yy , ε s = ε xx ε xy ε yx ε yy .
Finally we set t z = -σ zz and we define the two vectors

X = t s u z and Y = u s t z ,
which are hybrid in the sense that they mix both displacement field and stress field components.

The modes of the waveguide are solutions of (1) with separated variables of the standard form

u(x s , z) = U(x s ) e iβz ,
where the function U and the complex number β are respectively the field profile and the propagation constant of the mode. It can be proved that, if ω is not a cut-off frequency, there is a discrete family of rightgoing modes (β k , U + k (x s )), k = 1, 2, • • • , and a corresponding family of leftgoing modes (-β k , U - k (x s )). Among these modes, one can distinguish a finite number of propagative modes such that β k ∈ R and an infinite set of (possibly oscillating) evanescent modes such that β k ∈ C\R and Im(β k ) tends to +∞ [START_REF] Nasarov | The Mandelstam Energy Radiation Conditions and the UmovPoynting Vector in Elastic Waveguides[END_REF]. For the propagative modes, the direction of propagation is given by the sign of the group velocity, and for the evanescent modes by the sign of the imaginary part of β k : then, the propagative rightgoing modes satisfy ∂β k ∂ω > 0 (for a time dependance of the form e -iωt ), and the evanescent rightgoing modes Im(β k ) > 0.

In terms of X and Y variables, for simple symmetry reasons, the rightgoing and leftgoing modes are of the form (

β k , X k (x s ), Y k (x s )) and (-β k , -X k (x s ), Y k (x s )).
As mentioned in the introduction, a difficulty comes from the fact that the modal transverse displacement fields U ± k are not orthogonal in L 2 (S). However, the modes satisfy a bi-orthogonality relation, known as Fraser's biorthogonality relation (see [START_REF] Fraser | Orthogonality relation for the Rayleigh-Lamb modes of vibration of a plate[END_REF][START_REF] Gregory | A note on bi-orthogonality relations for elastic cylinders of general cross section[END_REF]), that takes the simple following expression in terms of X and Y variables:

S X (x s ) Y k (x s ) dx s = 0 for = k.
Except at the cutoff frequencies, the modes can be normalized so that we have finally:

S X (x s ) Y k (x s ) dx s = δ k . (3)
In the following, we suppose that ω is not a cutoff frequency.

The half-guide solution

Let us denote by Ω + a the semi-infinite straight waveguide defined by Ω + a = S × {z > a}. It will be useful in what follows to have a simple representation of the solution of the equations [START_REF] Auld | Acoustic fields and waves in solids[END_REF] in Ω + a for some given data on the boundary Γ + a = S × {z = a}.

More precisely, we will look for a solution which is a superposition of rightgoing modes, that we call an outgoing solution (since rightgoing modes propagate from Γ + a to infinity). Such explicit solution is not available for classical boundary conditions, like Dirichlet condition on u for instance. Fraser relation (3) leads us to consider instead that X or Y is given on Γ + a .

Suppose for instance that Y is known on Γ + a -in other words u s and t z = -σ zz -and let us denote by Y + a this data. Then, assuming the completeness of the elastic modes, we can expand Y + a on the Y k :

Y + a = k>0 α k Y k ,
and the value of α k can be directly derived by using Fraser relation (3):

α k = Γ + a Y + a (x s ) X k (x s ) dx s .
Finally,

u + (x s , z) = k>0 Γ + a Y + a (x s ) X k (x s ) dx s U + k (x s ) e iβ k (z-a) , (4)
is the outgoing solution of

     -div σ(u) -ω 2 ρ u = 0 in Ω + a , σ(u) • ν = 0 on ∂Ω + a , Y = Y + a on Γ + a . (5)
Remark 2.1. The proper mathematical sense of the modal series expansion (4) is not clear, because the sequence (U + k ) k>0 is not a Hilbert basis of L 2 (S) 2 in 2D, or L 2 (S) 3 in 3D. Nevertheless, we can point out that for z > a, since Im(β k ) tends to +∞ with k, the sequence e iβ k (z-a) tends to 0 exponentially fast. We will take advantage of this property in Section 5 .

Proceeding in the same way and setting

Ω - a = S × {z < -a} and Γ - a = S × {z = -a},
we prove that

u -(x s , z) = k>0 Γ - a Y - a (x s ) X k (x s ) dx s U - k (x s ) e -iβ k (z+a) , (6)
is the outgoing solution of

     -div σ(u) -ω 2 ρ u = 0 in Ω - a , σ(u) • ν = 0 on ∂Ω - a , Y = Y - a on Γ - a . (7)
Remark 2.2. A similar expansion holds if we impose a condition of the type X = X ± a on Γ ± a , where X ± a is decomposed on the family Y k (x s ). All the results obtained in the sequel with Y data on Γ ± a can be easily transposed to the case of X data. For the sake of conciseness, we consider here only the case of a Y data. The notion of outgoing waves being introduced, we can define the scattering problem we will study in the sequel. Let us denote by Ω a locally perturbed waveguide defined by (see Figure 1) :

The scattering problem

⌦ a ⌦ + a ⌦ a z = a z = a z + a a S S
Ω = Ω - a ∪ Ω a ∪ Ω + a
where Ω a = Ω ∩ {|z| ≤ a}.

We suppose that Ω a is a bounded domain and Ω is a connected domain. The scattering problem consists in finding the outgoing solution u of the following equations

-div σ(u) -ω 2 ρ u = f in Ω, σ(u) • ν = 0 on ∂S × R, (8)
where f ∈ L 2 (Ω) 2 in 2D and f ∈ L 2 (Ω) 3 in 3D is a compactly supported source term whose support is included in Ω a . In accordance to Section 2.2, we say that u is outgoing if it expands on the rightgoing modes in Ω + a , and on the leftgoing modes in Ω - a :

u =    k>0 α + k U + k (x s ) e iβ k (z-a) in Ω + a , k>0 α - k U - k (x s ) e -iβ k (z+a) in Ω - a , (9)
where (α + k ) k and (α - k ) k are two sequences of complex numbers. We define the scattering problem by equations ( 8) and ( 9). We suppose that ω is such that the scattering problem is well-posed.

In the rest of the paper, the outgoing behaviour of the solution in Ω - a and Ω + a will be implicitly imposed for each problem set in these unbounded domains. So, we will refer to the scattering problem simply by (8).

A first domain decomposition approach

The Schwarz algorithm

The first domain decomposition method was introduced by Schwarz [START_REF] Schwarz | Ueber einige Abbildungsaufgaben[END_REF] to solve a Laplace problem with non homogeneous Dirichlet boundary conditions, in a domain Ω which is the union of a disk Ω D and a rectangle Ω R that overlap each other: Ω D ∩ Ω R = ∅. Since the Laplace equation is analytically solvable in a rectangle and in a disk, the key idea is to solve a Laplace problem iteratively in each subdomain. Yet, to implement this process, boundary conditions must be chosen on the part of the boundary of the subdomains (the disk and the rectangle) not included in the boundary of the whole domain ∂Ω. This choice defines the iterative algorithm and determines its rate of convergence. In the 80's, this approach was extended to general elliptic equations set in an arbitrary geometry: in particular, it is well known that the domain decomposition techniques have a geometrical rate of convergence which increases with the size of the overlap [START_REF] Lions | On the Schwarz alternating method. I[END_REF] for this class of problems.

The application to time harmonic scalar wave equations raises specific difficulties [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique[END_REF][START_REF] Collino | Domain decomposition method for harmonic wave propagation: a general presentation[END_REF][START_REF] Gander | Optimized Schwarz methods with overlap for the Helmholtz equation[END_REF]. In that case and for a bounded domain, it has been proved that the Schwarz algorithm converges if Robin type conditions are imposed on the boundary of the subdomains. Yet, if one of the subdomains is unbounded, ensuring the convergence for the Helmholtz equation is not possible [START_REF] Boudendir | Techniques de décomposition de domaine et méthode d'équations intégrales[END_REF][START_REF] Belgacem | Convergence bounds of GMRES with Schwarz preconditioner for the scattering problem[END_REF].

In our case, the domain Ω is decomposed into three subdomains Ω - a , Ω b and Ω + a , where Ω ± a are defined as previously and Ω b = Ω∩{|z| < b}, with b-a = l > 0 the size of the overlap. We denote by Γ ± b = S×{z = ±b} the artificial boundaries of Ω b , and by B ± l the overlapping areas between Ω b and Ω ± a (see Figure 2). Based on these three domains Ω - a , Ω b and Ω + a , we consider the following class of iterative algorithms where u b

(n) ⌦ a ⌦ a z a b + b ⌦ b B + l B l ⌦ + a + a
     -div σ(u - (n) ) -ω 2 ρ u - (n) = 0 in Ω - a , σ(u - (n) ) • ν = 0 on ∂Ω ∩ ∂Ω - a , B -(u - (n) ) = B -(u b (n-1) ) on Γ - a ,      -div σ(u + (n) ) -ω 2 ρ u + (n) = 0 in Ω + a , σ(u + (n) ) • ν = 0 on ∂Ω ∩ ∂Ω + a , B + (u + (n) ) = B + (u b (n-1) ) on Γ + a ,      -div σ(u b (n) ) -ω 2 ρ u b (n) = f in Ω b , σ(u b (n) ) • ν = 0 on ∂Ω ∩ ∂Ω b , B b (u b (n) ) = B b (u ± (n) ) on Γ ± b . (10)
Let us emphasize that the three subproblems have to be solved in the same order as they are written in [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique[END_REF]. At the iteration n, u b (n-1) is known and provides the boundary data for the subproblems in Ω ± a . Once these subproblems have been solved, we can use u ±

(n) as boundary data for the subproblem in Ω b and get u b (n) .

Here B b , B + and B -are the operators which define the boundary conditions of the subproblems and whose choices will be discussed below. These choices should ensure that algorithm [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique[END_REF] is well-posed and consistent in the following sense: Definition 3.1. The iterative process [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique[END_REF] with the associated boundary conditions defined by B b , B + and B -is said to be well-posed if the three subproblems are well-posed. It is said to be consistent if for any initial data, in case of convergence, the limits (u -, u b , u + ) which are solutions to

     -div σ(u -) -ω 2 ρ u -= 0 in Ω - a , σ(u -) • ν = 0 on ∂Ω ∩ ∂Ω - a , B -(u -) = B -(u b ) on Γ - a ,      -div σ(u + ) -ω 2 ρ u + = 0 in Ω + a , σ(u + ) • ν = 0 on ∂Ω ∩ ∂Ω + a , B + (u + ) = B + (u b ) on Γ + a ,      -div σ(u b ) -ω 2 ρ u b = f in Ω b , σ(u b ) • ν = 0 on ∂Ω ∩ ∂Ω b , B b (u b ) = B b (u ± ) on Γ ± b , (11)
satisfy the compatibility relations:

u b = u -in B - l , u + in B + l . (12)
Obviously, the convergence of algorithm ( 10) is also required but this issue will be discussed in Section 5.

Concerning the subproblems in the half-guides Ω ± a , it is 'natural', as explained in Section 2.2, to choose the boundary condition B ± (u ± ) = Y ± on Γ ± a . On the other hand, we will choose a boundary condition on Γ ± b which allows us to derive a variational formulation of the subproblem set in Ω b .

Let us discuss the two following typical possibilities:

   B b (u b ) = X b on Γ ± b , B ± (u ± ) = Y ± on Γ ± a , (13)
and

   B b (u b ) = σ(u b ) • ν + α u b on Γ ± b , B ± (u ± ) = Y ± on Γ ± a , (14)
where α is a complex constant.

For these two possibilities, we will now discuss the well-posedness and the consistency of the associated iterative algorithms ( 10) + ( 13), and ( 10) + ( 14) (here and in the sequel, we mean by ( 10) + ( 13) algorithm [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique[END_REF] with boundary conditions ( 13)).

The well-posedness of the algorithms

In both cases, the subproblems set in Ω ± a are well-posed: indeed, given Y b (n-1) on Γ ± a , we get by formulae ( 4) and ( 6):

u ± (n) (x s , z) = k>0 Γ + a Y b (n-1) (x s ) X k (x s ) dx s U ± k (x s ) e ±iβ k (z∓a) , (15)
Therefore, to prove that the algorithm is well-posed, it suffices to prove that the problem in Ω b is well-posed as well. This is the object of the following lemma:

Lemma 3.1. Let f ∈ L 2 (Ω b ) and g ± ∈ L 2 (Γ ± b ). • The problem      -div σ(u b ) -ω 2 ρ u b = f in Ω b , σ(u b ) • ν = 0 on ∂Ω ∩ ∂Ω b , B b (u b ) = g ± on Γ ± b , (16)
where B b is defined by B b (u b ) = X b , is well-posed except for a countable set of frequencies ω. As a consequence, except for these frequencies, the algorithm defined by (10) + ( 13) is well-posed.

• Problem [START_REF] Belgacem | Convergence bounds of GMRES with Schwarz preconditioner for the scattering problem[END_REF], with

B b defined by B b (u b ) = σ(u b ) • ν + αu b
, has a unique solution for any value of α such that Im (α) = 0. As a consequence, if Im (α) = 0, the algorithm defined by (10) + ( 14) is well-posed for all frequencies ω.

Proof: For both boundary conditions, one can derive a variational formulation of problem [START_REF] Belgacem | Convergence bounds of GMRES with Schwarz preconditioner for the scattering problem[END_REF]. Then using Korn's inequality, one can prove that ( 16) is of Fredholm type. Therefore, the uniqueness of the solution is equivalent to the existence of a solution.

For the boundary condition B b (u b ) = X b , the homogeneous problem has a sequence of eigenfrequencies, tending to infinity. As a consequence, the uniqueness holds except when ω is an eigenfrequency.

For the second boundary condition B b (u b ) = σ(u b ) • ν + αu b , let us prove that the homogeneous problem admits as unique solution u b = 0. If the data f and g ± are null, then testing our equation with the complex conjugate of u b and taking the imaginary part of the resulting identity leads to Im (α)

Γ ± b |u b | 2 = 0. ( 17
) If Im(α) = 0, we get u b = 0 on Γ ± b . But, one has also σ(u b ) • ν = 0 on Γ ± b (because σ(u b ) • ν + α u b = 0)
. Finally, by a unique continuation argument [START_REF] Nakamura | Unique Continuation for Elliptic Systems and Crack Determination in Anisotropic Elasticity Contemp[END_REF] we prove that u b = 0 in Ω b .

The consistency of the algorithms

Considering the definition of the consistency, see Definition 3.1, it seems natural to focus on the two differences v ± = u bu ± in B ± l which are solutions of the following problems set in B ± l :

         -div σ(v ± ) -ω 2 ρ v ± = 0 in B ± l , σ(v ± ) • ν = 0 on ∂Ω ∩ ∂B ± l , B b (v ± ) = 0 on Γ ± b , B ± (v ± ) = 0 on Γ ± a . (18)
If the above problems are well-posed, then v ± are necessarily equal to 0, and the algorithm is consistent since the compatibility relations ( 12) are satisfied. In particular, for the boundary conditions defined by ( 13) and ( 14), we prove the following results:

Lemma 3.2.
• Algorithm (10) + ( 13) is consistent if and only if ω is not an eigenfrequency of problem (13) + [START_REF] Gregory | A note on bi-orthogonality relations for elastic cylinders of general cross section[END_REF].

• Algorithm (10) + ( 14) is consistent for all frequencies ω.

Proof: We prove using the same arguments as in the proof of Lemma 3.1 that problem (13) + (18) admits a countable set of eigenfrequencies, and that problem ( 14) + ( 18) has 0 as unique solution. The difficulty is to prove that conversely, if ω is an eigenfrequency of the homogeneous problem ( 13) + ( 18), then algorithm (10) + ( 13) is indeed not well-posed (see Appendix A for the proof).

Let us sum up the properties of the two iterative algorithms we have described:

• the iterative algorithm (10) + ( 13) is well-posed, except for a discrete set of frequencies (see Lemma 3.1), and is consistent, except for a second discrete set of frequencies (see Lemma 3.2).

• the iterative algorithm (10) + ( 14) is well-posed and consistent for all frequencies (see Lemmas 3.1 and 3.2).

Remark 3.3. The eigenfrequencies of problems (14)+( 16) and ( 14)+( 18), associated with algorithm (10)+ ( 14), are complex when α is complex. However some eigenfrequencies can come close to the real axis. After discretization, the conditioning of the problem then gets worse near these frequencies. An appropriate choice of the parameter α can address this difficulty [START_REF] Tonnoir | Conditions transparentes pour la diffraction d'ondes en milieu élastique anisotrope[END_REF].

From multi-domain formulations to transparent boundary conditions

In the previous algorithms, solving the subproblems in the half-guides Ω ± a only requires to project the boundary data Y ± on Γ ± a on the family of outgoing modes (see formula [START_REF] Gmati | Numerical study of some iterative solvers for acoustics in unbounded domains[END_REF]). This provides modal expressions of u ± (n) in the half-guides Ω ± a , that can be plugged in the boundary conditions

B b (u b (n) ) = B b (u ± (n) ) on Γ ± b .
For instance, taking the definition (13) of B b , we get the following relation:

X b (n) (x s , ±b) = k>0 Γ ± a Y b (n-1) (x s , ±a) X k (x s ) dx s X k (x s ) e iβ k (b-a) .
Consequently, we can rewrite the iterative algorithm (10) + (13), without using u ± (n) , as follows:

     -div σ(u b (n) ) -ω 2 ρ u b (n) = f in Ω b , σ(u b (n) ) • ν = 0 on ∂Ω ∩ ∂Ω b , X b (n) ∓ T ±,l Y →X (Y b (n-1) Γ ± a ) = 0 on Γ ± b , (19)
where the operator T ±,l Y →X is defined by: Definition 3.2. The operator that maps a given data Y ± on Γ ± a to X ± = X(u ± ) on Γ ± b , where u + and u -are the outgoing solutions of (5) and [START_REF] Boudendir | Techniques de décomposition de domaine et méthode d'équations intégrales[END_REF], is denoted by T ±,l Y →X . It has the following explicit expression:

T ±,l Y →X (Y ± ) = ± k>0 Γ ± a Y ± (x s ) X k (x s ) dx s e iβ k l X k (x s ), ( 20
)
where l = ba > 0. Similarly, the operator that maps

Y ± on Γ ± a to σ(u ± ) • ν + α u ± on Γ ± b is denoted by T ±,l
Y →R . Algorithm [START_REF] Hayashi | Guided wave simulation and visualisation by a semi-analytical finite element method[END_REF] can be seen as an iterative algorithm to solve the following problem

     -div σ(u b ) -ω 2 ρu b = f in Ω b , σ(u b ) • ν = 0 on ∂Ω ∩ ∂Ω b , X b ∓ T ±,l Y →X (Y b Γ ± a ) = 0 on Γ ± b , (21)
where the last equation is a so-called transparent boundary condition. We can derive in the same way the following transparent boundary condition: 21) was studied in [START_REF] Baronian | Transparent boundary conditions for the harmonic diffraction problem in an elastic waveguide[END_REF] in the case l = 0 (that is without overlapping). In the present paper, we propose extensions of the condition considered in [START_REF] Baronian | Transparent boundary conditions for the harmonic diffraction problem in an elastic waveguide[END_REF]. Let us notice that the behaviour of the modal series expansion of T ±,l Y →X is drastically different for l > 0 and l = 0. Indeed, for l > 0, the term e iβ k l is exponentially decreasing with respect to k. For this reason, the modal series has a normal convergence and one can prove that the operator T ±,l Y →X is regularizing. Clearly, this is not true for l = 0.

σ(u b ) • ν + αu b ∓ T ±,l Y →R (Y b Γ ± a ) = 0 on Γ ± b . Problem (
Let us finally emphasize that the boundary conditions on Γ ± b in problem ( 21) are non-local, which is classical with transparent boundary conditions. This leads to a partially full matrix after discretization. The iterative algorithm [START_REF] Hayashi | Guided wave simulation and visualisation by a semi-analytical finite element method[END_REF] (or equivalently (10) + ( 13)) avoids this difficulty since, at each iteration, the non-local term appears in the right hand side. For this reason, it is not more expensive to consider non-local transmission conditions in the iterative algorithm [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique[END_REF], and this is what we propose in the next Section.

A new modal 'outgoing' transmission condition

In Section 3, we have explained how to design a consistent and well-posed iterative algorithm, thanks to appropriate boundary conditions on the boundaries Γ ± a and Γ ± b . Thanks to Fraser relation (3), the conditions B ± (u ± ) = Y ± provide the modal expansions of the solution in Ω ± a , as explained in Section 2.2. Moreover, the Robin condition B b (u b ) = σ(u b ) + αu b with a complex constant α makes the algorithm consistent since the associated problems in the boxes B ± l (cf problem ( 18)) are always well-posed (see Lemma 3.2).

We will now go further by changing the interior boundary conditions B ± (u ± ). We will see below that considering instead non-local boundary conditions also provides naturally the modal expansions and gathers the following advantages: 1. it avoids spurious effects from the boxes B ± l (in particular the resonances mentioned in Remark 3.3), 2. it has the advantage that it is generalizable to anisotropic waveguides, contrary to the previous algorithms, 3. it has nice convergence properties that will be shown in Section 5 and 6.

Description of the method

In the following, the explanations will be given for the transmission condition on Γ + a . Obviously, the same technique can be applied for the transmission condition on Γ - a .

Let us recall that the transmission conditions on Γ + a and Γ + b must ensure the compatibility u + = u b in B + l (see Definition 3.1). The idea now is to use the modal expansions of both u b and u + to express a non-local modal transmission condition. Due to the geometry of B + l , u b and u + can be decomposed a priori on the rightgoing and leftgoing modes

X(u) Y (u) = k≥0 A + k (u) e iβ k (z-a) X k (x s ) Y k (x s ) + A - k (u) e -iβ k (z-a) -X k (x s ) Y k (x s ) in B + l , (22)
where u = u b or u + , and A ± k (u) are the modal amplitudes of u given, thanks to Fraser's relation (3), by

A + k (u) = 1 2 Γ + a X Y k + Γ + a Y X k and A - k (u) = 1 2 Γ + a Y X k - Γ + a X Y k . (23)
Since u + must be outgoing in Ω + a , we propose to match:

A + k (u b ) = A + k (u + ). ( 24
)
Using formula [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF], this condition is nothing else but a non-local condition involving the data X and Y (or (u s , u z ) and (t s , t z )) taken on Γ + a . It means that rightgoing components of u b and u + coincide, so that u bu + propagates to the left direction. This is why we call it a modal outgoing transmission condition, since we match the outgoing modal amplitudes in the half-guide Ω + a .

Remark 4.1. Using the operator T -,0 Y→X defined in (20) (with l = 0), condition ( 24) can be rewritten as follows:

X b -T -,0 Y→X (Y b ) = X + -T -,0 Y→X (Y + ) on Γ + a (25)
Thus, doing similarly on Γ - a , the boundary conditions we propose to associate to the iterative algorithm (10) are:

• either, similarly to [START_REF] Givoli | Numerical Methods for Problems in Infinite Domains[END_REF]:

   B b (u b ) = X b on Γ ± b , B ± (u ± ) = A ± k (u ± ) k∈N on Γ ± a , (26)
• or, similarly to [START_REF] Belgacem | On the Schwarz algorithms for elliptic exterior boundary value problems[END_REF]:

   B b (u b ) = σ(u b ) • ν + α u b on Γ ± b , B ± (u ± ) = A ± k (u ± ) k∈N on Γ ± a . (27)

The well-posedness and the consistency of the method

Well-posedness of the method. The resolution of the subproblems in Ω ± a in algorithm [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique[END_REF], given the data

B ± (u ± ) = A ± k (u b (n-1)
) on Γ ± a , is obvious since we directly have the modal amplitudes of u ± (n) :

u ± (n) (x s , z) = k>0 A ± k (u b (n-1) ) U ± k (x s ) e ±iβ k (z∓a) .
For the subproblem in Ω b , the well-posedness has been already studied in Section 3.1.1: we know by Lemma 3.1 that for boundary conditions [START_REF] Halla | Hardy space infinite elements for time-harmonic two-dimensional elastic waveguide problems[END_REF] the algorithm is well-posed for almost all frequencies, and for boundary conditions [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF] it is well-posed for all frequencies.

Consistency of the method. Once again, to prove the consistency of the algorithm, if suffices to check that the unique solution of the following problems set in B ± l :

         -div σ(v ± ) -ω 2 ρ v ± = 0 in B ± l , σ(v ± ) • ν = 0 on ∂Ω ∩ ∂B ± l , B b (v ± ) = 0 on Γ ± b , A ± k (v ± ) = 0, ∀k ≥ 0 on Γ ± a , (28)
is the trivial solution. Let us show the result for v + . As we explained, due to the geometry of B + l , v + can be decomposed on the rightgoing and the leftgoing modes. However, the modal condition on Γ + a implies that the rightgoing modal components of v + are equal to zero, so we have

X(v + ) Y (v + ) = k≥0 A - k (v + ) e -iβ k (z-a) X k (x s ) Y k (x s ) in B + l .
This shows that v + is the restriction to B + l of the outgoing solution in the fictitious half-guide

Γ + b × {z < b} of      -div σ(w) -2 ρ w = 0 in Γ + b × {z < b}, σ(w) • ν = 0 on ∂Γ + b × {z < b}, B b (w) = 0 on Γ + b , (29)
Therefore, the consistency of algorithm ( 10) + ( 26) results from the following lemma:

Lemma 4.2. The homogeneous problem (29) with the boundary condition B b (w) = X(w) admits as unique solution w = 0. As a consequence, the algorithm defined by (10) + ( 26) is consistent.

Proof: As in Section 2.2, we get that the solution w of ( 29) is given by:

w(x s , z) = k>0 Γ + b X(w) Y k U - k (x s ) e -iβ k (z-b) ,
which vanishes when X(w) = 0 on Γ + b . Now, for boundary conditions [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF], we have the following result: Proof: Because w is an outgoing solution of this homogeneous problem, it can be proved that its propagative components are equal to 0. This implies that w is exponentially decaying. For α = 0, non trivial solution can exist for exceptional frequencies and are called edge modes [START_REF] Pagneux | Revisiting the edge resonance for Lamb waves in a semi-infinite plate[END_REF]. For Im(α) > 0, the only solution of the above problem is the trivial solution (edge mode never exists), the proof being similar to the proof of Lemma 3.1. We refer the reader to [START_REF] Tonnoir | Conditions transparentes pour la diffraction d'ondes en milieu élastique anisotrope[END_REF] for more details.

Let us sum up the properties of the two iterative algorithms we have described in this section.

• The iterative algorithm ( 10) + ( 26) is well-posed except for a discrete set of frequencies (see Lemma 3.1), and it is consistent for all frequencies (see Lemma 4.2).

• When Im(α) > 0, the iterative algorithm ( 10) + ( 27) is well-posed and consistent for all frequencies (see Lemma 3.1 and Proposition 4.3).

Let us emphasize that the condition set on Γ + a has completely eliminated the spurious effects of the boxes B ± l : the boxes have in some sense disappeared.

Generalization to the anisotropic case

The modal outgoing condition we propose to impose on the interior boundaries Γ ± a is required to be able to compute the amplitudes

A ± k (u b (n-1)
). In the case of a general anisotropy, Fraser's relation [START_REF] Baronian | Hybrid modal/FE simulation of guided waves inspections[END_REF] does not hold anymore, so we cannot use formula [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. However, the general bi-orthogonality relations hold (see [START_REF] Auld | Acoustic fields and waves in solids[END_REF]): for all k, j,

Γ + a U ± k S ± j -S ± k U ± j = 0, and if k = j Γ + a U ± k S ∓ j -S ± k U ∓ j = 0, ( 30
)
where S ± k denotes the normal stress of the mode:

e iβ k z S ± k (x s ) = σ e iβ k z U ± k (x s ) • e z .
These relations express orthogonality relations between the two families of rightgoing and leftgoing modes. Thanks to them, given a data [u b (n-1) , t b (n-1) ] on Γ + a for instance, we can get the modal amplitudes by projection over the two families of modes:

A ± k (u b (n-1) ) = 1 J ± k Γ + a u b (n-1) S ∓ k -t b (n-1) U ∓ k , where J ± k = Γ + a U ± k S ∓ k -S ± k U ∓ k . (31)
Let us mention that J ± k = 0 because we supposed that ω is not a cut-off frequency. Then we deduce that u +

(n) has the following expression in Ω + a :

u + (n) (x s , z) = k>0 1 J + k Γ + a u b (n-1) S - k -t b (n-1) U - k e iβ k (z-a) U + k (x s ).
Similarly to the operator T ±,l Y →R (see Definition 3.2), we can construct the operator T ±,l O→R for the transmission conditions [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF], and derive transparent boundary conditions for anisotropic waveguides. Definition 4.1. The operator that maps a given data [u, t] on Γ ± a to σ(u) • ν + αu on Γ ± b is denoted by T ±,l O→R and it has the following explicit expression:

T ±,l O→R ([u, t]) = k>0 1 J ± k Γ + a u S ∓ k -t U ∓ k e iβ k l S ± k + α U ± k Remark 4.4.
In the isotropic case, formulae [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF] or [START_REF] Royer | Ondes élastiques dans les solides[END_REF] for computing the modal amplitudes are exactly the same. Therefore, the definition of the operators T ±,l O→R could be rewritten for isotropic waveguides as follows:

T ±,l O→R ([u, t]) = k>0 1 2 Γ ± a X Y k + Γ ± a Y X k e iβ k l S ± k + α U ± k .
With these operators, we can formulate the diffraction problem as follows:

     -div σ(u b ) -ω 2 ρu b = f in Ω b , σ(u b ) • ν = 0 on ∂Ω ∩ ∂Ω b , σ(u b ) + αu b ∓ T ±,l O→R ([u b , t b ] Γ ± a ) = 0 on Γ ± b , (32)
where the last equation is a transparent boundary condition.

To illustrate these new conditions, we have solved problem [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving non-symmetric linear systems[END_REF] for artificial 2D and 3D unperturbed anisotropic waveguides. Let us recall the general Hooke's law for anisotropic media:

σ(u) = C ε(u), ( 33
)
where C is the stiffness tensor.

1. 2D example : the size of the section is h = 1mm, the stiffness tensor (using classical Voigt notations) is given by Using triangular Lagrange finite elements of order 4 (characteristic size of the mesh h/20), we have computed the solution for two sizes b ∈ {0.5mm, 1mm} of the computational area Ω b , with the same overlap l = 0.1mm. Figure 3 represents the modulus of the displacement field obtained in these two cases. We also represent the whole solution reconstructed in the half-guides Ω ± a thanks to the modal expansion. 2. 3D example : the section is a rectangle of width L x = 2mm and length L y = 3mm, the stiffness tensor (using classical Voigt notations) is given by the frequency ω = 1 MHz is such that there exists 24 propagative modes, the density ρ = 1.56 kg / m 3 and finally the source term is given by

C =   6
C =        
f = [0 0 exp(-10z 2 -5y 2 -7.5x 2 )] t .
Using quadrangular Lagrange finite elements of order 2 (characteristic size of the mesh L y /20), we have computed the solution for three sizes b ∈ {1.35mm, 1.95mm, 2.55mm} of the computational area Ω b , with the same overlap l = 0.3mm. Figure 4 represents the modulus of the displacement field obtained in these three cases. As we can see in the 2D and 3D configurations, the restrictions of the solutions to the smallest domain match very well. To our knowledge, it is the first time that a transparent boundary condition for anisotropic waveguides is proposed.

Convergence analysis in a simple case

In this section, we will study and compare the convergence of the iterative algorithms presented in Sections 3 and 4. We consider an isotropic waveguide, since the algorithms of Section 3 only apply to this case, and a particular geometry for which analytical calculations can be done. Though, as we will see, the convergence

To study the convergence, it suffices to consider the case of a null source term f . Then the algorithm is convergent if, for a non zero initial data, u b (n) and u + (n) tend to 0 when n goes to +∞. We can focus only on the boundaries and prove that B b (u + (n) ) on Γ b tend to zero. Let us consider for instance that B b (u + (n) ) is known (i.e X + (n) is known) on Γ b and compute the next iteration B b (u + (n+1) ).

Using boundary conditions (36)

This computation is done in three steps:

1. Thanks to the regular geometry of Ω b , using Fraser's relation [START_REF] Baronian | Hybrid modal/FE simulation of guided waves inspections[END_REF] and proceeding as in the proof of Lemma Appendix A.1, we can solve the problem in Ω b and we get:

u b (n) = k>0 Γ b X + (n) Y k 1 e iβ k b -e -iβ k b U + k (x s ) e iβ k z + U - k (x s ) e -iβ k z . (38)
Let us remark that in the above expression, the term e iβ k be -iβ k b may vanish for a particular k. This corresponds to the frequencies for which the problem set in Ω b is not well-posed (see Lemma 3.1), which are excluded here.

2. Then, using

B + (u + (n+1) ) = B + (u b (n)
) and thanks to the modal series expansion [START_REF] Nakamura | Unique Continuation for Elliptic Systems and Crack Determination in Anisotropic Elasticity Contemp[END_REF], we get that

B + (u + (n+1) ) on Γ a is given by Y + (n+1) = Y b (n) = k>0 Γ b X + (n) Y k e iβ k a + e -iβ k a e iβ k b -e -iβ k b Y k (x s ).
3. This provides the modal series expansion of u + (n+1) in the half-guide Ω + a (by formula (4)), from which we easily deduce B b (u + (n+1) ) on Γ + b :

X + (n+1) = k>0 Γa Y + (n+1) X k X k (x s )e iβ k l = k>0 λ k Γ b X + (n) Y k X k (x s ), (39)
where

λ k = 1 + e -2iβ k a 1 -e -2iβ k (a+l) , k = 1, 2, • • • (40)
As usual, the convergence relies on the localization of the λ k 's with respect to the unit disk. The specific difficulty here is that the families (X k ) k and (Y k ) k are not orthonormal.

Since lim

k→+∞

Im(β k ) = +∞, it is clear that |λ k | ∼ k→+∞ e -2Im(β k )l .
Therefore, |λ k | tends to 0 exponentially fast and there are at most a finite number of λ k with modulus greater than 1. Furthermore, by (39), we easily deduce that for all ε ≥ 0:

X + (n+1) L 2 (S) ≤ sup k∈N (|λ k |e εIm(β k )l ) k>0 e -εIm(β k )l X k L 2 (S) Γ b X + (n) Y k .
Then, applying the Cauchy-Schwarz inequality and noticing that

Y k 2 L 2 (S) + X k 2 L 2 (S) ≥ 2 Y k L 2 (S) X k L 2 (S)
, we get:

X + (n+1) L 2 (S) ≤ X + (n) L 2 (S) sup k∈N (|λ k |e εIm(β k )l ) 1 2 k∈N e -εIm(β k )l Y k 2 L 2 (S) + X k 2 L 2 (S) .
Note that for ε = 0, the series on the right hand side diverges. We conjecture that the norm of X k and Y k grow at most polynomially so the series converges for all ε > 0. If all the λ k 's are inside the unit disk, there exists ε ∈ (0, 2) such that sup

k∈N (|λ k |e εIm(β k )l ) < 1
and then the algorithm converges. But numerical examples below show that this condition is not always satisfied. It depends on the configuration and the convergence cannot be ensured systematically.

Remark 5.1. If we take l = 0 (i.e. the non-overlapping case) in formula (40), we can check that lim

k→+∞ λ k = -1.
In that case, even if all the λ k 's are by chance inside the unit circle, the algorithm can never converge geometrically. This shows that the overlap has greatly improved the convergence for high order modes.

Using boundary conditions (37)

The computation of B b (u + (n+1) ) as a function of B b (u + (n) ) can also be made in three steps. 1. The first step is the same as above, and we get [START_REF] Nakamura | Unique Continuation for Elliptic Systems and Crack Determination in Anisotropic Elasticity Contemp[END_REF].

Then, using B

+ (u + (n+1) ) = B + (u b (n)
) and thanks to the modal expansion [START_REF] Nakamura | Unique Continuation for Elliptic Systems and Crack Determination in Anisotropic Elasticity Contemp[END_REF], we get the modal amplitudes

A + k (u + (n+1) ) on Γ a by A + k (u + (n+1) ) = A + k (u b (n) ) = Γ b X + (n) Y k e iβ k a e iβ k b -e -iβ k b , ∀k ≥ 0.
3. Finally, from the modal decomposition of u + (n+1) in Ω + a , we easily deduce B b (u + (n+1) ) on Γ + b :

X + (n+1) = k≥0 A + k (u + (n+1) ) X k (x s )e iβ k l = k≥0 λk Γ b X + (n) Y k X k (x s ), where λk = 1 1 -e -2iβ k b k = 1, 2, • • • . (41)
Compared to (40), we notice a better behaviour of λk for large k. Indeed, λk tends to 0 exponentially fast as λk

∼ k→+∞ e -2Im(β k )b ,
instead of e -2Im(β k )l in the previous algorithm. Surprisingly, λk tends to 0 exponentially fast, at least as fast as e -2Im(β k )a independently of the size of the overlap !

To illustrate these results, we have represented on Figure 6 the localization of the coefficients λ k and λk in the complex plane for a half-waveguide of geometry Ω = [-0.5, 0.5] × R made of aluminium, with the following Lamé's coefficients and density (see [START_REF] Royer | Ondes élastiques dans les solides[END_REF]): λ = 6.08 × 10 10 Pa, µ = 2.83 × 10 10 Pa and ρ = 2.702 × 10 3 kg/m 3 , (42) at a given frequency ω ∼ 31.62 kHz (ω 2 = 1000), such that there are 6 propagative modes. The coefficients λ k and λk related to the propagative modes (resp. the evanescent modes) are represented by circles (resp. by stars). We can observe that:

• in both cases, we cannot control the coefficients λ k or λk associated with the propagative modes by the size l of the overlap,

• for the coefficients λ k , the larger is the overlap, the closer to 0 are the λ k 's associated with the evanescent modes,

• the coefficients λk related to the evanescent modes accumulate at 0, even for the smallest overlap. Note also that the coefficients λk related to the propagative modes have a constant real part which is equal to 1/2. Summing up, this example shows that Schwarz algorithms [START_REF] Tonnoir | Conditions transparentes pour la diffraction d'ondes en milieu élastique anisotrope[END_REF] do not converge in general. However, thanks to the overlap, the sequences (λ k ) k and ( λk ) k tend to 0 very fast, and only a finite number of coefficients may have a modulus greater than 1. This leads us to use the previous ideas in order to build a Schwarz preconditioner for the GMRES algorithm. This will be explained in the next section. Due to the previous analysis, we expect that the convergence of the GMRES will be improved by increasing the overlap, in case of Y type interior transmission condition. On the other hand, for the modal outgoing interior transmission conditions, we can expect a rapid convergence of the associated GMRES algorithm quasi independently of the size of the overlap. In practice, this property enables us to use a small overlap and improve the convergence of the algorithm without increasing too much the computational cost.

Implementation and numerical results

Since the convergence of the iterative algorithms that we have studied in this paper is not ensured, we propose to use a preconditioned GMRES algorithm (see [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving non-symmetric linear systems[END_REF]), where the preconditioner is based on the same tools (we give more details on that point below). It is well known that the GMRES method is a Krylov subspace method (as the conjugate gradient) whose convergence is always ensured, contrary to classical iterative methods like Gauss-Seidel and Jacobi. Let us emphasize that the linear system we are led to solve is not symmetric and we cannot apply the conjugate gradient. It is widely recognized that the choice of the preconditioner is the most critical ingredient to improve the convergence of the GMRES algorithm.

Description of the method

At the discrete level, thanks to Section 3.2, we can show that algorithm (10) + ( 14) can be rewritten as

AU b (n+1) = F -TU b (n) . (43)
where

• F is related to the source term f ;

• U b is expected to be an approximation of the solution in Ω b ;

• A is the sparse matrix associated with the interior problem (inverting the matrix A corresponds to solving the problem set in Ω b ), see Appendix B;

• T is a partially dense matrix related to the non-local operator T ±,l Y→R , see Appendix B.

This is an iterative algorithm to solve the following linear system:

(A + T)U b = F. ( 44
)
Because the matrix T is partially dense, it is prohibitive for realistic problems to make the direct inversion of the matrix A + T. An iterative resolution of this problem is then required.

Algorithm (43) presents several advantages. First A is a sparse matrix, and its inversion is relatively cheap. Secondly, the matrix-vector product TU b (n) can be achieved efficiently without assembling the matrix T (the number of operations is of order N S × M where N S is the number of d.o.f. in S and M is the number of significant modes instead of N 2 S , if T is assembled, see Appendix B).

But again, this iterative algorithm does not always converge, as we have seen in Section 5. The idea then is to use GMRES with A -1 as a preconditioner. Exactly the same idea applies for algorithm [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique[END_REF] + [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF].

The convergence properties of the GMRES algorithm preconditioned by A -1 are linked to the position of the eigenvalues of the matrix M = (I+A -1 T) in the complex plane (see [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving non-symmetric linear systems[END_REF] for instance). More precisely, let us recall that if M a diagonalizable matrix, then the residual norm provided at the m th step of the GMRES algorithm satisfies

r (m+1) ≤ κ(M) r (0) min p∈Pm max λ∈σ(M) |p(λ)| (45)
where κ(M) depends on the spectral decomposition of M, σ(M) is the set of its eigenvalues and P m is the set of polynomial functions p of degree at most m such that p(0) = 1. This means that the convergence of the GMRES depends on the number of the so-called cluster points of the set σ(M).

Here, as explained in Section 3.2, due to the overlap, the operators T ±,l Y→R and T ±,l O→R are regularizing. As a consequence, their eigenvalues tend to zero (the decay is in fact exponential). Thus, the eigenvalues of the associated matrix T concentrate around zero. Therefore, the eigenvalues of M are clustered around 1. This is why our choice of preconditioner is relevant.

2D numerical results

We consider a 2D semi-infinite waveguide made of aluminium (see parameters (42)) of geometry [-0.5, 0.5]× R + . It contains a square cavity centered in (0.5, 0). The length of the diagonal of the square is equal to 0.2. The frequency ω ∼ 31.62 kHz (ω 2 = 1000) is such that there are 6 propagative modes. We consider as incident mode the first propagative mode (that is to say the propagative mode with the largest wavenumber) by imposing a Dirichlet condition on the left boundary of the half-guide. We keep all propagative modes and 20 evanescent modes in the modal series expansions. Figure 7 represents the computed solution. For this configuration, we compare the efficiency of the GMRES algorithm for two choices of transmission conditions and different sizes of the overlap. Figure 8 represents the evolution of the residual error M U bb versus the iteration of GMRES for various sizes of overlap l. As we can see, increasing l improves the convergence if we use the boundary conditions [START_REF] Belgacem | On the Schwarz algorithms for elliptic exterior boundary value problems[END_REF] (see left figure). Now using the modal outgoing transmission condition [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF] (see right figure), as expected, the effect of a small overlap is the same as the effect of larger overlaps when considering (14).

3D numerical results

We consider now a more realistic problem: a tube of diameter 115 mm and thickness 7.5 mm made of steel which is locally perturbed because of corrosion, see Figure 9. The material characteristic coefficients are given by: λ = 8 × 10 10 Pa, µ = 8 × 10 10 Pa, ρ = 7.8 × 10 3 kg/m 3 .

The frequency ω = 25 kHz is such that there exists 21 propagative modes and we keep 79 evanescent modes in our computations. Let us point out that a pre-processing step based on a SAFE method [START_REF] Hayashi | Guided wave simulation and visualisation by a semi-analytical finite element method[END_REF] enables us to compute numerically the modes of the problem. We consider the diffraction problem for an incident field which is the first compressional mode, represented in Figure 10. Figures 11 represents the 3 components of the displacement around the perturbation.

In Figure 12, we compare, as in the 2D case, the efficiency of the GMRES algorithm for the same two different boundary conditions and different sizes of the overlap. The conclusions are the same as in the 2D case.

We compare the computational time of the method using boundary conditions [START_REF] Belgacem | On the Schwarz algorithms for elliptic exterior boundary value problems[END_REF] or [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF] for different values of the overlap (see Figure 13). We represent in red the computational time for the factorization of the preconditioner (which is the same for the two boundary conditions). It increases with the overlap, due to the fact that the size of the interior domain Ω b (and then the size of the matrix) increases. Then we compare for two boundary conditions the remaining computational time. It corresponds to the matrixvector product, multiplied by the number of iterations. For a given overlap, the number of iterations of the GMRES algorithm is smaller for the second choice of transmission condition: it seems to be the best option. Moreover, for the first choice of transmission condition, the computational time first decreases with the overlap, because the convergence is faster. However for large overlaps, it increases with the overlap. The size of the matrix being bigger and bigger, the matrix-vector product spends more time at each iteration of the algorithm. For the second choice of transmission condition, the advantage of the overlap is obtained for the smallest overlap. Considering bigger and bigger overlaps then increases the computational time.

Conclusion

In this paper, we have presented new formulations for scattering problems in isotropic and anisotropic elastic waveguides. A natural splitting of these formulations allows us to define efficient preconditioners for a GMRES algorithm. The formulation based on the modal outgoing condition [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF] seems particularly attractive: indeed, the rate of convergence of the associated iterative algorithm is quasi independent of the size of the overlap l, and the convergence is even faster than for other algorithms. Since a small overlap is used, the storage and the computational cost of the method are also reduced. We prove here that if ω is an eigenfrequency of the homogeneous problem defined by ( 13) + (18), then the limit problem (11) + (13) admits solutions that do not satisfy the compatibility relations [START_REF] Gander | Optimized Schwarz methods with overlap for the Helmholtz equation[END_REF].

First, let us give the following characterization of the eigenfrequencies: Lemma Appendix A.1. The frequency ω is an eigenfrequency of the homogeneous problem (13) + ( 18) if and only if there is k 0 such that β k0 l ≡ π 2 mod (π), (A.1)

where β k is the axial wave number of the k-th mode at the frequency ω (see Section 2.2) and l = ba. Moreover, the associated eigenfunctions in B ± l are proportional to e iβ k 0 (z∓a) U + k0e -iβ k 0 (z∓a) U - k0 .

(A.2)

Proof: Since the calculations are the same in B - l and B + l , we will only make them in B + l . Let v + be a solution of the problem defined by ( 13) + [START_REF] Gregory | A note on bi-orthogonality relations for elastic cylinders of general cross section[END_REF]. Due to the regular geometry of B + l , v + can be decomposed on the rightgoing and the leftgoing modes:

X(v + ) Y (v + ) = k≥0 A + k e iβ k (z-a) X k (x s ) Y k (x s ) + A - k e -iβ k (z-a) -X k (x s ) Y k (x s ) ,
which is associated with the subproblem in Ω b in the algorithm (10) + ( 14).

Let us now explain how to compute the matrix vector product TU b (n) . The i th component of TU b (n) is defined by:

TU b (n) i = M k=0 Γ ± a Y b (n) (x s )X k (x s )dx s e iβ k l Γ ± b S ± k (x s ) + αU ± k (x s ) v i (x s )dx s ,
where M is the number of modes considered in the modal series expansion, u b (n) the function associated with the vector U b (n) and v i is the i th test function of the discrete variational space. The computation of the integrals on Γ ± b are classical when considering non local boundary conditions and are carried out by interpolating S ± k and U ± k . The integrals on Γ ± a are given by definition by:

Γ ± a Y b (n) X k = Γ ± a (u b (n) ) S (S ± k ) S + (t b (n) ) z (U ± k ) z , for k = 1, 2, • • • , M,
The only difficulty is to compute the second part of the integral in the right hand side. To do so, we use the fact that U b (n) was computed by AU b (n) = F -TU b (n-1) which implies that u b (n) satisfies the variational problem in B ± l :

B ± l σ(u b (n) ) : ε(v) -ω 2 ρu b (n) v - Γ ± a t b (n) v - Γ ± b t b (n) v = 0,
for all test function v. Taking a particular test function v that satisfies v z = (U ± k ) z on Γ ± a , v S = 0 on Γ ± a and v = 0 on Γ ± b , we get:

Γ ± a (t b (n) ) z (U ± k ) z = B ± l σ(u b (n) ) : ε(v) -ω 2 ρu b (n) v.
The computational cost of matrix vector product TU b (n) is then of order N S × M if we use a test function v whose support is localised along Γ ± a .

Figure 1 :

 1 Figure 1: Geometry of Ω and notations

Figure 2 :

 2 Figure 2: Decomposition of Ω and notations

Proposition 4 . 3 .

 43 The homogeneous problem[START_REF] Pagneux | Lamb wave propagation in elastic waveguides with variable thickness[END_REF] with the boundary condition defined by B b (w) = σ(w) • ν + α w admits as unique solution w = 0 if Im(α) > 0. As a consequence, if Im(α) > 0, the algorithm defined by (10) + (27) is consistent.

8 2 21 10 30  × 10 10

 3010 Pa, the frequency ω ∼ 31.62 kHz (ω 2 = 1000) is such that there exists 4 propagative modes, the density ρ = 1 kg / m 3 and finally the source term is given by f = [e -200(xy) 2 0] t .

Figure 3 :

 3 Figure 3: Modulus of the computed solution in the 2D anisotropic waveguide for different Ω b : b = 0.5mm (top) and b = 1mm (center) for the same overlap l = 0.1mm. Solution reconstructed in the half-guides Ω ± a (bottom).

Figure 4 :

 4 Figure 4: Modulus of the computed solution in the 3D anisotropic waveguide for different Ω b : b = 1.35mm (top), b = 1.95mm (center) and b = 2.55mm (bottom) for the same overlap l = 0.3mm.

Figure 6 :

 6 Figure 6: Localization of the coefficients λ k (on the left) and λk (on the right) in the complex plane: the color scale indicates the size l of the overlap.

Figure 7 :

 7 Figure 7: Modulus the reconstructed solution in a semi-infinite waveguide made of aluminium of width h = 1.

Figure 8 :

 8 Figure 8: Evolution of the residual error versus the iteration of GMRES for various values of l where l the size of the overlap : using the transmission condition (14) (left figure) and (27) (right figure).

Figure 9 :

 9 Figure 9: Mesh of the regular part (left figure) and the rusted part (right figure) of the tube.

Figure 10 :

 10 Figure 10: The three components of the incident field (ux, uy, uz).

Figure 11 :

 11 Figure 11: The three components of the displacement around the perturbation for the problem (i).

Figure 12 :

 12 Figure 12: Evolution of the residual error versus the iteration of GMRES for various values of l/h where h is the width of the waveguide and l the size of the overlap, using the transmission condition (14) (left figure) and (27) (right figure)

cannot be ensured, the analysis helps to explain the benefit of the overlapping and to present a second important advantage of the modal outgoing boundary conditions. We consider the following problem set in the unperturbed half-guide Ω = S×]0, ∞[ (see Figure 5):

The domain is decomposed into two subdomains Ω + a and Ω b = S×]0, b[ with l = ba ≥ 0 and we consider the following algorithms: z

where B b and B + are defined either by:

or by:

where we recall that the modal amplitudes are given by formula [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. With these choices, we will be able to compute analytically the solution of the previous algorithms (35) + ( 36) and ( 35) + (37) at each step. Let us emphasize that this is possible because the boundary conditions do not couple the modes, contrary to boundary conditions ( 14) and [START_REF] Pagneux | Lamb wave propagation in inhomogeneous elastic waveguides[END_REF]. Besides, let us point out that the iterative processes [START_REF] Tonnoir | Conditions transparentes pour la diffraction d'ondes en milieu élastique anisotrope[END_REF] with boundary conditions [START_REF] Benmeddour | Numerical modeling of guided wave interaction with non-axisymmetric cracks in elastic cylinders[END_REF] or [START_REF] Treyssède | Mode propagation in curved waveguides and scattering by inhomogeneities: application to the elastodynamics of helical structures[END_REF] are well-posed, except for a discrete set of frequencies (see Lemma 3.1) that we will recover in the sequel. Moreover, the algorithm (35) + ( 36) is similar to (10) + (13), and is consistent, except for a discrete set of frequencies linked to the overlapping area (see Lemma 3.2). Algorithm ( 35) + ( 37) is similar to (10) + [START_REF] Halla | Hardy space infinite elements for time-harmonic two-dimensional elastic waveguide problems[END_REF], and is consistent for all frequencies (see Lemma 4.2).

and from the boundary condition X(v + ) on Γ + b , we get

which gives the result. Now, let us consider k 0 such that

is a non-zero solution of ( 13) + ( 18) in B + l .

We want to find (u -, u b , u + ) solution of ( 11) + ( 13) such that u bu -= 0 in B - l and u b -

The idea is to consider the diffraction problem in Ω taking

as incident field. Then the total field u is such that u = u inc + u dif where the diffracted field u dif is outgoing. Finally, setting

, and u + = u dif + e ıβ k 0 (z-a) U + k0 , we get a solution of (11) + (13) which does not satisfy compatibility relations [START_REF] Gander | Optimized Schwarz methods with overlap for the Helmholtz equation[END_REF].

Appendix B. Computational aspects

In (43), the matrix A is related to the following bilinear form: