# Local on-surface radiation condition for multiple scattering of waves from convex obstacles

We propose a novel on-surface radiation condition to approximate the outgoing solution to the Helmholtz equation in the exterior of several impenetrable convex obstacles. Based on a local approximation of the Dirichlet-to-Neumann operator and a local formula for wave propagation, this new method simultaneously accounts for the outgoing behavior of the solution as well as the reflections arising from the multiple obstacles. The method involves tangential derivatives only, avoiding the use of integration over the surfaces of the obstacles. As a consequence, the method leads to sparse matrices and O(N) complexity. Numerical results are presented to illustrate the performance of the proposed method. Possible improvements and extensions are also discussed.

## Authors

• 3 publications
• ### Gravity-capillary flows over obstacles for the fifth-order forced Korteweg-de Vries equation

The aim of this work is to investigate gravity-capillary waves resonantl...
03/06/2021 ∙ by Marcelo V. Flamarion, et al. ∙ 0

• ### A continuation method for building invisible obstacles in waveguides

We consider the propagation of acoustic waves at a given wavenumber in a...
05/14/2020 ∙ by Antoine Bera, et al. ∙ 0

• ### The Whitham Equation with Surface Tension

The viability of the Whitham equation as a nonlocal model for capillary-...
02/20/2020 ∙ by Evgueni Dinvay, et al. ∙ 0

• ### An approximation method for electromagnetic wave models based on fractional partial derivative

The present article is devoting a numerical approach for solving a fract...
02/12/2021 ∙ by Vijay Kumar Patel, et al. ∙ 0

• ### Source Seeking in Unknown Environments with Convex Obstacles

Navigation tasks often cannot be defined in terms of a target, either be...
09/16/2019 ∙ by Bruno A. Angelico, et al. ∙ 0

• ### Reduced Sum Implementation of the BURA Method for Spectral Fractional Diffusion Problems

The numerical solution of spectral fractional diffusion problems in the ...
05/19/2021 ∙ by Stanislav Harizanov, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The on-surface radiation conditions (OSRCs), originally developed by Kriegsmann, Taflove and Umashankar Kriegsmann1987, are semi-analytical methods to approximate the solution of wave scattering problems. In general, OSRCs lead to rough approximations, meant to sacrifice accuracy in favor of tremendous computational speed. This property is useful to explore parameter spaces for optimization problems, to produce a good initial guess for iterative methods, and to design inexpensive preconditioners to solve boundary integral equations (BIE) numerically using Krylov subspace methods Antoine2004; Antoine2005; Antoine2007; Darbas2013. The OSRCs have been studied by many researchers including Antoine and collaborators who formulated them in a differential geometric setting in order to handle non-canonical surfaces Antoine1999; Antoine2001; Antoine2001b; Antoine2006; Antoine2008. Other improvements have been accomplished by Jones Jones1988; Jones1990; Jones1992, Ammari Ammari1998; Ammari1998b, Calvo et al. Calvo2004; Calvo2003, Barucq et al. Barucq2003; Barucq2010a; Barucq2012, Chaillat et al. Chaillat2015; Chaillat2017 and Darbas et al. Antoine2004; Antoine2005; Antoine2007; Darbas2013; Darbas2006. See also Atle2007; Acosta2015f; Stupfel1994; Murch1993; Teymur1996; Medvinsky2010; Chniti2016; Alzubaidi2016; Acosta2015f; Acosta2017c.

The underlying assumption in most of the aforementioned works is that the solution radiates outwardly at every point on the surface of the scatterer. This assumption may be violated when the scatterer is not convex since the wave field may exhibit reflections from the scatterer to itself. Hence, one of the main drawbacks of the OSRC method is that its formulation and accuracy are limited to convex scatterers. This issue was addressed by the author in Acosta2015f and by Alzubaidi, Antoine and Chniti in Alzubaidi2016 for multiple scattering problems where the scatterer consists of several disjoint obstacles. The formulations in Acosta2015f; Alzubaidi2016 successfully account for the propagation of waves from one obstacle to another. However, the wave propagation formulas are non-local. Due to the use of boundary integral operators, integration of the wave fields over the surface of the obstacles is required. As a consequence, matrices with dense blocks are obtained upon discretization, which leads to high memory demands and costly matrix inversion, especially in the three-dimensional setting and also for high frequencies.

Our main goal is to modify the formulation of the OSRC for multiple obstacles in such a way that its discretization leads to sparse matrices. This is accomplished by defining a local wave propagation formula to transmit the influence of one obstacle on another using differential rather than integral operators. In Section 2 we provided the mathematical formulation of the multiple scattering problem and decompose it into a system of single-scattering problems. Based on Antoine1999; Acosta2015f; Acosta2017c, we derive a local formula for the propagation of waves which allows us to account for the wave reflections between obstacles. This formula is developed in Section 3. In Section 4 we use the propagation formula to evaluate the far-field pattern using local operations as well. In Section 5, we propose how to numerically implement the OSRC based on a triangulation of the surfaces. A few numerical results and comparisons with exact solutions are shown in Sections 6. Limitations, future work and conclusions are discussed in Section 7.

## 2 Mathematical formulation

We seek to approximate an outgoing solution to the Helmholtz equation in the exterior to a set of impenetrable convex obstacles. Each obstacle is enclosed by a closed smooth surface that separates the space into a simply connected bounded domain and an exterior domain . We also define

 Ω−=J⋃j=1Ω−j,Ω+=J⋂j=1Ω+jand∂Ω=J⋃j=1∂Ωj. (1)

The wave field is assumed to satisfy the following boundary value problem,

 Δu+k2u=0in Ω+, (2a) u=fon ∂Ω, (2b) limr→∞r(∂ru−ιku)=0, (2c)

where the wavenumber is constant, is the imaginary unit, , and is the boundary of . The proposed approach to approximate the solution of the boundary value problem (2) is based on the following theorem whose objective is to express the solution as the sum of purely-outgoing wave fields each radiating from a single surface for . A proof is found in Balabane2004.

###### Theorem 1.

Let the closures of and of be mutually disjoint for , and let be a radiating solution to the Helmholtz equation in . Then can be uniquely decomposed into purely-outgoing wave fields for such that

 u=J∑j=1uj,in Ω+, (3)

where radiates only from , that is,

 Δuj+k2uj=0 in Ω+j,% andlimr→∞r(∂ruj−ιkuj)=0. (4)

Using Green’s theorem for each of the purely-outgoing wave fields , we obtain a representation of the solution

 u(\bf x)=J∑j=1∫∂Ωj[∂ν(\bf y)Φ(\bf x,\bf y)uj(\bf y)−Φ(\bf x,\bf y)∂νuj(\bf y)]dS(\bf y),\bf x∈Ω+, (5)

where is the fundamental solution to the Helmholtz equation. If the exact Dirichlet-to-Neumann (DtN) map on each surface is available, then and we obtain

 u(\bf x)=J∑j=1∫∂Ωj[∂ν(\bf y)Φ(\bf x,\bf y)−Φ(\bf x,\bf y)Λj]uj(\bf y)dS(\bf y),%x∈Ω+. (6)

The representation (6) renders the solution in provided that the Dirichlet data of each purely-outgoing wave field is known. This data can be found by imposing the boundary condition (2b) on each boundary which leads to the following system of equations,

 ui+∑j≠iPi,juj=fion ∂Ωi,i=1,2,...,J (7)

where is the evaluation of on the surface , and where the propagation of the wave field from the surface to the surface for is given by the following operator,

 (Pi,juj)(\bf x)=∫∂Ωj[∂ν(\bf y)Φ(\bf x,\bf y)−Φ(\bf x,\bf y)Λj]uj(\bf y)dS(% \bf y),\bf x∈∂Ωi. (8)

For practical purposes, the evaluation of the propagation operator in (8) has two main challenges. First, the DtN map must be evaluated or at least approximated which is the objective of OSRC in general. Second, the propagation operator involves integration over the boundary of each obstacle. This leads to the appearance of dense blocks once the governing system (7) is discretized. Our objective in the following section is to mitigate this latter problem.

## 3 Local formula for propagation of wave fields

We seek to derive an analytical formula for the propagation of a wave field away from a generic convex surface . This formula will allow us to account for the interaction between obstacles in order to approximate the system (7) using sparse matrices. This formula is for the propagation of waves on a expansive foliation of parallel surfaces generated by the smooth surface . We follow closely the approach from Antoine1999. The domain exterior to is denoted . If is a convex surface, then it generates a family of parallel surfaces parametrized by defined as follows

 Γs={\bf y=\bf x+s\bf n(% \bf x):\bf x∈Γ} (9)

where

is the outward normal vector of the surface

at the point . Notice . Since is convex, the family of surfaces foliates . For any there exists unique and unique such that . Moreover, the outward normal vector of the surface at a point coincides with the normal vector . See details in (Kress-Book-1999, §6.2), (Gil-Tru-2001, §14.6), (DoCarmo1976, Probl. 11 §3.5) and Antoine1999. The point can be regarded as the projection of y on the surface , and is the distance from x to y (Deu-2001, Ch. 2). The pair satisfies

 \bf x=argmin\bf z∈Γ|\bf z−\bf y|ands=min% \bf z∈Γ|\bf z−\bf y|. (10)

Now we proceed to write the Helmholtz differential operator in terms of a tangential system of coordinates on the surface

 L(s)w=Δw+k2w=∂2sw+2Hs∂sw+ΔΓsw+k2w. (11)

Here is the mean curvature of the surface and is the Laplace–Beltrami operator on the same surface. The differential represents the derivative in the outward normal direction on . See Antoine1999 for a concise review of the differential geometry of surfaces, including the definition of curvatures and the Laplace–Beltrami operator. The mean curvature and Gauss curvature of the surface satisfy the following relations,

 Hs(\bf x)=H(\bf x)+sK(\bf x)μs(\bf x)andKs(\bf x)=K(\bf x)μs(\bf x) (12)

where the symbol stands for

 μs(\bf x)=det(I+sR(% \bf x))=1+2sH(\bf x)+s2K(\bf x). (13)

Here and are mean and Gauss curvature of the surface , respectively. Also,

is the curvature tensor and

is the identity on the tangent plane. The curvature tensor

can be represented as a self-adjoint matrix with eigenvalues

and , called the principal curvatures, such that

 H=κ1+κ22andK=κ1κ2.

The Laplace-Beltrami operator satisfies

where this latter approximation is valid when so that .

Now we seek to decompose the wave field propagating across a surface into incoming and outgoing components using Nirenberg’s factorization theorem Nirenberg1973; Antoine1999. There are two pseudo–differential operators and of order , such that

 L(s)w=(∂s−Λ−(s))(∂s−Λ+(s))w. (15)

We refer to and as the outgoing and incoming DtN operators for the radiating boundary value problem defined in the exterior of . The operator admits the following second order approximation derived in Antoine1999,

 Λ+(s)=ιk−Hs−12ιk(ΔΓs+H2s−Ks)+O(k−2). (16)

If is purely-outgoing from the surface , then . Using the approximation (16) for the outgoing DtN operator, we obtain that satisfies the following differential equation in ,

 ∂w∂s≈ιkw−Hsw−12ιk(ΔΓs+H2s−Ks)w. (17)

This is a separable differential equation. In order to solve it approximately in closed-form, we need the following expressions

 ∫s0Hzdz=12lnμs ∫s0[H2z−Kz]dz≈−12[sK1+sH+Hs−H] ∫s0ΔΓzdz≈s1+sHΔΓ

where we have employed the approximation in (14) for the third integral. Hence, we obtain

 ln(wwo) ≈ιks−12lnμs−12ιk(s1+sHΔΓ−12[sK1+sH+Hs−H])

 (18)

where we have neglected the tangential derivatives of and in order to factor the exponential. Here represents the trace of on the boundary . In order to maintain the locality of the evaluation, we further approximate the exponential of the Laplace–Beltrami operator using an implicit linear approximation valid for large frequency ,

 exp(−sΔΓ2ιk(1+sH))≈(1+sΔΓ2ιk(1+sH))−1. (19)

This implicit approximation is chosen to preserve the stability (boundedness) of the operations. Notice that the implementation of (19

) corresponds to solving an elliptic partial differential equation on the surface

.

Now we go back to the multi-scattering problem with its notation formulated in Section 2. In order to propagate the wave field from the surface onto the surface , we approximate the propagator (8) as follows,

 (Pi,juj)(\bf y)≈eιksμs(\bf x)1/2exp(14ιk[sK(\bf x)1+sH(\bf x)+Hs(\bf x)−H(\bf x)])(1+sΔ∂Ωj2ιk(1+sH(\bf x)))−1uj(\bf x) (20)

where . Due to the assumed convexity of , the point and are determined uniquely by satisfying (10). The curvatures and , and the Laplace–Beltrami operator correspond to the surface .

## 4 Far-Field Pattern

It is commonly of interest to evaluate the radiating field away from the obstacles. The asymptotic behavior is characterized by the far-field pattern of the radiating field which is given by

 u(\bf y)=eikrru∞(^\bf y% )+O(r−1),as r→∞, (21)

where and . From the superposition of the purely-outgoing fields, we find that

 u∞=J∑j=1u∞j (22)

so that it only remains to find the far-field pattern of each purely-outgoing field . We proceed with the following asymptotic expansions valid as ,

 s2 =r2−2r\bf x⋅^\bf y+|\bf x|2 (23a) s =r−\bf x⋅^\bf y+O(r−1) (23b) 1μ1/2s =1rK1/2+O(r−2) (23c) sK1+sH+Hs−H (23d) s1+sH =1H+O(r−1) (23e)

Using (23) and (18), we find that the purely-outgoing wave field has the following asymptotic behavior

 uj(r^\bf y)=eιkrr⎡⎣e−ιk\bf x⋅^\bf yK(\bf x)1/2exp(K(\bf x)−H(\bf x)24ιkH(\bf x))(1+Δ∂Ωj2ιkH(\bf x))−1uj(\bf x)+O(r−1)⎤⎦

which implies that the far-field pattern corresponding to is given by

 u∞j(^\bf y)=e−ιk\bf x⋅^\bf yK(\bf x)1/2exp(K(\bf x)−H(\bf x)24ιkH(\bf x))(1+Δ∂Ωj2ιkH(\bf x))−1uj(\bf x) (24)

where is uniquely determined by such that . Hence, once the Dirichlet data of the field is found by solving the system (7), then plugging (24) into (22) renders the far-field pattern .

## 5 Discrete Implementation

In this section, we propose a numerical implementation of the on-surface radiation condition defined by the system (7) where the propagation operator is approximated by (20). The approach is based on approximating the Laplace–Beltrami operator, the mean and Gauss curvatures using a triangulation of the surfaces for . Our guiding references for approximating geometrical properties and operators on triangulated surfaces are Xu2006a; Xu2004; Xu2004a; Meyer2003; Magid2007; Bobenko2007.

### 5.1 Discrete Laplace–Beltrami operator and curvatures

The discrete Laplace–Beltrami operator is described as follows. Let be the collection of vertices of the triangulation of the surface . For a fixed vertex , let be the neighboring vertices of , and let be the edges of the triangulation that connect and . Now for each edge , let and be the angles opposing the edge in the two triangles that share that edge. For a smooth function defined on the surface , we use the following discrete Laplace–Beltrami operator,

 (ΔΓu)(\bf xj)=12∑icotαj(i)+cotβj(i)(AjAj(i))1/2(u(\bf xj(i))−u(\bf xj)), (25)

where is the area associated with the vertex defined as 13 of the total area of the triangular elements sharing the vertex .

For the discrete mean curvature, we first define the discrete mean curvature vector as

 H(\bf xj)=−14∑icotαj(i)+cotβj(i)(AjAj(i))1/2(\bf xj(i)−\bf xj). (26)

Then the mean curvature is given by

 H(\bf xj)=H(\bf xj)⋅n(\bf xj) (27)

where is the normal vector at the vertex .

The discrete Gauss curvature at each vertex of the triangulation, is defined as follows Xu2006a; Magid2007; Surazhsky2003. Let be the angle between two successive edges sharing the vertex . Then, the approximate Gauss curvature is given by

 K(\bf xj)=2π−∑iθiAj (28)

where is the area associated with the vertex as defined above.

### 5.2 Discrete wave propagator

Given the triangulation of the surface , we have the necessary definitions in Subsection 5.1 to implement a discrete version of the propagator (20) to pose the system (7) at the discrete level. The discrete operator propagates the waves from the triangulation of the surface to the triangulation of the surface . Let the triangulations and have and vertices, respectively. Then the propagation operator can be represented as a matrix . We proceed to describe the construction of this matrix. Let where is the index in the list of vertices of . Then we can find and as the discrete version of (10), that is,

 m=argmin1≤p≤Mj|\bf xp−\bf yn|andsm=|\bf xm−\bf yn|. (29)

Then the matrix can be written as

 Pi,j=Ai,j(I+ΔTj2ιkH)−1 (30)

where is an matrix defined by (25) which approximates the Laplace–Beltrami operator, and is an matrix. Each row of this matrix has exactly one non-zero entry. For the row, this non-zero entry is at the column and it is given by

 (Ai,j)n,m (31)

where depends on by satisfying (29). With this notation, we can write the discrete version of the system (7) in block-matrix form as follows,

 (32)

where the total number of degrees of freedom is

.

### 5.3 Complexity O(N)

The standard numerical techniques for BIE, such as Nystrom, collocation and boundary element methods, lead to fully populated matrices. Similarly, for the OSRCs developed in Acosta2015f; Alzubaidi2016, matrices with fully populated off-diagonal blocks are obtained. In these cases, high computational costs and large memory demands are needed for matrix-vector multiplication in any iterative solver (e.g. fixed-point or GMRES). By contrast, the proposed OSRC leads to complexity. This follows from the propagator matrix (30) that defines the governing system (32). The implementation of the propagator (30) requires the multiplication by the matrix after the inversion of a system governed by . The former is a sparse matrix, with exactly one non-zero entry per row, leading operations per matrix-vector multiplication. The latter is the solution for the discrete version of a two-dimensional elliptic equation. With proper sorting of the triangulation nodes, this system is both sparse and narrowly banded. Therefore, a direct solver such as LU decomposition requires operations. As a consequence, any iterative solver for the system (32) will require operations per iteration. Figure 1 displays (a) the sparsity pattern of the matrix for the configuration of two spheres from subsection 6.1, (b) the sparsity pattern of the discrete Laplace-Beltrami operator , and (c) CPU time for the application of the propagator (30) as a function of DOF.

For the numerical results presented in the next section, we apply a fixed-point iteration to solve the system (32). Due to the convexity of the surfaces for all , the norm of the propagation operator decreases as the distance between the and obstacles increases. Hence, the system (32) can be solved using the following fixed-point iteration,

 (33)

for vanishing initial guess at .

## 6 Numerical Results

In this section we present a few numerical results obtained from the implementation of the discrete formulation defined in Section 5. We consider three convex surfaces for the numerical experiments presented here. These surfaces are the unit sphere, a rounded cube, and a marshmallow–like surface. Their defining equations are as follows,

 Sphere: |x|2+|y|2+|z|2=1. (34a) Rounded Cube: |x|3+|y|3+|z|3=1. (34b) Marshmallow: |x|2+|y|2+|z|4=1. (34c)

Coarse triangulations of these three surfaces are shown in Figure 2. The approximate mean and Gauss curvatures of the rounded cube and of the marshmallow, defined by (26)-(27), are shown in Figure 3.

For comparison, we manufacture an exact solution to the boundary values problem (2), valid for any geometry of the surfaces in . This is accomplished by defining boundary data as the boundary trace of a radiating wave field. We consider the field,

 F(\bf x)=J∑j=1Φ(\bf x−\bf cj) (35)

where is the outgoing fundamental solution to the Helmholtz equation with frequency . Hence, represents the superposition of point-sources with respective locations at . If each point is enclosed by the respective surface , for , then the exact solution to the boundary values problem (2) with Dirichlet data is given by .

### 6.1 Example 1: Two spheres

Now we present some numerical results for the wave radiating problem in the exterior of unit-spheres. These spheres are centered at the follow points:

 \bf c1=(2,0,0)and\bf c2=(−2,0,0). (36)

Table 1 displays the -norm relative error for the far-field pattern for various wavenumbers and mesh refinements. The DOF approximately quadruples with each mesh refinement, which corresponds to halving the number of points per wavelength.

Figures 4 and 5 show the exact and numerical far-field patterns and the error between them, for wavenumbers and , respectively.

For visual comparison, cross-sections of the far-field pattern are shown Figures 6 and 7 for wavenumbers and , respectively.

### 6.2 Example 2: Three different shapes

Now we choose surfaces shown in Figure 2, translated to the respective points

 \bf c1=4(1,0,0),\bf c2=4(cos\sfrac2π3,sin\sfrac2π3,0),and\bf c3=4(cos\sfrac4π3,sin\sfrac4π3,0). (37)

Table 2 displays the -norm relative error for the far-field pattern for various wavenumbers and mesh refinements of these three shapes. The DOF approximately quadruples with each mesh refinement, which corresponds to halving the number of points per wavelength in any tangential direction.

For this example, we also display some cross-sections of the far-field pattern as shown Figures 10 and 11 for wavenumbers and , respectively.

## 7 Final Remarks

We have formulated an OSRC for multiple scattering problems that involves local operations only. As opposed to Acosta2015f; Alzubaidi2016, integration over the surfaces is avoided which leads to sparse matrices upon discretization and implementation with complexity. The proposed method is a local OSRC that simultaneously accounts for the outgoing behavior of the solution as well as the wave reflections between the multiple obstacles. We expect that this OSRC will provide an inexpensive, physics-based preconditioner to solve multiple scattering problems using boundary integral equations (BIE) via Krylov subspace methods. The derivation of adequate preconditioners for BIE is a mandatory requirement to solve large scale wave scattering problems Antoine2004; Antoine2005; Antoine2007; Darbas2013; Darbas2015; Chaillat2015; Chaillat2017.

The major limitations of the proposed method are similar to those of other OSRCs. The formulation is valid for convex surfaces and the OSRC provides a rough approximation of the solution without a-priori error estimates. Therefore, the OSRC should not be employed if a high degree of accuracy is needed. However, as seen from Table

1 and 2, and Figures 6-7 and 10-11, the proposed OSRC can produce qualitatively good approximations to the solution of the boundary value problem (2). Possible future work includes the extension to electromagnetic and elastic waves that govern important engineering applications. It is also possible to increase the order of approximations of the Dirichlet–to–Neumann map to improve the accuracy of the OSRC Acosta2017c. Here we have considered the second order differential approximation (16). The major challenge that still remains to be addressed is the formulation of an OSRC and the associated propagation formula for non-convex obstacles. The same is true for piecewise smooth surfaces in order to handle edges and corners.

## Acknowledgements

This work was partially supported by NSF grant DMS-1712725. The author would like to thank Texas Children’s Hospital for its support and for the research-oriented environment provided by the Predictive Analytics Laboratory.