Towards high frequency boundary element methods for multiple scattering

Oliver Phillips 1

University of Reading Department of Mathematics and Statistics, University of Reading, Reading, UK

Simon Chandler-Wilde 2

University of Reading Department of Mathematics and Statistics, University of Reading, Reading, UK

Stephen Langdon 3

Brunel University London Department of Mathematics, Brunel University London, London, UK

ABSTRACT Standard Boundary Element Methods (BEM) for time-harmonic acoustics, using piecewise polynomial finite-element type approximation spaces, have a computational cost that grows rapidly with frequency, to ensure at least a fixed number of degrees of freedom per wavelength. Hybrid Numerical-Asymptotic (HNA) BEMs, based on enriched approximation spaces consisting of the products of piecewise polynomials with carefully chosen oscillatory functions, have a computational cost that is almost frequency-independent for some problem classes, but the technology is largely undeveloped for problems where multiple scattering is important. In this paper we present a computational method, supported by mathematical analysis, which suggests that multiple scattering configurations may be within reach. Specifically, we propose an algorithm to solve, by a HNA BEM, scattering by a pair of screens in an arbitrary configuration, which we anticipate may serve as a building block towards general multiple scattering problems with computational cost independent of frequency. The specific configuration considered, as we discuss, is relevant to the simulation of multiple outdoor noise barriers.

1. INTRODUCTION

It has been shown that noise pollution in highly populated areas can have a negative e ff ect on people’s health [1]. One potential method to reduce noise pollution from motorways or high-speed train lines is the use of noise barriers designed to di ff ract, reflect and / or absorb unwanted noise, leading to

1 o.j.phillips@pgr.reading.ac.uk

2 s.n.chandler-wilde@reading.ac.uk

3 stephen.langdon@brunel.ac.uk

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

Figure 1: A cross-section through an infinite line source between parallel noise barriers on a hard ground surface (left). By simple method of images arguments this is equivalent to the di ff raction of the original source plus an image source in the ground surface by two noise barriers in free space (right).

significant sound level reductions in the shielded areas [2,3]. In this paper we describe a method that e ffi ciently computes the scattered field in 2D due to the interaction of an incident coherent line source of sound interacting with two screens in free space, which (see Figure 1 and [3]) can be used to model parallel noise barriers sitting on a ground surface. Our time-harmonic acoustics approach is to use a so-called hybrid-numerical-asymptotic (HNA) boundary element method (see, e.g., [4,5]), i.e. a BEM with enriched basis functions derived from consideration of the high-frequency asymptotic behaviour of the solution. These enriched basis functions are designed to capture e ffi ciently the oscillatory acoustic field with a number of degrees of freedom which is essentially independent of the frequency. This contrasts with the need, in standard BEM based on piecewise-polynomial approximation, to increase the number of degrees of freedom in proportion to the frequency, e.g. using “10 degrees of freedom per wavelength”, to maintain accuracy as the frequency increases. For simplicity we focus in our description on the idealised case when the screen surfaces are sound soft (i.e., the acoustic pressure vanishes), but the method should extend straightforwardly to sound hard or impedance boundary conditions (cf., e.g., [6]).

2. PROBLEM STATEMENT

We consider the scattering of a time-harmonic incident wave u i (e − i ω t time dependence), that could be either one or m ore poi nt sources or a plane wave, by the union of two disjoint 1 D screens, Γ = Γ 1 S Γ 2 , in D = R 2 \ Γ , where Γ denotes the closure of Γ . The two screens can be in any orientation so long as they are not touching (examples are Figure 2 or the right hand side of Figure 1). The scattering problem we are looking to solve is to find the total field u ∈ C 2 ( D ) ∩ W 1 loc ( D ) (our standard function space notations are those of, e.g., [7]) such that

∆ u + k 2 u 0 in D , (1) u 0 on Γ , (2)

where k > 0 is the wavenumber and the scattered field u s = u − u i satisfies the standard Sommerfeld radiation condition. It follows from Green’s 2nd identity (see, e.g., [4,5]) that

Z

u ( x ) = u i ( x ) − i

H (1) 0 ( k | x − y | ) φ ( y ) ds ( y ) , x ∈ D , (3)

4

where φ ∈ e H − 1 / 2 ( Γ ) is the jump in the normal derivative of u across the boundary and H (1) 0 is the Hankel function of the first kind of order zero, see [5] for details. Further, φ satisfies the boundary

integral equation

Z

u i ( x ) = i

H (1) 0 ( k | x − y | ) φ ( y ) ds ( y ) , x ∈ Γ . (4)

4

Our solution method will be to solve Equation 4 for the unknown function φ on Γ by a HNA BEM, following which the solution in the domain D is given as an integral over Γ by Equation 3.

3. MULTIPLE SCATTERING ITERATIVE METHOD

For simplicity we consider the case where u i is the incident plane wave u i ( x ) : = e i kx · d , where d is a unit vector in the direction of the plane wave. For ease of notation, define φ j : = φ | Γ j ∈ e H − 1 / 2 ( Γ j ), and let S ℓ, j : e H − 1 / 2 ( Γ j ) → H 1 / 2 ( Γ ℓ ) be defined by

S ℓ, j ψ ( x ) : = Z

Γ j H (1) 0 ( k | x − y | ) ψ ( y ) ds ( y ) , x ∈ Γ ℓ , (5)

for ψ ∈ e H − 1 / 2 ( Γ j ) and ℓ, j = 1 , 2. Equation 4 can then be written as

S 11 φ 1 + S 12 φ 2 u i | Γ 1 , (6) S 21 φ 1 + S 22 φ 2 u i | Γ 2 . (7)

Figure 2: Re( u ) in D , with the screens Γ 1 on the left and Γ 2 on the right. The incident wave is a plane wave with direction d indicated by the arrow, and k = 5.

We will solve the above 2 × 2 system of equations by an iterative method. The first step (iteration zero) is to ignore the e ff ect of Γ 2 so Equation 6 becomes

S 11 φ (0) 1 = u i | Γ 1 , (8)

where the 0 in the superscript refers to the number of the iteration considered. We next solve Equation 7, replacing φ 1 by φ (0) 1 , so that we solve

S 22 φ (1) 2 = u i | Γ 2 − S 21 φ (0) 1 . (9)

Physically, φ (1) 2 is an approximation to φ 2 that includes the e ff ect of the first reflection from Γ 1 incident onto Γ 2 . We next solve Equation 6 with φ 2 replaced by φ (1) 2 , etc., so that our iterative scheme is to solve for r = 0 , 1 , 2 , ...,

S 11 φ (2 r ) 1 u i | Γ 1 − S 12 φ (2 r − 1) 2 , (10)

S 22 φ (2 r + 1) 2 u i | Γ 2 − S 21 φ (2 r ) 1 . (11)

As we will see below φ (2 r ) 1 → φ 1 and φ (2 r + 1) 2 → φ 2 , i.e. the true solution is approached, as r →∞ , indeed accurate results are obtained for rather small r .

4. HIGH FREQUENCY APPROXIMATION SPACE

To solve Equations 10 and 11 for a given r with a small number of degrees of freedom, essentially independent of the wavenumber k , we propose a HNA BEM approximation space adapting the construction for the case of a single screen in [5]. The solution φ (2 r ) 1 to Equation 10 can be decomposed as

φ (2 r ) 1 ( s ) = Ψ (2 r ) 1 ( s ) + v + , 2 r 1 ( s )e i ks + v − , 2 r 1 ( s )e − i ks , (12)

for s ∈ [0 , L 1 ], where L 1 is the length of Γ 1 and s denotes the distance from one of the end points. (Essentially the same decomposition, with obvious changes, holds also for the solution φ (2 r + 1) 2 ( s ) to Equation 11.) The term Ψ (2 r ) 1 represents the leading order high frequency behaviour, the physical optics or Kirchho ff approximation, defined as twice the normal derivative of the total field incident on Γ 1 , including in this total both the original incident field u i and the wave incident on Γ 1 from the screen Γ 2 . Precisely, at this iteration,

Ψ (2 r ) 1 = 2 ∂

 u i −S 2 [ φ (2 r − 1) 2   Γ 1 ,

∂ n

where, for ψ ∈ e H − 1 / 2 ( Γ j ), S j ψ ∈ C 2 ( D ) ∩ W 1 loc ( D ) is given, for j = 1 , 2, by

Z

S j ψ ( x ) : = i

Γ j H (1) 0 ( k | x − y | ) ψ ( y ) ds ( y ) , x ∈ D ,

4

and b ψ ( x ) : = ψ ( x ) if a point source at x is incident on the same side of Γ 1 as u i , otherwise b ψ ( x ) : = − ψ ( x ). The remaining part of the right hand side of Equation 12, i.e.

ϕ (2 r ) 1 ( s ) : = φ (2 r ) 1 ( s ) − Ψ (2 r ) 1 ( s ) = v + , 2 r 1 ( s )e i ks + v − , 2 r 1 ( s )e − i ks , (13)

captures di ff raction from the corners. The first term in the sum on the right hand side captures waves di ff racting from the corner at s = 0 and travelling along the screen in the positive s direction. The other term captures waves di ff racting from the other corner at s = L 1 and travelling along the screen in the negative s direction. As in [5], it can be shown that the functions v ± , 2 r 1 in Equation 12 are not oscillatory and hence can be approximated using standard piecewise polynomials with a number of degrees of freedom essentially independent of the wavenumber k . Therefore we can approximate ϕ (2 r ) 1 by a sum of products of piecewise polynomials (i.e. standard BEM basis functions) and e ± i ks (this is our HNA BEM approximation space). Substituting the decomposition of Equation 13 into Equation 10 means we are solving, for r = 0 , 1 , 2 , . . . ,

S 11 ϕ (2 r ) 1 = u i | Γ 1 − S 12 φ (2 r − 1) 2 − S 11 Ψ (2 r ) 1 . (14)

These equations can each be solved by either the Galerkin method, as in [5], or the least squares collocation method of [8], using (whichever solution method we choose) the HNA BEM approximation space we have just described. For each r the Galerkin or collocation matrix that is the discretisation of the operator S 11 on the left-hand side of Equation 14 remains the same and only the right-hand side depends on r .

5. RESULTS

Figure 3: The real parts of the iterates φ ( r ) 1 ( s ) on Γ 1 (top, r even) and φ ( r ) 2 ( s ) on Γ 2 (bottom, r odd) for the configuration of Figure 2 with k = 5.

In this section we show results which test the first component of our proposed algorithm, the multiple scattering iterative method, showing that φ (2 r ) 1 → φ 1 and φ (2 r + 1) 2 → φ 2 , i.e. the true solution is approached, as r →∞ , indeed that convergence is exponentially fast. For this purpose we compute φ (2 r ) 1 and φ (2 r + 1) 2 by a standard piecewise constant collocation BEM, using a uniform mesh on each of Γ 1 and Γ 2 , with a number of degrees of freedom per wavelength (dof) that we indicate below. For the configuration of Figure 2 the real parts of the solutions φ ( r ) j , for r between 0 and 9 and j = 1 , 2, are plotted in Figure 3, together with plots of φ ( r ) j for r large ( r = 58 in the upper figure, r = 59 in the lower figure) as approximations to the limits φ 1 and φ 2 . (The piecewise constant BEM is used with 320 degrees of freedom per wavelength, which is expected to give high accuracy. Plots for the imaginary parts are similar.) As we can see from the upper part of Figure 3, φ (2 r ) 1 converges to φ (58) 1 ≈ φ 1 within a very few iterations, indeed φ (2 r ) 1 is indistinguishable from φ (58) 1 for r ≥ 1. (The initial estimate φ (0) 1 is not a good approximation as it contains no contribution from the reflections and di ff ractions due to Γ 2 .) In the lower part of Figure 3, φ (2 r + 1) 2 is indistinguishable from φ (59) 2 ≈ φ 2 already for r = 0. The convergence is more clearly shown in Figure 4 where we plot the normalised ℓ 1 error defined

—r=0 s/Ly s/Ly

by

err ( r ) j , N : = ∥ φ ( R ) j , N − φ ( r ) j , N ∥ ℓ 1

∥ φ ( R ) j , N ∥ ℓ 1 , (15)

where φ ( r ) j , N denotes the approximation to φ ( r ) j obtained by the piecewise constant collocation BEM when N degrees of freedom per wavelength are used, and R = 58 when j = 1, R = 59 when j = 2, so that φ ( R ) j , N ≈ lim r →∞ φ ( r ) j , N . In Figure 4 we see that as r is increased the error decreases exponentially until r ≈ 10, after which it remains constant at a value between 10 − 12 and 10 − 15 depending on N .

Figure 4: The relative error err ( r ) j , N in the approximation φ ( r ) j , N , plotted against iteration number r , for screen 1 (top, j = 1) and screen 2 (bottom, j = 2), in each case for a range of values of N , the number of degrees of freedom per wavelength (dof).

6. CONCLUSION

In this paper we have described an e ffi cient method for computing scattering of an incident wave by a pair of screens in any orientation, which captures all reflections and interactions between the two screens. The method has two components, a multiple scattering iteration, described in §3, which reduces the solution to a sequence of solves; at each iteration we solve a problem of scattering by one of the screens of the original incident field plus a field incident from the other screen. The second component to our method is an e ffi cient HNA BEM representation for the solution at each iteration, an extension of the methodology used for scattering by single screens in [5] and [8], described in §4. We have presented, in §5, results that suggest that the iterative scheme converges exponentially fast, achieving high accuracy after a few iterations.

Normalised ¢; error Normalised ¢; error 19-10 | ~~ dof =5 ~~ dof = 10 dof = 20 ~~ dof = 40 ~~ dof = 80 ~— dof = 160 49-10 | ~~ dof = 5 ~~ dof = 10

ACKNOWLEDGEMENTS

OP is supported by a PhD studentship from the Mathematics of Planet Earth Centre for Doctoral Training at the University of Reading, provided through financial support from the Engineering and Physical Sciences Research Council (EPSRC), grant number EP / L016613 / 1. SNCW is supported by EPSRC grant EP / V007866 / 1. We thank Andrew Gibbs (University College London) for many helpful discussions.

REFERENCES

[1] L. Fritschi, A. L. Brown, R. Kim, D. Schwela, and S. Kephalopoulos. Burden of disease from environmental noise: Quantification of healthy life years lost in Europe . World Health Organization Regional O ffi ce for Europe, 2011. https: // apps.who.int / iris / handle / 10665 / 326424. [2] Benz Kotzen and Colin English. Environmental Noise Barriers: A Guide To Their Acoustic and Visual Design, 2nd Ed. CRC Press, 2019. [3] Qiutong Li, Denis Duhamel, Yanyun Luo, and Honore Yin. Analysing the acoustic performance of a nearly-enclosed noise barrier using scale model experiments and a 2.5-D BEM approach. Applied Acoustics , 158:107079, 2020. [4] S. N. Chandler-Wilde, I. G. Graham, S. Langdon, and E. A. Spence. Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering. Acta Numerica , 21:89–305, 2012. [5] D. P. Hewett, S. Langdon, and S. N. Chandler-Wilde. A frequency-independent boundary element method for scattering by two-dimensional screens and apertures. IMA Journal of Numerical Analysis , 35(4):1698–1728, 2015. [6] S. N. Chandler-Wilde, S. Langdon, and M. Mokgolele. A high frequency boundary element method for scattering by convex polygons with impedance boundary conditions. Communications in Computational Physics , 11:573–593, 2012. [7] W. C. H. McLean. Strongly elliptic systems and boundary integral equations . Cambridge University Press, 2000. [8] A. Gibbs, D. P. Hewett, D. Huybrechs, and E. Parolin. Fast hybrid numerical-asymptotic boundary element methods for high frequency screen and aperture problems based on least- squares collocation. SN Partial Di ff erential Equations and Applications , 1:21, 2020.