R Kechroud
email: r.kechroud@yahoo.fr

A Soulaimani

X Antoine
email: xavier.antoine@mines.inpl-nancy.fr

Performance study of plane wave finite element methods with a Padé-type artificial boundary condition in acoustic scattering

Performance study of plane wave finite element methods with a Padé-type artificial boundary condition in acoustic scattering

Ryiad Kerchroud, Azzedine Soulaimani, Xavier Antoine

Introduction

Solving exterior scattering problems at high-frequency is an attractive and challenging computational problem to solve [START_REF] Zienkiewicz | Achievements and some unsolved problems of the finite element method[END_REF]. Significative efforts and improvements have been achieved over the last two decades but severe problems remain unsolved. Various numerical approaches have been developed until now. To this end, let us cite e.g. the integral equation method [START_REF] Burton | The application of integral equation methods to the numerical solution of some exterior boundary-value problems. A discussion on numerical analysis of partial differential equations[END_REF][START_REF] Colton | Integral Equations in Scattering Theory[END_REF][START_REF] Nédélec | Acoustic and Electromagnetic Equations, Integral Representations for Harmonic Problems[END_REF][START_REF] Bendali | Boundary Integral Equations in Acoustic Scattering[END_REF] coupled to an iterative Fast Multipole Method (FMM) or other accelerated solvers [START_REF] Chew | Fast and Efficient Algorithms in Computational Electromagnetics[END_REF][START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Rokhlin | Rapid solution of integral equations of scattering theory in two dimensions[END_REF][START_REF] Bruno | New high-order integral methods in computational electromagnetism[END_REF], the volume Finite Element Method (FEM) [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF][START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF] with truncation boundary (artificial or transparent boundary condition, DtN map, Perfectly Matched Layer,...) [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF][START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF][START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF][START_REF] Bayliss | Boundary conditions for the numerical solution of elliptic equations in exterior regions[END_REF][START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF][START_REF] Givoli | Exact representation on artificial interfaces and applications in mechanics[END_REF][START_REF] Hagstrom | Radiation boundary conditions for the numerical simulation of waves[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains, A review[END_REF][START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Turkel | Boundary Conditions and Iterative Schemes for the Helmholtz Equation in Unbounded Regions[END_REF], the infinite element method [START_REF] Astley | Infinite element for wave problems: a review of current formulations and an assessment of accuracy[END_REF][START_REF] Gerdes | A review of infinite element methods for exterior Helmholtz problems[END_REF] or also the asymptotic techniques [START_REF] Molinet | Asymptotic and Hybrid Methods in Electromagnetics[END_REF]. Despite their respective efficiency, important difficulties still need to be solved. In this paper, we particularly focus on the approach based on FEM.

The Finite Element Method have received much attention because of its flexibility, the availability now to get efficient basic finite element solvers which can be parallelized using for example domain decomposition methods [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF][START_REF] Djellouli | Finite element solution of twodimensional acoustic scattering problems using arbitrarily-shaped convex artificial boundaries[END_REF][START_REF] Tezaur | Three-dimensional finite element calculations in acoustic scattering using arbitrarily-shaped convex artificial boundaries[END_REF][START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized Newton method[END_REF] and its ability to handle complex targets and media [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF]. However, when the frequency increases, some difficulties are met. Indeed, it is well-known that pollution effects [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF][START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF][START_REF] Babuška | Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers[END_REF] arise in standard FEM limiting hence the possibility of solving large scale scattering problems. For this reason, new finite element methods have been explored to assess their ability to solve the problem. These include e.g. the Galerkin Least Squares methods, the element-free methods, the partition of unity method, the discontinuous enrichment method. Rather than citing the numerous contributions, we refer to the complete recent review paper by Thompson [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF] describing most of the available finite element methods with references. Another problem is related to truncating the exterior propagation medium in a suitable way to get the smallest finite domain but at the same time also trying to minimize the spurious reflection at the fictitious boundary. This problem is well-known as building a transparent, an artificial or an absorbing boundary condition. It has received a lot of contributions since the pioneering works of Engquist and Majda [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF] and Bayliss, Gunzburger and Turkel [START_REF] Bayliss | Boundary conditions for the numerical solution of elliptic equations in exterior regions[END_REF][START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF] at the end of the seventies. Since then, many major improvements have been realized [START_REF] Givoli | Exact representation on artificial interfaces and applications in mechanics[END_REF][START_REF] Hagstrom | Radiation boundary conditions for the numerical simulation of waves[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains, A review[END_REF] like for example with the Perfectly Matched Layer method [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Turkel | Absorbing PML boundary layers for wave-like equations[END_REF][START_REF] Turkel | Boundary Conditions and Iterative Schemes for the Helmholtz Equation in Unbounded Regions[END_REF]. Another competitive and related approach which have received much attention these last years is artificial boundary conditions based on rational approximants [START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF][START_REF] Guddati | Continued fraction absorbing boundary conditions for convex polygonal domains[END_REF]. It is shown for example in [START_REF] Guddati | Continued fraction absorbing boundary conditions for convex polygonal domains[END_REF] that similar accuracy can be expected for these methods and the PML technique. In the case of elongated scatterers which interest us in this paper, a general high-order accurate Padé-type On-Surface Radiation Condition (OSRC) [START_REF] Kriegsmann | A new formulation of electromagnetic wave scattering using the on-surface radiation condition method[END_REF][START_REF] Antoine | Advances in the On-Surface Radiation Condition Method: Theory, Numerics and Applications[END_REF] has been derived in [START_REF] Antoine | An improved On-Surface Radiation Condition for acoustic scattering in the high-frequency spectrum[END_REF][START_REF] Antoine | An improved surface radiation condition for highfrequency acoustics scattering problems[END_REF] and its performance as an Artificial Boundary Condition (ABC) has been provided in [START_REF] Kechroud | Numerical accuracy of a Padé-type nonreflecting boundary condition for the finite element solution of acoustic scattering problems at high-frequency[END_REF][START_REF] Kerchroud | Padé absorbing boundary conditions for finite element solution of high frequency scattering problems[END_REF] using standard finite element methods. The conclusion is that the method leads to a drastic diminution of the size of the computational domain leading hence to the possibility of prospecting high-frequency problems. However, despite this noticeable improvement, the pollution problem still remains present into the finite element method. The aim of this paper is to show that this problem can be significantly relaxed by using more adapted finite element techniques. To this aim, we propose here to analyze the performance of the Plane Wave FEM coupled to the Padé-type ABC for two-dimensional scattering problems in order to strongly decrease the number of degrees of freedom of the final linear system.

The structure of the paper is the following. In Section 2, we introduce the scattering problem and describe the Padé-type ABC and the associated variational formulation. In Section 3, the standard and plane wave-based FEM are described and discussed. Section 4 provides a complete numerical study of the performance of the Padé-type ABC coupled to the Plane Wave FEM. In particular, Section 4.3 reports some results obtained in the case of a submarine-like shaped scatterer. Finally, a conclusion is given in Section 5.

A Padé-type Artificial Boundary Condition

The two-dimensional scattering problem

Define Ω -⊂ R 2 as a two-dimensional impenetrable bounded domain with boundary Γ := ∂Ω -. The associated homogeneous exterior domain of propagation, which is the complemetary set of the scatterer Ω -in R 2 , is denoted by Ω e . Then, the scattering of an incident time-harmonic acoustic wavefield u inc by Ω -can be formulated as the following exterior Boundary Value Problem (BVP): find the scattered field u solution to

   ∆u + k 2 u = 0, in Ω e , ∂ n Γ u = -∂ n Γ u inc or u = -u inc , on Γ, lim r→∞ √ r(∂ r u -iku) = 0. (1)
If a and b are two complex-valued vector fields (and z denotes the complex conjugate of a complex number z ∈ C), their inner product is defined by a • b = 2 j=1 a j b j , and the associated norm • is:

a 2 = a • a. Let x = (x 1 , x 2 ) ∈ R 2 ,
then the gradient ∇ of a complex-valued scalar field f and the divergence div of a are defined respectively by: ∇f = (∂ x 1 f, ∂ x 2 f ) T and div a = 2 j=1 ∂ x j a j , designating by a T the transposed of a. Under these notations, the Laplace operator ∆ is classically defined by: ∆ = div ∇. We consider that the incident wave u inc is plane: u inc (x) = e ikd•x . The wave number k is related to the wavelength λ by the relation k = 2π/λ. The direction of incidence d is given through the relation: d = (cos(θ inc ), sin(θ inc )) T , where θ inc is the scattering angle. If we define by n Γ the outwardly directed unit normal to Ω e at the boundary Γ, then, the sound-hard or Neumann (respectively sound-soft or Dirichlet) boundary condition on Γ corresponds to the second (respectively third) equation of [START_REF] Zienkiewicz | Achievements and some unsolved problems of the finite element method[END_REF]. Finally, the last equation is the Sommerfeld radiation condition which allows only outgoing waves at infinity, setting r = x . This thereby guarantees the uniqueness of the solution to the BVP (1).

Bounding the domain by using a Padé-type ABC

It is well-known that the BVP (1) cannot be solved by usual domain based methods like the Finite Element Method (FEM) or the Finite Difference Method (FDM) if the infinite domain Ω e is not truncated via a fictitious boundary Σ enclosing Ω -. This implies that the considered domain of computation is now the one, denoted by Ω, delimited by Γ and Σ. To avoid in the best case or at least to minimize the reflection at the nonphysical boundary Σ, one must impose a suitable boundary condition at Σ. This condition takes different denominations like non-reflecting, artificial or absorbing boundary condition, according to its properties and aims [START_REF] Givoli | Exact representation on artificial interfaces and applications in mechanics[END_REF], but its goal is always to try to give a good compromise between flexibility in terms of implementation into an existing code and minimization of the reflection or/and of the size of the computational domain. Generally, this condition is given through the Dirichlet-to-Neumann (DtN) operator or an approximation of this operator. More precisely, if M is an approximation of the DtN operator, we get the following new approximate BVP with an Artificial Boundary Condition (ABC)

   ∆u + k 2 u = 0, in Ω, ∂ n Γ u = -∂ n Γ u inc or u = -u inc , on Γ ∂ n Σ u = -M u, on Σ. (2)
In a series of recent papers [START_REF] Djellouli | Finite element solution of twodimensional acoustic scattering problems using arbitrarily-shaped convex artificial boundaries[END_REF][START_REF] Tezaur | Three-dimensional finite element calculations in acoustic scattering using arbitrarily-shaped convex artificial boundaries[END_REF][START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized Newton method[END_REF], Farhat et al. show that the generalized secondorder Bayliss-Gunzburger-Turkel (BGT2) like ABC derived in [START_REF] Antoine | Bayliss-Turkel-like Radiation Condition on Surfaces of Arbitrary Shape[END_REF] yields accurate finite element solutions in the mid-frequency regime for two-and three-dimensional direct and inverse acoustic scattering problems using arbitrarily-shaped convex fictitious boundaries Σ. Even if this condition is accurate in the medium frequency range, its precision decreases in the high-frequency regime if the fictitious boundary is too close to the scatterer. To obtain a much better accuracy, a larger computational domain must be considered resulting in the need to solve a very large size sparse linear system of equations. This is then computationally expensive both in memory and time even with sophisticated solvers like for instance the FETI-H method [START_REF] Djellouli | Finite element solution of twodimensional acoustic scattering problems using arbitrarily-shaped convex artificial boundaries[END_REF][START_REF] Tezaur | Three-dimensional finite element calculations in acoustic scattering using arbitrarily-shaped convex artificial boundaries[END_REF][START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized Newton method[END_REF].

An alternative to the BGT2-like ABC is to rather consider the Padé-type ABC derived in and validated in [START_REF] Kechroud | Numerical accuracy of a Padé-type nonreflecting boundary condition for the finite element solution of acoustic scattering problems at high-frequency[END_REF][START_REF] Kerchroud | Padé absorbing boundary conditions for finite element solution of high frequency scattering problems[END_REF] for high-frequency acoustic scattering. This condition is expressed through a square-root operator which can be efficiently simulated by paraxial approximation techniques. More specifically, this operator is given by the relation

-M u = ik 1 + ∂ s ( 1 k 2 ∂ s )u - κ 2 u + κ 2 8(κ -ik) u -∂ s ( κ 2k 2 ∂ s )u. (3)
In this expression, s is the counterclockwise directed arclenght along Σ, ∂ s is the curvilinear derivative along Σ and κ is the curvature at a point of the surface. The notation √ z designates the principal determination of the square-root of a complex number z with branch-cut along the negative real axis. This operator which is in fact defined by a non-local pseudo-differential square-root operator must be efficiently represented through differential operators to have a sparse matricial representation. Following [START_REF] Antoine | An improved On-Surface Radiation Condition for acoustic scattering in the high-frequency spectrum[END_REF][START_REF] Antoine | An improved surface radiation condition for highfrequency acoustics scattering problems[END_REF], this can be accurately done by a rotating branch-cut approximation of the square-root and Padé approximants [START_REF] Milinazzo | Rational square-root approximations for parabolic equation algorithms[END_REF]. Moreover, to get an accurate representation of the tangential rays into the ABC (3), a local regularization procedure must be applied to the square-root operator. This is done by the introduction of a complex wave number k ǫ with small dissipation ǫ (much more details are given in [START_REF] Antoine | An improved surface radiation condition for highfrequency acoustics scattering problems[END_REF]). Finally, the modified regularized square-root operator used in ( 3) is replaced and approximated by

1 + ∂ s ( 1 k 2 ǫ ∂ s )u ≈ C 0 u + N j=1 A j ∂ s (k -2 ǫ ∂ s )(1 + B j ∂ s (k -2 ǫ ∂ s )) -1 u (4)
where k ǫ = k + iǫ. An optimized choice of the damping parameter ǫ stated in [START_REF] Antoine | An improved surface radiation condition for highfrequency acoustics scattering problems[END_REF] is: ǫ = 0.4k 1/3 κ 2/3 . In ( 4), the complex coefficients C 0 , A j and B j are given by

C 0 = e i α 2 R N (e -iα -1), A j = e -iα 2 a j (1 + b j (e -iα -1)) 2 , B j = e -iα b j (1 + b j (e -iα -1)) 2 .
(5)

The angle of rotation is α and (a j , b j ), for j = 1, ..., N , are the standard real Padé coefficients given by

a j = 2 2N + 1 sin 2 ( jπ 2N + 1 ) , b j = cos 2 ( jπ 2N + 1 ), (6)
denoting by R N the Padé approximant of order N

√ 1 + z ≈ R N (z) = 1 + N j=1 a j z 1 + b j z . (7)
In view of an efficient numerical treatment, the approximation of the Padé-type ABC (3)-( 7) is represented by using Lindman's [START_REF] Lindmann | Free-space boundary conditions for the time dependent wave equation[END_REF] auxiliary coupled functions trick

-M u = ik(u + N j=1 A j ϕ j ) - κ 2 u + κ 2 8(κ -ik) u -∂ s ( κ 2k 2 ∂ s )u, on Σ, (8)
where the functions ϕ j , j = 1, .., N , defined on Σ, are solutions of the following differential equations

(1 + B j ∂ s ( 1 k 2 ǫ ∂ s ))ϕ j = ∂ s ( 1 k 2 ǫ ∂ s )u. ( 9
)

Variational formulation with Padé-type ABC

Let us introduce V as the Sobolev space H 1 (Ω) for a Neumann boundary condition or H 1 0 (Ω) for a Dirichlet boundary condition. Let us define the product spaces: W N := H 1 (Σ) × ... × H 1 (Σ), N times, and X N := V × W N . Under these notations, a variational formulation of (2) with the Padé-type ABC (3)-( 7) is: find [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF] with j = 1, ..., N , and for test-functions (v, ψ 1 , ..., ψ N ) ∈ X N . The symmetrical bilinear form A acts from V × V into C and is defined by

(u, ϕ 1 , ..., ϕ N ) in X N such that      A(u, v) + N j=1 B j (ϕ j , v) = b(v), C(u, ψ j ) + D j (ϕ j , ψ j ) = 0,
A(u, v) = Ω ∇u • ∇v -k 2 uv dΩ -ik Σ uvdΣ + 1 2 Σ κuvdΣ - 1 8 Σ κ 2 κ -ik uvdΣ + 1 2k 2 Σ κ∂ s u∂ s vdΣ. (11)
Moreover, the bilinear forms {B} j=1,...,N , C and {D} j=1,...,N are defined from

H 1 (Σ) × H 1 (Σ) into C by B j (ϕ j , v) = -ikA j Σ ϕ j vdΣ, C(u, ψ j ) = 1 k 2 ǫ Σ ∂ s u∂ s ψ j dΣ, D j (ϕ j , ψ j ) = Σ ϕ j ψ j dΣ - B j k 2 ǫ Σ ∂ s ϕ j ∂ s ψ j dΣ. (12)
We do not precise the definition of the linear form b which follows from the nature of the boundary condition on Γ. The variational formulation ( 10) is an unconjugated formulation with bilinear operators. The following alternative conjugated formulation with sesquilinear operators can be obtained if the weighting functions v and ψ j are substituted respectively by v and ψ j ,

     A(u, v) + N j=1 B j (ϕ j , v) = b(v), C(u, ψ j ) + D j (ϕ j , ψ j ) = 0, (13)
for j = 1, ..., N , and ∀(v, ψ 1 , ..., ψ N ) ∈ X N .

Finite element approximation

To discretize the variational formulations, the computational domain Ω is partitioned into quadratic triangular finite elements resulting in a covering Ω h . The unknown scattered field u within each finite element can be approximated by using standard polynomial shape functions N j , j = 1, .., 6, as follows

u = 6 j=1 N j u j , (14)
denoting by u j , j = 1, .., 6, the nodal values of u at the interpolation points. This finite element is designated in the sequel as the T6 finite element (triangular finite element with six nodes) Another way of building interpolating finite element functions, called plane wave finite element functions, is to approximate the unknown field u within each triangular finite element by using the standard polynomial shape functions N j , j = 1, .., 6, each function being enriched by n q radiating plane waves e ikdq•(x-x j ) centered at the j-th nodes and for n q equally spaced directions of propagation

d q = (cos( 2π(q -1) n q ), sin( 2π(q -1) n q )) T , q = 1, .., n q . ( 15
)
It results that an approximation of the wavefield, solution to [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF], is seek through the following expansion on each triangle

u = 6 j=1 nq q=1
N j e ikdq•(x-x j ) u jq [START_REF] Givoli | Exact representation on artificial interfaces and applications in mechanics[END_REF] setting u jq as its nodal values. This approximation is designated in the sequel of the paper as PWT6 finite element (Plane Wave with T6). One can also choose to discretize the alternative variational formulation [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF]. In this case, the finite element method is referred to as CPWT6 (Conjugate PWT6). The aim of these new basis functions is to partially reproduce the oscillations of the scattered field in Ω. The approximation [START_REF] Givoli | Exact representation on artificial interfaces and applications in mechanics[END_REF] writes in a certain way the restriction of u on each triangle as a superposition of elementary plane waves of amplitude given by the shape functions N j and with known phase through the idea of enrichment. The elementary matrices can be computed over each element analytically using symbolic calculation or numerically by using adapted integration procedures [START_REF] Bettess | A numerical integration scheme for special finite elements for the Helmholtz equation[END_REF] or high-order Gauss-Legendre quadratures. Depending on the frequency, we have used up to 27 integration points during the numerical experiments presented in the next section.

The finite element discretization of the variational problems [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF] and ( 13) and using basis functions [START_REF] Bayliss | Boundary conditions for the numerical solution of elliptic equations in exterior regions[END_REF] or ( 16) leads to a linear system of coupled equations of the generic form

A h B h C h D h u h ϕ h = b h 0 . ( 17
)
The solution (u h , ϕ h ) is composed from the approximate wavefield u h in C n h , where n h is the number of degrees of freedom of the finite element method (T6, PWT6 or CPWT6) for the covered domain Ω h , h being the largest finite element size, and ϕ h is in C N n Σ h , designating by n Σ h the number of discretization points of the interpolated boundary Σ h := ∂Ω h . We classically introduce n λ := h/λ as the density of discretization points per wavelength which is often used in scattering problems to measure the thickness of the mesh. The complex linear system ( 17) is sparse, globally non-symmetric, non-Hermitian and non-diagonally dominant. Its size is equal to #dof, where #dof designates the total number of degrees of freedom involved in the numerical solution. In the case of the unconjugated formulation, all sub-matrices A h , B h , C h and D h are however symmetric.

In this paper, we are mainly interested in studying the accuracy of the plane wave finite element method in conjunction with the Padé-type ABC. For this reason, we solve the linear system with a direct solver. It is known that plane wave finite element methods lead to ill-conditioning. This issue is not addressed here since only two-dimensional problems are studied and linear systems can be handled by a direct solver.

Numerical study

Based on the numerical experiments presented by Kechroud et al. in [START_REF] Kechroud | Numerical accuracy of a Padé-type nonreflecting boundary condition for the finite element solution of acoustic scattering problems at high-frequency[END_REF], we set the parameters of the Padé-type ABC to N = 2 and θ = π/6 which yield a minimal spurious reflection at the boundary when solving exterior problems. The ABC placement is measured by a parameter m which is the distance in terms of wavelengths between the boundary Γ h of the scatterer Ω h and the fictitious boundary Σ h .

To measure the accuracy of the different finite element methods and formulations, we compute the relative Root-Mean Square (RMS) error (in percents %) between the reference solution and the approximate solution onto the computational domain Ω h (respectively, Γ h ) in the L 2 (Ω h )-norm (respectively, L 2 (Γ h )-norm). The reference solution is computed analytically if available or numerically with the CHIEF integral equation method [START_REF] Kirkup | The boundary element method in acoustics. A development in Fortran[END_REF] for a high density n λ . We also represent the far-field pattern given by the scattering amplitude

u ∞ (θ) = e i π 4 √ 8πik Γ (∂ n Γ u + ikn Γ • d ′ u)e -ikx•d ′ dΓ ( 18
)
through the bistatic Radar Cross Section (RCS) (also called the target strength)

RCS(θ) = 10 log 10 (2π|u ∞ (θ)| 2 ) (dB) (19)
setting d ′ = (cos(θ), sin(θ)) T as the vector of observation in the polar coordinates system (r, θ).

The circular cylinder

The first test-case is the scattering problem of a plane wave of incidence d = (1, 0) T by a circular cylinder of radius a = 1. The Padé-type condition is placed at a distance b = 2. It is known that the exact analytic solution to the exterior scattering problem (1) expands in Mie series (see Equation 23). However, comparing this solution to the one numerically computed with an artificial boundary condition would include some errors coming from the truncation by the fictitious boundary Σ. To avoid this problem and only observe the pollution involved into the finite element method, we begin by considering the solution to (2) with the Padétype artificial boundary condition. Since this reference wavefield is solution to a bounded problem set in the crown Ω, then we can expand it as

u ex (r, θ) = m∈N ǫ m (-i) m (A m H (1) m (kr) + B m H (2) m (kr)) cos(mθ), (20)
where ǫ m is the Neumann function which is equal to 1 for m = 0 and 2 otherwise. Functions H

m and H

m are respectively the first-and second-kind Hankel functions of order m. In the sequel, the prime ′ denotes the derivative of a function with respect to the radial variable r and J m is the Bessel function of order m. The complex-valued coefficients A m and B m are computed by imposing the Neumann boundary condition at Γ and the Padé-type condition at Σ. This leads to find (A m , B m ) as the solution to the following linear system

H (1) ′ m (ka) H (2) ′ m (ka) H (1) ′ m (kb) -ik m H (1) m (kb) H (2) ′ m (kb) -ik m H (2) m (kb) A m B m = -J ′ m (ka) 0 , (21)
setting

ik m = ikC 0 -ik N j=1 A j m 2 /(b 2 k 2 ε ) 1 -B j m 2 /(b 2 k 2 ε ) - 1 2b + 1 8b 2 (1/b -ik) + m 2 2k 2 b 3 .

The solution to (21) is given by

A m = -J ′ m (ka) ∆ m (H (2) ′ m (kb) -ik m H (2) m (kb)), B m = J ′ m (ka) ∆ m (H (1) ′ m (kb) -ik m H (1) m (kb)). ( 22
)
where

∆ m = H (1) ′ m (ka)(H (2) ′ m (kb) -ik m H 2 m (kb)) -H (2) ′ m (ka)(H (1) ′ m (kb) -ik m H (1)
m (kb)). The analytical solution [START_REF] Turkel | Absorbing PML boundary layers for wave-like equations[END_REF] is computed by using a finite number of Mie coefficients equal to 2kr + 5 to ensure the convergence.

We consider a structured mesh built using using n r elements in the radial direction and n θ elements in the angular direction θ. We report on Table 1 the RMS error of the solution in the computational domain Ω h and on the scatterer boundary Γ h using the T6 FEM. We remark that the mesh refinement leads to a convergence of the numerical solution. Moreover, the error increases classically with ka. Most specifically, for ka = 20, a mesh resolution less than n r = 12 and n θ = 60 does not give a satisfactory error. The situation is clearly much better if one considers now the PWT6 and CPWT6 FEM as seen on Table 2. Indeed, we can clearly observe an accuracy improvement of the error for the two coarse meshes n θ × n r = 3 × 15 and n θ × n r = 6 × 30 when the number of directions n q increases. The CPWT6 always yields better results while it is necessary to consider an even number of directions n q for the PWT6 FEM. For example, an error on the domain solution equal to 0.56% is obtained by the CPWT6 using 9360 dof while an error of 0.4% with 47040 dof is achieved for the T6 FEM. If we now fix an error of 5% at ka = 20, then we see that it can obtained using 12000 dof for the T6 FEM while only 3360 dof are required for the CPWT6 FEM using a mesh with 3 × 15 elements and n q = 4 directions.

From the previous numerical experiments, it appears that the error control due to the pollution with high wave numbers could be achieved by increasing the mesh resolution while maintaining the number of directions constant in the CPWT6 and PWT6 schemes, or on the contrary, increasing the number of directions while keeping the mesh resolution constant. Furthermore, the CPWT6 and PWT6 schemes deliver globally the same accuracy for an even number of directions n q . This is expected since for every wave going in one direction there is another one in the opposite direction. In this case, the complex representation of plane waves can be replaced by real-valued trigonometric basis functions. A consequence of this alternative representation is that the unconjugated variational formulation with sesquilinear operators [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF] and the conjugated variational formulation with bilinear operators [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF] become equivalent.

We now consider again the same scatterer for a wave number ka = 60. The direction of incidence of the plane wave is imposed such that d = (1, 0) T . The exterior circular boundary Σ has a radius b equal to a + mλ. To take into account now the effect of truncating the exterior domain, the reference solution solves the boundary value problem (1) in the polar coordinates system (r, θ) 3: Sound-hard and sound-soft circular cylinder: RMS error in the computational domain Ω h (and on Γ h for the sound-hard case) at ka = 60 and θ inc = 0 degree for the T6 FEM .

u ex (r, θ) = m∈N ǫ m (-i) m J ′ m (ka)(kr)H (1) m (kr) H (1) ′ m (ka) cos(mθ), ( 23
)
n θ × n r #dof ka = 1 ka = 5 ka = 10 ka = 20
for the sound-hard scatterer and

u ex (r, θ) = m∈N ǫ m (-i) m J m (ka)(kr)H (1) m (kr) H (1) m (ka) cos(mθ), ( 24
)
for sound-soft problem. We report in Tables 3 and4 the relative errors with respect to the analytical solution of the scattered field in Ω h and on Γ h for the sound-hard disk and only in Ω h for the sound-soft case. We can observe different important points

• using the PWT6 or CPWT6 FEM lead to a significant gain of accuracy in the solution for a low density of discretization points n λ .

• increasing n q yields an interesting accuracy improvement. Again, an even value of n q is required for stability reasons linked to the PWT6 FEM.

• an excellent accuracy is given even for a close fictitious boundary Σ leading to a drastic diminution of the total number of degrees of freedom.

• the accuracy is generally better for the sound-soft problem.

Another quantity of interest for practitioners is the RCS. To this end, we report on Figures 1 and2 the RCS for a position of the boundary m = 0.15 and n λ = 1, for n q = 2 directions in the CPWT6 FEM for respectively the sound-hard and sound-soft problems. We observe an excellent agreement between the reference and the numerical solutions.

Figure 1: Radar Cross Section of the sound-hard circular cylinder at ka = 60, θ inc = 0 degree, m = 0.15, n λ = 1 and n q = 2 using the Padé-type ABC.

m n λ n q #dof hard, CPWT6 soft, CPWT6 hard, PWT6 Soft, PWT6 err. in % err. in % err. in % err. in % 0. [START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF] 4: Sound-hard and sound-soft circular cylinder: RMS error in the computational domain Ω h (and on Γ h for the sound-hard case) at ka = 60 and θ inc = 0 degree for the PWT6 and CPWT6 FEM.

Figure 2: Radar Cross Section of the sound-soft circular cylinder at ka = 60, θ inc = 0 degree, m = 0.15, n λ = 1 and n q = 2 using the Padé-type ABC.

The sound-hard elliptical cylinder

We consider now a sound-hard elliptical cylinder with a major and minor semi-axis respectively equal to a = 1 and b = 0.25 along the x-and y-directions. The obstacle is centered at the origin. The wave number is set to ka = 60 and the direction of the incident plane wave is d = (

√ 2 2 , √ 2
2 ) T . The reference solution is computed by the (CHIEF) integral equation method [START_REF] Kirkup | The boundary element method in acoustics. A development in Fortran[END_REF] using a mesh resolution of n λ = 80 elements per wavelength to ensure a high accuracy. The fictitious boundary Σ is an ellipse with major (respectively minor) semi-axis equal to a ′ = a + mλ (respectively b ′ = b + mλ).

We report in Tables 5 and6 the relative error with respect to the reference solution of the trace on Γ h of the scattered field. On this new example different points must be noticed

• an important accuracy improvement is again obtained when using the CPWT6 and PWT6 FEM. Moreover, the CPWT6 method generally leads to a better accuracy. Using the PWT6 FEM requires an even value of n q to get a satisfactory precision. Finally, increasing n q gives a higher accuracy.

• the fictitious boundary can be set very close to the scatterer but, to obtain a good accuracy for low densities of discretization points, the CPWT6 or PWT6 FEM must be used with n q ≥ 2. Of course, a limited accuracy is obtained for a fixed distance m and, as expected, increasing m yields an improved solution.

We see on this example that we can get a low error approximation if we combine the CPWT6 (or PWT6) FEM with the Padé-type ABC. We must notice at this point that using a low-order ABC for high wave number would have the effect of leading to large computational domains and much more pollution. This is why the two following points must be used in conjunction: 1) to use a high-order accurate ABC for large wave numbers to have a small size computational domain and therefore also less pollution and 2) to consider a suitable approximation method like the CPWT6 FEM to reduce the pollution effect.

To see the impact on the RCS computation, we present some results on Figures 3 and4 obtained for different mesh resolutions and values of n q for the CPWT6 and T6 FEM. We see that a much better accuracy is obtained for the CPWT6 approach using for example the discretization corresponding to 3 × 120 dof . This is again confirmed on Figure 5 where the fictitious boundary is placed at a larger distance to improve the accuracy.

m n λ #dof Err. on Γ h (%) 0. [START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF] Figure 4: Comparison of the computed RCS of the sound-hard elliptical cylinder at ka = 60, for θ inc = 45 degrees and m = 0.15. We again increase the mesh resolution and compare it to the T6 FEM.

The submarine-like shaped scatterer

To consider a more realistic case, we choose the submarine-like shaped scatterer shown on Fig. 6. We set its dimensions as follows: the total length is L = 11m, its thickness is D = 1m, the height and length of the tower are respectively l = 1m and 2l = 2m. We perform several acoustic scattering computations at for kD = 20 for an incidence vector

d = -( √ 2 2 , √ 2
2 ) T . We place an elliptical fictitious boundary Σ (Fig. 6) with a major semiaxis (resp. a minor semi-axis) equal to a ′ = L/2 + mλ (respectively b ′ = D/2 + l + mλ) which circumscribe the submarine. The reference solution is obtained again with the CHIEF integral equation method [START_REF] Kirkup | The boundary element method in acoustics. A development in Fortran[END_REF] for a mesh resolution corresponding to n λ = 40 elements per wavelength to ensure a high accuracy.

We report the RMS error on Γ h for the computed scattered field using the different FEM in Tables 7 and8. For completeness, we plot on Figures 7, 8 and 9 the bistatic RCS computed by using the reference CHIEF solution and the solutions obtained using the T6 and CPWT6 FEM for different mesh resolutions and positions of the ABC. From the reported results, we can notice the following different points to attain an suitable accuracy for engineering calculations

• We have to increase the number of directions in the CPWT6 and PWT6 FEM when the mesh resolution is low in order to control the pollution error. Indeed, a number of directions n q = 6 seems to be sufficient for a mesh resolution of n λ = 1. This can be further decreased to only two directions per node if the mesh resolution consists of four elements per wavelength but the cost is an increasing number of degrees of freedom.

• A better accuracy is obtained when the fictitious boundary is placed at m = 2 as expected.

• For a prescribed a priori engineering accuracy (Table 9), the CPWT6 and PWT6 FEM require only coarse meshes and leads to small systems of equations in comparison with the classical T6 finite element.

Another important point is that the predicted RCS (see Figures 7,8 and 9) with the CPWT6 finite element are in good agreement with the reference solution even with a small number of degrees of freedom. The situation is quite different for the RCS computed by the T6 finite element solution since a correct agreement with the reference RCS is lost when the mesh resolution is less than three elements per wavelength. These last computations on a more realistic test-case shows that it is finally important for efficiency and accuracy reasons to combine both a high-order ABC and a suitable finite element method like the CPWT6 FEM for large frequencies computations. Table 9: Sound-hard submarine-like scatterer: mesh resolution and #dof needed to achieve a prescribed accuracy using the T6, CPWT6 and PWT6 finite elements for kD = 20 and θ inc = 225 degrees. element methods for solving two-dimensional high-frequency scattering problems. It is shown that small computational domains in terms of wavelength are needed to achieve a satisfactory accuracy on the trace of the wavefield and radar cross section. In particular, the method leads to a reduction of the effect of pollution into the global numerical method. The superiority of the conjugate finite element method is shown compared to a classical quadratic Lagrange finite element method and the unconjugate plane wave version. The performance is analyzed in detail for three numerical examples and in particular for a submarine-like scatterer. In a forthcoming work, we will focus on the extension to three-dimensional scattering problems.

Figure 3 :

 3 Figure 3: Comparison of the computed RCS of the sound-hard elliptical cylinder at ka = 60, for θ inc = 45 degrees and m = 0.15. We use the CPWT6 FEM for different values of n q .

Figure 5 :

 5 Figure 5: Comparison of the computed RCS of the sound-hard elliptical cylinder at ka = 60, for θ inc = 45 degrees and m = 1.2.

Figure 6 :

 6 Figure 6: Configuration for the computations: the submarine-like shaped scatterer is enclosed by an elliptical fictitious boundary Σ.

Figure 7 :

 7 Figure 7: Comparison of the RCS of the sound-hard submarine-like scatterer for kD = 20, m = 1.0 and θ inc = 225 degrees using the Padé-type ABC with the T6 and the CPWT6 finite elements (setting n λ = 1 and n q = 5 for the CPWT6 FEM and n λ = 2.5 for the T6 FEM).

Figure 8 :

 8 Figure 8: Comparison of the RCS of the sound-hard submarine-like scatterer for kD = 20, m = 1.0 and θ inc = 225 degrees using the Padé-type ABC with the T6 and CPWT6 finite elements (setting now n λ = 2 and n q = 3 for the CPWT6 FEM and n λ = 3.5 for the T6 FEM).

Figure 9 :

 9 Figure 9: Comparison of the RCS of the sound-hard submarine-like scatterer for kD = 20, m = 2.0 and θ inc = 225 degrees using the Padé-type ABC with the T6 and CPWT6 FEM (n λ = 1 and n q = 5 for the CPWT6 and n λ = 2.3 for the T6 FEM).

Table 1 :

 1 Sound-hard circular cylinder: relative RMS error (in %) in the computational domain Ω h (respectively on Γ h ) of the T6 finite element for different meshes.

	840 0.04 (0.02) 1.45 (0.87) 43.48 (40.68) 138.58 (48.55) 3120 0.01 (0.01) 0.10 (0.06) 2.39 (0.39) 67.87 (64.86) 5.34 (1.47) 60 × 12 12000 0.01 (0.01) 0.04 (0.03) 0.19 (0.06) 15 × 3 30 × 6 120 × 24 47040 0.01 (0.01) 0.03 (0.02) 0.03 (0.02) 0.40 (0.13)
	FEM	#dof	n q	ka = 1	ka = 5	ka = 10	ka = 20
	74.89 (5.94) 6.81 (4.60) 0.29 (0.28) 0.05 (0.05) 1 × 840 1 0.05 (0.04) 23.97 (26.42) 135.36 (94.35) 226.23 (179.65) 84.13 (16.30) 1 × 840 1 0.05 (0.04) 20.87 (1.15) 54.02 (14.53) 2 × 840 2 0.05 (0.04) 0.23 (0.23) 19.74 (5.76) CPWT6 3 × 840 3 0.05 (0.04) 0.02 (0.03) 2.70 (1.48) 4 × 840 4 0.05 (0.04) 0.02 (0.03) 6.81 (4.60) 54.02 (14.53) 2 × 840 2 0.05 (0.04) 0.23 (0.23) PWT6 2.74 (2.62) 259.63 (225.30) 3 × 840 3 0.05 (0.04) 0.12 (0.10) 0.05 (0.05) 2.70 (1.48) 4 × 840 4 0.05 (0.04) 0.02 (0.03) 37.91 (0.78) 80.40 (4.03) 1 × 3120 1 0.01 (0.01) 1.67 (0.12) 0.37 (0.19) 12.85 (7.03) 2 × 3120 2 0.01 (0.01) 0.18 (0.18) 0.10 (0.08) 0.56 (0.42) CPWT6 3 × 3120 3 0.01 (0.01) 0.02 (0.03) 0.04 (0.02) 0.06 (0.06) 4 × 3120 4 0.01 (0.01) 0.02 (0.03) 19.45 (12.59) 175.94 (120.43) 1 × 3120 1 0.01 (0.01) 2.27 (1.22) 0.37 (0.19) 12.85 (7.03) 2 × 3120 2 0.01 (0.01) 0.18 (0.18) 0.21 (0.11) 8.25 (4.29) PWT6 3 × 3120 3 0.01 (0.01) 0.10 (0.10) 4 × 3120 4 0.01 (0.01) 0.02 (0.03) 0.04 (0.02) 0.06 (0.06)

Table 2 :

 2 Sound-hard circular cylinder: relative RMS error (in %) in the computational domain Ω h (respectively on Γ h ) of the CPWT6 and PWT6 finite elements for two coarse mesh n θ × n r corresponding respectively to 3 × 15 and 6 × 30.

	m n λ # dof	hard	soft
				err. in %	err. in %
	0.15 1	600	76.56(79.17)	47.07
		2	752	21.65(21.72)	5.71
		4	1452	4.48(4.79)	2.67
		8	2904	2.59(2.59)	0.63
		16 11548	2.59(2.59)	0.57
	1.2	2	1344 38.60(38.20)	27.36
		4	5896	3.74(1.85)	3.74
		8 22676	0.89(1.11)	0.88
		16 89376	0.89(1.11)	0.88
	Table

Table 5 :

 5 Sound-hard elliptical cylinder: RMS error on Γ h for the T6 FEM at ka = 60, θ inc = 45 degrees.

		1	272	66.89
		2	524	18.58
		4	996	4.43
		8	1992	4.43
	1.2	2	944	29.56
		3	2172	5.27
		4	3824	2.33
		8 15336	1.49
		16 18352	1.47

Table 6 :

 6 Sound-hard elliptical cylinder: RMS error on Γ h for the CPWT6 and PWT6 FEM at ka = 60, θ inc = 45 degrees.

	m	n λ n q	#dof	CPWT6 Err. PWT6 Err.
				on Γ h in %	on Γ h in %
	0.15 0.5 1 1 × 120 2 2 × 120 3 3 × 120 4 4 × 120 0.15 1 1 1 × 240 2 2 × 240 3 3 × 240 4 4 × 240 0.15 2 1 1 × 524 2 2 × 524 3 3 × 524 4 4 × 524 0.15 4 1 1 × 996 2 2 × 996 3 3 × 996 4 4 × 996 1.2 1 3 3 × 264 4 4 × 264 5 5 × 264 6 6 × 264 7 7 × 264 8 8 × 264 1.2 2 1 1 × 944 2 2 × 944 3 3 × 944 4 4 × 944 5 5 × 944 1.2 4 1 1 × 3824 2 2 × 3824 3 3 × 3824	14.84 12.55 7.54 7.51 8.59 5.30 5.07 5.03 6.08 4.74 4.74 4.74 6.05 4.64 4.64 4.64 20.22 4.66 2.44 1.69 1.64 1.63 12.00 11.94 2.42 1.52 1.50 3.36 1.94 1.47	71.13 12.55 28.37 7.51 23.95 5.30 7.00 5.03 11.82 4.74 4.88 4.74 4.85 4.64 4.80 4.64 97.66 4.66 15.28 1.69 4.69 1.63 122.26 12.62 7.01 1.52 1.69 24.34 1.94 1.49

Table 7 :

 7 Sound-hard submarine-like scatterer: RMS error on the computational domain Γ h for the T6 finite element method for kD = 20 and θ inc = 225 degrees.We have conducted in this paper a numerical study to analyze the performance of coupling a high-order artificial boundary condition and conjugate and unconjugate plane wave finite m n λ n q

	m n λ #dof Err. on Γ h in %
	1 2	2582	26.7
	4 10102	5.25
	6 22410	3.82
	8 38544	3.13
	2 2	3942	27.95
	4 15493	4.04
	6 34037	2.46
	8 60622	1.85
	5 Conclusion

Table 8 :

 8 Sound-hard submarine-like scatterer: RMS error on the computational domain Γ h for CPWT6 and PWT6 finite elements for kD = 20 and θ inc = 225 degrees.

	element prescribed accuracy in % m n λ n q #dof
	T6 CPWT6	≤ 5.5%	1 4 1	10102 5 3505
	PWT6		1	6 4206
	T6 CPWT6	≤ 3.2%	1 8 2 1	38544 5 5205
	PWT6		2 1	6 6246
	T6 CPWT6	≤ 2.1%	2 8 2	60622 4 15768
	PWT6		2	4 15768