On the Numerical Solution of Diffraction Problem by Random Spheres Using Electric Field Integral Equation

Sahar Malwi Al Malky^1*, Sharefa Eisa Ali Alhazmi²

¹Department of Mathematics, University College Adham, Makkah, Kingdom of Saudi Arabia.
²Department of Mathematics, University College in Qunfudah, Umm Al-Qura University, Makkah, Kingdom of Saudi Arabia.

DOI: 10.4236/am.2023.141003                

Abstract

The aim of this paper is to solve the two-dimensional acoustic scattering problems by random sphere using Electric field integral equation. Some approximations for the two-dimensional case are derived. These various approximations are next numerically validated in the case of high-frequency.

Keywords

Multiple Scattering, Integral Equation, High Frequency, Acoustic, Iterative Solver

1. Introduction

Wave scattering in the time harmonic regime is known to be an extremely active and difficult area of research with many possible applications in civil and military domains. The applications can require different physics related to wave scattering: acoustics, electromagnetic, elasticity. Applications can be of different natures: radar, sonar, aerospace and aeronautics design, medical imaging, underwater acoustics, diffraction gratings, optics. Therefore, this is a huge area with very exciting applications that need top level computational developments and emerging new ideas to solve the most advanced challenging problems. In particular, high-frequency scattering remains an unsolved numerical problem with many applications. An example of application requires very powerful numerical methods for the prediction of scattering phenomena. We can indeed see the surface electromagnetic field created on an aircraft by an incident wave field. This can be useful when for example you investigate Electro-Magnetic Compatibility (EMC) problems where electromagnetic waves can penetrate inside the aircraft and badly damage its electronic components. This is of first importance for aerospace companies for security reasons and prevention of accidents (but also for cost reduction). Considering such a problem is clearly extremely difficult and then numerical methods provide a good understanding of the physical problems at low cost. Another example in aeronautics is related to the control of noise generation when the aircraft lands near a town, therefore creating noise pollution and problems near airports since (hopefully) restrictive administration rules must (are expected to) be followed by the aircraft companies. One common difficult point to all these problems is that they are set in an unbounded domain making the design of a suitable, and accurate numerical method very challenging. In this paper we consider the solution of the waves equation in an unbounded domain. The homogeneous and isotropic exterior domain of propagation is ${\Omega }^{+}={ℝ}^2\backslash \overline{{\Omega }^{-}}$. For the sake of conciseness in the presentation, we first assume that the scatterers are sound-soft (Dirichlet boundary condition). We now consider a time-harmonic incident acoustic plane wave $u^{\text{inc}}\left(x\right)={\text{e}}^{ik\beta \cdot x}$ (with $x=\left(x_1\mathrm{,}{x}_2\right)\in {ℝ}^2$ ) illuminating ${\Omega }^{-}$, with an incidence direction $\beta =\left(\cos \left({\theta }^{\text{inc}}\right)\mathrm{,}\sin \left({\theta }^{\text{inc}}\right)\right)$ and a time dependence ${\text{e}}^{-i\omega t}$, where $\omega$ is the wave pulsation and $k=2\pi /\lambda$ is the wavenumber, $\lambda$ being the wavelength. It’s natural that there many scatterers in the domain of computation, this case is more difficult, see [1]. In this paper we consider M random, bounded and disjoint scatterers ${\Omega }^p$, $p=\mathrm{1,}\cdots \mathrm{,}M$, distributed in ${ℝ}^2$, with boundary ${\Gamma }_p:=\partial {\Omega }_p^{-}$. The scatterer ${\Omega }^{-}$ is defined as the collection of the M separate obstacles, i.e. ${\Omega }^{-}={\displaystyle {\cup }_{p=1}^M{\Omega }^p}$, with boundary $\Gamma ={\displaystyle {\cup }_{p=1}^M{\Gamma }_p}$. The sound-soft multiple scattering problem of $u^{\text{inc}}$ by ${\Omega }^{-}$ consists in computing the scattered wavefield u as the solution to the boundary-value problem [2] [3].

$\{\begin{array}{l}\left(\Delta +k^2\right)u=\mathrm{0,}\text{in}{\Omega }^{+}\mathrm{,}\\ u=-u^{\text{inc}}\mathrm{,}\text{on}\Gamma \mathrm{,}\\ \underset{\Vert x\Vert \to +\infty }{\lim }{\Vert x\Vert }^{1/2}\left(\nabla u\cdot \frac{x}{\Vert x\Vert }-iku\right)=0.\end{array}$ (1)

The last equation of (1) is the well-known Sommerfeld’s radiation condition at infinity that ensures the uniqueness of u [4] [5].

Multiple scattering is known to be a very complicated and challenging problem in terms of computational method [1] [2] [6] - [15] since the incident wave is multiply diffracted by all the single scatterers involved in the geometrical configuration. As a consequence, the scattered wave-field has a highly complicated structure and exhibits some particular physics properties. We limit our study to the case of scattering by circular cylinder. The structure of the paper is the following. In Section 2, we introduce some definitions that will be used in our paper. In Section 3, we compute an exact solution of the single scattering problem using Mie series expansion. In Section 4, a general introduction will be proposed to explain the scattering problem. In Section 5, we present some problems arising when solving numerically the multiple scattering problem at high-frequency. In Section 6, we introduce a spectral method used to approximate the integral operator, so a linear system will be presented and solved. Then we conclude.

2. Preliminary

2.1. Plane Waves

We seek first solutions depending only on the coordinated $x_1$ of the Helmoltz equation. We have

$u^{″}\left(x_1\right)+k^{2}u\left(x_1\right)=0$

This is a linear differential equation of the second order; its solutions thus form a vector space of dimension two; they are given by

$u\left(x\right)=A{\text{e}}^{ik{x}_1}+B{\text{e}}^{-ik{x}_1}$

It already appears in this very simple situation the existence of two types of waves, they are distinguished by their propagation direction. To fix the idea, suppose for example that we are interested to pressure; we have

$P\left(x,t\right)=\mathcal{R}\left(u\left(x\right){\text{e}}^{-i\omega t}\right)=\mathcal{R}\left(A{\text{e}}^{i\left(k{x}_1-\omega t\right)}\right)+\mathcal{R}\left(B{\text{e}}^{-i\left(k{x}_1+\omega t\right)}\right)$

The first type of waves $A{\text{e}}^{i\left(k{x}_1-\omega t\right)}$ are a progressive wave in the direction of increasing $x_1$.

More generally, a plane wave propagating in the direction of the $\beta$ will be as follows

$u\left(x\right)={\text{e}}^{ikx\cdot \beta }\mathrm{,}\Vert \beta \Vert =1$.

2.2. Hankel Function

The polar coordinates are defined as

$x_1=r\cos \theta$

$x_2=r\sin \theta$

where $r\ge 0$ and $\theta \in \left]-\pi \mathrm{,}\pi \right]$.

Let $\tilde{u}\left(r\mathrm{,}\theta \right)=u\left(x_1\mathrm{,}{x}_2\right)$, the Laplace operator in two dimensions is given by

$\Delta u=\frac{{\partial }^{2}u}{\partial x_1^2}+\frac{{\partial }^{2}u}{\partial x_2^2}$

where $x_1$ and $x_2$ are the standard Cartesian coordinates of the $x_1{x}_2$ plane. In polar coordinates

$\Delta u\left(x_1,x_2\right)=\left(\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \tilde{u}}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2\tilde{u}}{\partial {\theta }^2}\right)\left(r,\theta \right)$ (2)

$=\left(\frac{{\partial }^2\tilde{u}}{\partial r^2}+\frac{1}{r}\frac{\partial \tilde{u}}{\partial r}+\frac{1}{r^2}\frac{{\partial }^2\tilde{u}}{\partial {\theta }^2}\right)\left(r,\theta \right)$ (3)

If we seek for a solution of the Helmoltz equation independent of $\theta$ we are going to solve

$\left(\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial }{\partial r}\right)+k^2\right)\tilde{u}=0$.

Let $z=rk$ and $\tilde{u}\left(z\right)=\tilde{u}\left(r\mathrm{,}\theta \right)$ then $\tilde{u}\left(z\right)$ is solution of

$z^2\frac{{\partial }^2\tilde{u}\left(z\right)}{\partial z^2}+z\frac{\partial \tilde{u}\left(z\right)}{\partial z}+z^2\tilde{u}\left(z\right)=0$ (4)

The Equation (4) is a linear equation of the second order with variable coefficients whose solutions have been extensively studied. They are of the form

$\tilde{u}\left(z\right)=c_1{H}_0^1\left(z\right)+c_2{H}_0^2(\; z\; )$

and then

$u\left(r\right)=c_1{H}_0^1\left(rk\right)+c_2{H}_0^2(\; r\; k\; )$

where $H_0^1$ and $H_0^2$ are respectively the Hankel function of first and second kind of order zero.

We recall that $J_n$ and $Y_n$ are the Bessel functions of the first and second kinds, respectively. They satisfy Bessel’s Equation (5):

$x^2\frac{{\partial }^{2}y\left(x\right)}{\partial x^2}+x\frac{\partial y\left(x\right)}{\partial x}+\left(x^2-n^2\right)y\left(x\right)=0$ (5)

and the $H_n^1\left(x\right)\mathrm{,}{H}_n^2\left(x\right)$ also known as the Bessel functions of the third kind and the satisfy Bessel’s equation, and are related to $J_n$ and $Y_n$ by

$H_n^1\left(x\right)=J_n\left(x\right)+i{Y}_n\left(x\right),H_n^2\left(x\right)=J_n\left(x\right)-i{Y}_n\left(x\right).$

Therefore, if one looks for the Fourier series expansion of the solution as

$u\left(r\mathrm{,}\theta \right)={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{c}_m\left(r\right){\text{e}}^{im\theta }$

and the coefficient $c_m\left(r\right)$ are solution of

$\left(r^2\frac{{\partial }^2}{\partial r^2}+r\frac{\partial }{\partial r}+\left(k^2-m^2\right)\right)c_m=0$.

2.3. Sommerfeld Condition

The Radial solutions we have just calculated form a two-dimensional space; we will see that is possible to distinguish two categories of solutions that correspond to outgoing waves (or divergent) and incoming waves (or converged). This distinction is based on their asymptotic behavior in the neighborhood of infinity (as $r\to +\infty$ ). The first term of the asymptotic expansion of Hankel functions for large real arguments is given

$H_0^1\left(z\right)=\sqrt{\frac{2}{\pi z}}{\text{e}}^{i\left(z-\pi /4\right)}+O\left(\frac{1}{z^{3/2}}\right)$

$H_0^2\left(z\right)=\sqrt{\frac{2}{\pi z}}{\text{e}}^{-i\left(z-\pi /4\right)}+O\left(\frac{1}{z^{3/2}}\right)$.

Its obvious that the outgoing solution is given by $H_0^1$, the same conclusion can be obtained from the following hypothesis, called outgoing radiation condition or Sommerfeld condition

$\frac{\partial u\left(r\right)}{\partial r}-iku=o\left(\frac{1}{r^{1/2}}\right),r\to +\infty$.

Proposition 1. The outgoing radial waves are given by

$u\left(r\right)=c_1{H}_0^1\left(rk\right)$ (6)

Proof. Using

$\frac{\text{d}{H}_0^j\left(z\right)}{\text{d}z}=-H_1^j\left(z\right),j=1,2$

$H_1^1\left(z\right)=\sqrt{\frac{2}{\pi z}}{\text{e}}^{i\left(z-3\pi /4\right)}+O\left(\frac{1}{z^{3/2}}\right)$

$H_1^2\left(z\right)=\sqrt{\frac{2}{\pi z}}{\text{e}}^{-i\left(z-3\pi /4\right)}+O\left(\frac{1}{z^{3/2}}\right)$.

Then

$\begin{array}{c}\frac{\partial u\left(r\right)}{\partial r}=c_1\frac{\partial H_0^1\left(rk\right)}{\partial r}+c_2\frac{\partial H_0^2\left(rk\right)}{\partial r}\\ =-k\left(c_1{H}_1^1\left(rk\right)+c_2{H}_1^2\left(rk\right)\right)\\ =-k\sqrt{\frac{2}{\pi kr}}\left(c_1{\text{e}}^{i\left(kr-3\pi /4\right)}+c_2{\text{e}}^{-i\left(z-3\pi /4\right)}\right)+o(\; 1\; r\; )\end{array}$

on the other hand

$iku\left(r\right)=k\sqrt{\frac{2}{\pi kr}}\left(c_1{\text{e}}^{i\left(kr-\pi /4\right)}+c_2{\text{e}}^{-i\left(z-\pi /4\right)}\right)+o(\; 1\; r\; )$

then we have

$\frac{\partial u\left(r\right)}{\partial r}-iku\left(r\right)=-2i{c}_{2}k\sqrt{\frac{2}{\pi kr}}{\text{e}}^{-i\left(kr-\pi /4\right)}+o\left(\frac{1}{\sqrt{r}}\right)$.

Sommerfeld condition gives $c_2=0$. Then the outgoing solution are given by

$u\left(r\right)=c_1{H}_0^1\left(rk\right)$.

2.4. Integral Formulation

Let G be the two-dimensional free-space Green’s function defined by

$\forall x\mathrm{,}y\in {ℝ}^2\mathrm{,}x\ne y\mathrm{,}G\left(x\mathrm{,}y\right)=\frac{i}{4}{H}_0^{\left(1\right)}\left(k\Vert x-y\Vert \right)\mathrm{.}$

The function $H_0^{\left(1\right)}$ is the first-kind Hankel function of order zero. Integral equations are essentially based upon the Helmholtz integral representation formula ( [4], Theorems 3.1 and 3.3).

Proposition 2. If v is a solution to the Helmholtz equation in the unbounded connected domain ${\Omega }^{+}$ and satisfies the Sommerfeld radiation condition, then we have

${\displaystyle {\int }_{\Gamma }}-G\left(x\mathrm{,}y\right){\partial }_{n}v\left(y\right)+{\partial }_{n_y}G\left(x\mathrm{,}y\right)v\left(y\right)\text{d}\Gamma \left(y\right)=(\begin{array}{ll}v\left(x\right)\hfill & \text{if}x\in {\Omega }^{+}\mathrm{,}\hfill \\ 0\hfill & \text{otherwise}\text{.}\hfill \end{array}$ (7)

The unit normal vector $n$ is outwardly directed to ${\Omega }^{-}$. The integrals on $\Gamma$ must be understood as duality brackets between the Sobolev space $H^{1/2}\left(\Gamma \right)$ and its dual space $H^{-1/2}\left(\Gamma \right)$. Nevertheless, when the incident wavefield $u^{\text{inc}}$ and the curve $\Gamma$ are sufficiently smooth, the scattered field is then regular and the duality bracket can be identified (this is systematically the case in the presentation) to the (non-Hermitian) inner product in $L^2(\; \Gamma \; )$

${\langle f\mathrm{,}g\rangle }_{H^{-1/2}\mathrm{,}{H}^{1/2}}={\displaystyle {\int }_{\Gamma }}fg\text{d}\Gamma \mathrm{.}$

Let us now introduce the volume single- and double-layer integral operators, respectively denoted by $\mathcal{L}$ and $\mathcal{M}$, and defined by: $\forall x\in {ℝ}^2\backslash \Gamma$

$\begin{array}{l}\mathcal{L}\mathrm{<mo>\; :\; </mo>}\rho \mapsto \mathcal{L}\rho \left(x\right)={\displaystyle {\int }_{\Gamma }}G\left(x\mathrm{,}y\right)\rho \left(y\right)\text{d}\Gamma \left(y\right)\mathrm{,}\\ \mathcal{M}\mathrm{<mo>\; :\; </mo>}\lambda \mapsto \mathcal{M}\lambda \left(x\right)=-{\displaystyle {\int }_{\Gamma }}{\partial }_{n_y}G\left(x\mathrm{,}y\right)\lambda \left(y\right)\text{d}\Gamma \left(y\right)\mathrm{.}\end{array}$

We can then express the wavefields v as

$v\left(x\right)=-\mathcal{L}\left({{\partial }_{n}v|}_{\Gamma }\right)\left(x\right)-\mathcal{M}\left({v|}_{\Gamma }\right)\left(x\right)\mathrm{,}\forall x\in {\Omega }^{+}\mathrm{,}$

Furthermore, the single- and double-layer integral operators provide some outgoing solutions to the Helmholtz equation [16].

Proposition 3. For any densities $\rho \in H^{-1/2}\left(\Gamma \right)$ and $\lambda \in H^{1/2}\left(\Gamma \right)$ , the functions $\mathcal{L}\rho$ , $\mathcal{M}\lambda$ and any linear combination of them are some outgoing solutions to the Helmholtz equation in ${ℝ}^2\backslash \Gamma$ for some boundary conditions.

We now recall the expressions of the trace and normal derivative trace of the volume single- and double-layer potentials which are commonly called jump relations ( [16], Theorem 3.1).

Proposition 4. For any $x$ in $\Gamma$ , the trace and normal derivative traces of the operators $\mathcal{L}$ and $\mathcal{M}$ are given by the following relations (the signs indicate that z tends towards x from the exterior or the interior of $\Gamma$ )

$\underset{z\in {\Omega }^{\pm }\to x}{lim}\mathcal{L}\rho \left(z\right)=L\rho \left(x\right)\mathrm{,}\underset{z\in {\Omega }^{\pm }\to x}{lim}\mathcal{M}\lambda \left(z\right)=\left(\mp \frac{1}{2}I+M\right)\lambda \left(x\right)\mathrm{,}$ (8)

where I is the identity operator, for $x\in \Gamma$,

$L\rho \left(x\right)={\displaystyle {\int }_{\Gamma }}G\left(x\mathrm{,}y\right)\rho \left(y\right)\text{d}\Gamma \left(y\right)\mathrm{,}M\lambda \left(x\right)=-{\displaystyle {\int }_{\Gamma }}{\partial }_{n_y}G\left(x\mathrm{,}y\right)\lambda \left(y\right)\text{d}\Gamma \left(y\right)\mathrm{.}$

Throughout the paper, the boundary integral operators are denoted by a roman letter (e.g. L) while the volume integral operators use a calligraphic letter (e.g. $\mathcal{L}$ ). The operator $M^{*}=-N$ is the adjoint operator of M, that is

${\langle g\mathrm{,}Mf\rangle }_{H^{-1/2}\mathrm{,}{H}^{1/2}}={\langle -Ng\mathrm{,}f\rangle }_{H^{-1/2}\mathrm{,}{H}^{1/2}}\mathrm{,}\forall \left(f\mathrm{,}g\right)\in H^{1/2}\left(\Gamma \right)\times H^{-1/2}\left(\Gamma \right)\mathrm{.}$

Other properties like compactness or invertibility of integral operators can also be stated [2] [4] [5].

3. Single Scattering Using Mie Series

We consider simple scattering geometries in the two-dimensional domain. The goal of this section is to precise the case of one single disk, multiple scattering configurations being harder to treat. Let us consider a single disk centered at the origin and with radius R. The polar coordinates system is designated by $\left(r\mathrm{,}\theta \right)$. We consider an incident plane wave $u^{\text{inc}}$ defined by $u^{\text{inc}}\left(x\right)={\text{e}}^{ik\beta \cdot x}$ where $\beta ={\left(\cos \left({\theta }^{\text{inc}}\right)\mathrm{,}\sin \left({\theta }^{\text{inc}}\right)\right)}^{\text{T}}$, where T indicate the transpose of the vector. The angle of incidence ${\theta }^{\text{inc}}$ is given in the polar system. Plane waves are standard physically admissible incident wave fields that are used in practice. In the polar coordinates system, we then have

$\begin{array}{c}{u}^{inc}={\text{e}}^{ik\left(x_1\cos \left({\theta }^{\text{inc}}\right)+x_2\sin \left({\theta }^{\text{inc}}\right)\right)}\\ ={\text{e}}^{ikr\left(\cos \left(\theta \right)\cos \left({\theta }^{\text{inc}}\right)+\sin \left(\theta \right)+\sin \left({\theta }^{\text{inc}}\right)\right)}\\ ={\text{e}}^{ikr\cos \left(\theta -{\theta }^{\text{inc}}\right)}\end{array}$

Carl Gustav Jacob Jacobi (1804-1851) was a German mathematician by using standard trigonometric formulae prove the following proposition.

Proposition 5. The Jacobi expansion

${\text{e}}^{i\omega \cos \left(\varphi \right)}={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{\left(-i\right)}^m{J}_m\left(\omega \right){\text{e}}^{im\varphi }$

where $J_m$ are the Bessel special functions of the first kind.

Then, one wants to compute the solution to the Helmholtz equation in the polar coordinates system. The scattered field is given by

$u\left(r,\theta \right)={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{c}_n^1{H}_n^1\left(kr\right){\text{e}}^{in\theta }+c_n^2{H}_n^2\left(kr\right){\text{e}}^{in\theta }$

where we consider that ${\theta }^{\text{inc}}=0$.

The new special functions $H_n^1$ and $H_n^2$ are respectively called the Hankel function of the first and second kind of order n. They are defined by

$H_n^1\left(x\right)=J_n\left(x\right)+i{Y}_n(\; x\; )$

$H_n^2\left(x\right)=J_n\left(x\right)-i{Y}_n(\; x\; )$

In our case we solve the Dirichlet problem, so the boundary condition applied on the scatterer surface $r=R$ (R the radii of the disc), we get

$c_n^1{H}_n^1\left(kR\right)+c_n^2{H}_n^2\left(kR\right)=-{\left(-i\right)}^n{J}_n(\; k\; R\; )$

Proposition 6. The outgoing solution of (9) is given by

$u\left(r,\theta \right)=-{\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{\left(-i\right)}^n\frac{J_n\left(kR\right)}{H_n^1\left(KR\right)}{H}_n^1\left(kr\right){\text{e}}^{in\theta };$

Proof. The derivative of the Hankel function is given as follows:

$\frac{\text{d}{H}_n^j\left(z\right)}{\text{d}z}=-H_{n+1}^1\left(z\right)+\frac{n}{z}{H}_n^1\left(z\right),j=1,2$

The asymptotic expansions of the Hankel function are given by

$H_n^1\left(z\right)=\sqrt{\frac{2}{\pi z}}{\text{e}}^{i\left(z-n\pi /2-\pi /4\right)}+O\left(\frac{1}{z^{3/2}}\right)$

$H_n^2\left(z\right)=\sqrt{\frac{2}{\pi z}}{\text{e}}^{-i\left(z-n\pi /2-\pi /4\right)}+O\left(\frac{1}{z^{3/2}}\right)$

$\begin{array}{l}\frac{\partial u\left(r\mathrm{,}\theta \right)}{\partial r}-iku\left(r\mathrm{,}\theta \right)\\ =k\left(c_n^1\frac{\left(1-i\right)n{\text{e}}^{-i\pi n/2}{\text{e}}^{ikr}}{kr\sqrt{\pi }\sqrt{kr}}+c_n^2\frac{\left(1-i\right){\text{e}}^{i\pi n/2}\left(2kr+in\right)}{kr\sqrt{\pi }{\text{e}}^{ikr}\sqrt{kr}}\right)+O\left(\frac{1}{r^{3/2}}\right)\end{array}$

Sommerfeld’s radiation condition at infinity:

$\underset{r\to +\infty }{lim}{r}^{1/2}\left(\frac{\partial u\left(r\mathrm{,}\theta \right)}{\partial r}-iku\left(r\mathrm{,}\theta \right)\right)=0$

lead to

$\underset{r\to +\infty }{lim}\frac{\left(1-i\right)n{\text{e}}^{-1/2i\pi n}{\text{e}}^{ikr}}{\sqrt{r}k\sqrt{\pi }\sqrt{kr}}{c}_n^1+\underset{r\to +\infty }{lim}\frac{\left(1-i\right){\text{e}}^{1/2i\pi n}\left(2kr+in\right)}{\sqrt{r}k\sqrt{\pi }{\text{e}}^{ikr}\sqrt{kr}}{c}_n^2=0$

then we have

$\underset{r\to +\infty }{lim}{c}_n^2\frac{\left(2-2i\right){\text{e}}^{1/2i\pi n}}{\sqrt{\pi }{\text{e}}^{ikr}\sqrt{k}}=0\Rightarrow c_n^2=0$

then the solution is given by

$u\left(r\mathrm{,}\theta \right)={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{c}_n^1{H}_n^1\left(kr\right){\text{e}}^{in\theta }\mathrm{<mo>\; ;\; </mo>}$

To find $c_n^1$ we apply the boundary condition at $r=R$, we get

$c_n^1=-{\left(-i\right)}^n\frac{J_n\left(kR\right)}{H_n^1(\; K\; R\; )}$

Considering the case of single scattering of an incident plane wave ( $k=10$ and $\beta =0$ (rad.) (coming from the right)) by sound-soft obstacle, we study diffracted field and total field produced, this is depicted in Figure 1.

4. Multiple Scattering Problem

Multiple diffraction describes the phenomenon of interaction between a wave and a medium of propagation containing two or more obstacles (see Figure 2 for multiple diffraction by obstacles ${\Omega }_j,j=1,\cdots ,6$ ) as opposed to simple diffraction where only an obstacle is present in the propagation medium. A similar definition of single diffraction, in which the effects of multiple scattering are completely ignored: the total scattered field is simply the sum of the fields scattered by the various obstacles, each of which is illuminated by an incident field in isolation from other obstacles. This approximation is widely used; its validity requires a very large spacing between the obstacles compared to the size of the obstacles and the length of the incident waves. Several mathematical and numerical difficulties arise in the typical problem of multiple diffraction which arises in several applications in acoustics, electromagnetism, elasticity or hydrodynamics.

Figure 1. Single scattering of an incident plane wave ( $k=10$ and $\beta =0$ (rad.) (coming from the right)) by sound-soft obstacle. (a) Real part of the scattered field; (b) Imaginary part of the scattered field; (c) Absolute value of the scattered field; (d) Real part of the total field; (e) Imaginary part of the total field; (f) Absolute value of the total field.

Figure 2. Diffraction by circular cylinders.

Consider a homogeneous acoustic medium filling the whole space ${ℝ}^2$ and containing M disjoint scatterers ${\Omega }_1\mathrm{,}\cdots \mathrm{,}{\Omega }_M$. We assume that each scatterer ${\Omega }_p$, for $p=1,\cdots ,M$, is a bounded subdomain of ${ℝ}^2$ of boundary $\partial {\Omega }_p$. We denote by ${\Omega }_{-}={\displaystyle \underset{p=1}{\overset{M}{\cup }}{\Omega }_p}$ the domain occupied by the scattering obstacles. We consider the scattering problem of an incident plane wave $u^{inc}$ by ${\Omega }_{-}$. In other words, we want to determine the scattered wave u, which is the solution of the exterior boundary value problem

$\left(\mathcal{E}\right)\{\begin{array}{l}\left(\Delta +k^2\right)u=\mathrm{0,}\text{in}{\Omega }^{+}\mathrm{,}\\ u=-u^{\text{inc}}\mathrm{,}\text{on}\Gamma \mathrm{,}\\ \underset{\Vert x\Vert \to +\infty }{lim}{\Vert x\Vert }^{1/2}\left(\nabla u\cdot \frac{x}{\Vert x\Vert }-iku\right)=0.\end{array}$ (9)

In order to solve numerically the multiple scattering problem $\left(\mathcal{E}\right)$, a natural idea is to reduce it to a family of single scattering problems. By linearity, this can be achieved by introducing M fictitious scattered waves $u^1\mathrm{,}\cdots \mathrm{,}{u}^M$, where each field $u^p$ corresponds to the wave reflected only by the scatterer p when it is illuminated by the incident wave and the scattered waves $u^q$, for $q=1,\cdots ,M$ with $q\ne p$. More precisely, we have the following result.

Theorem 1. Let u be the solution of the multiple scattering problem $\left(\mathcal{E}\right)$ . Then, the family of M coupled single scattering problems for $p=1,\cdots ,M$ :

$\left({\mathcal{E}}^p\right)\{\begin{array}{ll}\Delta u^p+k^2{u}^p=0\hfill & \left({ℝ}^2\backslash {\bar{\Omega }}_p\right)\hfill \\ u^p=-\left(u^{\text{inc}}+{\displaystyle \underset{q=1,q\ne p}{\overset{M}{\sum }}}{u}^q\right)\hfill & \left(\partial {\Omega }_p\right)\hfill \\ \underset{\left|x\right|\to +\infty }{lim}{\left|x\right|}^{1/2}\left(\nabla u\cdot \frac{x}{\left|x\right|}-iku\right)=0.\hfill & \hfill \end{array}$

admitsa unique solution $\left(u^1\mathrm{,}\cdots \mathrm{,}{u}^p\right)$ .Furthermore,the following decomposition holds true:

$u={\displaystyle \underset{p=1}{\overset{M}{\sum }}}{u}^p.$

Proof. This theorem was proved in [17]. Here we give an alternative proof for the case of Dirichlet boundary conditions.

We start by proving the uniqueness of the solution of the coupled problems ${\left({\mathcal{E}}^p\right)}_{p=1,M}$. Let then $u^1\mathrm{,}\cdots \mathrm{,}{u}^M$, with $u^p\in H_{\mathcal{l}oc}^1\left({ℝ}^2\backslash \overline{{\Omega }_p}\right)$, satisfy

$\{\begin{array}{ll}\Delta u^p+k^2{u}^p=0\hfill & \left({ℝ}^2\backslash {\bar{\Omega }}_p\right)\hfill \\ u^p=-{\displaystyle \underset{q=1,q\ne p}{\overset{M}{\sum }}}{u}^q\hfill & \left(\partial {\Omega }_p\right)\hfill \\ \underset{\left|x\right|\to +\infty }{lim}{\left|x\right|}^{1/2}\left(\nabla u\cdot \frac{x}{\left|x\right|}-iku\right)=0.\hfill & \hfill \end{array}$

Assume by contradiction that there exists $p\in \left\{\mathrm{1,}\cdots \mathrm{,}M\right\}$ such that $u^p$ does not vanish identically. We first note that the function

$v={\displaystyle \underset{q=1}{\overset{M}{\sum }}}{u}^q$ (10)

is an outgoing solution of Helmholtz equation in ${ℝ}^2\backslash \overline{{\Omega }_{-}}$, which satisfies in addition the boundary condition $\Lambda v=0$ on $\partial {\Omega }_{-}={\displaystyle {\cup }_{q=1}^M}\partial {\Omega }_q$. By classical uniqueness results for Helmholtz equation, this implies that $v\equiv 0$ in ${ℝ}^2\backslash \overline{{\Omega }_{-}}$. In particular, we have

$v={\partial }_{n}v=0\text{on}\partial {\Omega }_p.$ (11)

Let us define then the function

$w\left(x\right)=\{\begin{array}{ll}{u}^p\left(x\right)\hfill & \text{for}x\in {ℝ}^2\backslash {\bar{\Omega }}_p\hfill \\ -{\displaystyle \underset{q=\mathrm{1,}q\ne p}{\overset{M}{\sum }}}{u}^q\left(x\right)\hfill & \text{for}x\in {\Omega }_p\mathrm{,}\hfill \end{array}$

Clearly, w is an outgoing solution of Helmholtz equation in ${ℝ}^2\backslash {\bar{\Omega }}_p$ and ${\Omega }_p$. Moreover, according to (10) and (11), w has continuous trace and normal derivative through $\partial {\Omega }_p$. Consequently, $w\equiv 0$ in ${ℝ}^2$, and thus $u^p\equiv 0$ in ${ℝ}^2\backslash {\bar{\Omega }}_p$, which provides the desired contradiction.

In order to prove the existence of the solution $\left(u^1\mathrm{,}\cdots \mathrm{,}{u}^p\right)$ of problems ${\left({\mathcal{E}}^p\right)}_{p=1,M}$, we show that these coupled scattering problems are equivalent to a system of M integral equations of Fredholm type. The announced existence follows then from the uniqueness result proved above. Let us seek the solution $u^p$ in the form of a single-layer potential

$u^p\left(x\right)={\displaystyle {\int }_{\partial {\Omega }_p}}G\left(x\mathrm{,}y\right){\psi }^p\left(y\right)\text{d}y\forall x\in {ℝ}^2\backslash {\Omega }_p\mathrm{,}$ (12)

where ${\psi }^p\in H^{-1/2}\left(\partial {\Omega }_p\right)$ and $G\left(x\mathrm{,}y\right)=\left(1/4i\right)H_0^1\left(k\left|x-y\right|\right)$ denotes the Helmholtz Green’s function in ${ℝ}^2$. Let us define for all $p,q=1,\cdots ,M$, the following integral operators:

$\begin{array}{l}{\mathcal{L}}^{pq}:H^{-1/2}\left(\partial {\Omega }_q\right)\to H^{1/2}\left(\partial {\Omega }_p\right)\\ \psi \mapsto {\mathcal{L}}^{pq}\psi \left(x_p\right)={\displaystyle {\int }_{\partial {\Omega }_q}}G\left(x_p\mathrm{,}{y}_q\right)\psi \left(y_q\right)\text{d}{y}_q\end{array}$

It is well-known from classical results of potential theory that for all $p=\mathrm{1,}\cdots \mathrm{,}M$, the operator ${\mathcal{L}}^{pp}\mathrm{<mo>\; :\; </mo>}{H}^{-1/2}\left(\partial {\Omega }_p\right)\to H^{1/2}\left(\partial {\Omega }_p\right)$ defines an isomorphism, while for $p\ne q$, the operators ${\mathcal{L}}^{pq}\mathrm{<mo>\; :\; </mo>}{H}^{-1/2}\left(\partial {\Omega }_q\right)\to H^{1/2}\left(\partial {\Omega }_q\right)$ are compact (since their kernels are analytic). With these notations, the single-layer potential (12) solves $\left({\mathcal{E}}^p\right)$ if and only if the unknown density ${\psi }^p\in H^{-1/2}\left(\partial {\Omega }_q\right)$ solves the integral equation

${\mathcal{L}}^{pp}{\psi }^p=f_p-{\displaystyle \underset{q=1,q\ne p}{\overset{M}{\sum }}}{\mathcal{L}}^{pq}{\psi }^q,$ (13)

where we have set $f_p=-u_{\mathrm{<mo>\; |\; </mo>}\partial {\Omega }_p}^{inc}$. The M coupled integral Equations (13), for $p=1,\cdots ,M$, can be rewritten in the abstract form

$\left(\mathcal{L}+\mathcal{K}\right)\Psi =F,$ (14)

provided we set $\Psi \mathrm{<mo>\; :\; </mo>}=\left({\psi }^1\mathrm{,}\cdots \mathrm{,}{\psi }^M\right)\in {\mathcal{H}}^{-1/2}\mathrm{<mo>\; :\; </mo>}=H^{-1/2}\left(\partial {\Omega }_1\right)\times \cdots \times H^{-1/2}\left(\partial {\Omega }_M\right)$ and

$\mathcal{L}=\left(\begin{array}{ccc}{\mathcal{L}}^{11}& \cdots & 0\\ ⋮& & ⋮\\ 0& \cdots & {\mathcal{L}}^{MM}\end{array}\right),\mathcal{K}=\left(\begin{array}{cccc}0& {\mathcal{L}}^{12}& \cdots & {\mathcal{L}}^{1M}\\ ⋮& & ⋮& \\ {\mathcal{L}}^{M1}& {\mathcal{L}}^{M2}& \cdots & 0\end{array}\right),F=\left(\begin{array}{c}{f}_1\\ ⋮\\ f_M\end{array}\right).$

If we introduce ${\mathcal{H}}^{1/2}:=H^{1/2}\left(\partial {\Omega }_1\right)\times \cdots \times H^{1/2}\left(\partial {\Omega }_M\right)$, then the operator $\mathcal{L}\mathrm{<mo>\; :\; </mo>}{\mathcal{H}}^{-1/2}\to {\mathcal{H}}^{1/2}$ defines an isomorphism while $\mathcal{K}\mathrm{<mo>\; :\; </mo>}{\mathcal{H}}^{-1/2}\to {\mathcal{H}}^{1/2}$ is a compact operator. Therefore, Equation (14) is of Fredholm type and the proof is now complete.

5. Difficulties

Therefore, the design of efficient and robust numerical resolution methods for multiple diffraction by a large number of spheres and for a wide band of frequencies has a large importance. Concerning these numerical methods, the literature existing one essentially deals with low frequencies. The low or higher frequencies problem has been much less studied because it leads to very specific numerical difficulties, related on the one hand to the effects of strong coupling between the obstacles and, on the other hand, to the undefined nature of the matrix linear system complex. In general, the numerical methods classically used in low frequency become ineffective at higher frequency. It is extremely important but also difficult to understand the properties of light scattered by many particles, in the complete spectrum, from low- to high-frequencies. We propose here to study the case of multiple diffraction by circular cylinders and we will propose solutions allowing to deal with both a high number of obstacles and a wide frequency band. Problems arising: infinite domain: ( ${ℝ}^2\backslash \bar{\Omega }$ ), multiple scattering, high-frequency (wavelength $\lambda =2\pi /k$ small compared to the characteristic size of ${\Omega }_p$ ). For a recent and comprehensive overview of these issues with extensive bibliographic references, we refer the reader to the manual [3] by P.A. Martin.

6. Numerical Simulation

6.1. Spectral Fourier Basis

We consider the case where the scatter are spheres in the two dimensional domain, see Figure 3. One can compute the boundary integral equations in a Fourier basis, leading therefore to an efficient computational spectral method when used direct or iterative solvers. Let us consider an orthonormal system $\left(O\mathrm{,}\stackrel{\leftarrow }{O{x}_1}\mathrm{,}\stackrel{\leftarrow }{O{x}_2}\right)$. We assume that the scattering obstacle ${\Omega }^{-}$ is the union of M disks ${\Omega }_p^{-}$, for $p=1,\cdots ,M$, of radius $a_p$ and center $O_p$. We define ${\Gamma }_p$ as the boundary of ${\Omega }_p^{-}$ and by $\Gamma ={\displaystyle {\cup }_{p=1,\cdots ,M}{\Gamma }_p}$ the boundary of ${\Omega }^{-}$. The unit normal vector $n$ to ${\Omega }^{-}$ is outgoing.

For any $p=1,\cdots ,M$, we introduce $b_p$ as the vector between the center $O_p$ and the origin $O$

$b_p=O{O}_p\mathrm{,}{b}_p=\Vert b_p\Vert \mathrm{,}{\alpha }_p=Angle\left(\stackrel{\leftarrow }{O{x}_1}\mathrm{,}{b}_p\right)\mathrm{,}$

Figure 3. Sketch two typical cylinders.

and, for $q=\mathrm{1,}\cdots \mathrm{,}M$, with $q\ne p$, $b_{pq}$ as the vector between the centers $O_q$ and $O_p$

$b_{pq}=O_q{O}_p\mathrm{,}{b}_{pq}=\Vert b_{pq}\Vert ,{\alpha }_{pq}\mathrm{<mo>\; =\; </mo>}Angle\left(\stackrel{\leftarrow }{O{x}_1}\mathrm{,}{b}_{pq}\right)\mathrm{.}$

Furthermore, any point $x$ is described by its global polar coordinates

$r\left(x\right)=Ox\mathrm{,}r\left(x\right)=\Vert r\left(x\right)\Vert \mathrm{,}\theta \left(x\right)=Angle\left(\stackrel{\leftarrow }{O{x}_1}\mathrm{,}r\left(x\right)\right)\mathrm{,}$

or by its polar coordinates in the orthonormal system associated with the obstacle ${\Omega }_p^{-}$, with $p=1,\cdots ,M$,

$r_p\left(x\right)=O_{p}x\mathrm{,}{r}_p\left(x\right)=\Vert r_p\left(x\right)\Vert \mathrm{,}{\theta }_p\left(x\right)=Angle\left(\stackrel{\leftarrow }{O_p{x}_1}\mathrm{,}{r}_p\left(x\right)\right)\mathrm{.}$

Let us now build a basis of $L^2\left(\Gamma \right)$ to approximate the integral operators. To this end, we first construct a basis of $L^2\left({\Gamma }_p\right)$ associated with ${\Omega }_p^{-}$, for $p=1,\cdots ,M$. If the circle ${\Gamma }_p$ has a radius one and is centered at the origin, then a suitable basis of $L^2\left({\Gamma }_p\right)$ is the spectral Fourier basis of functions ${\left({\text{e}}^{im\theta }\right)}_{m\in \mathbb{Z}}$. We adapt this basis to the general case where $a_p\ne 1$ by introducing, on one hand, the functions ${\left({\phi }_m\right)}_{m\in \mathbb{Z}}$ defined on ${ℝ}^2$ by: $\forall m\in \mathbb{Z}$, $\forall x\in {ℝ}^2$, ${\phi }_m\left(x\right)={\text{e}}^{im\theta \left(x\right)}$, and, on the other hand, the functions ${\left({\phi }_m^p\right)}_{1\le p\le M\mathrm{,}m\in \mathbb{Z}}$ given by

$\forall p=\mathrm{1,}\cdots \mathrm{,}M\mathrm{,}\forall m\in \mathbb{Z}\mathrm{,}\forall x\in {\Gamma }_p\mathrm{,}{\phi }_m^p\left(x\right)=\frac{{\phi }_m\left(r_p\left(x\right)\right)}{\sqrt{2\pi a_p}}=\frac{{\text{e}}^{im{\theta }_p\left(x\right)}}{\sqrt{2\pi a_p}}\mathrm{.}$

For $p=1,\cdots ,M$, the family ${\left({\phi }_m^p\right)}_{m\in \mathbb{Z}}$ forms an orthonormal basis of $L^2\left({\Gamma }_p\right)$ for the Hermitian inner product ${\left(\cdot \mathrm{,}\cdot \right)}_{L^2(\; \Gamma \; p\; )}$

$\forall f\mathrm{,}g\in L^2\left({\Gamma }_p\right)\mathrm{,}{\left(f\mathrm{,}g\right)}_{L^2\left({\Gamma }_p\right)}={\displaystyle {\int }_{{\Gamma }_p}}f\left(x\right)\overline{g\left(x\right)}\text{d}{\Gamma }_p\left(x\right)\mathrm{.}$

To build a basis of $L^2\left(\Gamma \right)$, we introduce the functions ${\Phi }_m^p$ of $L^2\left(\Gamma \right)$ as the union of these M families

$\forall p\mathrm{,}q=\mathrm{1,}\cdots \mathrm{,}M\mathrm{,}\forall m\in \mathbb{Z}\mathrm{,}{{\Phi }_m^p|}_{{\Gamma }_q}=(\begin{array}{ll}0\hfill & \text{if}q\ne p\mathrm{,}\hfill \\ {\phi }_m^p\hfill & \text{if}q=p\mathrm{.}\hfill \end{array}$

The family $\Omega =\left\{{\Phi }_m^p\mathrm{,}m\in \mathbb{Z}\mathrm{,}p=\mathrm{1,}\cdots \mathrm{,}M\right\}$, also called Fourier or spectral basis, is a Hilbert basis of $L^2\left(\Gamma \right)$ for the usual scalar product ${\left(\cdot \mathrm{,}\cdot \right)}_{L^2\left(\Gamma \right)}$. This family is used to approximate all integral equation used to solve the exterior scattering problem.

6.2. System of Algebraic Equations

We consider here the single-layer representation of the scattered field

$u=\mathcal{L}\rho ,$

(even if various integral equations can also be written for solving the Dirichlet problem (like the Combined Field Integral Equation (CFIE) or the Brakhage-Werner integral equation). One can prove that the surface density $\rho$ is equal to ${\left(-{\partial }_{n}u-{\partial }_n{u}^{\text{inc}}\right)|}_{\Gamma }$. A first-kind integral equation, which is usually called Electric Field Integral Equation (EFIE), is based on the trace of the single-layer operator

$L\rho =-{u^{\text{inc}}|}_{\Gamma }\mathrm{.}$ (15)

when ${\Omega }^{-}={\displaystyle {\cup }_{p=1}^M}{\Omega }_p$ is multiply connected, all the integral operators can be written by blocks. For example, the single-layer potential $\mathcal{L}\rho$ can be expressed as the sum of elementary potentials

$\mathcal{L}\rho ={\displaystyle \underset{p=1}{\overset{M}{\sum }}}{\mathcal{L}}_p{\rho }_p,$

where ${\rho }_p={\rho |}_{{\Gamma }_p}$ and

${\mathcal{L}}_p{\rho }_p\left(x\right)={\displaystyle {\int }_{{\Gamma }_p}}G\left(x\mathrm{,}y\right){\rho }_p\left(y\right)\text{d}x\mathrm{,}\forall x\in {ℝ}^2\backslash {\bar{\Omega }}_p\mathrm{.}$

Introducing some matrix operator notations, we also have the compact form of the EFIE for the acoustic multiple scattering problem

$L\rho =-U^{\text{inc}}\mathrm{,}$ (16)

where

$L\mathrm{<mo>\; :\; </mo>}=\left(\begin{array}{cccc}{L}_{\mathrm{1,1}}& L_{\mathrm{1,2}}& \cdots & L_{\mathrm{1,}M}\\ L_{\mathrm{2,1}}& L_{\mathrm{2,2}}& \cdots & L_{\mathrm{2,}M}\\ ⋮& ⋮& \ddots & ⋮\\ L_{M\mathrm{,1}}& L_{M\mathrm{,2}}& \cdots & L_{M\mathrm{,}M}\end{array}\right)\mathrm{,}\rho \mathrm{<mo>\; :\; </mo>}=\left(\begin{array}{c}{\rho }_1\\ {\rho }_2\\ ⋮\\ {\rho }_M\end{array}\right)\mathrm{,}{U}^{\text{inc}}\mathrm{<mo>\; :\; </mo>}=\left(\begin{array}{c}{u^{\text{inc}}|}_{{\Gamma }_1}\\ {u^{\text{inc}}|}_{{\Gamma }_2}\\ ⋮\\ {u^{\text{inc}}|}_{{\Gamma }_M}\end{array}\right)$

Remark 1.

● The integral formulation (16)provides an exact representation of the wavefields;

● Surface field $\rho$ +far-field pattern can be computed;

● However,computationally extensive calculation =dense linear system (16)+size of the system increases thanks to M when a discretization by a boundary element is applied;

● Huge memory requirement as well as a costly numerical solution for the linear system;

● Clearly,a numerical solution by Krylov subspace iterative solvers can be expected +Multilevel Fast Multipole Methods;

● Nevertheless,the numerical solution remains expensive to obtain since M can be very large.

6.3. Approximation

Let us for example consider the EFIE: find the density $\rho$ in $H^{-1/2}\left(\Gamma \right)$ ( $L^2\left(\Gamma \right)$ by smoothness) such that

$\forall x\in \Gamma \mathrm{,}L\rho \left(x\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_{\Gamma }}G\left(x\mathrm{,}y\right)\rho \left(y\right)\text{d}y=-u^{\text{inc}}(\; x\; )$

Proposition 7. Weak formulation of the EFIE in $L^2(\; \Gamma \; )$

${\langle L\rho \mathrm{,}\phi \rangle }_{L^2\left(\Gamma \right)}={\langle -{u^{\text{inc}}|}_{\Gamma }\mathrm{,}\phi \rangle }_{L^2\left(\Gamma \right)}\mathrm{,}\forall \phi \in L^2\left(\Gamma \right)\mathrm{,}$

where ${\langle f\mathrm{,}g\rangle }_{L^2\left(\Gamma \right)}={\displaystyle {\int }_{\Gamma }}f\bar{g}\text{d}\Gamma$.

6.3.1. Case of One Disk

Considering the case with only one spherical scatter, see Figure 4. First of all, we need the following proposition that will be used in our calculus.

Proposition 8. Case of one disk ${\Omega }_1$ : The functions ${\left({\text{e}}^{im{\theta }_1}\right)}_{m\in \mathbb{Z}}$ form an orthogonal basis of $L^2\left(\Gamma \right)=L^2\left({\Gamma }_1\right)$ .

Decomposition in the Fourier Basis:

$\text{}\text{}\text{}\rho ={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{\rho }_m^1{\text{e}}^{im{\theta }_1},\text{}\text{}\text{}-{u^{\text{inc}}|}_{{\Gamma }_1}={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{d}_m^1{\text{e}}^{im{\theta }_1}$

● the coefficients ${\left({\rho }_m^1\right)}_{m\in \mathbb{Z}}$ are unknown;

● the ${\left(d_m^1\right)}_{m\in \mathbb{Z}}$ are known.

Weak formulation in the Fourier basis:

Figure 4. Case of one disk.

${\displaystyle \underset{n\in \mathbb{Z}}{\sum }}{\mathbb{L}}_{m\mathrm{,}n}^{\mathrm{1,1}}{\rho }_n^1=d_m^1\mathrm{,}\forall m\in \mathbb{Z}\mathrm{,}$

where ${\mathbb{L}}_{m\mathrm{,}n}^{\mathrm{1,1}}={\langle L{\text{e}}^{in{\theta }_1}\mathrm{,}{\text{e}}^{im{\theta }_1}\rangle }_{L^2\left(\Gamma \right)}$.

Theorem 2 (Kress, 1985). The functions ${\text{e}}^{im{\theta }_1}$ are eigenvectors of the operator L, and we have, for all $m\mathrm{,}n\in \mathbb{Z}$ :

${\mathbb{L}}_{m\mathrm{,}n}^{\mathrm{1,1}}=(\begin{array}{ll}i{\pi }^2{a}_1^2{J}_m\left(k{a}_1\right)H_m^{\left(1\right)}\left(k{a}_1\right)\hfill & \text{if}m=n\mathrm{,}\hfill \\ 0\hfill & \text{if}m\ne n\mathrm{,}\hfill \end{array}$

$J_m$ and $H_m^{\left(1\right)}$ :Bessel and Hankel functions of first kind of order m.

The matrix of the resulting system is diagonal. Now, What about 2 disks?

6.3.2. Case of Two Disks

Considering the case with two cylindrical scatters, see Figure 5. Use one Fourier basis per scatterer.

${\Omega }_1\mathrm{<mo>\; :\; </mo>}(\begin{array}{l}{{\Phi }_m^1|}_{{\Gamma }_1}={\text{e}}^{im{\theta }_1}\\ {{\Phi }_m^1|}_{{\Gamma }_2}=0\end{array}$

$-{u^{\text{inc}}|}_{{\Gamma }_1}={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{d}_m^1{\text{e}}^{im{\theta }_1}$

${\Omega }_2\mathrm{<mo>\; :\; </mo>}(\begin{array}{l}{{\Phi }_m^2|}_{{\Gamma }_1}=0\\ {{\Phi }_m^2|}_{{\Gamma }_2}={\text{e}}^{im{\theta }_2}\end{array}$

$-{u^{\text{inc}}|}_{{\Gamma }_2}={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{d}_m^2{\text{e}}^{im{\theta }_2}$

$\rho ={\displaystyle \underset{n\in \mathbb{Z}}{\sum }}{\rho }_n^1{\Phi }_n^1+{\displaystyle \underset{n\in \mathbb{Z}}{\sum }}{\rho }_n^2{\Phi }_n^2\mathrm{.}$

The weak formulation of the EFIE is presented in the following theorem.

Proposition 9.

$\forall p=\mathrm{1,2,}\forall m\in \mathbb{Z}\mathrm{,}{\displaystyle \underset{q=1}{\overset{2}{\sum }}}{\displaystyle \underset{n\in \mathbb{Z}}{\sum }}{\mathbb{L}}_{m\mathrm{,}n}^{p\mathrm{,}q}{\rho }_n^q=d_m^p\mathrm{,}$

with ${\mathbb{L}}_{m\mathrm{,}n}^{p\mathrm{,}q}={\langle L{\text{e}}^{in{\theta }_q}\mathrm{,}{\text{e}}^{im{\theta }_p}\rangle }_{L^2(\; \Gamma \; p\; )}$

Remark 2. Calculation of the coefficients ${\mathbb{L}}_{m\mathrm{,}n}^{p\mathrm{,}q}$ .

● if $q=p$ ,self-interaction of the $p^{th}$ scatterer with itself,already computed;

● if $q\ne p$ ,interaction from ${\Gamma }_q$ to ${\Gamma }_p$ ,then we use the Graf’s addition theorem to compute analytically the coefficients

Figure 5. Two cylindrical scatters.

${\mathbb{L}}_{m\mathrm{,}n}^{p\mathrm{,}q}=\frac{i}{4}{\displaystyle {\int }_{{\Gamma }_p}}{\displaystyle {\int }_{{\Gamma }_q}}{H}_0^{\left(1\right)}\left(k\left|x-y\right|\right){\text{e}}^{in{\theta }_q\left(y\right)}{\text{e}}^{-im{\theta }_p\left(x\right)}\text{d}{\Gamma }_q\left(y\right)\text{d}{\Gamma }_p\left(x\right)\mathrm{.}$

Theorem 3 (Two-center expansion). Let $x\in {\Gamma }_p$ , $y\in {\Gamma }_q$ , $b_{pq}=\stackrel{\leftarrow }{O_q{O}_p}$ , ${\alpha }_{pq}=Angle\left(\stackrel{\leftarrow }{Ox}\mathrm{,}{b}_{pq}\right)$ , and $b_{pq}=\left|b_{pq}\right|$ ,

$H_0^{\left(1\right)}\left(k\left|x-y\right|\right)={\displaystyle \underset{m\in \mathbb{Z}}{\sum }}{\displaystyle \underset{n\in \mathbb{Z}}{\sum }}\text{}{J}_m\left(k{a}_p\right){\text{e}}^{im{\theta }_p\left(x\right)}{S}_{nm}\left(k{b}_{pq}\right)J_n\left(k{a}_q\right){\text{e}}^{-in{\theta }_q(\; y\; )}$

with $S_{nm}\left(k{b}_{pq}\right)=H_{n-m}^{\left(1\right)}\left(k{b}_{pq}\right){\text{e}}^{i\left(n-m\right){\alpha }_{pq}}$.

Using this and the orthogonality of the functions ( ${\text{e}}^{im{\theta }_p}$ ):

${\mathbb{L}}_{m\mathrm{,}n}^{p\mathrm{,}q}=i{\pi }^2{a}_p{a}_q{J}_m\left(k{a}_p\right)S_{nm}\left(k{b}_{pq}\right)J_n\left(k{a}_q\right)\mathrm{.}$

$\left[\begin{array}{cc}{\mathbb{J}}^1{ℍ}^1& {\mathbb{J}}^1{\left({\mathbb{S}}^{\mathrm{1,2}}\right)}^{\text{T}}{\mathbb{J}}^2\\ {\mathbb{J}}^2{\left({\mathbb{S}}^{\mathrm{2,1}}\right)}^{\text{T}}{\mathbb{J}}^1& {\mathbb{J}}^2{ℍ}^2\end{array}\right]\left[\begin{array}{c}{\rho }^1\\ {\rho }^2\end{array}\right]=\left[\begin{array}{c}{B}^1\\ B^2\end{array}\right]$

● ${\rho }^p={\left({\rho }_m^p\right)}_{m\in \mathbb{Z}}$ vector of the unknown Fourier coefficients;

● $B^p={\left(d_m^p\right)}_{m\in \mathbb{Z}}$ vector of the Fourier coefficients for the incident field;

● Self-interaction (diagonal blocks): ${\mathbb{J}}^p=\text{diag}\left(J_m\left(k{a}_p\right)\right)$ and ${ℍ}^p=\text{diag}\left(H_m^{\left(1\right)}\left(k{a}_p\right)\right)$ diagonal operators which represent the interaction of the scatterer p with itself;

● Interaction between ${\Gamma }_p$ and ${\Gamma }_q$ (off-diagonal blocks): full matrices ${\left({\left({\mathbb{S}}^{p\mathrm{,}q}\right)}^{\text{T}}\right)}_{mn}=S_{nm}\left(b_{pq}\right)=H_{n-m}^{\left(1\right)}\left(k{b}_{pq}\right){\text{e}}^{i\left(n-m\right){\alpha }_{pq}}$.

6.3.3. General Configuration with M Scatterers

Using the same approach we get the following linear system of algebraic equations.

$\mathbb{L}\rho =B$

setting

$\mathbb{L}=\left[\begin{array}{cccc}{\mathbb{J}}^1{ℍ}^1& {\mathbb{J}}^1{\left({\mathbb{S}}^{\mathrm{1,2}}\right)}^{\text{T}}{\mathbb{J}}^2& \cdots & {\mathbb{J}}^1{\left({\mathbb{S}}^{\mathrm{1,}M}\right)}^{\text{T}}{\mathbb{J}}^M\\ {\mathbb{J}}^2{\left({\mathbb{S}}^{\mathrm{2,1}}\right)}^{\text{T}}{\mathbb{J}}^1& {\mathbb{J}}^2{ℍ}^2& \cdots & {\mathbb{J}}^2{\left({\mathbb{S}}^{\mathrm{2,}M}\right)}^{\text{T}}{\mathbb{J}}^M\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbb{J}}^M{\left({\mathbb{S}}^{M\mathrm{,1}}\right)}^{\text{T}}{\mathbb{J}}^1& {\mathbb{J}}^M{\left({\mathbb{S}}^{M\mathrm{,2}}\right)}^{\text{T}}& \cdots & {\mathbb{J}}^M{ℍ}^M\end{array}\right],\rho =\left[\begin{array}{c}{\rho }^1\\ ⋮\\ {\rho }^M\end{array}\right],B=\left[\begin{array}{c}{B}^1\\ ⋮\\ B^M\end{array}\right]$

6.4. Examples

Radar cross-section (RCS) is a measure of how detectable an object is with a radar. A larger RCS indicates that an object is more easily detected. Radar cross-section is used to detect planes in a wide variation of ranges. In the polar coordinates system $\left(r\mathrm{,}\theta \right)$ and by using an asymptotic expansion of u when $r\to +\infty$, the following relation holds [4]

$\forall \theta \in \left[\mathrm{0,2}\pi \right]\mathrm{,}u\left(r\mathrm{,}\theta \right)=\frac{{\text{e}}^{ikr}}{r^{1/2}}\left[a_{\mathcal{L}}\left(\theta \right)+a_{\mathcal{M}}\left(\theta \right)\right]+O\left(\frac{1}{r^{3/2}}\right)\mathrm{,}$

where $a_{\mathcal{L}}$ and $a_{\mathcal{M}}$ are the radiated far-fields for the single- and double-layer potentials, respectively, defined for any angle $\theta$ of $\left[\mathrm{0,2}\pi \right]$ by

(a)(b)

Figure 6. Convergence study and Radar cross section. Multiple scattering of an incident plane wave ( $k=10$ and $\beta =0$ (rad.) (coming from the right)) by $M=10$ sound-soft obstacles, randomly distributed in ${\left[-10,10\right]}^2$. (a) History of convergence of the GMRES without restart; (b) Radar cross section.

$(\begin{array}{l}{a}_L\left(\theta \right)=\frac{1}{\sqrt{8k\pi }}{\text{e}}^{i\pi /4}{\displaystyle {\int }_{\Gamma }}{\text{e}}^{-ik\theta \cdot y}\rho \left(y\right)\text{d}\Gamma \left(y\right)\mathrm{,}\\ a_M\left(\theta \right)=\frac{1}{\sqrt{8k\pi }}{\text{e}}^{i\pi /4}{\displaystyle {\int }_{\Gamma }}-\frac{ik}{\Vert y\Vert }\theta \cdot y{\text{e}}^{-ik\theta \cdot y}\lambda \left(y\right)\text{d}\Gamma \left(y\right)\mathrm{,}\end{array}$

with $\theta \mathrm{<mo>\; :\; </mo>}=\left(\cos \left(\theta \right)\mathrm{,}\sin \left(\theta \right)\right)$. In addition, the Radar Cross Section (RCS) is defined by

$\forall \theta \in \left[\mathrm{0,2}\pi \right]\mathrm{,}\text{RCS}\left(\theta \right)=10{\log }_{10}\left(2\pi {\left|a_{\mathcal{L}}\left(\theta \right)+a_{\mathcal{M}}\left(\theta \right)\right|}^2\right)\left(\text{dB}\right)\mathrm{.}$

Figure 7. Absolute part scattred and total fields. Multiple scattering of an incident plane wave ( $k=10$ and $\beta =0$ (rad.) (coming from the right)) by $M=10$ sound-soft obstacles, randomly distributed in ${\left[-10,10\right]}^2$. (a) Absolute part scattred field; (b) Absolute part total field.

Figure 8. Properties of total field. Multiple scattering of an incident plane wave ( $k=10$ and $\beta =0$ (rad.) (coming from the right)) by $M=10$ sound-soft obstacles, randomly distributed in ${\left[-10,10\right]}^2$. (a) Real part of the total field; (b) Imaginary part the total field; (c) Absolute value of the total field.

The coefficients ${\rho }_m^p$ allow us to calculate any physical interesting quantity as e.g. the Radar Cross Section and other quantities.

In our numerical simulation, we consider the scattering of an incident plane wave with $k=10$ and $\beta =0$ (rad.) coming from the right by $M=10$ sound-soft obstacles, randomly distributed in the square domain ${\left[-10,10\right]}^2$. In Figure 6, we present different field (scatted and total) to characterize the scattering problem. In Figure 7, Absolute part scattred and total field is depicted. In Figure 8, we have reported the history of convergence of the GMRES with $k=10$ without restart and plotted the behavior the RCS. Numerical results suggest that an accurate approximation of the RCS can be obtained with two obstacles, leading to a relative error of the order 10⁻⁴. All these result are obtained using Matlab software.

7. Conclusion

The single scattering waves was computed using Mie series expansion and EFIE technique. This is the first step to understand the problem and open a new huge way for studying multiple scattering by different shapes. Fourier series expansion can be used for multiple scattering and then a linear system must be solved. Spectral technique can be used also in the case of Integral representation of the solution. This reliable numerical solution would be able to provide a fast (iterative) algorithm for multiple scattering but could also lead to efficient preconditioners for the full solution by integral equations.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Antoine, X., Chniti, C. and Ramdani, K. (2008) On the Numerical Approximation of High-Frequency Acoustic Multiple Scattering Problems by Circular Cylinders. Journal of Computational Physics, 227, 1754-1771. https://doi.org/10.1016/j.jcp.2007.09.030
[2]	Antoine, X., Geuzaine, C. and Ramdani, K. (2010) Computational Methods for Multiple Scattering at High Frequency with Applications to Periodic Structures Calculations. In: Ehrhardt, M., Ed., Wave Propagation in Periodic Media—Analysis, Numerical Techniques and Practical Applications, Vol. 1, Bentham Science Publishers Ltd., UAE, 73-107. https://doi.org/10.2174/978160805150211001010073
[3]	Martin, P.A. (2006) Multiple Scattering. Interaction of Time-Harmonic Waves with N Obstacles, Vol. 107 of Encyclopedia of Mathematics and Its Applications. Cambridge University Press, Cambridge.
[4]	Colton, D.L. and Kress, R. (1983) Integral Equation Methods in Scattering Theory. Pure and Applied Mathematics. John Wiley & Sons Inc., New York.
[5]	Nédélec, J.-C. (2001) Acoustic and Electromagnetic Equations. Integral Representations for Harmonic Problems, Vol. 144 of Applied Mathematical Sciences. Springer-Verlag, New York.
[6]	Acosta, S. (2015) On-Surface Radiation Condition for Multiple Scattering of Waves. Computer Methods in Applied Mechanics and Engineering, 283, 1296-1309. https://doi.org/10.1016/j.cma.2014.08.022
[7]	Acosta, S. and Villamizar, V. (2010) Coupling of Dirichlet-to-Neumann Boundary Condition and Finite Difference Methods in Curvilinear Coordinates for Multiple Scattering. Journal of Computational Physics, 229, 5498-5517. https://doi.org/10.1016/j.jcp.2010.04.011
[8]	Antoine, X., Ramdani, K. and Thierry, B. (2012) Wide Frequency Band Numerical Approaches for Multiple Scattering Problems by Disks. Journal of Algorithms & Computational Technology, 6, 241-259. https://doi.org/10.1260/1748-3018.6.2.241
[9]	Chen, J., Lee, Y., Lin, Y., Chen, I. and Lee, J. (2011) Scattering of Sound from Point Sources by Multiple Circular Cylinders Using Addition Theorem and Superposition Technique. Numerical Methods for Partial Differential Equations, 27, 1365-1383. https://doi.org/10.1002/num.20583
[10]	Ehrhardt, M. (2010) Wave Propagation in Periodic Media Analysis, Numerical Techniques and Practical Applications, E-Book Series Progress in Computational Physics (PiCP), Volume 1. Bentham Science Publishers, UAE.
[11]	Ehrhardt, M., Han, H. and Zheng, C. (2009) Numerical Simulation of Waves in Periodic Structures. Communications in Computational Physics, 5, 849-870.
[12]	Geuzaine, C., Bruno, O. and Reitich, F. (2005) On the O(1) Solution of Multiple-Scattering Problems. IEEE Transactions on Magnetics, 41, 1488-1491. https://doi.org/10.1109/TMAG.2005.844567
[13]	Grote, M. and Kirsch, C. (2004) Dirichlet-to-Neumann Boundary Conditions for Multiple Scattering Problems. Journal of Computational Physics, 201, 630-650. https://doi.org/10.1016/j.jcp.2004.06.012
[14]	Kharlamov, A. and Filip, P. (2012) Generalisation of the Method of Images for the Calculation of Inviscid Potential Flow past Several Arbitrarily Moving Parallel Circular Cylinders. Journal of Engineering Mathematics, 77, 77-85. https://doi.org/10.1007/s10665-012-9532-6
[15]	Van Genechten, B., Bergen, B., Vandepitte, D. and Desmet, W. (2010) A Trefftz-Based Numerical Modelling Framework for Helmholtz Problems with Complex Multiple-Scatterer Configurations. Journal of Computational Physics, 229, 6623-6643. https://doi.org/10.1016/j.jcp.2010.05.016
[16]	Colton, D. and Kress, R. (1998) Inverse Acoustic and Electromagnetic Scattering Theory, Vol. 93 of Applied Mathematical Sciences. 2nd Edition, Springer-Verlag, Berlin. https://doi.org/10.1007/978-3-662-03537-5
[17]	Balabane, M. (2004) Boundary Decomposition for Helmholtz and Maxwell Equations. I. Disjoint Subscatterers. Asymptotic Analysis, 38, 1-10.