Hélène Barucq

Morgane Bergot

Juliette Chabassier

Elodie Estecahandy

Derivation of high order absorbing boundary conditions for the Helmholtz equation in 2D

Keywords: absorbing boundary conditions, Helmholtz equation, high-order approximation, finite elements

We present high order absorbing boundary conditions (ABC) for the Helmholtz equation in 2D, that can adapt to any regular shaped surfaces. The new ABCs are derived by using the technique of micro-diagonalisation to approximate the Dirichlet-to-Neumann map. Numerical results on different shapes illustrate the behavior of the new ABCs along with high-order finite elements.

Derivation of high order absorbing boundary conditions for the Helmholtz equation in 2D.

Résumé : Nous présentons des conditions aux limites absorbantes (CLA) d'ordre élevé pour l'équation de Helmholtz en 2D, s'adaptant à une surface régulière quelconque. Les nouvelles CLAs sont dérivées grace à la technique de micro-diagonalisation afin d'approcher l'opérateur Dirichlet-to-Neumann. Des résultats numériques sur plusieurs formes d'obstacle et de frontière illustrent le comportement des nouvelles CLAs dans un contexte d'éléments finis.

Mots-clés : conditions aux limites absorbantes, équation de Helmholtz, approximation d'ordre élevé, éléments finis

Numerical simulation of wave propagation raises the issue of dealing with outgoing waves. In most of the applications, the physical domain is unbounded and an artificial truncation needs indeed to be carried out for applying numerical methods like finite element approximations. Adapted boundary conditions that avoid the reflection of outgoing waves and provide a well-posed mathematical problem must then be derived. With ideal boundary conditions, the solution on the new mixed boundary valued problem in the truncated domain would actually be equal to the restriction of the mathematical solution in the unbounded domain. However, such ideal boundary conditions, called "transparent boundary conditions", can be shown to be nonlocal, which leads to dramatic computational overcosts. The seek of local boundary conditions, called "absorbing boundary conditions" (ABC), has been the object of numerous works trying to perform efficient conditions based on different techniques of derivation. Among them, the technique of micro-diagonalisation has been employed to the wave equation and more generally to hyperbolic systems in [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF], leading to a hierarchy of absorbing local boundary conditions based on the approximation of the Dirichlet-to-Neumann map. A comprehensive review of different used strategies and higher order conditions can be found in [START_REF] Tsynkov | Numerical solution of problems on unbounded domains[END_REF]. One desirable property of ABCs is that the reflection of the waves on the artificial boundary generates an error of the same order as the one generated by the spatial discretization inside the domain. The computational effort is thus optimized in terms of modeling and numerical inaccuracies. Moreover, the ABC must fit the artificial boundary chosen by the user of the method.

In the context of high order spatial discretization (spectral finite elements [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF], Interior Penalized Discontinuous Galerkin [START_REF] Agut | Stability analysis of the Interior Penalty Discontinuous Galerkin method for the wave equation ESAIM: Mathematical Modelling and Numerical Analysis[END_REF]), there is nowadays a need for high order ABCs that can adapt on non flat geometries since these methods prove very efficient for capturing arbitrary shaped domains.

The aim of the present work is to develop high order ABCs for the Helmholtz equation, that can adapt to regular shaped surfaces. A classical way of designing ABCs is to use Nirenberg theorem [START_REF] Nirenberg | Lectures on linear partial differential equations[END_REF] on the second order formulation of the Helmholtz equation, which enables us to decompose the operator as a product of two first order operators. Here our approach is to rewrite the Helmholtz equation as a first order system of equations before developing ABCs using M.E. Taylor's micro-diagonalisation method [START_REF] Taylor | Reflection of singularities of solutions to systems of differential equations Comm[END_REF]. Then an asymptotic truncation must be performed in order to make the ABC local, and we will see that the high frequency approximation will lead to more usable ABCs than the one stating that the angle of incidence is small. During the process, while increasing the degree of the pseudo differential operator decomposition along with the order of asymptotic truncation, we retrieve classical ABCs that have been found with other techniques by other authors. For now, we have restricted ourselves to two dimensions of space, but despite the fact that 3D generalization should obviously generate more calculation, no further theoretical difficulties are expected. In the following section, we recall the general approach and steps needed to apply the micro-diagonalization technique on the Helmholtz equation. Then a section is devoted to the derivation of a first family of ABCs based on a rough truncation of the pseudo-differential symbol of the operator. Next, a section aims at designing higher order ABCs using a more complete truncation of the symbol. Throughout the report, each ABC will be illustrated with examples. When not specifically said, theses examples will be modeling the scattering of a plane wave field on a circle obstacle of radius R 1 = 2 (homogeneous Neumann boundary conditions). In this simple configuration, analytic solutions can be computed via Hankel expansion of the solution in the open domain {R ≥ R 1 }, or in the truncated domain {R 1 ≤ R ≤ R 2 } with a concentric and circular ABC imposed on the artificial boundary of radius R 2 . The value of R 2 will be set to several values, testing the sensibility of the ABC with respect to the distance between the obstacle and the artificial boundary. Two frequencies are considered : 0.5 Hz and 5 Hz, in order to test the sensibility of the ABCs with respect to the frequency. A last section will present numerical simulations using the newly designed ABCs in realistic and difficult contexts where analytical computations are no longer available. Since a great variety of parameters will be introduced in this report, we do not aim at providing a thorough test cases zoology, but rather at illustrating the performance and behavior of each ABC. As a matter of fact, a definitive conclusion will not be reached as for the supremacy of the newly designed high order ABCs over other possible techniques.

General approach

The first step of the micro-local approximation is to rewrite the equation in a local coordinate system near the artificial boundary that we will call Σ in the following. The Helmholtz equation written in the same local coordinate system as [START_REF] Antoine | Bayliss-Turkel like radiation conditions on surfaces of arbitrary shape[END_REF] (tangent coordinate s and radial coordinate r to the boundary Σ) reads:

(iω) 2 c 2 u -∂ 2 r u -κ r ∂ r u -h -1 ∂ s h -1 ∂ s u = 0 (2.1)
where κ(s) is the curvature of the boundary Σ, h = 1 + r κ(s) and κ r = h -1 κ(s). In order to involve only first order derivatives, and apply Taylor's theorem, we introduce an auxiliary unknown v such that

iω v + ∂ r u = 0 (2.2)
If we introduce the vectorial unknown U = t (u, v), the system reads

∂ r U = LU (2.3)
where the symbol of the pseudo-differential operator L is given by

L = σ(L) =   0 -iω - iω c 2 -h -3 ∂ s (h) ξ ω - h -2 ξ 2 iω -κ r   (2.4)
This symbol can be decomposed as a sum of two symbol L 1 + L 0 where L 1 is homogeneous of degree 1 and L 0 of degree 0 1 . The principal symbol is the less regularizing part of the operator, here we have:

σ p (L) = L 1 =   0 -iω - iω c 2 - h -2 ξ 2 iω 0   , L 0 = 0 0 - ∂ s h ξ h 3 ω - κ h (2.5)
where ξ is the dual variables associated to s.

Remark 2.1

The symbol is, in general, different from the principal symbol. The order p of a pseudodifferential operator is the greatest homogeneity degree of its symbol. The principal operator is exactly of order p. We will call σ p (A) the principal symbol of an operator A, σ r (A) the rest of order p -1, and σ j (A) the part of the symbol being of order j. Finally, σ * j (A) will denote the symbol of the part of the operator of exactly order j. For commodity, we will denote OP p (A) := OP σ p (A) as well as OP j (A) := OP σ j (A) .

Remark 2.2

Notice that on the boundary (r = 0), we have

∂ 2 r u = iω c 2 u -κ ∂ r u + ξ 2 u (2.6)
Our aim is to find:

1 The lower the degree, the more regularizing the operator.

Inria

a diagonal pseudo-differential operator Λ such that Λ = Λ 1 + Λ 0 + Λ -1 + . . .
where σ(Λ j ) = D j is homogeneous of degree exactly j.

a pseudo-differential operator P such that P = P 0 + P -1 + P -2 + . . .

where σ(P j ) = P j is homogeneous of degree j.

such that V = P U and (2.3) is equivalent to

∂ r V = ΛV (2.7)
We will then end up with a diagonal system, in which the first component of the auxiliary unknown V corresponds to the ingoing wave, while the second component stands for the outgoing wave. Hence, the exact absorbing boundary condition will be obtained by ordering the first component to vanish on Σ. As this condition is non explicit, non local and therefore inconvenient, an approximation process will be done step by step : our approach will consist in describing P as a factor of correcting terms of decreasing order. Each step will allow us to write an "approximate ABC" involving only first order derivatives.

Remark 2.3

Recall that the symbol of a product of operators (see [START_REF] Hörmander | Pseudodifferential operators and hypoelliptic equations[END_REF]) is computed as:

σ(AB) = k∈N (-i) k k! ∂ k s σ(A)∂ k ξ σ(B) = σ(A)σ(B) -i∂ s σ(A)∂ ξ σ(B) - 1 2 ∂ 2 s σ(A)∂ 2 ξ σ(B) + . . . (2.8)
We follow a three-step approach that we describe now.

Step 1. The symbol L 1 admits two eigenvalues ±λ 1 with

λ 1 = h -2 ξ 2 -c -2 ω 2 .
(2.9)

The frequencies (ω, ξ) are in R 2 and the square root is defined according to the sign of h -2 ξ 2 -c -2 ω 2 . If (ω, ξ) belongs to the hyperbolic region, h -2 ξ 2 -c -2 ω 2 is negative and λ 1 matches a propagating wave. In that case, λ 1 should be written as where

λ 1 = i c -2 ω 2 -h -2 ξ 2 (2.10)
In the event of h -2 ξ 2 -c -2 ω 2 being positive, the frequencies (ω, ξ) follow the elliptic region and λ 1 matches evanescent waves. In the following, we restrict our study to the case where (ω, ξ) belong to the hyperbolic region or "propagating cone", letting λ 1 being purely imaginary.

We introduce the change of basis operator P 0 such that

σ(P 0 ) = P 0 = 1 √ 2    - λ 1 iω 1 1 iω λ 1    , P -1 0 = 1 √ 2    -iω λ 1 1 1 λ 1 iω    (2.11)
We have:

P 0 L 1 P -1 0 = D 1 = λ 1 0 0 -λ 1 ⇒ L 1 = P -1 0 D 1 P 0 (2.12)
We denote

V 0 = P 0 U (2.13) Then ∂ r V 0 = ∂ r (P 0 U ) (2.14) = (∂ r P 0 ) U + P 0 ∂ r U (2.15) = (∂ r P 0 )P -1 0 V 0 + P 0 LP -1 0 P 0 L 1 P -1 0 + P 0 L 0 P -1 0 V 0 (2.16) = D 1 V 0 + R 0 V 0 (2.17)
where

RR n°8632
D 1 is the principal symbol of the operator P 0 L 1 P -1 0 , R 0 is the operator of order 0 having symbol σ r (P 0 L 1 P -1 0 ) + σ((∂ r P 0 )P -1 0 ) + σ(P 0 L 0 P -1 0 )

Step 2. We seek

V 1 such that ∂ r V 1 = (D 1 + D 0 ) V 1 + R -1 V 1 (2.18)
We set V 1 under the following form:

V 1 = (1 + K -1 )P 0 U (2.19)
and we are going to construct K -1 such that D 0 is a diagonal operator of order 0 exactly. Operator K -1 is of order -1 and can be determined so that (2.18) is true.

A simplest ABC can be found (as in [START_REF] Barucq | Mathematical Models and Methods in Applied Sciences[END_REF]) by stating that the first component of V 1 must vanish on the boundary (when r = 0) and taking the symbol of order -1.

Step 3. We now seek V 2 such that

∂ r V 2 = (D 1 + D 0 + D -1 ) V 2 + R -2 V 2 (2.20)
We set V 2 under the following form:

V 2 = (1 + K -2 )(1 + K -1 )P 0 U (2.21)
such that D -1 is a diagonal operator of order -1 and K -2 is an operator of order -2 that we will determine so that (2.20) is true.

When all the operators K j and D j have an explicit expression, we will state that the first component of V 2 must vanish on the boundary, which will lead us to a second approximate boundary condition after taking the symbol of order -2. Obviously, the process can be iterated as long as necessary to design higher order conditions (see [START_REF] Antoine | Microlocal Diagonalization of Strictly Hyperbolic Pseudodifferential Systems and Application to the Design of Radiation Conditions in Electromagnetism[END_REF]).

A first family of ABCs

Derivation of K -1

By definition,

∂ r V 0 = D 1 V 0 + R 0 V 0 where R 0 = (∂ r P 0 )P -1 0 + OP σ r (P 0 L 1 P -1 0 ) + σ(P 0 L 0 P -1 0 ) (3.1)
We will need the commutator notation :

{A ; B} = AB -BA. ∂ r V 1 = ∂ r (I + K -1 )V 0 (make V 0 appear) = ∂ r (I + K -1 )V 0 + (I + K -1 )∂ r V 0 = ∂ r (K -1 )V 0 + (I + K -1 ) D 1 + R 0 V 0 (use the knowledge on V 0 ) = ∂ r (K -1 )(I + K -1 ) -1 V 1 + D 1 + R 0 V 1 + {(I + K -1 ) ; D 1 + R 0 } (I + K -1 ) -1 V 1 (return to V 1 to identify terms) = D 1 V 1 + D 0 V 1 + R -1 V 1
where by identification,

D 0 = OP p ∂ r (K -1 )
order -1

(I + K -1 ) -1
order 0 order -1

+R 0 + {I + K -1 ; D 1 + R 0 } {K -1 ; D 1 } order 0 + {K -1 ; R 0 } order -1 (I + K -1 ) -1 (3.2) = OP p (R 0 ) + OP p (K -1 D 1 -D 1 K -1 ) (3.3)
by writing that (I + K -1 ) -1 = I -K -1 + . . . At this point, we can give a necessary condition on K -1 so that D 0 is diagonal. Indeed, we can use formula 2.8 and only the first product will be of order 0:

σ p (K -1 D 1 -D 1 K -1 ) = K -1 D 1 -D 1 K -1 = 0 -2λ 1 (K -1 ) 1,2 2λ 1 (K -1 ) 2,1 0 (3.4)

Inria

The operator D 0 will be diagonal if and only if the extradiagonal terms of the sum of this quantity with σ p (R 0 ) = σ * 0 (R 0 ) are zero. This will give a "simple" relation between (K -1 ) 1,2 , (K -1 ) 2,1 and (ω, ξ), involving the curvature and its derivatives. This calculation will be done hereafter.

Determination of a first ABC

We will write that the first component of V 1 must vanish on the boundary (when r = 0):

V 1 = (1 + K -1 )P 0 t (u, v) (3.5)
Since only the extradiagonal terms of K -1 are constrained by the fact that D 0 is diagonal, we can take the diagonal terms we want. We choose to consider γ and ζ such that

K -1 =     γ(s) λ 1 σ * 0 (R 0 ) 1,2 2λ 1 - σ * 0 (R 0 ) 2,1 2λ 1 ζ(s) λ 1     (3.6)
where R 0 was given in equation (3.1) and P 0 in equation (2.11). It gives:

1 + γ λ 1 - λ 1 iω u + v + σ * 0 (R 0 ) 1,2 2λ 1 u + iω λ 1 v = 0 (3.7)
We use result (A.18) of appendix A to evaluate σ * 0 (R 0 ) 1,2 for r = 0:

σ * 0 (R 0 ) 1,2 (r = 0) = κω 2iλ 1 c 2 (3.8)
We get the condition:

1 + γ λ 1 - λ 1 iω u + v + κω 4iλ 2 1 c 2 u + iω λ 1 v = 0 (3.9)
Moreover, we recall that, from (2.3)

∂ r u = -iω v
We obtain the following non local absorbing boundary condition:

1 + γ λ 1 - λ 1 iω u - ∂ r u iω + κω 4iλ 2 1 c 2 u - ∂ r u λ 1 = 0 (3.10)

Asymptotic approximation: small "Angle of incidence"

A first approach is to introduce the "angle of incidence" δ = ξ/ω and suppose it is small. Retaining the first terms of the approximation gives:

u - 1 c - γ(s) iω - κ 4iω + c 2 - κ 4iω δ 2 + ∂ r u - 1 iω + cγ(s) ω 2 - κc 4ω 2 + c 3 γ(s) 2ω 2 - 3κc 3 8ω 2 δ 2 + O(δ 4 ) = 0 (3.11)
A zeroth & first order condition. By neglecting the terms in O(δ) and in O(δ 2 ), we obtain the same condition:

u - 1 c - γ(s) iω - κ 4iω + ∂ r u - 1 iω + cγ(s) ω 2 - κc 4ω 2 = 0 (3.12)
We can now multiply by -iω to get:

u iω c + γ(s) + κ 4 + ∂ r u 1 + cγ(s) iω - κc 4iω = 0 (3.13)
The first order ABC reads in frequency domain:

k = ω c , (∂ r u + ik u) + γ + κ 4 u + γ - κ 4 ∂ r u ik = 0 (3.14)
Fig. 2 shows the relative L2 error of the trace of the solution on the obstacle with respect to the parameter γ/κ. We have chosen to evaluate this quantity because it is independent of the size of the computational domain. We indeed modify it by increasing the radius of the external boundary. If we had considered the L 2 -norm of the numerical solution inside the whole computational box, we would have compared incomparable things. As said in the introduction, analytical computations are made for this circular configuration, in an open domain and in the truncated domain with ABC (3.14) on the artificial boundary. Left and right figures correspond to frequencies 0.5 and 5 Hz. Different values of the artificial boundary radius R 2 are considered and displayed in different color curves. The legend indicates the value of R 2 /R 1 . It is clear from these figures that the ABC performs better when the artificial boundary is far from the obstacle. No specific value of the parameter γ seems to provide better results. The ABC performs as well for both frequencies.

Remark 3.1

Notice that the specific value γ(s) = κ(s)/4 leads to the following simplifications:

∂ r u + ik u + κ 2 u = 0 (3.15)
which is (the sign of k is sometimes inverted in the definition) the so called "curvature-ABC" or "C-ABC" well known in the literature ( [START_REF] Stupfel | Absorbing boundary conditions on arbitrary boundaries for the scalar and vector wave equations[END_REF][START_REF] Jones | Surface radiation conditions[END_REF][START_REF] Hanouzet | Absorbing boundary conditions for Maxwell's equations[END_REF][START_REF] Bayliss | Boundary conditions for the numerical solution of elliptic equations in exterior regions[END_REF]). Notice also that a coefficient 2 appears in 3d because of the 2d/3d difference of definition of κ. Since the variable r is in fact the normal oriented coordinate, we obtain after a invert Fourier transform in time the following ABC:

1 c ∂ t ∂ n u + 1 c ∂ t u = κ(x) 4 -γ(x) ∂ n u - 1 c κ(x) 4 + γ(x) ∂ t u , x ∈ Σ (3.16)
This condition was found with the same reasoning in [START_REF] Barucq | Mathematical Models and Methods in Applied Sciences[END_REF].

A second order condition. We retain all terms up to O(δ 2 ), we multiply by -iω and we replace δ with its value. We get:

k = ω c ,   1 + γ - κ 4 ik   ∂ r u -    γ - 3κ 4 2(ik) 3    ξ 2 ∂ r u + ik + γ + κ 4 u + 1 2ik - κ 4(ik) 2 ξ 2 u = 0 (3.17) Inria
This equation involves ξ 2 u and ξ 2 ∂ r u terms, which are in the primal variables the Laplace-Beltrami operator applied to u and ∂ r u. This means that to implement this condition, one has to evaluate the Laplace-Beltrami of u, which implies more computational costs, but also to invert a surfacic operator of second order to retrieve ∂ r u. Fig. 3 shows the relative L2 error for ABC (3.17) in the same conditions as Fig. 2. As expected by the theory, the performance of the ABC improves as the frequency increases. We can also observe that here the specific value γ = κ/2 seems to give much better results, especially when the boundary is far from the obstacle.

Asymptotic approximation : "High frequency"

In this second approach, we want to truncate equation (3.10) when ω tends to infinity, ξ being considered as a parameter (not involved in the asymptotic). Let us recall equation (3.10):

u κω 4iλ 2 1 c 2 - λ 1 + γ iω -∂ r u λ 1 + γ λ 1 iω + κω 4iλ 3 1 c 2 = 0 (3.18)
We get

u - 1 c + c 2 ξ 2 2ω 2 + O 1 ω 4 - γ iω - κ 4iω + c 2 κξ 2 4iω 3 + O 1 ω 4 -∂ r u 1 iω + γc (iω) 2 + O 1 ω 4 + κωc 3 4ic 2 (iω) 3 + O 1 ω 4 = 0 (3.19)
We can now multiply by -iω to get:

u iω c + γ + κ 4 + c 2 ξ 2 2iω + c 2 κ 4ω 2 + ∂ r u 1 + γc iω - κc 4iω = O 1 ω 3 (3.20)
A zeroth order condition. Neglecting terms in O 1 ω gives the following zeroth order ABC:

k = ω c , ∂ r u + ik u + γ + κ 4 u = 0 (3.21) Remark 3.3
The specific value γ(s) = κ(s) 4 also leads to the "C-ABC" found in [START_REF] Stupfel | Absorbing boundary conditions on arbitrary boundaries for the scalar and vector wave equations[END_REF][START_REF] Jones | Surface radiation conditions[END_REF][START_REF] Hanouzet | Absorbing boundary conditions for Maxwell's equations[END_REF][START_REF] Bayliss | Boundary conditions for the numerical solution of elliptic equations in exterior regions[END_REF].

A first order condition. Neglecting terms in O 1 ω 2 gives the following first order ABC:

k = ω c , ∂ r u + ik u + γ(s) + κ(s) 4 u + ξ 2 2 ik u + γ(s) - κ(s) 4 ∂ r u ik = 0 (3.22)
which can be written in the following form

k = ω c ,   1 + γ - κ 4 ik   ∂ r u + ik + γ + κ 4 u + ξ 2 2 ik u = 0 (3.23)
Fig 4 shows the relative L2 error for ABC (3.23) in the same conditions as Fig. 2. The performance of the ABC seems to improve as the frequency increases. Here again the specific value γ = κ/2 seems to give much better results, especially when the boundary is far from the obstacle.

Remark 3.4

Notice that the specific value γ(s) = κ(s) 4 leads to the well known ABC formulated in the primal variables (ξ 2 ↔ -∆ ⊥ ): It is not equivalent to multiply (3.10) by λ 1 and then perform the high frequency asymptotic. Indeed, in this case, the equation becomes:

∂ n u + iku + κ(s) 2 u - 1 2ik ∆ ⊥ u = 0 (3.
-∂ n u λ 1 + γ iω + κω 4iλ 2 1 c 2 + u κω 4iλ 1 c 2 -(λ 1 + γ) λ 1 iω = 0 (3.25) -∂ n u 1 c + c iω γ(s) - κ(s) 4 + O 1 ω 2 + u - γ(s) c - κ(s) 4c - iω c 2 - ξ 2 iω + O 1 ω 2 = 0 (3.26)
Multiplying by c and neglecting all terms in O( 1 ω 2 ) we get the following ABC:

u iω c + γ(s) + κ(s) 4 + cξ 2 iω + ∂ n u 1 + c iω γ(s) - κ(s) 4 = 0 (3.27)
which differs from (3.23) by a factor 1 2 in the last term in ξ 2 .

Inria

A second order condition. By neglecting terms in O 1 ω 3 , we obtain the following second order ABC: The obtained ABCs (3.23) and (3.28) involve the Laplace-Beltrami operator on u but not on ∂ r u (as opposed to the ABC (3.17)), which avoids the inversion of a second order operator on the boundary.

k = ω c ,   1 + γ - κ 4 ik   ∂ r u + ik + γ + κ 4 u + 1 2 ik - κ 4 (ik) 2 ξ 2 u = 0 (3.28) - 5

Remark 3.6

The reader can notice the fact that the ABC (3.17), (3.23) and (3.28) show a super convergence behavior when the parameter γ is chosen equal to κ/2. The last one can be explained by the remark 4.2.

A second ABC

Calculation of K -2

With the same approach as in 3.1, we find that

∂ r V 2 = (∂ r K -2 )(I + K -1 ) -1 V 2 + (D 1 + D 0 + D -1 ) V 2 + [{K -2 ; D 1 } + {K -2 ; D 0 } + {K -2 ; D -1 }] (I + K -2 ) -1 V 2 (4.1)
and by definition, on the other hand,

∂ r V 2 = (D 1 + D 0 + D -1 ) V 2 + R -2 V 2 (4.2)
Identifying terms gives

D -1 = σ * -1 (R -1 ) + K -2 D 1 -D 1 K -2 (4.3)
As before,

K -2 D 1 -D 1 K -2 = 0 -2λ 1 (K -2 ) 1,2 2λ 1 (K -2 ) 2,1 0 (4.4)
We can give a necessary condition on K -2 so that D -1 is diagonal : the sum of this quantity with σ * -1 (R -1 ) must have null extradiagonal terms.

Determination of a second ABC

We deduce from above considerations that

K -2 =    θ(s) λ 1 σ -1 (R -1 ) 1,2 2λ 1 - σ -1 (R -1 ) 2,1 2λ 1 0    (4.5)
We write that the first component of V 2 must vanish on the boundary (when r = 0):

V 2 = (1 + K -2 )(1 + K -1 )P 0 t (u, v) (4.6)
We truncate the operators of order lower than -3, and (4.6) becomes

Ṽ2 =OP -2 ((1 + K -2 )(1 + K -1 )P 0 ) t (u, v) = ((1 + K -2 )(1 + K -1 )P 0 -i∂ s K -1 ∂ ξ P 0 ) t (u, v)
which reads (see equation (C.5)):

1 + θ(s) λ 1 1 + γ(s) λ 1 - λ 1 u iω + v + σ 12 2λ 1 u + iω λ 1 v + σ 12 2λ 1 - σ 21 2λ 1 - λ 1 iω u + v + 1 + ζ(s) λ 1 u + iω λ 1 v + γ (s)ξ ωλ 2 1 u + κ (s)ξiω 2 4λ 5 1 v = 0 (4.7)
where

σ 12 = σ * 0 (R 0 ) 12 , σ 21 = σ * 0 (R 0 ) 21 and σ 12 = σ * -1 (R -1 ) 12 .
According to (A.18) and (B.16), we have:

               σ 12 = κω 2iλ 1 c 2 (4.8a) σ 21 = iκω 3 2λ 3 1 c 2 (4.8b) σ 12 = κ 2 ω 4iλ 4 1 c 2 ω 2 c 2 + 2ξ 2 + κ ξω 2λ 3 1 c 2 + κωγ(s) 2iλ 2 1 c 2 + κ ξ 2ωλ 1 (4.8c)
We recall that, from (2.3)

∂ r u = -iω v
We obtain the following non local absorbing boundary condition

1 + θ(s) λ 1 1 + γ(s) λ 1 - λ 1 u iω - ∂ r u iω + σ 12 2λ 1 u - ∂ r u λ 1 + σ 12 2λ 1 σ 21 2λ 1 λ 1 iω u + ∂ r u iω + 1 + ζ(s) λ 1 u - ∂ r u λ 1 + γ (s)ξ ωλ 2 1 u - κ (s)ξω 4λ 5 1 ∂ r u = 0 (4.9)

Asymptotic approximation: small "Angle of incidence"

In that section we choose the "angle of incidence" δ = ξ/ω to be the small parameter. We then get conditions by applying Taylor expansions to each of the symbols.

A zeroth order condition.

If we neglect the terms in O(δ), we get A first order condition.

k = ω c ,   1 + θ + γ - κ 4 ik + θγ -θ + γ - κ 2 κ 4 (ik) 2 - γ - κ 2 ζ + κ 4 κ 4 (ik) 3   ∂ r u +   ik + γ + θ + κ 4 + θγ + θ + γ - κ 2 κ 4 ik + γ - κ 2 ζ - κ 4 κ 4 (ik) 2   u = 0 (4.
If we neglect the terms in O(δ 2 ), we get

k = ω c ,   1 + θ + γ - κ 4 ik + θγ -θ + γ - κ 2 κ 4 (ik) 2 - γ - κ 2 ζ + κ 4 κ 4 (ik) 3   ∂ r u - κ 4(ik) 3 iξ ∂ r u +   ik + γ + θ + κ 4 + θγ + θ + γ - κ 2 κ 4 ik + γ - κ 2 ζ - κ 4 κ 4 (ik) 2   u - γ (ik) 2 iξ u = 0 (4.11)
We can notice that for κ (s) = 0, (4.11) is equal to (4.10). Therefore the two ABCs are equal for a circular artificial boundary. For this test case, we will thus only present the L2 relative errors corresponding to the condition (4.10) in Fig. 4.11. Since the ABC depends on three parameters, we present the error on a grid (radius ratio, ζ) where each subfigure is a small (θ,γ) map. The error, which is displayed in color scale, highly diminishes as the radius ratio increases. However, all the subfigures are very different, meaning that the ABC's performance depends very strongly on the parameter's values. The main drawback is that this dependency is not consistent when the radius ratio growth. In practice, the obstacle can be more complex than a circle, thus the distance between the obstacle and the boundary is not a constant value throughout the domain. In Fig. 6 we can see that a fixed parameter set cannot be a good choice everywhere if this distance varies in the domain. Moreover, the presence of the term ξ 2 ∂ r u leads to the necessary inversion of a Laplace-Beltrami operator on the boundary which generates additional computational burden.

A second order condition.

If we neglect the terms in O(δ 3 ), we get

k = ω c ,   1 + θ + γ - κ 4 ik + θγ -θ + γ - κ 2 κ 4 (ik) 2 - γ - κ 2 ζ + κ 4 κ 4 (ik) 3   ∂ r u - κ (s) 4(ik) 3 iξ ∂ r u +     - θ + γ - 3κ 4 2(ik) 3 + κ 2 γ - 5κ 4 -θ γ - κ 2 (ik) 4 + κ 4 κ 4 7γ - 13κ 2 + ζ 5γ - 11κ 2 2(ik) 5     ξ 2 ∂ r u +   ik + γ + θ + κ 4 + θγ + θ + γ - κ 2 κ 4 ik + γ - κ 2 ζ - κ 4 κ 4 (ik) 2   u - γ (s) (ik) 2 iξ u +     ik - κ 2 2(ik) 2 - 3κ 4 γ - 3κ 2 + θ γ + 3κ 4 2(ik) 3 + 3κ 2 8 (γ -κ) -κζ γ - 5κ 4 2(ik) 4     ξ 2 u = 0 (4.12)
In Fig. 7 the L2 relative errors for the ABC (4.12) are displayed in the same conditions as Fig. (4.10). Condition (4.12) is a complete condition in the sense that the truncation error of the Taylor expansion is the same as the order of the pseudo-differential operator K -2 . That explains the very low errors obtained with this condition, especially for high frequency. As for ABC (4.11), the subfigures are very different, meaning that it will be difficult to choose a parameter set which will give good results for any situation. Moreover, the ABC involves ξ 2 ∂ r u, meaning that a second order operator must be inverted on the boundary.

Remark 4.1

Recall that ξ is the dual variable associated to the curvilinear coordinate s. If we denote ↔ the corresponding between primal (physical) and dual domains, we can remind the reader that: 4 -

-iξ ↔ ∂ s and -ξ 2 ↔ ∂ 2 s Inria Let us call                                            a(s) = γ (s) (ik) 2 (4.13a) A(s) =     - 1 2ik + κ 4(ik) 2 + 3κ 4 γ - 3κ 2 + θ γ + 3κ 4 2(ik) 3 - 3κ 2 8 (γ -κ) -κζ γ - 5κ 4 2(ik) 4     (4.13b) b(s) = κ (s) 4(ik) 3 (4.13c) B(s) =     θ + γ - 3κ 4 2(ik) 3 - κ 2 γ - 5κ 4 -θ γ - κ 2 (ik)
κ 4 κ 4 7γ - 13κ 2 + ζ 5γ - 11κ 2 2(ik) 5     (4.13d) Therefore,          -a(s) iξ u ←→ a(s)∂ s u (4.14a) -A(s)ξ 2 u ←→ A(s)∂ 2 s u (4.14b) -b(s) iξ∂ r u ←→ b(s)∂ s ∂ r u (4.14c) -B(s)ξ 2 ∂ r u ←→ B(s)∂ 2 s ∂ r u (4.14d)
In the prospect of a variational numerical method, we must be able to write these terms in a symmetric way:

-a(s) iξu -A(s)ξ 2 u ←→ ∂ s (A(s) ∂ s u) (4.15)
It is then necessary to have ∂ s A(s) = a(s) and ∂ s B(s) = b(s). Since these equalities must be true for any k, we must adequate the terms in increasing power of (ik). Consider the part in u (A and a): we must have

∂ s κ 4(ik) 2 = γ (s) (ik) 2 , ∂ s     3κ 4 γ - 3κ 2 + θ γ + 3κ 4 2(ik) 3     = 0, ∂ s     3κ 2 8 (γ -κ) -κζ γ - 5κ 4 2(ik) 4     = 0
Therefore there exists constants α 2 , α 3 and α 4 independent of s such that:

γ(s) = κ(s) 4 + α 2 , 3κ 4 γ - 3κ 2 + θ γ + 3κ 4 = α 3 , 3κ 2 8 (γ -κ) -κζ γ - 5κ 4 = α 4 (4.16) γ(s) = κ(s) 4 + α 2 , θ(s) = α 3 - 3κ 4 α 2 - 5κ 4 κ + α 2 , ζ(s) = -α 4 + 3κ 2 8 α 2 - 3κ 4 κ (α 2 -κ) (4.17)
Using these values, let us consider now the part in ∂ r u (B and b): we must have

∂ s    θ + γ - 3κ 4 2(ik) 3    = κ (s) 4(ik) 3 , ∂ s     κ 2 γ - 5κ 4 -θ γ - κ 2 (ik) 4     = 0, ∂ s     κ 4 κ 4 7γ - 13κ 2 + ζ 5γ - 11κ 2 2(ik) 5     = 0
The first equation leads to the existence of a constant β 3 such that

α 3 -3κ 4 α 2 -5κ 4 κ + α 2 + κ 4 + α 2 - 3κ 4 = κ 2 + β 3 ⇒ α 3 + α 2 2 -β 3 α 2 = κ 3α 2 4 + β 3 - κ 2 16
The left-hand size of this equation is then equal to a quantity that depends on s, which is in contradiction with the hypothesis. This ABC is thus not usable in a variational context as soon as the artificial boundary is not circular. We now perform a high frequency approximation in the sense that now 1/ω is the small parameter. We get the following conditions

A zeroth order condition.

Neglecting terms in O 1 ω gives the following zeroth order ABC:

k = ω c , ∂ r u + ik u + γ + θ + κ 4 u = 0 (4.18)
The obtained ABC is the same as (3.21) and leads to the "curvature-ABC" when θ + γ = κ/4. Its behavior is hence the same and we do not plot the error curve again.

A first order condition.

Neglecting terms in O 1 ω 2 gives the following first order ABC:

k = ω c ,   1 + θ + γ - κ 4 ik   ∂ r u +   ik + γ + θ + κ 4 + θγ + γ + θ - κ 2 κ 4 ik   u + ξ 2 2 ik u = 0 (4.19)
Fig. 8 shows the L2 relative errors on the obstacle in color scale, in (θ,γ) planes. The dependency on the radius ratio between the obstacle and the artificial boundary is displayed in abscissa. The ABC gives better results in high frequency.

A second order condition.

Neglecting terms in O 1 ω 3 gives the following second order ABC:

k = ω c ,   1 + θ + γ - κ 4 ik + θ γ - κ 4 + κ 4 κ 2 -γ (ik) 2   ∂ r u +   ik + γ + θ + κ 4 + θ γ - κ 4 + κ 4 γ - κ 2 ik + γ - κ 2 ζ - κ 4 κ 4 (ik) 2   u - γ (ik) 2 iξ u + 1 2 ik - κ 4 (ik) 2 ξ 2 u = 0 (4.20)
Fig. 9 shows the L2 errors obtained with the ABC (4.20) in the same conditions as Fig. 6. The first noticeable property is the very low dependance of the ABC's performance on the parameter ζ. Moreover, all subfigures offer the same kind of behavior, meaning that a given parameter set will provide results of the same quality for different situations. This robustness is a very nice feature to this ABC. The obtained condition (4.20) is to be compared by the reader with (4.12) since it results from a different asymptotic truncation, at the same order, of the same non local ABC (4.9). It is however much nicer since no Laplace-Beltrami operator needs to be inverted on the boundary (no term ξ 2 in front of ∂ r u). Moreover it involves lower orders of ik, meaning that the ABC may also be nicer in time domain. We observe that the ratio performance/computational burden is very interesting with the ABC (4.20). In the prospect of using this ABC in the time domain, one can wonder if the parameters γ, θ and ζ could be optimized in order to simplify its expression. Especially, we know that the terms involving ∂ r u divided by ik once or more times will require the solution of auxiliary systems on the artificial boundary. A natural idea is then to seek parameters values that make the coefficients before these terms vanish. The specific choice γ = 0 and θ = κ 2 , ζ = κ 4 (4.21) leads to the following ABC :

∂ r u + ik u + 3κ 4 u + κ 4 ∂ r u ik + ξ 2 u 2 ik 1 - κ 2 ik = 0 (4.22)
which is a complete second order ABC. Notice that the latter can be seen as the first order ABC in equation (3.28) with the specific choice γ = κ/2. This explains why this value leads to better results than the other ones.

A second order symmetric condition.

Recall that ξ is the dual variable associated to the curvilinear coordinate s. If we denote ↔ the corresponding between primal (physical) and dual domains, we can remind the reader that:

-iξ ↔ ∂ s -ξ 2 ↔ ∂ 2 s Therefore, - γ (s) (ik) 2 iξ u ←→ γ (s) (ik) 2 ∂ s u (4.23) - κ(s) 4(ik) 2 ξ 2 u ←→ κ(s) 4(ik) 2 ∂ 2 s u (4.24)
In the prospect of a variational numerical method, the only suitable choice is to take

γ = K + κ 4
, where K is constant with respect to s, so that

- iξκ (s) 4(ik) 2 u -K + κ(s) 4 1 (ik) 2 ξ 2 u ←→ ∂ s K + κ(s) 4 1 (ik) 2 ∂ s u (4.25)
which corresponds to a symmetric term in a variational context. The resulting ABC reads:

k = ω c ,   1 + θ + K ik + θK + κ 4 κ 4 -K (ik) 2   ∂ r u +   ik + θ + K + κ 2 + θ K + κ 2 + κ 4 K - κ 4 ik - K - κ 4 κ 4 ζ - κ 4 (ik) 2   u - κ (s) 4(ik) 2 iξu + 1 2 ik -K + κ 4 1 (ik) 2 ξ 2 u = 0 (4.26)
This ABC thus writes under the following variational form, where ϕ is a test function: (4.26). The color scale represents the radius ratio between the artificial boundary and the obstacle (from 1.05 to 500). It can be seen that the ABC's performance does not really depend on the parameter ζ, which can therefore be taken equal to κ/4 in order to make the 1/(ik) 2 contribution disappear in the u ϕ term. The dependance in θ is less pronounced than in previously found ABCs but the value θ = κ/2 seems to be slightly better, especially when the artificial boundary is far from the obstacle. 5 To sum up...

  1 + θ + K ik + θK + κ 4 κ 4 -K (ik) 2   ∂ r u ϕ +   ik + θ + K + κ 2 + θ K + κ 2 + κ 4 K - κ 4 ik - K - κ 4 κ 4 ζ - κ 4 (ik) 2   u ϕ + 1 2 ik -K + κ 4 1 (ik) 2 ∂ s u ∂ s ϕ = 0 (4.
All ABC can be written under the following form

a 0 (ik) + a 1 (ik) iξ + a 2 (ik) ξ 2 ∂ r u + b 0 (ik) + b 1 (ik) iξ + b 2 (ik) ξ 2 u = 0 (5.1)
The coefficients for each ABC are summarized in Table 1 for the ABC of section 3 (with a 1 (ik) = b 1 (ik) = 0) and in Table 2 for the ABC of section 4 (with a 2 (ik) = 0).

a 0 (ik) b 0 (ik) a 2 (ik) b 2 (ik) (3.14) Km1 δ 0 , δ 1 1 + γ -κ/4 ik ik + γ + κ 4 0 0 (3.15) Km1 δ 0 , δ 1 (γ = κ/4) 1 ik + κ 2 0 0 (3.17) Km1 δ 2 1 + γ -κ/4 ik ik + γ + κ 4 - γ -3κ/4 2(ik) 3 1 2ik - κ/4 (ik) 2 Km1 δ 2 (γ = κ/2) 1 + κ 4ik ik + 3κ 4 κ/4 (ik) 3 1 2ik - κ/4 (ik) 2 (3.21) Km1 ω 0 1 ik + γ + κ 4 0 0 Km1 ω 0 (γ = κ/4) 1 ik + κ 2 0 0 (3.23) Km1 ω 1 1 + γ -κ/4 ik ik + γ + κ 4 0 1 2ik (3.24) Km1 ω 1 (γ = κ/4) 1 ik + κ 2 0 1 2ik Km1 ω 1 (γ = κ/2) 1 + κ 4ik ik + 3κ 4 0 1 2ik (3.28) Km1 ω 2 1 + γ -κ/4 ik ik + γ + κ 4 0 1 2ik - κ/4 (ik) 2 (4.22) Km1 ω 2 (γ = κ/2) 1 + κ 4ik ik + 3κ 4 0 1 2ik - κ/4 (ik) 2
Table 1: Summary of the coefficients for the ABCs of section 3 Frequency is taken equal to 0.5 Hz and 5 Hz. Solutions are computed with the same numerical tool as all previously displayed figures i.e. with Hankel expansion. Each curve corresponds to the use of a different ABC (see the legend). It is clear on the figures that only complete ABCs converge to order 4 : Km2 delta 2 and Km2 omega 2. All other ABCs eventually exhibit a second order rate of convergence. This means that performing a second order Taylor expansion on conditions coming from the first term of the pseudo-differential operator decomposition does not lead to complete ABCs. However, the convergence regime is sometimes reached only for very large radius ratio, especially for the case of a high frequency. This illustrates the fact that the newly designed ABCs will be competitive with the usual ones for low frequencies, and this is consistent with the fact that the Taylor expansion has been pushed one term ahead. In the legend, the indicated ABCs are followed by the relative L2 error percentage of the first point of the curve.

Inria a 0 (ik) b 0 (ik) a 1 (ik) b 1 (ik) b 2 (ik) (4.10) Km2 δ 0 1 + θ + γ - κ 4 ik ik + γ + θ + κ 4 0 0 0 + θγ -θ + γ - κ 2 κ 4 (ik) 2 + θγ + θ + γ - κ 2 κ 4 ik - γ - κ 2 ζ + κ 4 κ 4 (ik) 3 + γ - κ 2 ζ - κ 4 κ 4 (ik) 2 Km2 δ 0 1 + κ 4ik ik + 3κ 4 0 0 0 (γ = κ 2 , θ = 0) (4.11) Km2 δ 1 1 + θ + γ - κ 4 ik ik + γ + θ + κ 4 - κ 4(ik) 3 - γ (ik) 2 0 + θγ -θ + γ - κ 2 κ 4 (ik) 2 + θγ + θ + γ - κ 2 κ 4 ik - γ - κ 2 ζ + κ 4 κ 4 (ik) 3 + γ - κ 2 ζ - κ 4 κ 4 (ik) 2 Km2 δ 1 1 + κ 4ik ik + 3κ 4 - κ 4(ik) 3 - κ 2(ik) 2 0 (γ = κ 2 , θ = 0, ζ = κ 4 ) (4.18) Km2 ω 0 1 ik + γ + θ + κ 4 0 0 0 Km2 ω 0 1 ik + 3κ 4 0 0 0 (γ = κ 2 , θ = 0) (4.19) Km2 ω 1 1 + θ + γ - κ 4 ik ik + γ + θ + κ 4 0 0 1 2ik + θγ + θ + γ - κ 2 κ 4 ik Km2 ω 1 1 + κ 4ik ik + 3κ 4 0 0 1 2ik (γ = κ 2 , θ = 0, ζ = κ 4 ) (4.20) Km2 ω 2 1 + θ + γ - κ 4 ik ik + γ + θ + κ 4 0 - γ (ik) 2 1 2ik - κ 4(ik) 2 + θγ -θ + γ - κ 2 κ 4 (ik) 2 + θγ + θ + γ - κ 2 κ 4 ik + γ - κ 2 ζ - κ 4 κ 4 (ik) 2 (4.22) Km2 ω 2 1 + κ 4ik ik + 3κ 4 0 0 1 2ik - κ 4(ik) 2 (γ = 0, θ = κ 2 , ζ = κ 4 ) (4.26) Km2 ω 2 sym 1 + θ ik + κ 16(ik) 2 ik + θ + κ 2 0 - κ 4(ik) 2 1 2ik - κ 4(ik) 2 (γ = κ/4, ζ = κ/4) + κ θ - κ 8 2(ik)

Inria

This point is obtained with the artificial boundary very close to the obstacle (R 2 /R 1 = 1.005). It illustrates that the newly designed ABCs really improve the numerical error only when the wavelength is comparable to the size of the obstacle. Indeed, if the frequency is equal to 5 Hz, the wavelength is small compared to the size of the obstacle. In this case, Km2 omega2 performs as well as Km1 omega but not better (error = 25 %). But if the frequency is smaller, equal to 0.5 Hz, the wavelength is comparable to the size of the obstacle and Km2 omega2 performs better (32 %) than Km1 omega (64 %) or Km1 delta (37 %).

Finite Element simulations

Realistic simulations on arbitrary domains can only be performed by using efficient numerical methods such as finite elements methods. In the following we present numerical results obtained with Montjoie code (http: //montjoie.gforge.inria.fr) using Galerkin finite elements Q 8 . The meshes are generated using Gmsh software (http://geuz.org/gmsh). Quantifying the efficiency and accuracy of the ABCs on these arbitrary shapes requires a reference solution. This is done using a transparent condition based on an integral representation (see [START_REF] Duruflé | Intégration numérique et éléments finis d'ordre élevé appliqués aux équations de Maxwell en régime harmonique[END_REF]). To avoid integration accuracy loss, and obtain a good reference solution, the artificial boundary carrying the transparent condition will be put far away from the obstacle (several wavelengths away).

In the context of finite elements, a variational form of the ABCs is welcome. Unfortunately some of the ABCs do not own this property : Km1 delta 2 and Km1 omega 2. Another difficulty arises from the ABCs that include a term of the form a 2 (ik)ξ 2 ∂ r u because such a term requires the inversion of a second order operator on the artificial boundary. This is technically possible but has not been done yet in the finite elements code Montjoie. For this reason, Km1 delta 2 and Km2 delta 2 ABCs are not considered in that work. We have thus implemented the following ABCs : Km1 delta, Km1 omega, Km2 delta, Km2 omega and Km2 omega 2. We recall that they read as:

(∂ r u + ik u) + γ + κ 4 u + γ - κ 4 ∂ r u ik = 0 (Km1 delta)   1 + γ - κ 4   ∂ r u + ik + γ + κ 4 u + ξ 2 2 ik u = 0 (Km1 omega)   1 + θ + γ - κ 4 ik + θγ -θ - κ 2 κ 4 (ik) 2 - γ - κ 2 ζ + κ 4 κ 4 (ik) 3   ∂ r u (Km2 delta) +   ik + γ + θ + κ 4 + θγ + θ + γ - κ 2 κ 4 ik + γ - κ 2 ζ - κ 4 κ 4 (ik) 2   u = 0   1 + θ + γ - κ 4 ik   ∂ r u +   ik + γ + θ + κ 4 + θγ + γ + θ - κ 2 κ 4 ik   u + ξ 2 2 ik u = 0 (Km2 omega) 1 + θ ik + κ 16(ik) 2 ∂ r u +   ik + θ + κ 2 + κ 2 θ - κ 8 ik - κ 16 ζ - κ 4 (ik) 2   u (Km2 omega 2) - κ (s) 4(ik) 2 iξu + 1 2 ik - κ 4 (ik) 2 ξ 2 u = 0 6.2.

Validation of finite element ABC implementation on concentric circles

The first numerical test is done on concentric circles, in order to compare the numerical solutions to the analytical ones, which are available for this specific geometry. The exterior radius is chosen very close to the obstacle radius (R 1 = 2.0, R 2 = 2.01), and the frequency is equal to 0.5 Hz. The reader will notice that this is a very difficult test case, hence the high values of the error. The transparent condition is computed with an artificial boundary with radius R 2 = 10.0 (see the total and diffracted field's real part on figure 12). The computed error is the L2 relative error of the Dirichlet trace of the solution on the obstacle (R = 2.0) between the numerical solution and the analytic solution obtained as described earlier (see section 6.1). The first result is that the transparent condition leads to an error of 3.7 × 10 -9 % compared to the analytical solution, and can therefore serve as a reference solution.

Figure 12: Reference solution obtained with eighth order finite elements associated with a transparent condition.

In Fig. 13 we display the mesh of the domain {r, 2.00 ≤ r ≤ 2.01} with ABC on the exterior boundary, and the real part of the Dirichlet trace of the solution on the obstacle for several simulations. The L2 relative errors are displayed as a legend are listed in table 3 where they are compared to the analytic errors obtained in section 6.1 (no spatial discretization). The mesh can hardly be seen because the domain is very thin compared to its size. The results are very comparable (the 1 % discrepancy can be attributed to spatial discretisation, round off and projection errors) and show that the finite element code leads to the same conclusions as the analytical one. In Fig. 14 we represent the L2 relative error on all the domain of each ABC, as a function of parameter values. When several parameters exist for the ABC (as for Km2 delta and Km2 omega), the legend specifies at which value the non-varying parameter is set. We can notice on this curve that when the boundary is very close to the obstacle, the ABCs coming from the 1/ω asymptotic present singularities at some parameter values Inria other types of dependency. We can observe that Km1 delta and Km2 delta do not depend strongly on the parameter (the curves are rather flat) while Km1 omega, Km2 omega and Km2 omega2 show a stronger dependency on their parameter θ or γ. With the frequency 0.1 Hz, the two first ABCs provide higher errors than the three last, for which the best value of θ or γ seems to be close to 1. The value of the incidence angle does not seem to change the tendency of the results, and whatever value chosen for θ, the ABCs coming from the 1/ω asymptotic seem to perform better than the other ones (see the three first sub-figures). Recall that only the ABCs coming from the 1/ω asymptotic posses a Laplace-Beltrami term which could therefore explain their better behavior. Another set of experiments has been lead with and incident angle of 90 • and a frequency of 0.05 Hz (twice smaller). The results are displayed in the fourth sub-figure and are graphically displayed in Fig. 18. In this case, the θ or γ-dependency of Km1 omega, Km2 omega and Km2 omega2 is stronger, leading to crossing curves : for some values of the parameter θ or γ (close to 1), the 1/ω ABCs perform better than the other ones, but for other values of the parameter θ or γ (close to -2), the 1/ω ABCs perform less well than the other ones. This raises the question of how to choose the parameter, whether there is a universal answer or if the chosen value must be case-dependent. Notice for instance that for the ABC Km1 omega, the value γ = κ/4 coincides with the classical C-ABC Laplace Beltrami, but this is not the best possible value in these cases.

In the second example, the non convex shape is a regularized star, which can be seen in Fig. 19 It has been suggested in [START_REF] Barucq | Mathematical Models and Methods in Applied Sciences[END_REF] that the lack of performance of classical ABCs in such cases comes from the fact that the grazing waves are not controlled by the usual ABCs. The grazing waves are waves for which the couple (ω, ξ) satisfies λ 1 = 0. Following the approach done in page 140 of [START_REF] Barucq | Mathematical Models and Methods in Applied Sciences[END_REF], it is possible to write an enriched condition that controls grazing waves, and still falls in our family (5.1), by combining an ABC that controls grazing waves with any ABC of our family respecting a 1 = 0 and a 2 = 0:

(∂ r + c) a 0 (ik)∂ r + b 0 (ik) + b 1 (ik) iξ + b 2 (ik) ξ 2 u = 0 (6.1)
where

c ≡ c(ik, κ) = (6κ) 1/3 Γ(2/3) Γ(1/3) (ik) 2/3 (6.2)
Using the identity of remark 2.2, it gives a new ABC of the form (5.1) with tilded coefficients :

  (a0 (c -κ) + b 0 ) a0 + b 1 a1 iξ + b 2 a2 ξ 2    ∂ r u +    (ik) 2 a 0 + cb 0 b0 + cb 1 b1 iξ + (cb 2 + a 0 ) b2 ξ 2    u = 0 (6.3)
For instance, for the classical C-ABC, we get:

1 + c -κ/2 ik ∂ r u + ik + c + κ c 2ik u + ξ 2 ik u = 0 (6.4)
It is also possible to use it on Km2 delta since no term in ξ 2 appears in it (a 1 = 0 and a 2 = 0). Let us make computations for the ellipse where the grazing ABC is combined with the C-ABC (Km1 delta with γ = κ/4) or Km2 delta (the value of parameters will be chosen equal to (γ, θ, ζ) = (0.25, 0.5, 0.25) arbitrarily) and see in Fig. 25 which errors are obtained. The two ABCs obtained with this method are plotted in the two last subfigures and their errors indicated as subtitles. In the table below are indicated the errors obtained for the three frequencies (recall that the displayed error for Km2 delta is the best possible for all values of θ while here we are interested in the value θ = 0.5, which explains the discrepancy between the errors in Tab. 4 and in Fig. 25).

It is clear that using the combined ABC is always better than the original ABC (the column C-ABC+grazing is always better than the column C-ABC and the column Km2 delta + grazing than the column Km2 delta). Notice that the obtained ABC is better than the high order ABC Km2 omega2 for the frequencies 1.00 and 0.10 Hz but is not better when the frequency is equal to 0.05 Hz.

freq (Hz) C-ABC C-ABC + grazing Km2 delta Km2 delta + grazing Km2 omega2 (0.25, 0.5, 0.25) (0.25, 0. In the following of the report (and in previous figures as well), the error obtained with the control of grazing waves for Km1 delta or Km2 delta will be specified in parenthesis in the subtitle corresponding to the associated ABC. Notice that the value corresponding to Km2 delta is the best possible value and can therefore hardly be compared to the value in parenthesis. The purpose is to compare the value in parenthesis to the newly designed ABC Km2 omega2. In the prospect of using the ABCs in an "On-Surface Radiation Conditions" (OSRC) context, let us put the artificial boundary very close to the obstacle, and respecting its shape. This leads to non-convex artificial boundaries, and very hard numerical problems.

Inria

The first example is the peanut, as shown in Fig. 27 and 26. The frequency is 0.05, 0.10 and 0.50 Hz. When the frequency is high, the ABCs all perform quasi equivalently (error around 70 %), while the grazing-controlling ABCs offer a 52 % error. When the frequency diminishes, the results are very disappointing since the ABCs coming from the 1/ω asymptotic perform less well than the 1/δ ones and anyway the grazing-controlling ABCs give mush better results. The conclusion of this set of experiments would therefore be that it is essential to control grazing waves in the context of OSRC. The second example is the regularized star, which is even more non-convex, for frequencies 0.01, 0.05, 0.10 and 0.50 Hz as shown in Fig. 28 and29. In high frequency the conclusions are the same as before (the control of grazing waves seems to be the best option), but the results are very surprising in low frequency where it can happen that the best option, from far, is to use Km2 omega2. We think this is linked to the very strong dependency of the ABCs performance on the parameters values (see Fig. 29). This last set of experiment instructs us to continue our investigation before drawing further conclusions. Considering another type of dependency on s for the parameters might be one possibility. We compute at the same time the principal symbol of R 0 and its symbol of degree -1, which will be useful in the following.

Inria

σ * 0 (R 0 ) = σ * 0 (∂ r P 0 )P -1 0 + σ * 0 P 0 L 1 P -1 0 + σ * 0 P 0 L 0 P -1 0 (A.1a) σ * -1 (R 0 ) = σ * -1 (∂ r P 0 )P -1 0 + σ * -1 P 0 L 1 P -1 0 + σ * -1 P 0 L 0 P -1 0 (A.1b)
A.1 Symbol of P -1 0 .

Let Q = P -1 0 , and

Q = σ(Q) = Q 0 + Q -1 + Q -2 + . . . We have: σ(P Q) = σ(I) = I 2 ⇒ P 0 Q -i∂ s P 0 ∂ ξ Q - 1 2 ∂ 2 s P 0 ∂ 2 ξ Q = I 2
, which is of order 0 : we can identify the terms :

⇒              P 0 Q 0 = I 2 ⇒ Q 0 = P -1 0 P 0 Q -1 -i∂ s P 0 ∂ ξ Q 0 = 0 ⇒ Q -1 = iP -1 0 ∂ s P 0 ∂ ξ Q 0 P 0 Q -2 -i∂ s P 0 ∂ ξ Q -1 - 1 2 ∂ 2 s P 0 ∂ 2 ξ Q 0 = 0 ⇒ Q -2 = iP -1 0 ∂ s P 0 ∂ ξ Q -1 + 1 2 P -1 0 ∂ 2 s P 0 ∂ 2 ξ Q 0
We have :

∂ ξ Q 0 = ∂ ξ λ 1 √ 2    iω λ 2 1 0 0 1 iω    , where ∂ ξ λ 1 = ξ h 2 λ 1 (A.2) ∂ s P 0 = ∂ s λ 1 √ 2    -1 iω 0 0 -iω λ 2 1    , where ∂ s λ 1 = -ξ 2 rκ (s) λ 1 h 3 (A.3) hence Q -1 = -i(∂ ξ λ 1 )(∂ s λ 1 ) 2λ 2 1 P -1 0 (A.4) A.2 First part. σ (∂ r P 0 )P -1 0 = ∂ r (P 0 )Q -i∂ s ∂ r P 0 ∂ ξ Q + σ -2 (. . .) ⇒ σ * 0 (∂ r P 0 )P -1 0 = ∂ r P 0 Q 0 σ * -1 (∂ r P 0 )P -1 0 = ∂ r P 0 Q -1 -i∂ s ∂ r P 0 ∂ ξ Q 0 (A.5) Calculation gives: ∂ r P 0 = ∂ r λ 1 √ 2    -1 iω 0 0 -iω λ 2 1    , where ∂ r λ 1 = -κ(s)ξ 2 h 3 λ 1 (A.6) hence, σ * 0 (∂ r P 0 )P -1 0 = ∂ r λ 1 2     1 λ 1 -1 iω -iω λ 2 1 -1 λ 1     , where ∂ r λ 1 = -κ(s)ξ 2 h 3 λ 1 (A.7) Inria A.3 Second part. σ P 0 L 1 P -1 0 = σ(P 0 )σ(L 1 P -1 0 ) -i∂ s σ(P 0 )∂ ξ σ(L 1 P -1 0 ) - 1 2 ∂ 2 s P 0 ∂ 2 ξ σ(L 1 P -1 0 ) + σ -2 (. . .) = P 0 L 1 Q -i∂ s L 1 ∂ ξ Q - 1 2 ∂ 2 s L 1 ∂ 2 ξ Q + σ -2 (. . .) -i∂ s P 0 ∂ ξ L 1 Q -i∂ s L 1 ∂ ξ Q + σ -2 (. . .) - 1 2 ∂ 2 s P 0 ∂ 2 ξ L 1 Q + σ -2 (. . .) σ P 0 L 1 P -1 0 = P 0 L 1 Q 0 + Q -1 + Q -2 -i∂ s L 1 ∂ ξ Q 0 + Q -1 - 1 2 ∂ 2 s L 1 ∂ 2 ξ Q 0 -i∂ s P 0 ∂ ξ L 1 Q 0 + Q -1 -i∂ s L 1 ∂ ξ Q 0 - 1 2 ∂ 2 s P 0 ∂ 2 ξ L 1 Q 0 + σ -2 (. . .) σ * 1 P 0 L 1 P -1 0 = P 0 L 1 Q 0 (A.8) σ * 0 P 0 L 1 P -1 0 = P 0 L 1 Q -1 -iP 0 ∂ s L 1 ∂ ξ Q 0 -i∂ s P 0 ∂ ξ (L 1 Q 0 ) (A.9) σ * -1 P 0 L 1 P -1 0 = P 0 L 1 Q -2 -iP 0 ∂ s L 1 ∂ ξ Q -1 - 1 2 P 0 ∂ 2 s L 1 ∂ 2 ξ Q 0 -i∂ s P 0 ∂ ξ (L 1 Q -1 ) -∂ s P 0 ∂ ξ (∂ s L 1 ∂ ξ Q 0 ) - 1 2 ∂ 2 s P 0 ∂ 2 ξ L 1 Q 0 (A.10) P 0 L 1 Q -1 = - i(∂ s λ 1 )(∂ ξ λ 1 ) 2λ 2 1 P 0 L 1 P -1 0 = - i(∂ s λ 1 )(∂ ξ λ 1 ) 2λ 2 1 λ 1 0 0 -λ 1 (A.11) iP 0 ∂ s L 1 ∂ ξ Q 0 = -i(∂ ξ λ 1 )(∂ s λ 1 ) λ 1   1 0 iω λ 1 0   = iξ 3 rκ (s) h 5 λ 3 1   1 0 iω λ 1 0   (A.12) i∂ s P 0 ∂ ξ (L 1 Q 0 ) = i(∂ s λ 1 )(∂ ξ λ 1 ) 2     0 1 iω -iω λ 2 1 2 λ 1     ' (A.13) A.4 Third part. σ P 0 L 0 P -1 0 = P 0 σ(L 0 Q) -i∂ s P 0 ∂ ξ σ(L 0 Q) + σ -2 (. . .) = P 0 L 0 Q -i∂ s L 0 ∂ ξ Q + σ -2 (. . .) -i∂ s P 0 ∂ ξ L 0 Q + σ -2 (. . .) + σ -2 (. . .) = P 0 L 0 Q 0 + Q -1 -i∂ s L 0 ∂ ξ Q 0 -i∂ s P 0 ∂ ξ L 0 Q 0 + σ -2 (. . .) σ * 0 P 0 L 0 P -1 0 = P 0 L 0 Q 0 (A.14) σ * -1 P 0 L 0 P -1 0 = P 0 L 0 Q -1 -iP 0 ∂ s L 0 ∂ ξ Q 0 -i∂ s P 0 ∂ ξ L 0 Q 0 (A.15) P 0 L 0 Q 0 = 1 2     i ∂ s h ξ h 3 λ 1 -κ r - ∂ s h ξ ωh 3 - κ r λ 1 iω - ∂ s h ωξ h 3 λ 2 1 - iω κ r λ 1 -i ∂ s h ξ h 3 λ 1 -κ r     (A.16) RR n°8632
A.5 It gives...

σ * 0 (R 0 ) = 1 2     ∂ r λ 1 λ 1 + i (∂ s λ 1 )(∂ ξ λ 1 ) λ 1 + i ∂ s h ξ h 3 λ 1 -κ r - ∂ r λ 1 iω - (∂ s λ 1 )(∂ ξ λ 1 ) ω - ∂ s h ξ ωh 3 + κ r λ 1 iω - ∂ r λ 1 iω λ 2 1 - 3ω(∂ s λ 1 )(∂ ξ λ 1 ) λ 2 1 - ∂ s h ωξ h 3 λ 2 1 + iω κ r λ 1 - ∂ r λ 1 λ 1 - i(∂ s λ 1 )(∂ ξ λ 1 ) λ 1 -i ∂ s h ξ h 3 λ 1 + κ r    
(A.17)

Taking r = 0 in this matrix reduces the number of terms:

σ * 0 (R 0 )(r = 0) =     - κ 2λ 2 1 (2ξ 2 - ω 2 c 2 ) κω 2iλ 1 c 2 iκω 3 2λ 3 1 c 2 κω 2 2λ 2 1 c 2     (A.18) Inria B Calculation of σ * -1 (R -1 ) 1,2
We need to know the term σ * -1 (R -1 ) 1,2 in order to establish the CLA. We know by identifying terms in 3.1 that

R -1 = ∂ r (K -1 )(I + K -1 ) -1
First term

+ OP -1 R 0 Second term + OP -1 {K -1 ; D 1 + R 0 } (I + K -1 ) -1
Last term W 1 (B.1)

B.1 First term.

We have:

σ * -1 ∂ r (K -1 )(I + K -1 ) -1 = σ * -1 (∂ r K -1 ) = ∂ r K -1
The term (∂ r K -1 ) 1,2 writes, since we have (3.6):

(∂ r K -1 ) 1,2 = ∂ r σ * 0 (R 0 ) 1,2 2λ 1 = ∂ r (σ * 0 (R 0 ) 1,2 ) 2λ 1 - σ * 0 (R 0 ) 1,2 (∂ r λ 1 ) 2λ 2 1 (B.2)
The second term of (B.2) will considerably diminish when we will evaluate r = 0. Indeed, We need to calculate the first term before taking r = 0. We have, ) and (∂ s h) will vanish when r = 0, it is unnecessary to calculate their multiplicating terms. Only four terms are to be considered.

∂ r (σ * 0 (R 0 ) 1,2 ) = - ∂ 2 r λ 1 2iω -(∂
                       ∂ 2 r λ 1 = κ ξ 2 3 κ h 4 λ 1 + ∂ r λ 1 λ 2 1 h 3 (B.5a) ∂ 2 r,s λ 1 = -∂ r ξ 2 rκ λ 1 h 3 = - ξ 2 κ λ 1 h 3 + r (...) (B.5b) ∂ 2 r,s h = κ (B.5c) ∂ r κ λ 1 h iω = κ∂ r λ 1 h iω - κλ 1 (∂ r h) h 2 iω (B.5d)
When r = 0, we get

                       (∂ 2 r λ 1 )(r = 0) = - κ 2 ξ 2 λ 1 -3 + ξ 2 λ 2 1 (B.6a) (∂ 2 r,s λ 1 )(r = 0) = - ξ 2 κ λ 1 (B.6b) (∂ 2 r,s h)(r = 0) = κ (B.6c) ∂ r κ λ 1 h iω (r = 0) = - κ 2 iω ξ 2 λ 1 + λ 1 (B.6d)
Consequently, equation (B.4) becomes, when r = 0, The symbol of OP -1 R 0 was given in (A.1b).

∂ r (σ * 0 (R 0 ) 1,2 ) = κ 2 ξ 2 2iωλ 1 -3 + ξ 2 λ 2 1 + ξ 3 κ 2ωλ 2 1 - 1 2 ξ κ ω - κ 2 iω ξ 2 λ 1 + λ 1 = κ 2 ξ 2 2iωλ 1 -3 + 1 + ξ 2 λ 2 1 + κ 2 λ 1 2iω + κ ξ 3 2ωλ 2 1 - ξ 2ω = κ 2 2iω -2ξ 2 λ 1 + ξ 4 λ 3 1 + λ 1 + κ ξ 2ω ξ 2 -λ 2 1 λ 2 1 = κ 2 2iω -2ξ 2 λ 2 1 + ξ 4 + λ 4 1 λ 2 1 + κ ξ 2ω ω 2 λ 2 1 c 2 = κ 2 2iωλ 3 1 ξ 2 -λ 2
The only contributing term comes from -iP 0 ∂ s L 0 ∂ ξ Q 0 and in the end, We use formula 2.8 to calculate the symbol of exactly order -1 of the last term:

σ * -1 (W 1 ) = σ -1 K -1 D 1 -D 1 K -1 + K -1 R 0 -R 0 K -1 (I -K -1 + OP -2 ) = σ -1 (K -1 D 1 ) -σ -1 (D 1 K -1 ) -σ -1 (K -1 D 1 K -1 ) + σ -1 (D 1 K -1 K -1 ) + σ -1 (K -1 R 0 ) -σ -1 (R 0 K -1 ) = $ $ $ $ K -1 D 1 order 0 -i∂ s K -1 ∂ ξ D 1 - $ $ $ $ D 1 K -1 order 0 +i∂ s D 1 ∂ ξ K -1 -K -1 D 1 K -1 + D 1 K -1 K -1 + K -1 σ p (R 0 ) -σ p (R 0 )K -1 σ * -1 (W 1 ) = -i∂ s K -1 ∂ ξ D 1 + i∂ s D 1 ∂ ξ K -1 -K -1 D 1 K -1 + D 1 K -1 K -1 + K -1 σ p (R 0 ) -σ p (R 0 )K -1
where σ p (R 0 ) is given by (A.1a). Taking r = 0, we get :

∂ s K -1 (r = 0)∂ ξ P 0 (r = 0) t (u, v) 1 =   γ (s) λ 1 κ (s)ω 4iλ 2 1 × ×   ∂ ξ λ 1 √ 2    -u iω -iωv λ 2 1    (C.4) ∂ s K -1 (r = 0)∂ ξ P 0 (r = 0) t (u, v) 1 = - γ (s) ξ λ 2 1 iω u - κ (s) ξω 2 4λ 5 1 v (C.5) RR n°8632

Figure 1 :

 1 Figure 1: Diffracted filed obtained for the scattering of a plane wave on a circular obstacle of radius R 1 = 2 using an absorbing boundary condition on the circle of radius R 2 = 5.

Figure 2 :

 2 Figure 2: L2 relative error using ABC (3.14), dependency on the parameter γ. Legend : radius ratio R2/R1.

Figure 3 :

 3 Figure 3: L2 relative error using ABC (3.17), dependency on the parameter γ. Legend : radius ratio R2/R1.

Figure 4 :

 4 Figure 4: L2 relative error using ABC (3.23), dependency on the parameter γ. Legend : radius ratio R2/R1.

Figure 5 :

 5 Figure 5: L2 relative error using ABC (3.28), dependency on the parameter γ. Legend : radius ratio R2/R1.

Fig 5

 5 Fig 5 shows the relative L2 error for ABC (3.28) in the same conditions as Fig.2. The performance of the ABC seems to improve as the frequency increases. Here again the specific value γ = κ/2 seems to give much better results, especially when the boundary is far from the obstacle.

Figure 6 :

 6 Figure 6: L2 relative error using ABC (4.10), dependency on the parameters γ, θ and ζ.

Figure 7 :

 7 Figure 7: L2 relative error using ABC (4.12), dependency on the parameters γ, θ and ζ.

Inria

4. 2 . 2

 22 Asymptotic approximation: "High frequency"

Figure 8 :Figure 9 :

 89 Figure 8: L2 relative error using ABC (4.19), dependency on the parameters γ and θ.

  Fig 10 shows the relative L2 errors with respect to the parameters θ and ζ of ABC (4.26). The color scale represents the radius ratio between the artificial boundary and the obstacle (from 1.05 to 500). It can be seen that the ABC's performance does not really depend on the parameter ζ, which can therefore be taken equal

Figure 10 :

 10 Figure 10: L2 relative error using ABC (4.26), dependency on the parameters θ and ζ.

6 Numerical illustration 6 . 1

 61 Convergence curvesIn Fig.11we depict the relative L2 error of the Dirichlet trace of the solution on the circular obstacle of radius R 1 = 2, with different ABCs on the concentric circular artificial boundary of radius R 2 ∈ [2.01, 2560].

= 5 4 Figure 11 :

 5411 Figure 11: Convergence curves for all the ABCs. Remind that with these values of the parameters, Km1 delta corresponds to the C-ABC and Km1 omega to the C-ABC with Laplace-Beltrami term. Following the legend are indicated the errors for the first point of the curve where the absorbing boundary is very close to the obstacle (R 2 /R 1 = 1.005).

  Real part of the Dirichlet trace of the solution on the obstacle. L2 relative error compared to the analytical solution is indicated in percent.

Figure 13 :

 13 Figure 13: Finite element simulations for a circular obstacle of radius R 1 = 2 and an artificial circular boundary of radius R 2 = 2.01 with frequency 0.5 Hz. The errors are very high because the artificial condition is set very close to the obstacle.

  . The frequency is set to 0.5 Hz and the incident wave has a 0 • angle with the horizontal axis. The numerical results are displayed in Fig 20 and Fig. 21, in the same configuration as before. The newly designed ABCs provide good results, although they depend strongly on the parameter's value. The best value of γ for Km1 omega, and of θ for Km2 omega and Km2 omega2 seems close to 1.

  (a) Mesh for the transparent condition (b) Mesh for ABCs

Figure 15 :

 15 Figure 15: Meshes used in the simulations for Fig. 16 and Fig. 17.

  Figure 15: Meshes used in the simulations for Fig. 16 and Fig. 17.

Figure 16 :

 16 Figure 16: Peanut-shaped obstacle in a circular artificial boundary. Freq = 0.1 Hz. The analytical solution is computed in a larger domain than the actually displayed square area. Km1 delta and Km1 omega are computed with γ = κ/4 while for Km2 delta, Km2 omega and Km2 omega2 we display the simulation obtained with the best value of θ obtained in Fig. 17. Titles : relative L2 error of the Dirichlet trace of the solution on the obstacle.

  gamma=kappa/4, theta, zeta=kappa/4) Km2 omega (gamma=kappa/4, theta) Km2 omega2 (theta) (c) Incident angle = 90 • , Freq = 0.1 Hz gamma=kappa/4, theta, zeta=kappa/4) Km2 omega (gamma=kappa/4, theta) Km2 omega2 (theta) (d) Incident angle = 90 • , Freq = 0.05 Hz

Figure 17 :

 17 Figure 17: Peanut-shaped obstacle in a circular artificial boundary. Freq = 0.1 or 0.05 Hz

  Figure 18: Peanut-shaped obstacle in a circular artificial boundary. Freq = 0.05 Hz. The analytical solution is computed in a larger domain than the actually displayed square area. Km1 delta and Km1 omega are computed with γ = κ/4 while for Km2 delta, Km2 omega and Km2 omega2 we display the simulation obtained with the best value of θ obtained in Fig. 17. Titles : relative L2 error of the Dirichlet trace of the solution on the obstacle.

Figure 19 :

 19 Figure 19: Meshes used in the simulations for Fig. 20 and Fig. 21.

Figure 20 :Figure 21 :

 2021 Figure 20: Star-shaped obstacle in a circular artificial boundary. Freq = 0.5 Hz. The analytical solution is computed in a larger domain than the actually displayed square area. Km1 delta and Km1 omega are computed with γ = κ/4 while for Km2 delta, Km2 omega and Km2 omega2 we display the simulation obtained with the best value of θ obtained in Fig. 21. Titles : relative L2 error of the Dirichlet trace of the solution on the obstacle.

  Figure 25: Comparison of the ABCs with and without the grazing condition on an elliptic domain. The analytical solution is computed in a larger domain than the actually displayed rectangular area. Km1 delta and Km1 omega are computed with γ = κ/4 while for Km2 delta, Km2 omega and Km2 omega2 we display the simulation obtained with the best value of θ obtained in Fig. 23. The last two figures concern the C-ABC + grazing and the Km2 delta + grazing with (γ, θ, ζ) = (0.25, 0.5, 0.25). Titles : relative L2 error of the Dirichlet trace of the solution on the obstacle.

Figure 26 :

 26 Figure 26: Very close peanut-shaped boundary. Freq = 0.05, 0.1 and 0.5 Hz

  (a) Freq = 0.05 Hz (b) Freq = 0.1 Hz (c) Freq = 0.5 Hz

Figure 27 :

 27 Figure 27: Very close peanut-shaped boundary. Freq = 0.05, 0.1 and 0.5 Hz. The analytical solution is computed in a larger domain than the actually displayed rectangular area. Km1 delta and Km1 omega are computed with γ = κ/4 while for Km2 delta, Km2 omega and Km2 omega2 we display the simulation obtained with the best value of θ obtained in Fig. 23. Titles : relative L2 error of the Dirichlet trace of the solution on the obstacle.

  (a) Freq = 0.01 Hz (b) Freq = 0.05 Hz (c) Freq = 0.1 Hz (d) Freq = 0.5 Hz

Figure 28 :

 28 Figure 28: Very close star-shaped boundary. Freq = 0.01, 0.05, 0.1 and 0.5 Hz. The analytical solution is computed in a larger domain than the actually displayed rectangular area. Km1 delta and Km1 omega are computed with γ = κ/4 while for Km2 delta, Km2 omega and Km2 omega2 we display the simulation obtained with the best value of θ obtained in Fig. 29. Titles : relative L2 error of the Dirichlet trace of the solution on the obstacle.

Inria

Figure 29 :

 29 Figure 29: Very close star-shaped boundary. Freq = 0.01, 0.05, 0.1 and 0.5 Hz

κ 2 ω 4iλ 4 1 c 2 ω 2 c 2 + ξ 2 + κ ξω 4λ 3 1 c 2 (B. 7 )B. 2

 2272 r K -1 ) 1,2 (r = 0) = Second term.

B. 3 . 1 2 λ 2 1 ∂ 1 ξ 2 κ ξω 4λ 3 1 c 2 (B. 10 )B. 3 . 2 B. 3 . 6 c 2 + κ 2 ξ 2 ω 4iλ 4 1 c 2 (B. 15 )B. 4 C∂

 31112210323622154 First term.(-i∂ s K -1 ∂ ξ D 1 ) 1,2 = i(∂ ξ λ 1 )∂ s σ * 0 (R 0 ) 1,2 2λ 1 = i(∂ ξ λ 1 ) ∂ s (σ * 0 (R 0 ) 1,2 ) 2λ 1 + (∂ s λ 1 )(...) (-i∂ s K -1 ∂ ξ D 1 ) 1,2 (r = 0) = iξ 2h s (σ * 0 (R 0 ) 1,2 ) (r = 0) (B.9)Calculation gives:∂ s (σ * 0 (R 0 ) 1,2 ) = -∂ 2 r,s λ 1 2iω -(∂ 2 r,s λ 1 )(∂ ξ λ 1 ) 2ω + (∂ s λ 1 )(...) -(∂ 2 r,s h)ξ 2h 3 ω + (∂ s h)(...) + ∂ s -κλ 1 2iωh ∂ s (σ * 0 (R 0 ) 1,2 ) (r = 0) = κ ξ 2 2iωλ 1 -κ λ 1 2iω + (∂ s λ 1 )(...) + (∂ s h)(...) = κ 2ωiλ s K -1 ∂ ξ D 1 ) 1,2 (r = 0) = Second term. Since ∂ s D 1 = (∂ s λ 1 )(...), it vanishes when r = 0. B.3.3 Third term. (-K -1 D 1 K -1 ) 1,2 = -σ * 0 (R 0 ) 1,2 2λ 1 (γ(s) -ζ(s)) (-K -1 D 1 K -1 ) 1,2 (r = 0) = κω 4iλ 2 1 c 2 (ζ(s) -γ(s)) (B.11) B.3.4 Fourth term. (D 1 K -1 K -1 ) 1,2 = σ * 0 (R 0 ) 1,2 2λ 1 (γ(s) + ζ(s)) (D 1 K -1 K -1 ) 1,2 (r = 0) = κω 4iλ 2 1 c 2 (γ(s) + ζ(s)) Sixth term. (-σ * 0 (R 0 )K -1 ) 1,2 = -σ * 0 (R 0 ) 1,2 2λ 1 (σ * 0 (R 0 ) 1,1 + 2ζ(s)) (-K -1 σ * 0 (R 0 )) 1,2 (r = 0) = -It gives...If we consider the exact expressions (B.7), (B.8) and (B.15) we have Calculation of ∂ s K -1 ∂ ξ P 0

Table 2 :

 2 Summary of the coefficients for the ABCs of section 4 (ABC (4.12) is omitted because of its complexity)

Table 4 :

 4 Comparison of the L2 relative error of the Dirichlet trace of the solution on the obstacle with and without the control of grazing waves.

	5, 0.25)

RR n°8632

Inria

Avenue de la Vieille Tour, 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

(which are negative). However, there exists some values of the parameters for which the error is very low compared to the other ABCs.

Non convex obstacle and circular ABC

Non convex obstacles are considered in the following numerical examples. The reference solution is computed using the transparent condition on the boundary of a circle put far away from the obstacle, while the computations done with the ABCs are done with a circular artificial boundary very close to the obstacle.

In the first example, the non convex shape is a peanut, which can be seen on the displayed meshes in Fig. 15. The frequency is set to 0.1 Hz, while the angle of the incident wave is 0, 45 and 90 • with the horizontal axis. The numerical results for each angle of incidence are displayed in Fig. 16. The reference solution, computed with a transparent boundary far from the obstacle, is displayed in the left-top corner in a square domain. The reader can notice that the artificial boundary hence does not appear, the solution is drawn in all the displayed square. In the same line, different solutions are obtained by using ABCs on the artificial boundary. The first two are Km1 delta with γ = κ/4, Km1 omega with γ = κ/4, because they coincide with very classical ABCs (the C-ABC and the C-ABC with Laplace-Beltrami term). The three others are Km2 delta, Km2 omega and Km2 omega2, and the chosen values of γ and ζ are κ/4 but the value of θ can be different. We have chosen to represent graphically in Fig. 16 the simulation obtained with the best value depending on θ as shown in Fig. 17 for the three ABCs Km2 delta, Km2 omega and Km2 omega2, while Km1 delta and Km1 omega are displayed for γ = κ/4 since they therefore coincide with classical ABCs. This choice will be followed throughout this report hence we will not specify it again. The second line shows the error as a function of space, and the L2 relative error of the Dirichlet trace on the obstacle is specified in the subtitles. The meaning of the value under parenthesis will be explained later (see 6.2.3). It is then possible to see where the solution is wrongly computed : it appears to be in the regions of shadow, which differ with the incidence angle.

In Fig. 17, we present the value of the L2 relative error of the Dirichlet trace on the obstacle obtained for the same experiments, as a function of the parameters value. Each parameter is chosen under the form c × κ where c is a constant, and appears as abscissa of Fig. 17. The ABCs Km1 delta and Km1 omega are depending on γ while we present the dependency of Km2 delta, Km2 omega and Km2 omega2 on the parameter θ (γ and ζ are fixed to κ/4 for the experiments). Recall that a more general form could be chosen for the parameters, for instance a polynomial function of κ or any s-depending function. In this report we choose to focus on the linear case, but interesting results shall be obtained by considering

Elliptic domain and elliptic artificial boundary

An elliptic obstacle is considered in the following numerical examples. The reference solution is computed using the transparent condition on the boundary of a scaling of the same ellipse put far away from the obstacle, while the computations done with the ABCs are done with an elliptic artificial boundary very close to the obstacle. The associated meshes can be seen in Fig. 22. The frequency is set to 0.1 and 1 Hz, while the angle of the incident wave is 0 • with the horizontal axis. The numerical results are shown in Fig. 23 and Fig. 24. When the frequency is 1 Hz, the wavelength is very small compared to the obstacle. In this case, the three ABCs coming from 1/ω asymptotic perform clearly better (error around 51 %) than the ones coming from the 1/δ asymptotic (error around 66 %), showing that the Laplace Beltrami term is very important. Moreover, the ABCs performance are about the same for all values of parameters. But when the frequency decreases, the ABCs depend more and more on the parameters. When the frequency is 0.05 Hz, it is possible to find parameters values such that the ABCs coming from the 1/ω asymptotic lead to errors close to 27 %, but for other values of the parameters, the error can grow to more than 100 % which is clearly unacceptable.