On the impact of the shape of the artiﬁcial boundary in exterior Helmholtz
problems
Nick WULBUSCH(1) , Reinhild RODEN(2) , Alexey CHERNOV(1) , Matthias BLAU(2,3) , Andrea MOIOLA(4)(1) Institut für Mathematik, Universität Oldenburg, Germany, nick.wulbusch@uni-oldenburg.de(2) Institut für Hörtechnik und Audiologie, Jade Hochschule Oldenburg, Germany(3) Cluster of Excellence Hearing4All, Germany(4) Department of Mathematics, University of Pavia, Italy

Abstract
Acoustic problems can be investigated numerically in the frequency domain by solving the Helmholtz equation
with the ﬁnite element method (FEM). In exterior problems we approximate the radiation condition by truncating the considered domain with an artiﬁcial absorbing boundary. Popular choices include absorbing boundary
conditions (ABC) or perfectly matched layers (PML), where the implementation of the latter is usually more
involved. Both methods have the disadvantage (although for PML this is less signiﬁcant) that they still reﬂect
certain parts of the wave. This reﬂection error is typically small for normal (0 degrees) incidence but increases
with larger incidence angles. In this numerical study we consider a model problem in two dimensions. For various ABCs we study the impact of changing the shape of the artiﬁcial boundary on the reﬂection error and theaccuracy of the numerical simulations. The outcome of numerical simulations for different ABCs and methodsof shape optimization are discussed in detail.
Keywords: FEM, Absorbing boundary conditions, Helmholtz equation

1INTRODUCTION

Solutions to exterior Helmholtz problems, which arise often in acoustics and electromagnetics, can be approximated in a truncated domain via the ﬁnite element method. The truncation leads to additional errors, which
depend largely on the method that is used to approximate the Sommerfeld radiation condition [8, 6, 5]. Many
tools have been proposed to deal with this problem, like perfectly matched layers [3, 9, 6], the Dirichlet-toNeumann-operator [6] and local absorbing boundary conditions [1, 2]. In this study we further investigate someof the latter, more precisely the impedance boundary condition and the curvature absorbing boundary condition(C-ABC). These boundary conditions are quite well understood [1, 6, 7, 2] and often used on a circular/spherical truncation domain in 2D/3D. Further they lead to large errors, if the "angle of incidence" of the wave islarge. "Angle of incidence" is meant to be the angle between a suitably deﬁned Poynting vector and the normalvector at the boundary. In this study we will consider a model problem in 2D and change the artiﬁcial boundary
in order to ﬁnd a truncation that leads to improved accuracy.

2

MODEL PROBLEMConsider the following exterior Dirichlet problem (1)1lim |x| 2|x|→∞−Δp − k2 pp
∂p
∂ |x| − ikp==0pD=0in Ωon ΓD(1)where Ω ⊂ R2 is the complement of a bounded domain and ΓD its boundary. p denotes the (complex) pressureand k is the wavenumber. In this study we will look at approximations of this problem by truncating the7504

boundary. To deal with the Sommerfeld radiation condition we use two different boundary conditions, which
are brieﬂy introduced in the next subsection.
2.1 Boundary Conditions
In this section we introduce the boundary conditions to approximate the Sommerfeld radiation condition usedin this study. To do so we truncate the inﬁnite domain by an artiﬁcial boundary Γ+ . We follow the derivationin [1] and compare two different absorbing boundary conditions.
2.1.1 Impedance boundary condition
The simplest boundary condition is the impedance boundary condition. Here we enforce∂p− ikp = 0 on Γ+ ,∂n(2)where ∂∂ np denotes the normal derivative (pointing outwards) of p. This boundary condition perfectly absorbsa planar wave that hits the boundary perpendicularly. In exterior problems the scattered ﬁeld is not a planarwave, so that in practice this boundary condition leads to errors. However, if we place the artiﬁcial boundary
far away from the obstacle, we actually get a good approximation to the Sommerfeld radiation condition. Thelarge distance required between the scatterer and the truncation boundary may lead to signiﬁcant computational
cost and should therefore be avoided.
2.1.2 C-ABC
Following [1, eq.(3.14)] we arrive at another approximation for the Sommerfeld radiation condition, assuming
the angle of incidence is small1 :∂pκκ 1 ∂p− ikp + γ += 0 on Γ+ ,(3)p+ γ −∂n
44 ik ∂ n
where κ is the curvature of Γ+ and γ is a parameter that can be chosen freely. In the numerical experiments in[1] it was shown that the value of γ did not have much inﬂuence, so (as also described in [1], see also [2]) for
the sake of simpliﬁcation we assume γ = κ4 and arrive at the so called curvature absorbing boundary condition(C-ABC):κ∂p− ikp + p = 0 on Γ+ .
(4)
∂n2
2.2 Shape deformation
Since the standard impedance boundary condition is exact for perpendicular incoming planar waves, we expectthat it approximates a non-planar wave better when this hits the boundary perpendicularly. We propose to
estimate the main direction of propagation of p by the Poynting vector ℜ( p̄∇pik ) or alternatively by the gradientof the complex argument arg(p), since a simple calculation revealsℜ(p̄∇p) = c(p, k)∇ arg(p),ik(5)for some scalar c, so the vectors in both expressions point into the same direction (except for the discontinuity
jumps of arg(p) from π to −π). Therefore the estimated main direction of propagation is perpendicular to acontour line of arg(p). For the sake of a fair comparison we want to ﬁx the area of the domain, such that thecomputational cost does not change signiﬁcantly (we will get roughly the same number of elements). Hence inthe numerical experiments we will consider a contour line that encloses nearly the same area as the nominal
circle and consider the interior as the new computational domain.
1 Notice that the sign in front of ikp is modiﬁed according to the sign convention in the time-harmonic ansatz from [1]7505

3NUMERICAL SIMULATIONSIn this section we present some numerical results that could be observed with shape deformation. In the following experiments we have wavenumber k = ωc with angular frequency ω = 2π f , frequency f and speed ofsound c = 343 ms . The spatial dimensions are given in meter (m).
3.1 Model problem I: Scattering
The ﬁrst model problem we are looking at is a standard scattering problem of the following form
−Δp − k2 pp
∂p
∂ |x| − ikp

1lim |x| 2|x|→∞=0= −pinc=in Ωon ΓD(6)0

where pinc , deﬁned in polar coordinates (r, θ ),pinc = p0 eikr sin(θinc −θ )(7)is the incoming (planar) wave hitting the obstacle from the right, i.e. (following the convention in [4]) θinc = 32 π,
and amplitude p0 = 1. Here ΓD describes the boundary of a circle, the scatterer, with radius a and Ω is theexterior of this scatterer. The analytic solution (see [4]), i.e. the scattered acoustical ﬁeld, of this problem isgiven (in polar coordinates) by∞p(r, θ ) = −p0 ∑Jn (ka)(1)
n=−∞ Hn (ka)(1)

Hn (kr)ein(θinc −θ ) , r ≥ a,(8)(1)

where p0 = 1, θinc = 32 π, Jn denotes the Bessel function of the ﬁrst kind of order n and Hnfunction of ﬁrst kind of order n.

denotes the Hankel

In our experiments we truncate the inﬁnite domain by different artiﬁcial boundaries and solve the correspondingproblem (9) with standard quadratic ﬁnite elements on that ﬁnite domain, where we use different approximations
for the Sommerfeld radiation condition:0−Δp − k2 p =p = −pinc
B(p) =0in Ωon ΓD
on Γ+ ,

(9)

where B(p) denotes one of the conditions (2) or (4). We choose a = 0.5 m for the radius of the scatterer.
3.1.1 Truncated domains
For the numerical simulations we used two different domains with the following boundaries:(1) Circle with radius 1.5 m, see ﬁgure 1(2) Contour line of arg(p), see ﬁgure 2
The boundary of the contour domain is generated with the help of Matlabs contour-function on the analytic
solution with f = 1000 Hz. The shape of the contour lines does not vary much with the frequency. For low
frequencies however, for which the wavelength is much larger than the diameter of the domain, there is nocontour line of arg(p) that could be used as boundary for the truncated domain, such that the area remains
approximately the same. Thus, for lower frequencies (< 150 Hz) the angles of incidence are not necessarilyclose to 0.

7506

Figure 1. |p| and arg(p) in circular domain, f = 1000 Hz

Figure 2. |p| and arg(p) in contour domain, f = 1000 Hz

Impedance BC on nominal domain
C-ABC on nominal domain
Impedance BC on contour domain
C-ABC on contour domain

10-1

10-2

102

103

104

Figure 3. Relative L2 -error as a function of frequency f for a ﬁxed mesh (h ≈ 0.0118, #(degrees of freedom)≈ 4.5 × 105 , 7 elements per wavelength for largest frequency) for the circular scattering problem.

7507

Figure 4. arg(p) in circular domain (left) and contour domain (right), f = 1000 Hz
3.1.2 Results
We compare the relative error in L2 -norm for the whole domainerr(h, f ) =pFEM − pexact L2 (Ω)
pexact L2 (Ω),(10)

where h denotes the mesh width, i.e. the maximal edge length of a triangular element of the mesh. Figure 3
shows the error for the impedance boundary condition and for the C-ABC in the nominal (circular) and contour
domain. The C-ABC performs better than the impedance boundary condition. For the impedance boundary
condition the error is smaller on the contour domain than on the nominal one. For the C-ABC this effectonly takes place for frequencies higher than 200 Hz. Notice that the mesh is ﬁxed, i.e. solutions for higher
frequencies might still improve from mesh reﬁnement, while for the solutions for lower frequencies the model
error due to the approximation of the Sommerfeld radiation dominates. From a practical stand point this is areasonable comparison, since for many problems we cannot reﬁne the mesh till the discretization error is below
a certain threshold due to the increasing computational cost and choose the mesh size according to some "ruleof thumb", such as 6 or 10 elements per wave length [6].
3.2 Different shape of scatterer
To validate the method, we look at a second scattering problem. This problem is identical to the previous
one, except that now the scatterer is polygonal, to be precise a square with vertices (−0.5, −0.5), (0.5, −0.5),(0.5, 0.5) and (−0.5, 0.5). No analytic solution is known for this problem. As reference solution we take aPML solution on a larger domain. The argument of the solution arg(p) can be seen in ﬁgure 4 for the nominal
domain and the contour domain. The error is shown in ﬁgure 5. We can observe a similar behavior as in theprevious problem.
3.3 Model problem II: Point sources
Finally we will take a look at another model problem, that generates a solution with a slightly more complicated(0)
contour line of arg(p). We consider equation (1) with pD as a sum of 50 Hankel functions H1 with centerat (0, 0.4 − 0.8
49 (n − 1)), n = 1, . . . , 50, i.e. an linear array of point sources on the y-axis from −0.4 to 0.4. The
modulus of the solution |p|, as well as arg(p) are shown in ﬁgure 6 for f = 2000 Hz. The resulting truncatedcontour domain is visualized in ﬁgure 7 (left). For this problem the shape of the contour lines of arg(p) largely
depends on the frequency. We consider f = 2000 Hz for generating the contour domain, but still compare theerror for other frequencies. Here the results do not improve for the contour domain and are only similar forhigh frequencies. Even for f = 2000 Hz the contour domains performs worse than the nominal one. This might
be due to the non-convexity, numerical inaccuracies in the computation of the contour line and in the computed7508

Impedance BC on nominal domain
C-ABC on nominal domain
Impedance BC on contour domain
C-ABC on contour domain

10-1

10-2
102

103

104

Figure 5. Relative L2 -error as a function of frequency f for a ﬁxed mesh (h ≈ 0.0118, #(degrees of freedom)≈ 4.5 × 105 , 7 elements per wavelength for largest frequency) for the polygonal scatterer problem.
curvature κ in case of the C-ABC. A simpler, convex domain is considered: An ellipse that is as close as
possible to the contour line in the L2 -sense. We compute the parameters of this elliptic domain by solving theoptimization problem (12) considered in the next subsection.
3.3.1 Ellipse optimization
The contour line will be approximated by an ellipse, parameterized via
f (x, y) = (x0 + a cos(θ ) cos(ϕ) − b sin(θ ) sin(ϕ), y0 + a sin(θ ) cos(ϕ) − b cos(θ ) sin(ϕ))(11)

where (x0 , y0 ) denotes the center of the ellipse, a and b the semiaxes’ lengths, ϕ the rotation of the ellipse
around its center and θ := θ (x, y) is the azimuth variable.
We solve the following optimization problem numericallymin γ(θ ) − f (xop (x0 , a, b, ϕ; θ ), yop (y0 , a, b, φ ; θ )) 2
x0 ,y0 ,a,b,ϕL ([0,2π])xop = x0 + a cos(θ ) cos(ϕ) − b sin(θ ) sin(ϕ)
yop = y0 + a sin(θ ) cos(ϕ) − b cos(θ ) sin(ϕ)(12)
Here γ denotes the parametrization of the contour line with respect to θ . x0 , y0 and ϕ are included, eventhough it is obvious that we will get x0 = y0 = ϕ = 0 up to machine precision due to symmetry. We solve thisproblem numerically by the gradient descent method, where we initialize a = 1.5 = b, x0 = 0 = y0 and ϕ = 0,i.e. a circle with radius 1.5 m centered at the origin (so it is equal to the nominal domain). The algorithm
returns a = 1.395 and b = 1.579. The resulting domain can be seen in ﬁgure 7 (right). For this elliptic domain
the error decreases for some frequencies, see ﬁgure 8. The impedance boundary condition behaves similar on
all domains for low frequencies, but performs better for the elliptic domain for high frequencies. For low
frequencies C-ABC performs best for the nominal domain, for high frequencies the elliptic domain is slightly
superior in accuracy.
such that

4CONCLUSIONSIn this numerical study we compared two simple absorbing boundary conditions on different shapes of domains.
We showed that deforming the shape of the artiﬁcial circular boundary, such that the angle of incidence of thewave is normal to the boundary, can improve numerical results. This was especially recognizable for higher
frequencies and for the impedance boundary condition. This also indicates that the shape deformation is more7509

Figure 6. |p| and arg(p) in circular domain, f = 2000 Hz

Figure 7. arg(p) in contour domain (left) and in optimized elliptic domain (right), f = 2000 Hz

Impedance BC on nominal domain
C-ABC on nominal domain
Impedance BC on contour domain
C-ABC on contour domain
Impedance BC on elliptic domain
C-ABC on elliptic domain

10-1

10-2102103104

Figure 8. Relative L2 (Ω)-error as a function of frequency f for a ﬁxed mesh (h ≈ 0.0120, #(degrees of freedom)≈ 4.5 × 105 , 7 elements per wavelength for largest frequency) for the point source problem

7510

useful for simpler boundary conditions, since more involved conditions already deal better with some range ofangles of incidence. Some problems, however, do not yield improved accuracy by using a contour domain.

5OUTLOOKIn this study we used contour lines of the exact solution as boundary for new domains. Usually the exact
solution is not known and hence a contour line can only be constructed by a numerical solution. Since solving
the problem on a larger domain ﬁrst and then truncating it is not necessarily reasonable, one could implement analgorithm that starts with a circular domain and changes it iteratively to a domain that is similar to the contour
domain discussed above. One could further consider higher order absorbing boundary conditions and see atwhich order the improvement of deforming the domain vanishes. The transition to 3D problems is possible, butleads to high computational cost.

ACKNOWLEDGEMENTS
Partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 352015383 – SFB 1330 C1.

REFERENCES[1] Barucq, H.; Bergot, M., Chabassier, J., Estecahandy, E. Derivation of high order absorbing boundary conditions for the Helmholtz equation in 2D. [Research Report] RR-8632, INRIA, Bordeaux; INRIA. 2014.
<hal-01085180>[2] Bayliss, A.; Gunzburger, M.; Turkel, E. Boundary conditions for the numerical solution of elliptic equations
in exterior regions. SIAM J. Appl. Math., vol 42, pp 430–451, 1982.
[3] Bérenger, J. A Perfectly Matched Layer for the Absorption of Electromagnetic Waves. J. Comput. Phys.
114 (2), pp. 185–200, 1994
[4] Bourlier, C.; Pinel, N.; Kubické, G. Method of Moments for 2D Scattering problems: Basic concepts andapplications. John Wiley & Sons 2013. 10.1002/9781118648674.
[5] Colton, D.; Kress, R., Integral Equation Methods in Scattering Theory. New York et al., John Wiley & Sons
1983
[6] Ihlenburg, F. Finite Element Analysis of Acousting Scattering. Springer, New York, 1998.
[7] Melenk, J. M.; Sauter, S. Wavenumber explicit convergence analysis for Galerkin discretizations of theHelmholtz equation. SIAM Journal on Numerical Analysis, vol. 49, no. 3/4, 2011, pp. 1210–1243.
[8] Sommerfeld, A. Partial Differential Equations in Physics, Academic Press, New York, 1949.
[9] Turkel, E.; Yefet, A. Absorbing PML boundary layers for wave-like equations. Appl. Numer. Math. 27 (4),pp. 533–557, 1998

7511