On the constant constitutive parameter (e.g., mass density) assumption in integral equation approaches to (acoustic) wave scattering

TL;DR

This paper evaluates the common assumption of constant mass density in acoustic wave scattering problems, demonstrating that a series expansion can correct errors from this assumption for various shapes, frequencies, and parameters.

Contribution

It introduces a series-based correction method to improve the accuracy of constant density assumptions in integral equation approaches for wave scattering.

Findings

Series expansion effectively corrects constant density approximation errors.

Two additional terms in the series suffice for high accuracy across parameters.

Method applies to non-canonical obstacle shapes.

Abstract

In 2D acoustic and elastodynamic problems the spatial variability of a constitutive parameter such as the mass density makes it difficult to employ boundary integral and domain integral techniques to solve the forward and inverse wave scattering problems. The oft-employed method for avoiding this problem is to assume this constitutive parameter (which is chosen herein to be the mass density) to be spatially-invariant throughout all space. The reliability of this assumption is evaluated both theoretically and numerically and it is shown, in the example of a canonical-shaped scattering obstacle, that the scattered field can be obtained in the form of a series of powers of the mass density contrast (the latter vanishing for constant mass density). The first term of this series is the solution for the scattered field corresponding to the constant density assumption and it is shown that taking into account only two more terms in the series enables to correct for practically all the errors incurred by the constant mass density assumption for a wide range of the other constitutive parameters and frequencies. It is shown how to apply this result for obstacles of non-canonical shape.

Equations

$${\sigma }_{kl,k}+\rho (f_l-{\dot{v}}_l)=0;l,k=1,2,3,$$

$${\sigma }_{kl,k}+\rho (f_l-u_{l,tt})=0;l,k=1,2,3,$$

$${\sigma }_{kl}=c_{klmn}{\epsilon }_{mn};k,l,m,n=1,2,3,$$

$${\epsilon }_{kl}=\frac{1}{2}(u_{k,l}+u_{l,k});k,l=1,2,3,$$

$${\sigma }_{kl}=\lambda {\delta }_{kl}{\epsilon }_{mm}+2\mu {\epsilon }_{kl}={\sigma }_{lk};k,l,m=1,2,3,$$

$$\begin{gathered} {\sigma }_{kl,k}=\lambda {\delta }_{kl}{\epsilon }_{mm,k}+{\lambda }_{,k}{\delta }_{kl}{\epsilon }_{mm}+2\mu {\epsilon }_{kl,k}+2{\mu }_{,k}{\epsilon }_{kl}= \\ \lambda {\epsilon }_{mm,l}+{\lambda }_{,l}{\epsilon }_{mm}+2\mu {\epsilon }_{kl,k}+2{\mu }_{,k}{\epsilon }_{kl};k,l,m=1,2,3. \end{gathered}$$

$$\begin{gathered} {\sigma }_{kl,k}=\lambda u_{m,ml}+{\lambda }_{,l}{u}_{m,m}+\mu (u_{k,lk}+u_{l,kk})+{\mu }_{,k}(u_{k,l}+u_{l,k})= \\ \lambda u_{m,lm}+{\lambda }_{,l}{u}_{m,m}+\mu (u_{k,lk}+u_{l,kk})+{\mu }_{,k}(u_{k,l}+u_{l,k})= \\ \lambda u_{k,lk}+{\lambda }_{,l}{u}_{m,m}+\mu (u_{k,lk}+u_{l,kk})+{\mu }_{,k}(u_{k,l}+u_{l,k});k,l,m=1,2,3, \end{gathered}$$

$${\sigma }_{kl,k}=(\lambda +\mu )u_{k,kl}+\mu u_{l,kk}+{\lambda }_{,l}{u}_{m,m}+{\mu }_{,k}(u_{k,l}+u_{l,k});k,l,m=1,2,3,$$

$$(\lambda +\mu )u_{k,kl}+\mu u_{l,kk}+{\lambda }_{,l}{u}_{m,m}+{\mu }_{,k}(u_{k,l}+u_{l,k})+\rho (f_l-{\dot{v}}_l)=0;k,l,m=1,2,3,$$

$$(\lambda +\mu )u_{k,kl}+\mu u_{l,kk}+{\lambda }_{,l}{u}_{m,m}+{\mu }_{,k}(u_{k,l}+u_{l,k})+\rho (f_l-u_{l,tt})=0;k,l,m=1,2,3.$$

$${\sigma }_{kl}=\lambda {\delta }_{kl}{\epsilon }_{mm}=\lambda u_{m,m}{\delta }_{kl}{\sigma }_{lk};k,l,m=1,2,3.$$

$$p:=-\lambda {\epsilon }_{mm}=-\lambda u_{m,m};m=1,2,3,$$

$${\sigma }_{kl}=-p{\delta }_{kl}={\sigma }_{lk};k,l=1,2,3.$$

$${\sigma }_{kl,k}=-p_{,k}{\delta }_{kl}=-p_{,l};k,l=1,2,3.$$

$$-p_{,l}+\rho (f_l-{\dot{v}}_l)=0;k,l,m=1,2,3,$$

$$-p_{,l}+\rho (f_l-u_{l,tt})=0;k,l,m=1,2,3.$$

$$p_{,l}=-\lambda u_{m,ml}-{\lambda }_{,l}{u}_{m,m};l,m=1,2,3,$$

$$\lambda u_{m,ml}+{\lambda }_{,l}{u}_{m,m}+\rho (f_l-u_{l,tt})=0;l,m=1,2,3.$$

$$-p_{,ll}+{\rho }_{,l}(f_l-u_{l,tt})+\rho (f_{l,l}-u_{l,ttl})=0;l=1,2,3.$$

$$p_{,tt}:=-\lambda u_{m,mtt}\Rightarrow u_{l,ltt}=u_{l,ttl}=-\frac{p_{,tt}}{\lambda }.$$

$$u_{l,tt}=-\frac{p_l}{\rho }+f_l;l=1,2,3,$$

$$-p_{,ll}+{\rho }_{,l}(f_l+\frac{p_l}{\rho }-f_l)+\rho (f_{l,l}+\frac{p_{,tt}}{\lambda })=0;l=1,2,3.$$

$$-p_{,ll}+\frac{\rho }{\lambda }{p}_{,tt}+\frac{{\rho }_{,l}}{\rho }{p}_{,l}+\rho f_{l,l}=0;l=1,2,3.$$

$$\frac{1}{c^2}:=\frac{\rho }{\lambda }=\frac{\rho }{K}=\rho \kappa ,$$

$$-p_{,ll}+\frac{1}{c^2}{p}_{,tt}+\frac{{\rho }_{,l}}{\rho }{p}_{,l}+\rho f_{l,l}=0;l=1,2,3.$$

$$-\nabla \cdot \nabla p(\mathbf{x},t)+\frac{1}{c(\mathbf{x}{)}^2}\frac{{\partial }^{2}p(\mathbf{x},t)}{\partial t^2}+\frac{\nabla \rho (\mathbf{x})}{\rho (\mathbf{x})}\cdot \nabla p(\mathbf{x},t)+\rho (\mathbf{x})\nabla \cdot \mathbf{f}(\mathbf{x},t)=0;\forall \mathbf{x}\in {\mathbb{R}}^3.$$

$$c=c(\mathbf{x},\omega )=c^{\prime }(\mathbf{x},\omega )+i{c}^{\prime \prime }(\mathbf{x},\omega ),\rho (\mathbf{x},\omega )={\rho }^{\prime }(\mathbf{x},\omega )+i{\rho }^{\prime \prime }(\mathbf{x},\omega ),$$

$$p(\mathbf{x},t)={\int }_{-\infty }^{\infty }p(\mathbf{x},\omega )\exp (-i\omega t)d\omega ,$$

$$\mathbf{f}(\mathbf{x},t)={\int }_{-\infty }^{\infty }\mathbf{f}(\mathbf{x},\omega )\exp (-i\omega t)d\omega ,$$

$$\nabla \cdot \nabla p(\mathbf{x},\omega )+k^2(\mathbf{x},\omega )p(\mathbf{x},\omega )-\frac{\nabla \rho (\mathbf{x},\omega )}{\rho (\mathbf{x},\omega )}\cdot \nabla p(\mathbf{x},\omega )=\rho (\mathbf{x},\omega )\nabla \cdot \mathbf{f}(\mathbf{x},\omega );\forall \mathbf{x}\in {\mathbb{R}}^3,$$

Full text

Keywords: 2D acoustics, 2D elastodynamics, scattering, domain integral equations, boundary integral equations, constant density

Abbreviated title: Constant mass density assumption in integral equation methods of wave scattering

Corresponding author: Armand Wirgin,

e-mail: wirgin@lma.cnrs-mrs.fr

1. 1 Introduction

2. 2 Governing equations

3. 2.1 Conservation of momentum

4. 2.2 Conservation of momentum for small deformations

5. 2.3 Constitutive relations in a linear material continuum

6. 2.4 Constitutive relations in a linear, isotropic, elastic continuum

7. 2.5 Wave equation in a linear, isotropic, elastic medium (Navier’s equation)

8. 2.6 Governing equations in an inviscid fluid

9. 2.7 Time domain wave equation in a fluid

10. 2.8 Frequency domain wave equation in an inviscid or viscous fluid

11. 2.9 Getting rid of the gradient of pressure term via the $p-q$ transformation

12. 2.10 Specifics of 2D elastic and acoustic wave mechanics

13. 3 Acoustic radiation from applied sources in the absence of an obstacle in unbounded space occupied by a homogeneous, lossless fluid

14. 3.1 Description of the problem

15. 3.2 Governing equations

16. 3.3 Field representation

17. 3.4 The 2D free-space Green’s function

18. 3.4.1 Demonstration that the 2D GF satisfies the 2D Helmholtz equation

19. 3.4.2 Demonstration that the 2D GF satisfies the radiation condition

20. 4 Scattering of the sound radiated from applied sources by a generally-heterogeneous fluid-like obstacle within an unbounded, homogeneous, background fluid medium

21. 4.1 Description of the problem

22. 4.2 The 2D scattering problem for a macroscopically-homogeneous, isotropic circular cylindrical obstacle by the domain decomposition separation of variables (DD-SOV) technique

23. 4.2.1 Governing equations via DD

24. 4.2.2 Field representations via SOV

25. 4.2.3 Determination of the unknown coefficients via the transmission conditions

26. 4.3 The coupled boundary integral (BI-BI) method applicable to the special case in which both the mass density of, and wavespeed in, the fluid-like obstacle do not depend on position

27. 4.3.1 BI-BI method: domain decomposition

28. 4.3.2 BI-BI method: iterative method for obtaining the solutions for $U$ and $V$

29. 4.3.3 BI-BI method: the 2D case (of cylindrical symmetry)

30. 4.3.4 BI-BI method: the case of a homogeneous, circular cylindrical obstacle

31. 4.4 The coupled domain-boundary (DI-BI) method applicable to the special case in which the mass density of the fluid-like obstacle does not depend on position

32. 4.4.1 Domain decomposition

33. 4.4.2 DI-BI method: iterative methods for obtaining the solutions for $U$ and $W$

34. 4.4.3 DI-BI method: the 2D (i.e., cylindrical symmetry) case in polar coordinates when $k^{[1]}(\mathbf{x})=k^{[1]}(r)$

35. 4.4.4 DI-BI method: the 2D (i.e., cylindrical symmetry) case in polar coordinates when $k^{[1]}(\mathbf{x})=k^{[1]}$

36. 4.5 The BI-DI method

37. 4.5.1 Governing equations

38. 4.5.2 BI-DI method: the integral equation and integral transform for a circular cylindrical obstacle

39. 4.5.3 BI-DI method: the integral equation and integral transform for a circular cylindrical obstacle in which the density and wavespeed are constants with respect to the space coordinates

40. 4.6 The Domain Integral (DI) method

41. 4.6.1 Governing equations

42. 4.6.2 DI field representation

43. 4.6.3 DI method: a two-step scheme deriving from the $p-q$ transformation

44. 4.6.4 DI method: an iterative scheme

45. 4.6.5 DI method: the 2D case of cylindrical symmetry

46. 4.6.6 The DI method: case in which the boundary of the ”obstacle” is a circle, $k^{[1]}=k^{[1]}(r)$ and $\rho^{[1]}$ is a constant as a function of $r$ and $\theta$

47. 4.6.7 The DI method: case in which the boundary of the ”obstacle” is a circle, and $k^{[1]}$ and $\rho^{[1]}$ are constants as a function of $r$ and $\theta$

48. 4.6.8 DI method: the constant density solution for the ”obstacle”

49. 4.6.9 DI method: an iterative solution for the ”obstacle” useful for small differences of density

50. 4.6.10 Reconciling various publications dealing with the DI method

51. 4.7 The reasons, in general 2D acoustic scattering problems, for choosing the configuration of two homogeneous, isotropic media separated by a boundary of canonical shape

52. 5 Numerical results

53. 5.1 The numerical procedure and choice of parameters

54. 5.2 Variation of $\rho^{[1]}$ for various orders $l$: case $\rho^{[1]}>\rho^{[0]}$ and low $f$

55. 5.3 Variation of $\rho^{[1]}$ for various orders $l$: case $\rho^{[1]}<\rho^{[0]}$ and low $f$

56. 5.4 Variation of $\rho^{[1]}$ for various orders $l$: case $\rho^{[1]}>\rho^{[0]}$ and higher $f$

57. 5.5 Variation of $\rho^{[1]}$ for various orders $l$: case $\rho^{[1]}<\rho^{[0]}$ and higher $f$

58. 5.6 Variation of (low-to-medium) frequency $f$ for various orders $l$: case of fairly-large ($\epsilon=0.231$) mass density contrast and large \big{|}\Im(c^{[1]})\big{|}

59. 5.7 Variation of frequency $f$ for various orders $l$: case of relatively-small ($\epsilon=0.091$) density contrast and large \big{|}\Im(c^{[1]})\big{|}

60. 5.8 Variation of (low-to-medium) frequency $f$ for various orders $l$: case of relatively-small ($\epsilon=0.091$) small density contrast and small \big{|}\Im(c^{[1]})\big{|}

61. 5.9 Variation of low-to-medium frequency $f$ for various orders $l$: case of fairly-large ($\varepsilon=0.231$) mass density contrast and large wavespeed contrast for $\Re(c^{[1]})>c^{[0]}$

62. 5.10 Variation of low-to-medium frequency $f$ for various orders $l$: case of fairly-large ($\varepsilon=0.231$) mass density contrast and small wavespeed contrast for $\Re(c^{[1]})>c^{[0]}$

63. 5.11 Variation of low-to-medium frequency $f$ for various orders $l$: case of fairly-large ($\varepsilon=0.231$) mass density contrast and small wavespeed velocity contrast for $\Re(c^{[1]})<c^{[0]}$

64. 5.12 Variation of high frequency $f$ for various orders $l$: case of fairly-large ($\varepsilon=0.231$) mass density contrast, small velocity contrast and large \big{|}\Im(c^{[1]})\big{|}

65. 5.13 Variation of frequency near the resonance frequency of the $l=20$ normal mode, for various $\Im(c^{[1]})$

66. 5.14 Variation of frequency $f$ near the former $l=20$ normal mode resonance frequency for various $\Re(c^{[1]})$

67. 5.15 Variation of $\Re(c^{[1]})$ for various $l$: case of $\epsilon=0.231$ density contrast at $f=20~{}KHz$

68. 5.16 Variation of $\Re(c^{[1]})$ for various $l$: case of $\epsilon=-0.429$ density contrast at $f=20~{}KHz$

69. 6 Conclusions

1 Introduction

When an oscillating or impulsive source is present in an unbounded, homogeneous, isotropic medium, it produces radiation (i.e., a radiated acoustic, elastic or electromagnetic wave). When the medium is locally-inhomogeneous (as for instance, due to the presence of a spatially-bounded or unbounded in one or two directions (e.g., half-space) region in which the constitutive properties are different from those of the previously-mentioned homogeneous medium), the wave produced by the source is said to be scattered. Interest in scattering is primarily of practical nature (see [11, 43, 99] for many references), but also of theoretical (mathematical) nature.

Scattering theory initially concerned forward (scattering) problems that deal with the prediction of the spatial and spectral distributions of the scattered wave, assuming known everything (boundaries and constitutive functions (in terms of space, time/frequency variables) of the source and scatterer (or scatterers which can be grouped into a single entity termed ’the obstacle’) as well as the partial differential equations (PDE’s) and boundary conditions (BC’s, the radiation condition at infinity being a sort of BC on the ’boundary’ at infinity) that describe the way the radiated field acts on the obstacle. At present, scattering theory concentrates mostly on inverse (scattering) problems [69] that are concerned with the recovery (from knowledge of the radiated and/or scattered fields) of something or everything (i.e., the boundaries and constitutive functions) about the source and obstacle, assuming once again that the PDE’s and BC’s are known.

Fundamentally, the PDE’s of the forward and inverse problems are the same [60], but in the forward case the unknown is a field (acoustic, elastic, electromagnetic) and usually linear in terms of this field, whereas in the inverse case the unknown(s) is(are) the coefficient(s) or driving term of the PDE (associated with the aforementioned source which we now qualify as ’applied’ because it is assumed to not depend on the scattered field) and are nonlinear in terms of these entities. Another difference between the two cases is that forward scattering problems require a single resolution of the PDE per (scattering) configuration and frequency, whereas inverse scattering problems usually require several (if not many) resolutions of the PDE per configuration and frequency [46]. For this reason, it has become increasingly-important to either modify (in the sense of simplification and enabling ease of resolution) the PDE and/or obtain the most appropriate method of resolution (notably to reduce the computational burden required for inversion).

The important point to remember is that both forward and inverse scattering problems require the resolution(s) of a PDE (or PDE’s) for a function of space and time (or space and frequency) that represents the wavefield of interest. This function can be scalar (representing pressure, as in acoustics) or vectorial (representing particle displacement, as in elastodynamics, or the electric or magnetic field, as in electrodynamics). In the vectorial case, the various field components are usually coupled so that the previously-mentioned PDE is, in fact, a set of coupled PDE’s for the scalar components of the vectorial field. In both cases, we refer to the set of PDE’s as ’the PDE’. Similarly, we shall refer to the multitude of boundary conditions of each problem as ’the boundary condition’.

The PDE can be solved in brute-force manner by some sort of discretization scheme (e.g., finite difference (FD; [60, 69, 92]), finite element (FE; [3, 36]), spectral element (SE; [45]), etc.). This is the approach in many of the so-called nonlinear full-wave inversion (FWI) schemes [4, 101] that are increasingly adopted [69] in the geophysics community. We shall not be further concerned with this purely-numerical approach other than to underline the fact that, at present, little or no efforts are made, nor required in this approach (assuming, as is often the case, that one has considerable computer power at one’s disposal) to simplify the governing PDE, notably by eliminating a term that one could estimate a priori to be small compared to the others. This, of course, was not the situation several decades ago, especially in the marine acoustics and geophysical (i.e., seismology in relation to the solid earth) contexts for large-scale scattering problems, when the tendency was rather to first reduce the set of coupled PDE’s to a single PDE for a scalar function (i.e., the so-called acoustic approximation) [87, 47, 62, 92, 12] and then perhaps assume a coefficient of the PDE to be either small or constant with respect to position (these coefficients were (and still are often) assumed to be independent of time when the time domain PDE was (is) solved, but this limitation was (is) artificially-overcome in the frequency domain PDE by making the velocity complex).

Of course, less numerically-oriented methods also exist to solve directly (i.e., without resorting to integral transformations) the PDE. The most widely-known is a high-frequency procedure [84] wherefrom the solution of the second-order scalar-wave PDE emerges iteratively via a series of simpler, first-order PDE’s, the first one of which is none other than the eikonal equation, enabling the establishment of the notions of wavefronts and rays. Kline [44] extends this method to the vectorial PDE of electrodynamics and shows that the successive solutions of the first-order equations combine to form an asymptotic (high-frequency) representation of the solution of the second-order equation. A less well-known procedure [85, 17] is low-frequency in nature, and consists in expanding the solution of the governing PDE in a series of powers of the wavenumber whereby the solution emerges via a set of simpler PDE’s, the first of which is the static PDE. This series of powers was first thought to be convergent in certain regions outside the obstacle, but it is more generally an asymptotic series. Last but not least, we mention the parabolic approximation procedure [15, 83, 101] which is particularly suited to problems in which the waves in the inhomogeneous medium travel preferentially in one (horizontally in the sea) direction. Using this fact, it is possible, using some approximations valid in the far-field zone, to reduce the hyberbolic (second-order) wave equation to a parabolic (first-order) equation which lends itself to the non-iterative so-called ’marching method’ (also related to the ’splitting method’ mentioned further on) for its solution. An interesting feature of this parabolic equation is that it resembles the Schrödinger equation of quantum mechanics and since the solution of the latter can be represented in path-integral form, it has been suggested that this feature might provide new insights on how to efficiently solve scalar wave scattering problems such as those which arise in acoustics [25, 27, 26, 76].

Perhaps the oldest method for solving the PDE is by separation of variables (SOV). This is usually preceded by domain decomposition (DD) (e.g., to separate the regions inside and outside the obstacle and/or various homogeneous subregions within the obstacle) in order to be confronted with a set of PDE’s having (spatially-) constant coefficients. SOV then enables (provided each decomposed domain of the obstacle is homogeneous and its boundary has certain shapes) the PDE(s) to be reduced to several ordinary differential equations (ODE), each of which can be solved (this term is inexact because the ODE’s actually only give rise to representations of the wavefield, and the actual solution is obtained by employing these representations in the boundary conditions) in terms of special functions (exponentials, ordinary and spherical Bessel functions, etc.). The final solution of each PDE is then a sum or integral of the products of these various functions [43]. This method (outlined herein in sect. 4.2 is very powerful (not just because computer libraries enable the computation of these special functions very efficiently), but is restricted to very special scattering problems (notably as concerns the shape of the boundary (planar, circular [91, 79, 99], spherical [41, 56, 43], etc.) of the scatterer). The reasons for paying attention to the DD-SOV method are that: (a) it provides solutions of so-called ’canonical scattering problems’ (with reference to the shape of the boundary) which are exact (to within numerical errors) and therefore useful to test the validity [36] of other methods of resolution of the forward scattering problem, (b) it leads to a functional (i.e., mathematical instead of purely-numerical) expression of the solution which can be further manipulated by perturbation or asymptotic schemes so as to furnish functionally-approximate approximations (of the solution), which can be very useful in the inverse problem context [86, 96] when one or several parameters of the obstacle and/or of the source is/are small [17, 80] or large [96], or when the obstacle has non-canonical shape but its response to the radiated wave can be approximated locally by that of obstacles with canonical shape (i.e., the ICBA [53, 78, 97, 98, 11, 28]), and (c) the canonical solutions can be incorporated in homogenization schemes [70, 82, 43, 68, 99] for composite or multi-body obstacles.

A useful variant of the DD-SOV method is its ’extension’, by Rayleigh [71, 9, 97], to obstacles of other-than-canonical shape. Another interesting variant of the DD-SOV method is the so-called null-field or T-matrix method, attributed to Waterman [93, 94, 9, 97]. Both the Rayleigh and Waterman methods will not be discussed any more in the present study although they too provide grounds for further approximations (via perturbations [71, 24, 97] for a size parameter of the obstacle that is small compared to the wavelength, and [72] for a size parameter that is large compared to the wavelength, as well as other asymptotic procedures [97]), and are, in any case, useful candidates for treating inverse scattering problems [78, 34]. A final line of attack, not discussed further on in our study, is that of the so-called wave-splitting method, which has proved to be useful too in inverse problems for continuously- or discontinuously-varying (in terms of position) constitutive parameters. The wave splitting technique is exact and easily-applicable only for layered one-dimensional slab-or layer-like obstacles [20, 19], but approximate and not easily-applicable for 2D and 3D obstacles [39].

This brings us to the very object of our study which concerns integral equation methods for solving forward, as well as inverse (for the aforementioned reasons) scattering problems.

In the preceding paragraphs, we followed the modern practice of discussing the governing equations of elastodynamic and electrodynamic phenomena in their (partial) differential forms (PDE’s). Actually, the procedure of going from PDE’s to integral equations is contrary to the way the said dynamic phenomena were first (physically-) formulated in integral terms and then reformulated in differential terms. For instance, Faraday’s law, linking the electric field $\mathbf{E}$ to temporal ($t$) changes of the magnetic field $\mathbf{B}$, is (spatially) in integral form, and from this, one can, with the help of Stoke’s theorem, derive the PDE known as (one of the) Maxwell equation(s) $\nabla\times\mathbf{E}=-\frac{\partial\mathbf{B}}{\partial t}$. Our procedure is perhaps justified by the fact that we are searching for ways, other than those mentioned previously, of actually solving the governing PDE, in the hope that these ways, based on integral representations of the wavefield, may offer useful alternatives to FD, FE, SE and DD-SOV methods, particularly in the inverse problem context.

Re-consider the radiation problem in a homogeneous, isotropic medium of infinite extent in all directions. In electrostatics, a radiated field reduces to a static electric field $\mathbf{E}$ and Coulomb’s law is expressed by a domain (volume for a 3D source, area for a 2D source) integral whose integrand is the product of the source density with the Green’s function of the medium (called the free-space Green’s function (FPGF) of this medium), the latter being the field radiated by a point or line source (and known by the law of action at a distance). By employing the divergence theorem, one can obtain a PDE relating the divergence of $\mathbf{E}$ to the source (i.e., charge) density divided by a (spatially-constant) constitutive parameter of the medium, and from this one can show [66] that $\nabla\times\mathbf{E}=\mathbf{0}$ which is another PDE. This illustrates once again that a physically-established (Coulomb’s) law, expressed in integral form, is, in fact, the starting point of the means by which the usual object of interest for mathematical theory (i.e., a PDE) is obtained. Of course, this does not exclude the reverse procedure, which is to obtain the integral representation from the PDE.

Consequently, we now address the issue of how, starting with the PDE of our scattering problem, to represent, by a boundary or domain integral, the scattered field. Henceforth, we shall be concerned mostly with elastic wave, and especially with, acoustic wave problems (because the latter are scalar in nature; note that for this reason, many publications dealing with geophysical elastic wave problems, or acoustic problems in which the obstacles are manifestly solid surrounded by a fluid host, make the so-called acoustic approximation which partially-enables the description of vectorial field phenomena by scalar field phenomena) [12, 14, 32, 48, 47]. Moreover, since the time domain expressions of the governing equations relative to the wavefield are related to their frequency domain expressions by a simple Fourier transform, we will be concerned (unless stated otherwise) with frequency domain expressions of the wavefield.

The PDE of interest is of the form $\mathcal{L}(\mathbf{\mathfrak{c}})u=-s$, where $\mathcal{L}$ is a linear partial differential operator, $u$ the field and $s$ the source density. $\mathcal{L}$, $u$ and $s$ are functions of the space and frequency) and $\mathcal{L}$ is second-order in terms of the space variables. $\mathcal{L}(\mathbf{\mathfrak{c}})$ contains certain coefficients (grouped into the vector $\mathbf{\mathfrak{c}}$) which define the constitutive properties of the medium filling the domain $\Omega$ in which the field $u$ is searched.

In the aforementioned radiation problem, $\Omega$ is all space, $s$ is usually of finite support in $\Omega$ and the coefficient vector $\mathbf{\mathfrak{c}}$ is constant (denoted by $\mathbf{\mathfrak{c}}^{[0]}$) with respect to the space, and also usually the frequency, variables. The free-space Green’s function (FPGF), denoted by $g$, is the solution of the radiation problem when $s=-\delta$, wherein $\delta$ is the Dirac delta distribution which is a function of the space variables of both the ’receiver’ (i.e., where the field is ’sensed’) and the ’emitter’ (i.e., the point or line applied source). Since $\Omega$ (for this problem) is unbounded, the radiation condition tells us that the radiated wave must be outgoing at points far-removed from the location of the source. Furthermore, the field is, by the very nature of $\mathcal{L}$, bounded everywhere except at the location of the source. It turns out that the field $g$ of this PDE can be represented by a sum of plane waves, and that the amplitudes of these waves can be fairly-easily obtained by projection with the help of the radiation condition [61]. Thus, the solution $g$ of $\mathcal{L}(\mathbf{\mathfrak{c}}^{[0]})g=-\delta$ can be considered to be known (a similar technique is employed in [23] to obtain the free-space Green’s function for the elastodynamic wave equation).

When a part of the previous $\Omega$ is occupied by the obstacle (assumed to be isotropic), we say that the latter is included in the domain $\Omega_{1}$ whereas the remaining portion of $\Omega$ is $\Omega_{0}$, both being separated by the boundary $\Gamma$. We already know that the medium in $\Omega_{0}$ is characterized by the constant vector $\mathbf{\mathfrak{c}}^{[0]}$. Generally speaking, the medium in $\Omega_{1}$ is inhomogeneous, which means that $\mathfrak{c}$ therein, denoted by $\mathbf{\mathfrak{c}}^{[1]}$, is not constant with respect to the space variables. This fact makes the problem of the representation of the field $u$ more complex.

The first trick is to write the governing PDE $\mathcal{L}(\mathbf{\mathfrak{c}})u=-s$ in $\Omega$ as $\mathcal{L}(\mathbf{\mathfrak{c}}^{[0]})u=-s-[\mathcal{L}(\mathbf{\mathfrak{c}})u-\mathcal{L}(\mathbf{\mathfrak{c}}^{[0]})u]$ and to consider the factor $[~{}]$ on the right-hand side of this equation as a new source density which is termed the induced source density because it is what induces (i.e., gives rise to) the scattered (=diffracted) field $u^{d}$, much as the the applied source density $s$ gives rise to the ’incident’ field $u^{i}$, their sum being none other than the total field $u=u^{i}+u^{d}$. The second trick is to perform the multiplications: $g\mathcal{L}(\mathbf{\mathfrak{c}}^{[0]})u=-gs-g[\mathcal{L}(\mathbf{\mathfrak{c}})u-\mathcal{L}(\mathbf{\mathfrak{c}}^{[0]})u]$ and $u\mathcal{L}(\mathbf{\mathfrak{c}}^{[0]})g=-u\delta$, subtract these two equations, and then integrate the result over $\Omega$. The third, and last, trick is to employ Green’s second identity, the sifting property of the Dirac delta, and the radiation condition satisfied   by the   field   on   the ’boundary at  infinity’,   to   obtain   the   domain  integral  (DI)  representation   $u=u^{i}+\int_{\Omega}g\left[\mathcal{L}(\mathbf{\mathfrak{c}})u-\mathcal{L}(\mathbf{\mathfrak{c}}^{[0]})u\right]d\Omega$ of $u$ in $\Omega$, wherein $u^{i}=\int_{\Omega}gsd\Omega$. We insist on the word ’representation’ because this DI relation does not solve the problem since what we are looking for, i.e., $u$, appears not only on the left-hand side of the relation, but also within the induced source density term $[~{}]$ on the right hand side. Solving the forward problem, i.e., actually determining $u$ throughout $\Omega$ by means of this DI representation, given $s$, $\mathbf{\mathfrak{c}}^{[0]}$, $\mathbf{\mathfrak{c}}^{[1]}$, and $\Gamma$, is another story that will be schematized further on in this introduction and treated in more detail in sect. 4.6.

The PDE for a complete (incident plus scattered) acoustic wavefield in a fluid-like medium with spatially-varying density and compressibility is generally-attributed to Bergman [6]. We call this the Bergman PDE (BPDE). A more detailed development of the way to obtain the BPDE from first mechanical continua principles can be found in [11] from which is borrowed the material in sect. 2 of the present study. A DI representation of the wavefield is also associated with Bergman’s name (see also [57]) and is named the Bergman domain integral representation (BDIR) hereafter. Actually there exist several versions of the BDIR [6, 86, 56, 49, 16] that are shown, in sect. 4.6.10 herein to be equivalent under certain conditions. In the same section, we also show that the 2D version of the elastodynamic domain integral representation of the antiplane component of the scattered displacement field in [55] is equivalent to the BDIR provided certain associations are made between the fields and constitutive parameters of the two problems. The 3D DI representation of the elastodynamic wavefield was proposed in [55, 51, 50, 67], but, since the stress herein will be on scalar fields, we shall not be interested in this DI representation any longer. Finally, an integral (actually, integrodifferential) representation similar to the BDIR, but for the electrodynamic wavefield, was first obtained by Shifrin [81]. Further discussions of Shifrin’s representation and how to employ it to solve electrodynamic scattering problems can be found in [2, 7, 8, 86] and we shall not discuss it any further for the reasons mentioned previously.

Now, let us focus on the scattered-wave component of the BDIR. The domain $\Omega$ of this representation is actually all space (i.e., $\mathbb{R}^{n}$ in one ($n=1$), two or three-dimensional problems). As we saw previously, the integrand of the DI is actually the product of the FSGF (i.e., g) with the induced source density, the latter involving: (1) the field within the obstacle and its gradient, and (2) the velocities and densities of the inner and outer media, or (in each of the two media) any two of the acoustic constitutive parameters (which may, or may not, depend on the spatial variables): mass density, wavespeed, bulk modulus and compressibility (we choose mass density and wavespeed hereafter). Since (1) depends on (2) and multiplies (2) in the integrand, the scattered field is a nonlinear function of the $\mathbf{\mathfrak{c}}$, this being the principal source of difficulty in the inverse problem (of determining one, several, or all of $s$, $\mathbf{\mathfrak{c}}^{[1]}$ and $\Gamma$ from the knowledge of $u$ in some subdomain of $\Omega$ and for a given set of frequencies). However, in the forward problem context we are not confronted with such a difficulty since we search for $u$ which appears only linearly in the integrand of the DI (although it is in a multiplicative relation with $g$ which is another field function, but the latter is known).

At present, we address the question of how to employ the BDIR to solve for $u$ in $\Omega$ (i.e., the forward problem). The purely-numerical method is schematically to first divide $\Omega$ into two subdomains $\Omega_{ext}$ and $\Omega_{int}$, with the properties of $\Omega_{ext}$ being the spatially-invariant $\mathbf{\mathfrak{c}}^{[0]}$. It then follows that the portion of the DI concerning $\Omega_{ext}$ vanishes so that we obtain a sort of integral equation of the second kind $u=u^{i}+\int_{\Omega_{int}}g[~{}]d\Omega~{};~{}\in\Omega_{int}$ which could be the means for obtaining $u\in\Omega_{int}$ if $[~{}]$ were only a function of $u$. Unfortunately, this is not the case when the mass density depends on position within $\Omega_{int}$ because then $[~{}]$ depends both on $u$ and $\nabla u$ (within $\Omega_{int}$).

Two strategies have been proposed to cope with this problem. The first [13, 49, 65], which we call the p-q transformation in sect. 2.9, enables the elimination of the field gradient in the induced source density term, but this is obtained at the expense of rendering more difficult (especially when there are discontinuites of density) the way the density component of $\mathbf{\mathfrak{c}}$ is taken into account. The second strategy is to suppose that the mass density does not depend on position within all of $\Omega$ [63, 62] and the majority of what is offered in the present investigation will deal with a critical examination of this supposition.

The constant-density supposition enables the DI to be reduced to what is often called the Lippmann-Schwinger DI (LSDI) and the aforementioned integral equation of the second kind is then called the Lippmann-Schwinger integral equation (LSIE) [54, 58, 74, 98, 102] for the sole unknown function $u\in\Omega_{int}$. Note that this unknown is of finite support when, as is often the case, the obstacle occupies a domain of finite extension.

The preference for the DI rather then for the PDE methods, in relation to solving scattering problems in brute-force numerical manner, is due to the fact that the support of the unknown function of the PDE is all space (i.e., $\Omega$) and must be fenced-in by artificial (usually of the absorption variety [36, 92]) boundary conditions whereas that of the LSDI is only a subdomain (i.e., $\Omega_{int}$) of $\Omega$. Thus, the LSIE lends itself quite easily to brute-force numerical resolution [74, 37, 38, 97], but also, and most often (in the scientific literature), to iterative resolution, starting with the so-called Born approximation [75, 38], which amounts to assuming $u=u^{i}~{};~{}\in\Omega_{int}$. The last step in the LSIE scheme is to obtain $u$ in $\Omega_{ext}$ by inserting the previously-obtained $u~{};~{}in\Omega_{int}$ into the LS field representation.

So much for the use of the LS method in conjunction with the forward scattering problem. The question now is whether there is any advantage in using such a DI method in the inverse scattering problem. For many years (actually, until powerful desktop computers became available), the scientific community responded almost unanimously in the positive manner for several reasons. The first has to do with what was written in the last lines of the preceding paragraph. But, quite early, it was realized that, even with the LS representation (which is based on the constant-density assumption, so that implicitly one renounces at retrieving the mass density in $\Omega_{int}$), the remaining to-be-determined parameter (more often a position-dependent function), i.e., the wavespeed in $\Omega_{int}$, intervenes in nonlinear manner in the induced source density so that the LSIE is nonlinear in terms of this parameter (or function). Recall that it was linear in terms of $u$. Thus, retrieving the wavespeed from the LSIE turns out to be a difficult task.

Let us examine the ingredients of this task. The problem, for a known $s$, $\Gamma$, $\mathbf{\mathfrak{c}}^{[ext]}$ and measured field (in $\Omega_{mes}\subset\Omega_{ext}$), is to determine the wavenumber contrast function $\chi(\mathbf{\mathfrak{c}}^{[ext]},\mathbf{\mathfrak{c}}^{[int]})$ via the LSIE $u=u^{i}+\int_{\Omega_{int}}g\chi ud\Omega~{};~{}\in\Omega_{mes}$. Since the $u$ in the integrand is unknown (i.e., not measured), there does not appear to exist any obvious way to ’invert’ this IE so as to retrieve $\chi$ [40].

This led to the idea of approximating $u$ within the induced source density term (the latter giving rise to what is termed the first-order Born approximation of the scattered field [75, 38, 42]) so that the LSIE becomes a linear function of the to-be-retrieved parameter. Above all, $u$ in the integrand is now a known function (just like $g$) so that the only unknown (in the integrand) is $\chi$. This, of course, makes it much easier to solve the LSIE for the wavespeed. To do this, various approximations (asymptotic far-field, high frequency, etc.) were made of the Green’s function [89] whereupon it was realized that $\chi$ emerges as a sort of (explicit) Fourier transform of $\chi$, which is the basis of what has become to be known as Diffraction Tomography ((DT)[73, 21, 22, 52, 75, 95, 59, 18, 38], widely used in geophysical and bioacoustical applications [42]. Other approximations of the Green’s function [58] enabled the explicit computation of the Fréchet derivative [60] in the full-wave inversion (FWI) methods of inversion [33], the latter based on retrieval by optimizing (i.e., seeking for the minimum) of a functional involving the difference between $u$ in $\Omega_{mes}$ and $u^{i}+\int_{\Omega_{int}}g\chi u^{i}d\Omega~{};~{}\in\Omega_{mes}$. Naturally, there exist refinements of the FWI method [97, 102, 101, 79, 99] which make use of the DI representation of the scattered field, but since our study is not focused on the inverse problem, we shall not delve further into this issue.

The main issue under our scrutiny is that of whether (or not) the constant mass-density assumption is justified, this, of course, being of some importance in the forward problem, and perhaps of greater importance in the inverse problem. Surprisingly, there exist practically no studies (perhaps [88, 5, 42]) that treat this issue in depth. In the bioacoustics community, which is concerned with imaging organs within the body [101], it is customary, and probably justified for in situ studies, to either consider the organ to be surrounded by a medium whose mass density is quite close to that of the organ [101], or in laboratory studies to surround the obstacle with a fluid medium (usually water) for the same purpose [38]. In the geophysics community, this assumption appears a priori to be less justified, but is nonetheless frequently made for reasons of ’simplicity’ or ’mathematical convenience’ [75]. In the porous media acoustics context (notably for applications dealing with the absorption of airborne sound) it has been found [100] that the effective (bulk) density of typical air-saturated porous materials can be complex, dispersive, and quite different (real part) from the density of air. Sometimes, theoretical studies of scattering are undertaken in a quite general setting (i.e., without assuming constant mass density), but the mass density is taken to be either continuously-varying [90] or constant in the final, numerical stage [32]. The reason for the latter assumption is (as was suggested in the preceding lines) that it is not easy to carry out the numerical calculations via the DI formulation when the mass density is not constant in $\Omega$.

This is the reason why a less-drastic assumption has been invoked in connection with the integral formulation. What surely comes to mind is to assume that the mass density of the obstacle be different from that of the host medium, but constant within $\Omega_{int}$. At first, one is not required to make such an assumption regarding the wavespeed, so that one can show, quite easily (this is done in sect. 4.4 of the present study and appeals to domain decomposition) that the scattered field can be represented by a sum of a DI and a boundary integral (BI), involving $u$ in $\Omega_{int}$ and $u$ on the boundary $\Gamma_{ext-int}$ separating $\Omega_{ext}$ from $\Omega_{int}$. In the forward problem, these two functions are unknown, and can be shown to be obtained via a system of two coupled equations involving both DI’s and BI’s [37, 30, 31], this being called the DI-BI scheme. This system lends itself to an iterative resolution scheme similar to that of the previously-evoked sole domain integral equation, or else, it can be solved in brute-force numerical manner. As regards the brute-force manner, the amount of required computations is less than in the DI method (before invocation of the constant mass assumption, which is more-properly called the DI-DI method because it involves two domain integrals, one involving variations of the mass density and the other variations of the wavespeed) due to the fact that one domain integral (involving a large amount of discretized unknowns) is replaced by a boundary integral (involving a much lesser amount of discretized unknowns). However, it should be recalled that this advantage is obtained at the expense of making a strong assumption regarding the mass density.

With the same assumptions on the mass density and wavespeed as in the preceding paragraph, it can be shown (as is done in sect. 4.5 herein), starting with domain decomposition, that the scattered field can be represented by a sum of a DI and a boundary integral (BI), involving $u$ in $\Omega_{int}$ and $\nabla u$ on the boundary $\Gamma_{ext-int}$. In the forward problem, these two functions are unknown, and are obtained via a system of two coupled equations involving both BI’s and DI’s, this being called the BI-DI scheme. This system again lends itself to an iterative resolution scheme similar to that of the previously-evoked sole domain integral equation, or else, it can be solved in brute-force numerical manner. For the same reasons as previously, the BI-DI brute-force scheme involves less computations than the DI-DI scheme. However, it should be recalled that this advantage is obtained at the expense of making a strong assumption regarding the mass density.

With more drastic assumptions: (a) mass densities spatially-constant, but different from each other, in the host and obstacle and (b) wavespeeds spatially-constant, but different from each other, in the host and obstacle, one can show (as is done in sect. 4.3 herein), starting with domain decomposition, that the scattered field in $\Omega_{ext}$ is represented by a BI involving $u$ and $\nabla u$ on $\Gamma_{ext-int}$ and the scattered field in $\Omega_{int}$ is represented by another BI involving the same two field functions on the boundary. The latter can then be shown to constitute the unknowns of a coupled system of boundary integral equations (this method is thus termed the BI-BI method [10, 97]) which again lends itself to an interative type of resolution scheme or else can be solved in brute-force manner. Since two domain integrals in the DI-DI scheme are replaced by two boundary integrals in the BI-BI scheme [77], the latter requires much less computations than the DI-DI scheme when brute force is employed. However, this enormous advantage is obtained at the expense of making the drastic assumptions regarding both the wavespeed and the density (see also sect.4.7).

At this point, one may wonder whether the DI-DI, DI-BI, BI-DI and BI-BI methods give rise to the same solution (i.e., the exact solution, as they should) when the assumptions of the BI-BI regarding mass density and wavespeed are made. We have been unable to make this demonstration in the general case of an obstacle of arbitrary shape, because, among other reasons, we do not dispose of an exact, reference solution for obstacles of arbitrary shape. However, we do dispose of such a solution for an obstacle of canonical shape, and since we are here interested essentially in 2D problems, this obstacle was chosen to be a circular cylinder. As mentioned before, the scattering problem for this homogeneous obstacle can be found by the DD-SOV method. Thus, we shall employ this solution (which is in explicit, functional form) as the reference (a similar procedure, using the DD-SOV solution for elastic waves impinging on a spherical obstacle as a reference, was followed in [41]) for determining whether the other four methods lead to the same (or different) solutions. As will be shown, it turns out that all these methods lead to the same solution for the homogeneous circular cylinder canonical obstacle even though the way of obtaining the solution is relatively simple only via the BI-BI method. Futhermore, and more importantly, only the DI-DI method leads to an intermediate form of the solution which readily separates the contributions of the mass density contrast from the wavespeed contrast. This enables us to generate the solution via an iterative scheme that is equivalent to its expansion as a geometrical series in terms of the mass density contrast parameter. Taking only the first term of this series is equivalent to approximating the solution by the one obtained from the constant mass density assumption. Finally, we show that retaining only three terms in this series is generally sufficient to obtain a quite accurate approximation of the solution for realistic mass density contrasts.

2 Governing equations

The material presented in this section is taken mostly from [11]. A solid or fluid (the medium of interest in acoustics) is here considered to be a linear, isotropic elastic continuum.

2.1 Conservation of momentum

The conservation of momentum relation in a mechanical material continuum is expressed by

[TABLE]

wherein: $\sigma_{kl}$ are the (cartesian) components of the stress tensor, $\rho$ the mass density, $f_{l}$ the components of the applied body force, and $\mathcal{F}_{,k}$ designates the partial derivative of $\mathcal{F}$ with respect to $x_{k}$ (with $x_{k}$ the $k$-th cartesian coordinate). By extension, $\mathcal{F}_{,t}$ designates the partial derivative of $\mathcal{F}$ with respect to $t$, with $t$ the time. Furthermore, $v_{l}$ are the components of particle velocity, $\dot{v}_{l}$ the components of particle acceleration (the dot above a variable meaning ”material derivative”, i.e., $\dot{\mathcal{F}}_{k}:=\mathcal{F}_{,kt}+\mathcal{F}_{k,l}\mathcal{F}_{l}$).

2.2 Conservation of momentum for small deformations

For small deformations, $\dot{v}_{l}\simeq u_{l,tt}$, wherein $u_{l}$ are the cartesian components of the particle displacement vector, and $\mathcal{F}_{,tt}$ designates the second partial derivative of $\mathcal{F}$ with respect to $t$, so that

[TABLE]

2.3 Constitutive relations in a linear material continuum

In a linear mechanical material continuum the stress tensor is related to the strain tensor (whose components are $\varepsilon_{mn}$) by

[TABLE]

with

[TABLE]

and $c_{klmn}$ the fourth-order material coefficient tensor.

2.4 Constitutive relations in a linear, isotropic, elastic continuum

In a linear, elastic continuum, (3) takes the form given by the Cauchy-Hooke relation

[TABLE]

wherein $\lambda,~{}\mu$ are the Lamé coefficients, $\delta_{kl}$ the Kronecker symbol, and repeated indices meaning summation over this index (Einstein convention), i.e., $\varepsilon_{mm}:=\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33}$.

Consequently,

[TABLE]

2.5 Wave equation in a linear, isotropic, elastic

medium (Navier’s equation)

Introducing (4) into (6) gives

[TABLE]

or, collecting terms

[TABLE]

so that the conservation of momentum relation (2) becomes

[TABLE]

or, in its linearized version

[TABLE]

2.6 Governing equations in an inviscid fluid

When the medium is a (heterogeneous) ideal (i.e., inviscid; see e.g., [36] for a method of handling the viscous fluid case) fluid, $\mu=0$, so that

[TABLE]

With the definition

[TABLE]

wherein $p$ denotes the pressure in the ideal fluid, it follows that

[TABLE]

which shows that the stress tensor is diagonal in an ideal fluid.

Then

[TABLE]

so that the conservation of momentum relation in an ideal fluid becomes

[TABLE]

or, in its linearized form

[TABLE]

2.7 Time domain wave equation in a fluid

Eq. (12) gives rise to

[TABLE]

so that (16) becomes

[TABLE]

which is the wave equation in an ideal fluid in terms of particle displacement.

An alternative (i.e., a partial differential equation exclusively in terms of pressure) can be found as follows. Taking derivatives of (16) gives

[TABLE]

However, (12) gives rise, under the assumption that $\lambda$ does not depend on $t$, to

[TABLE]

Another consequence of (16) is

[TABLE]

so that inserting (20) and (21) into (19) results in

[TABLE]

or

[TABLE]

which is the wave equation, in terms of pressure, in an ideal non-homogeneous fluid.

Letting

[TABLE]

with $c$, $K$ and $\kappa$ the phase velocity, isentropic bulk modulus and isentropic compressibility within/of the fluid medium, we obtain

[TABLE]

or, in vectorial notation (recalling that $\rho$ and $c$ do not depend on $t)$,

[TABLE]

wherein $\mathbf{x}=(x_{1},x_{2},x_{3})$. This is the time domain version of the Bergman equation (BPDE in sect. 1) [6].

Note that up till now we assumed the fluid to be lossless (inviscid) which means that $c$ and $\rho$ are real and independent of $t$. A convenient way to account for losses in the fluid is to consider $c$ and/or $\rho$ to be complex and dispersive, i.e.,

[TABLE]

wherein $c^{\prime\prime}<0$ and $\rho^{\prime\prime}>0$ for a passive (i.e., lossy ) fluid (the case considered hereafter) and $c^{\prime\prime}>0$ and $\rho^{\prime\prime}<0$ for an active fluid. This is what we shall do in the frequency-domain formulation of the wave equation.

2.8 Frequency domain wave equation in an inviscid or viscous fluid

By expanding the time domain pressure and force in the Fourier integrals

[TABLE]

[TABLE]

wherein $\omega=2\pi f$ is the angular frequency and $f$ the frequency, one obtains the frequency domain expression of the wave equation

[TABLE]

wherein

[TABLE]

is the wavenumber in the fluid, which is a generally-complex quantity due to the fact that $c$ is assumed to be generally-complex (to account for viscosity). Eq. (30) is the frequency domain version of the Bergman equation (BPDE).

2.9 Getting rid of the gradient of pressure term via the $p-q$ transformation

Eq. (30) was

[TABLE]

As we shall see later on, the presence of the gradient of pressure term (i.e., $p_{,ll}$) is a source of nuisances. To eliminate it, we attempt the change of variables [6, 57]

[TABLE]

whence

[TABLE]

so that (32) becomes

[TABLE]

Thus, the problem reduces to getting rid of $[~{}]$, i.e.,

[TABLE]

the only possible solution of which is $m=-1/2$. Thus, the proper transformation of variables is

[TABLE]

whence (35) becomes

[TABLE]

It can be shown that

[TABLE]

so that we finally obtain the inhomogeneous Helmholtz equation for $q$

[TABLE]

wherein

[TABLE]

Thus, the problem of finding the pressure $p$ reduces to that of first finding $q$ via (40) and then computing $p$ from $q$ via (37). This is apparently the method adopted in publications such as [49]. We shall return to this issue later on in sect. 4.6.3 of the present study.

2.10 Specifics of 2D elastic and acoustic wave mechanics

Since the remainder of this study will be concerned mostly with two-dimensional wave-mechanical problems, we now consider the particularities of such problems.

Let us return to the equation (10) for wave motion in an elastic solid in which we replace the density $\rho$ and two Lamé parameters $\lambda$ and $\mu$ by the letters $R$, $L$ and $M$ whereas $u$ and $f$ are replaced by $U$ and $F$.

[TABLE]

In a general 3D situation, $R$, $L$ and $M$ and $\mathbf{F}$ depend on the three cartesian coordinates $x_{1},x_{2},x_{3}$ so that the three components of $\mathbf{U}$ depend as well on these three spatial coordinates.

We now consider the so-called 2D situation in which

[TABLE]

and we name the $x_{1}-x_{2}$ plane the sagittal plane.

It is then easy to show that

[TABLE]

whence (42) becomes

[TABLE]

from which we obtain the three relations:

[TABLE]

[TABLE]

[TABLE]

The first two relations describe the wave motion (i.e., exclusively the $U_{1}$ and $U_{2}$ components of $\mathbf{U}$) in the sagittal plane so that this motion is termed ’in-plane’ and is due to in-plane forces. The third relation describes the wave motion (i.e., exclusively the $U_{3}$ component of $\mathbf{U}$) along the $x_{3}$ direction perpendicular to the sagittal plane so that this motion is termed ’out-of plane’ and is due to out-of-plane forces. Of utmost importance is to note that the out-of-plane motion is totally decoupled from the in-plane motion in this 2D situation.

We are especially interested here in this third relation (48) which can be re-written as

[TABLE]

wherein

[TABLE]

and $C$ is the out-of-plane (or shear-horizontal) phase velocity (i.e., of $U_{3}$ motion) in the solid.

Until further notice, we make the following definitions: $\mathbf{x}:=(x_{1},x_{2})$, $\nabla_{\mathbf{x}}:=(\frac{\partial}{\partial x_{1}},\frac{\partial}{\partial x_{2}})$, $U(\mathbf{x}):=U_{3}(\mathbf{x})$, $F(\mathbf{x}):=F_{3}(\mathbf{x})$, whence (49) can be written in the vectorial form

[TABLE]

The 2D version of the acoustic wave equation in inviscid fluids is, by virtue of (26)

[TABLE]

The important point is to notice that these two wave equations are exactly of the same form if the following associations are made: $U\leftrightarrow p$, $M\leftrightarrow 1/\rho$, $C\leftrightarrow c$ and $F/C^{2}\leftrightarrow-\rho\nabla\cdot\mathbf{f}$ which is consistent with the previous definitions $c^{2}=K/\rho$ and $C^{2}=M/R$ provided $1/R\leftrightarrow K$. This means that everything that will be found for the 2D pressure wave in fluids will be applicable as well for the 2D out-of-plane particle displacement wave in solids provided that the acoustic density $\rho$ is associated with the reciprocal of the elastic shear modulus $1/M$. Thus, the assumption of constant density in 2D acoustics is equivalent to the assumption of constant shear modulus in 2D elasticity (in which case $c$ is not assumed spatially-constant in acoustics and $C$ is not assumed spatially-constant in elasticity).

3 Acoustic radiation from applied sources in the absence of an obstacle in

unbounded space occupied by a homogeneous, lossless fluid

3.1 Description of the problem

The configuration is depicted in fig. 1. Until further notice, it represents a 3D situation. The source, of bounded support $\Omega^{s}$, is located in $\mathbb{R}^{3}$ which is filled with a homogeneous medium $M^{0}(\rho^{[0]},c^{[0]})$ such that $\rho^{[0]}$ and $c^{[0]}$ are position- and frequency-independent.

3.2 Governing equations

From now on, we place the analysis in the frequency domain, so that the governing equation (30)is, due to the supposed-homogeneous nature of the fluid medium,

[TABLE]

wherein $\omega$ is implicit in the arguments of the various functions, $k^{[0]}=\omega/c^{[0]}$ is position-independent, and $s:=-\rho\nabla\cdot\mathbf{f}$ represents the source density associated with applied forces, whose support is assumed to be the finite domain $\Omega^{s}$.

In addition, the pressure, which is bounded everywhere except in $\Omega^{s}$, satisfies the radiation condition whose mathematical translation is

[TABLE]

wherein $r^{\prime}$ is the radial coordinate in the $\mathbf{x}^{\prime}$ coordinate system.

3.3 Field representation

The so-called free-space Green’s function, which is bounded everywhere except for $\|\mathbf{x}-\mathbf{x^{\prime}}\|=0$, satisfies

[TABLE]

wherein $\delta(~{})$ is the Dirac delta distribution whose sifting property is expressed by

[TABLE]

Furthermore, this Green’s function satisfies the radiation condition

[TABLE]

It follows that

[TABLE]

Applying Green’s second identity in $\Omega=\mathbb{R}^{3}$, assumed to be ”closed” by a hypersphere of infinite radius whose boundary is $\Gamma_{\infty}$ (see fig. 1), leads to

[TABLE]

The second integral can be written as

[TABLE]

and since $p^{[0]}$ and $G^{[0]}$ are bounded on $\Gamma_{\infty}$ and these two functions satisfy the radiations conditions (54) and (57), both of the integrals on the right hand side of this expression vanish. Thus,

[TABLE]

or, using the fact that the support of the applied sources is $\Omega^{s}$

[TABLE]

Since $G^{[0]}$ is a function of $\mathbf{x}-\mathbf{x^{\prime}}$, one sees that (62) expresses the rather intuitive fact that the pressure field radiated by applied sources of density $s$ is expressed by a convolution between the latter and the pressure field radiated by a ”point” source.

It is readily verified by means of (62) and (56) that when $s$ corresponds to a ”point” source, the radiated field is the free-space Green’s function.

Remark

An important feature of (62) is that not only does it provide a domain integral representation of the radiated field, but also the solution to the forward problem of radiation by sources of known density $s$ and known support $\Omega^{s}$, the kernel of the integral by which this solution is computed being known as well.

The inverse radiation problem (also called the inverse source problem), which is to determine the sources from the radiated field, cannot be solved in such a simple manner. However, this inverse problem is linear in that the pressure depends linearly on the pressure field data.

3.4 The 2D free-space Green’s function

From this point on, everything will concern the 2D situation ($x_{3}$ being the ignorable coordinate due to the invariance of the constitutive parameters of the medium, impressed force and consequently the pressure field with respect to $x_{3}$) in which $\mathbf{x}=(x_{1},x_{2})$ and $\mathbf{x^{\prime}}=(x^{\prime}_{1},x^{\prime}_{2})$ (in cartesian coordinates).

As shown previously in this study, the free-space Green’s function (GF)constitutes a fundamental tool for the study of wave propagation. Rather than derive the 2D GF from the governing equations, we shall assume known the expression ([61], p. 822, 891) for this GF and show that it satisfies the governing equation and radiation condition.

The 2D GF is

[TABLE]

wherein $H_{0}^{(1)}$ is the zeroth-order Hankel function of the first kind [1]. This function admits the following representation in polar coordinates (([61], p. 827) in which $\mathbf{x}=(r,\theta)$ and $\mathbf{x^{\prime}}=(r^{\prime},\theta^{\prime})$:

[TABLE]

wherein $J_{n}(.)$ is the $n$-th order Bessel function and $H(\zeta)$ the Heaviside distribution ($H(\zeta>0)=1$, $H(\zeta<0)=0$).

3.4.1 Demonstration that the 2D GF satisfies the 2D Helmholtz equation

We first want to show that(64) satisfies the Helmholtz equation in polar coordinates ([61], p. 825; [29])

[TABLE]

From the facts that

[TABLE]

we obtain:

[TABLE]

wherein $\dot{H}_{n}(\zeta):=\frac{d}{d\zeta}H_{n}(\zeta)$, $\ddot{H}_{n}(\zeta):=\frac{d^{2}}{d\zeta^{2}}H_{n}(\zeta)$, $\dot{J}_{n}(\zeta):=\frac{d}{d\zeta}J_{n}(\zeta)$, $\ddot{J}_{n}(\zeta):=\frac{d^{2}}{d\zeta^{2}}J_{n}(\zeta)$.

[TABLE]

But ([1], p. 360)

[TABLE]

so that

[TABLE]

The Poisson sum formula ([61], p. 483) tells us that

[TABLE]

whence

[TABLE]

Consequently

[TABLE]

It ensues that

[TABLE]

But, by means of ([1], p. 358), and for $w_{n}=J_{n}$ or $w_{n}=H_{n}^{(1)}$:

[TABLE]

whence

[TABLE]

or

[TABLE]

which demonstrates (65).

3.4.2 Demonstration that the 2D GF satisfies the radiation condition

Our second task is to show that our 2D GF satisfies the radiation condition;

[TABLE]

assuming that $r^{\prime}$ is bounded. This means that we are interested in the case $r>>r^{\prime}$ for which $H(r-r^{\prime})=1$ and $H(r^{\prime}-r)=0$, so that

[TABLE]

[TABLE]

We make use of the asymptotic forms of the Hankel functions ([1], p. 364)

[TABLE]

whence

[TABLE]

it follows that

[TABLE]

so that (78) is satisfied.

4 Scattering of the sound radiated from applied sources

by a generally-heterogeneous fluid-like obstacle within an unbounded, homogeneous, background fluid medium

4.1 Description of the problem

The problem differs from the previous one in that the formerly-homogeneous fluid medium now contains a generally-heterogeneous, spatially-bounded or unbounded (in one direction only) obstacle (see fig. 2 for the case of a bounded obstacle). This obstacle can be composed of a fluid-like material with either continuously-varying or discontinuously-varying physical properties. Thus, the domain of the obstacle can be simply- or multiply-connected. In the latter case, the obstacle is actually a set of smaller obstacles, all contained within the virtual boundary of the larger finite-domain obstacle.

The beauty of the domain integral (DI) method, which will be discussed in detail further on, is that it can cope with a great variety of topologies, all with the same formalism. This feature is particularly-useful in the inverse problem context.

As will be shown, the obstacle, by its very presence, gives rise to a field outside (and inside) itself that is different from what it was in the the previous section. The difference between the two total fields in these two problems is termed the scattered or diffracted field and designated by $p^{d}$, and what was formerly called the field due to applied sources is here called the incident field and designated by $p^{i}$.

The unbounded portion of space outside of the obstacle and support of the sources is called $\Omega_{0}$ and is occupied by the homogeneous medium $M^{[0]}(\rho^{[0]},c^{[0]})$, with $\rho^{[0]},~{}c^{[0]}$ constants with respect to position. The obstacle forms the domain $\Omega_{1}$ and is occupied by the generally-heterogeneous fluid-like medium $M(\rho^{[1]},c^{[1]})$ such that $\rho^{[1]}$ and $c^{[1]}$ are generally position-dependent, complex and dispersive in the sense of (). The domain $\Omega_{s}$ of the sources is assumed to be wholly exterior to $\Omega_{1}$.

4.2 The 2D scattering problem for a macroscopically-homogeneous, isotropic circular cylindrical obstacle by the domain decomposition separation of variables (DD-SOV) technique

This scattering problem is canonical in that it serves to introduce many of the concepts that are employed in more general scattering problems, including those in which the constant density assumption is not made. However, it is important to stress that the DD-SOV technique that we shall employ to solve this problem is not applicable to obstacles with more general (than circular cylindrical, spherical, etc.) shapes in which the Helmholtz equation is not separable.

Let $(r,\theta,z)$ be the cylindrical coordinates. The obstacle is a circular cylinder, infinite in the $z$-direction and whose cross-section is the disk $\varpi_{1}=\{r<a~{};~{}\forall\theta\in[0,2\pi]\}$ the boundary of which is the (circular) curve $\gamma=\{r=a~{};~{}\forall\theta\in[0,2\pi]\}$. The incident field is assumed to not depend on $z$, so that, due to the homogeneous, isotropic nature of the obstacle and the shape-invariance of the latter with respect to $z$, it ensues that the scattered field is also invariant with respect to $z$. In the sagittal (orthogonal to the $z$ axis) plane, the region outside of $\varpi_{1}$ is designated by $\varpi_{0}$, so that the problem is to determine the scattered fields in $\varpi_{0}$, $\varpi_{1}$, and on $\gamma$, it being understood that the media in the two 2D regions $\varpi_{0}$, $\varpi_{1}$ have the same properties as those in the previous 3D regions $\Omega_{0}$, $\Omega_{1}$ respectively. The case addressed here assumes that both the mass density of, and wavespeed in, the obstacle are constants as a function of position, i.e.,

[TABLE]

the consequence of which (as concerns the wavespeed) is

[TABLE]

it being recalled that $\rho^{[0]}$ and $k^{[0]}$ are also (different) constants.

4.2.1 Governing equations via DD

We search for the solution of the following problem:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

wherein

[TABLE]

4.2.2 Field representations via SOV

We had:

[TABLE]

which, after the introduction of the SOV representation of the 2D GF becomes

[TABLE]

Suppose the support of the source is between $r=r^{-}$ and $r=r^{+}>r^{-}$ in terms of $r$ and between $\theta=\theta^{-}$ and $\theta=\theta^{+}>\theta{-}$ in terms of $\theta$, with the understanding that $r^{-}>a$, $a$ being the largest value of $r$ along the boundary $\Gamma$ of the obstacle. Then

[TABLE]

so that, from the definition of the Heaviside distribution,

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

Consequently, the representation of the incident wave due to applied sources in $\Omega^{s}$ that we have to take into account in the transmission conditions (88)-(89) is

[TABLE]

wherein we have replaced $B_{n}^{i[0]}$ simply by $B_{n}^{[0]}$.

Note that even in the case of a plane incident wave (incident angle $\theta^{i}$),

[TABLE]

the SOV expression of $p^{i}(r,\theta)$, now valid for all $r$, is the same ([1], p. 361) as in (98), but with

[TABLE]

Note that the $B_{n}^{[0]}$ are known since $s(\mathbf{x})$ was assumed known in the case of source wave incidence and $B^{i}$ and $\theta^{i}$ are assumed to be known in the case of plane wave incidence.

The SOV representations of the scattered fields satisfying the boundedness condition (notably within the cylinder) and the radiation condition (outside and far from the cylinder) are:

[TABLE]

[TABLE]

with the understanding that the coefficients $A_{n}^{[0]}$ and $B_{n}^{[1]}$ are unknown.

4.2.3 Determination of the unknown coefficients via the transmission conditions

A consequence of the transmission conditions is:

[TABLE]

[TABLE]

or

[TABLE]

[TABLE]

so that, from the identity

[TABLE]

(wherein $\delta_{mn}$ is the Kronecker delta symbol), it ensues

[TABLE]

[TABLE]

the solution of which, on account of ([1], p. 360), is

[TABLE]

[TABLE]

This solves the problem of acoustic scattering by a macroscopically-homogeneous, isotropic circular cylinder.

Note that (110)-(111) constitute the exact reference solutions for the cylindrical wave amplitudes inside and outside the homogeneous cylinder whose constitutive properties (i.e., mass density and wavespeed) are generally-different from those of the host medium.

Note further that it is not easy to discern in (110)-(111) how the constant density assumption (i.e.,$\rho^{[1]}=\rho^{[0]}$) affects the scattering amplitudes, nor, more generally, how the latter (and the scattered wavefield) depend on the contrast of mass densities $\epsilon=\frac{\rho^{[1]}-\rho^{[0]}}{\rho^{[1]}}$. We shall return to this issue in sect. 4.6.

4.3 The coupled boundary integral (BI-BI) method applicable to the special case in which both the mass density of, and wavespeed in, the fluid-like obstacle do not depend on position

This case, which is treated in publications such as [37], again assumes that both the mass density of, and wavespeed in, the obstacle are constants as a function of position, i.e.,

[TABLE]

the consequence of which (as concerns the wavespeed) is

[TABLE]

it being recalled that $\rho^{[0]}$ and $k^{[0]}$ are also (different) constants.

4.3.1 BI-BI method: domain decomposition

Instead of employing the frequency domain wave equation in $\mathbb{R}^{3}$ (26), we decompose the latter into two wave equations, one for $\Omega_{0}$ in which the total pressure is $p^{[0]}$ and the other for $\Omega_{1}$ in which the total pressure is $p^{[1]}$, these two domains being separated by $\Gamma$, through which the pressure is continuous and the velocity potential is continuous. Thus, assuming once again that the support of the source is in $\Omega_{0}$:

[TABLE]

[TABLE]

By employing Green’s second identity and the radiation condition we obtain, from (114) and the wave equation satisfied by the free space Green’s function:

[TABLE]

Similarly, by employing Green’s second identity, we obtain, from (115) and the wave equation satisfied by the free space Green’s function $G^{[1]}$ (satisfying the same wave equation as $G^{[0]}$ in which $k^{[0]}$ is replaced by $k^{[1]}$):

[TABLE]

wherein $\boldsymbol{\nu}$ is the unit vector normal to $\Gamma$ pointing into $\Omega_{1}$, $d\Gamma$ the differential area element of the boundary $\Gamma$, and we have made use of the fact that

[TABLE]

with $\mathcal{H}_{\Omega}$ the domain Heaviside distribution equal to 1 for $\mathbf{x}\in\Omega$ and to 0 for $\mathbf{x}\in(\mathbf{R}^{3}-\Omega)$.

If we assume that $\mathbf{x}$ is a point on $\Gamma$, $\mathbf{x^{+}}$ a point on $\Gamma^{+}$ (the surface just outside, and homothetic to, $\Gamma$), and $\mathbf{x^{-}}$ a point on $\Gamma^{-}$ (the surface just inside, and homothetic to, $\Gamma$), all three points lying along the line associated with $\boldsymbol{\nu(\mathbf{x})}$, then

[TABLE]

wherein $\varepsilon$ is a small, positive number. It follows that

[TABLE]

whence

[TABLE]

[TABLE]

With the notations:

[TABLE]

we obtain

[TABLE]

[TABLE]

Recall that the fields satisfy the continuity conditions:

[TABLE]

whence (124) takes the form

[TABLE]

so that, with the notations

[TABLE]

we are finally confronted with the following system of two coupled boundary integral equations in the two unknown functions $U(\mathbf{x})$ and $V(\mathbf{x})$:

[TABLE]

[TABLE]

4.3.2 BI-BI method: iterative method for obtaining the solutions for $U$ and $V$

Let $U^{(j)}$ and $V^{(j)}$ be the $j$-th order approximations of the coupled BI-BI equations (129)-(130). Then, after the intialization

[TABLE]

the iterative method proceeds (in the indicated order) as follows:

[TABLE]

[TABLE]

Note that this method requires the resolution of a single BIE of the second kind (132), followed by the computation of a single integral transform (133), for each approximation order $j\geq 1$.

4.3.3 BI-BI method: the 2D case (of cylindrical symmetry)

Let the origin $O$ be located within the obstacle, whence the cylindrical coordinates $r,\theta,z$ are appropriate, with $z$ the ignorable coordinate. Then, $\mathbf{x}=(r,\theta)$, $\mathbf{x^{\prime}}=(r^{\prime},\theta^{\prime})$, the closed boundary $\Gamma$ becomes the closed curve $\gamma$ and $d\Gamma$ becomes $d\gamma$. $\gamma$ is assumed to be describable by the single-valued function

[TABLE]

and

[TABLE]

The $2\pi$-periodicity (in terms of $\theta$) of the boundary makes it legitimate to assume

[TABLE]

Furthermore, the Greens’ functions admit the representations

[TABLE]

so that the two boundary integral equations (BIE) lead (with $r^{\prime}_{\gamma}=r_{\gamma}(\theta^{\prime})$) to

[TABLE]

[TABLE]

We project these relations so as to obtain the coupled system of two matrix integral equations (for the two unknown vector functions $\{U_{m};\forall m\in\mathbb{Z};\forall\theta\in[0,2\pi[\}$, $\{V_{m};\forall m\in\mathbb{Z};\forall\theta\in[0,2\pi[\}$)

[TABLE]

[TABLE]

4.3.4 BI-BI method: the case of a homogeneous, circular cylindrical obstacle

The origin $O$ is at the center of the circular boundary of radius $a$, so that

[TABLE]

whereupon the previous two matrix integral equations become

[TABLE]

[TABLE]

We now make use of the identity

[TABLE]

so as to find

[TABLE]

[TABLE]

or simply,

[TABLE]

[TABLE]

which represent two coupled systems of linear equations for the two unknown vectors $\{U_{l}(a);\forall l\in\mathbb{Z}\}$, $\{V_{l}(a);\forall l\in\mathbb{Z}\}$. Note that:

[TABLE]

To go a step further, we must recall that

[TABLE]

so that

[TABLE]

wherein $H(\zeta>0)=1$ and $H(\zeta<0)=0$. Consequently,

[TABLE]

[TABLE]

from which it follows that;

[TABLE]

[TABLE]

or

[TABLE]

in which

[TABLE]

We now make use of the Wronskian formula for the Bessel and Hankel functions ([1], p. 360)

[TABLE]

to obtain:

[TABLE]

We showed previously that

[TABLE]

so that the solution of the system of linear equations turns out to be:

[TABLE]

The last step is to obtain explicit expressions for the fields within and outside the obstacle from $U_{l}(a)~{},~{}V_{l}(a)$. We found previously (117)-(118) that:

[TABLE]

[TABLE]

wherein $\varpi^{[1]}$ is the disk $\{r<a~{};~{}\forall\theta\in[0,2\pi[\}$ and $\varpi^{[0]}=\mathbb{R}^{2}-\varpi^{[1]}$. Because of the $2\pi$-angular periodicity of the fields, we have

[TABLE]

so that, after evaluation of $\int_{0}^{2\pi}p^{[j]}(\mathbf{x})e^{-il\theta}d\theta$, we obtain the following explicit representations of the $r-$ dependent components of the fields:

[TABLE]

[TABLE]

or, more explicitly:

[TABLE]

[TABLE]

or, even more explicitly;

[TABLE]

[TABLE]

which, on account of (159) and (162), become:

[TABLE]

[TABLE]

these being identical to the exact, reference solutions (110)-(111) obtained previously by the DD-SOV method.

4.4 The coupled domain-boundary (DI-BI) method applicable to the special case in which the mass density of the fluid-like obstacle does not depend on position

This case, which is treated in publications such as [37], differs from the most general one only by the fact that the mass density of the obstacle is a constant as a function of position, i.e.,

[TABLE]

it being implicitly assumed that

[TABLE]

4.4.1 Domain decomposition

Instead of employing the frequency domain wave equation in $\mathbb{R}^{3}$, we decompose the latter into two wave equations, one for $\Omega_{0}$ in which the total pressure is $p^{[0]}$ and the other for $\Omega_{1}$ in which the total pressure is $p^{[1]}$, these two domains being separated by $\Gamma$, through which the pressure is continuous and the velocity potential is continuous. Thus:

[TABLE]

[TABLE]

By employing Green’s second identity and the radiation condition we obtain, from (176) and the wave equation satisfied by the free space Green’s function:

[TABLE]

wherein we have made use of the fact that

[TABLE]

with $\mathcal{H}_{\Omega}(\mathbf{x})$ the domain Heaviside distribution equal to 1 for $\mathbf{x}\in\Omega$ and to 0 for $\mathbf{x}\in(\mathbf{R}^{3}-\Omega)$.

Similarly, by employing Green’s second identity, we obtain, from (177) and the wave equation satisfied by the free space Green’s function:

[TABLE]

Eq. (178) yields

[TABLE]

whereas (180) gives rise to

[TABLE]

wherefrom it ensues (by linear superposition, and on account of the position-free nature of the two mass densities and the transmission condition on pressure and velocity potential):

[TABLE]

Eq. (178) gives rise to

[TABLE]

whereas (180) yields

[TABLE]

so that their linear superposition results (also on account of the position-free nature of the two mass densities and the transmission condition on pressure and velocity potential) in

[TABLE]

To retrieve $p^{[0]}(\mathbf{x}\in\Omega_{0})$ from these two functions, we proceed as follows. Eq. (178) yields

[TABLE]

whereas (180) gives rise to

[TABLE]

so that their linear superposition results (also on account of the position-free nature of the two mass densities and the transmission condition on pressure and velocity potential) in

[TABLE]

Let:

[TABLE]

Because of the continuity of $G^{[0]}(\mathbf{x};\mathbf{x^{\prime}})$ (illustrated hereafter), the term in $\boldsymbol{\nu}\cdot\nabla_{\mathbf{x^{\prime}}}p^{[1]}(\mathbf{x^{\prime}})$ disappears in (183), so that the latter becomes:

[TABLE]

whereas (186) becomes:

[TABLE]

Eqs. (191) and (192) constitute a system of two coupled (one BI, the other DI) integral equations for the two unknown functions $U(\mathbf{x}$ and $W(\mathbf{x})$. Once these two functions are determined, $p^{[0]}$ follows by simple integration via

[TABLE]

4.4.2 DI-BI method: iterative methods for obtaining the solutions for $U$ and $W$

Let $U^{(j)}$ and $W^{(j)}$ be the $j$-th order approximations of the coupled DI-BI equations (191)-(192).

A first iterative method proceeds, after the initialization

[TABLE]

as follows (in the indicated order):

[TABLE]

[TABLE]

Note that this iterative method requires the resolution of a single domain integral equation (DIE) (195) of the second kind for $W^{(j)}(\mathbf{x})$, followed by the computation of $U^{(j)}(\mathbf{x})$, via the simple integral transform (196), for each approximation order $j\geq 1$..

A perhaps more efficient iterative method proceeds, after the initialization

[TABLE]

as follows (in the indicated order):

[TABLE]

[TABLE]

Note that this iterative method requires the resolution of a single boundary integral equation (198) of the second kind for $U^{(j)}(\mathbf{x})$ (in principle, less computationally-intensive than a DIE), followed by the computation of $W^{(j)}(\mathbf{x})$, via the simple integral transform (199), for each approximation order $j\geq 1$..

4.4.3 DI-BI method: the 2D (i.e., cylindrical symmetry) case in polar coordinates when $k^{[1]}(\mathbf{x})=k^{[1]}(r)$

Let the origin $O$ be located within the obstacle. The constitutive parameters as well as the pressure are $2\pi-$periodic functions of $\theta$ and the boundary (formerly $\Gamma$, now $\gamma$) of the obstacle is a circle of radius $a$ in the sagittal plane, so that from (183) it ensues that

[TABLE]

whereas (186) gives rise to

[TABLE]

and () yields

[TABLE]

We make the following separation-of-variables (SOV) ansatzes (suggested by the $2\pi$-periodicity of the involved functions:

[TABLE]

[TABLE]

and make use of the identity

[TABLE]

wherein $\delta_{lm}$ is the Kronecker delta symbol, to find:

[TABLE]

[TABLE]

[TABLE]

Recall that

[TABLE]

[TABLE]

so that

[TABLE]

which fact demonstrates the previously-mentioned continuity of $G_{l}^{[0]}$, and

[TABLE]

Moreover,

[TABLE]

the consequences of which are:

[TABLE]

[TABLE]

[TABLE]

wherein

[TABLE]

4.4.4 DI-BI method: the 2D (i.e., cylindrical symmetry) case in polar coordinates when $k^{[1]}(\mathbf{x})=k^{[1]}$

In the case $k^{[1]}=$const., the relations (214)-(215) become:

[TABLE]

[TABLE]

[TABLE]

Let:

[TABLE]

Then (218)-(220) become:

[TABLE]

[TABLE]

[TABLE]

We first consider in more detail the boundary integral equation (222) and the domain integral equation (223). We make the following pressure field ansatz, dictated notably by the boundedness condition of the pressure field for $r\leq a$

[TABLE]

so that we are faced with the two integrals

[TABLE]

which can be evaluated with the help of ([1], p. 484) so as to yield;

[TABLE]

whence, with the help of the Wronskian ([1], p. 360)

[TABLE]

It follows that:

[TABLE]

[TABLE]

First, consider (229) in more detail. With the help of the previously-mentioned Wronskian, it is easily shown that:

[TABLE]

whence (229) becomes

[TABLE]

Again, employing the Wronskian, we find

[TABLE]

so that

[TABLE]

On account of (225) it seems reasonable to assume

[TABLE]

whence

[TABLE]

and after recalling that $U_{l}^{i}(a)=B_{l}^{[0]}J_{l}(k^{[0]}a)$, it finally ensues that:

[TABLE]

which agrees with the DD-SOV solution (111) for $B_{l}^{[1]}$.

To see if the ansatzes (225) and (235) are admissible we return to (231) which must yield the same expression for $B_{l}^{[1]}$ as previously. Thus, under the hypotheses (225) and (235), (231) yields

[TABLE]

whence

[TABLE]

the only possible solution of which is again (237) if it is recalled that $W_{l}^{i}(r)=B_{l}^{[0]}J_{l}(k^{[0]}r)$. Thus, (225) and (235) are admissible ansatzes in that the same expression for $B_{l}^{[1]}$ is obtained via both the first and second coupled BIE+DIE.

To obtain an explicit expression for the pressure field in $\varpi_{0}$, we return to (224) in which we make use of the now-known facts that

[TABLE]

together with the expression for $K_{l}(r)$ (227), to obtain

[TABLE]

which simplifies to:

[TABLE]

whence finally:

[TABLE]

This, of course, agrees with the previously-obtained exact, reference DD-SOV (110), as well as BI-BI, solutions of this problem.

4.5 The BI-DI method

4.5.1 Governing equations

We again make the following assumptions:

[TABLE]

wherein $\rho^{[1]}$ is generally different from $\rho^{[0]}$ and

[TABLE]

and appeal to the following results found in sect. 4.4:

[TABLE]

[TABLE]

Adding these two relations, and making use of the properties of the Heaviside distributions results in:

[TABLE]

or, after employment of the transmission conditions across $\Gamma$

[TABLE]

This relation forms the basis of what we term the coupled BI-DI method. We note that (249) gives rise to two relations:

[TABLE]

[TABLE]

so that we must first solve for $p^{[1]}$ via the integral equation (250) and then solve for $p^{[0]}$ via the integral transform (251). We now proceed to do this for a circular cylindrical obstacle.

4.5.2 BI-DI method: the integral equation and integral transform for a circular cylindrical obstacle

For a circular cylindrical obstacle we proceed as in sect. 4.4 to obtain the integral equation

[TABLE]

and the integral transform

[TABLE]

Employing the previously-obtained polar coordinate representations of the Green’s function results in:

[TABLE]

[TABLE]

4.5.3 BI-DI method: the integral equation and integral transform for a circular cylindrical obstacle in which the density and wavespeed are constants with respect to the space coordinates

Now $k^{[1]}$ is a constant so that the term \Big{(}\big{(}k^{[1]}(\mathbf{x^{\prime}})\big{)}^{2}-\big{(}k^{[0]}\big{)}^{2}\Big{)} can be placed outside the integrals. Moreover, we can adopt the ansatzes

[TABLE]

[TABLE]

and recall that

[TABLE]

so that (254)-(255) become

[TABLE]

[TABLE]

from which we obtain, via (226)-(228):

[TABLE]

[TABLE]

which agrees, once again, with the exact, reference DD-SOV solutions (110)-(111) as well as with (237) and (243).

4.6 The Domain Integral (DI) method

4.6.1 Governing equations

The governing (BPDE) equation for this problem is re-written here for convenience (in simplified form with the $\omega$-dependence being implicit in $p$, $s$, $\rho$ and $c$).

[TABLE]

In addition, the field is bounded everywhere except in the support of the sources and satisfies the radiation condition.

Recall that the field radiated by the sources, henceforth termed the ”incident field” and designated by $p^{i}$, satisfies

[TABLE]

Subtracting these two equations gives

[TABLE]

wherein

[TABLE]

is the diffracted pressure field and

[TABLE]

the induced source density.

4.6.2 DI field representation

The important feature of (265) is that it is exactly of the same form as the partial differential equation relative to radiation from a bounded source of sect. 3.3. Moreover, since $p$ and $p^{i}$ satisfy the radiation condition, their difference $p^{d}$ must also satisfy the radiation condition, which fact authorizes (with the appropriate substitutions) the use of the previous result (61), i.e.

[TABLE]

Remark a

Contrary to what was stated in the remark at the end of sect. 3.3, here the domain integral does not solve the forward problem at hand, which is to determine the pressure field, given the applied sources, spatial distributions of density and wavespeed, and geometry of the configuration. The reason for this is that the term $s^{d}$ under the integral sign in (268) is a function of $p$ and $\nabla p$, both of which are unknown functions. This is seen from:

[TABLE]

wherein

[TABLE]

We shall return to this expression, which is often termed the Bergman integral representation of the pressure, further on.

Remark b

If both the wavespeed $c$ (and therefore $k$) in the first term containing $k^{2}-(k^{[0]})^{2}$ is everywhere the same, i.e., equal to $k^{[0]}$, and the mass density $\rho$ in the term containing $\nabla\rho$ is everywhere the same, i.e., equal to $\rho^{[0]}$, then both these terms vanish so that the total pressure field $p$ reduces to the incident pressure field $p^{i}$, as one would expect in the case of absence of the obstacle.

Remark c

If the mass density $\rho$ in the term containing $\nabla\rho$ is everywhere the same, i.e., equal to $\rho^{[0]}$, then this term vanishes, and, as we shall see hereafter, the determination of the pressure field within (and outside) the obstacle is straightforward via the domain integral representation

[TABLE]

which is what is often termed the Lippmann-Schwinger integral representation of the pressure field.

4.6.3 DI method: a two-step scheme deriving from the $p-q$ transformation

The until-now-considered (standard) version of the DI method seeks to obtain the pressure $p$ of the equation

[TABLE]

whereas the $p-q$ transform version of the DI method first seeks to obtain the transformed pressure $q$ of the equation

[TABLE]

wherein

[TABLE]

followed by the obtention of $p$ from $q$ from

[TABLE]

In the same manner that we obtained the standard version DI equation (269) from (273) we find the $p-q$ version DI equation

[TABLE]

wherein

[TABLE]

This shows that the $p-q$ transform has enabled the elimination of the $\nabla_{\mathbf{x^{\prime}}}p$ term in the integral but this ’advantage’ is obtained at the expense of requiring the computation of $\nabla_{\mathbf{x^{\prime}}}\cdot\nabla_{\mathbf{x^{\prime}}}\rho^{-1/2}(\mathbf{x^{\prime}})$ which may be difficult if there exists a discontinuity of $\rho$ in the medium (this actually occurring on the boundary of an obstacle in which the density is constant and different from the density in the host medium). Our opinion on this matter is that the $p-q$ method cannot be employed in cases in which the density is discontinuous across the boundary of the obstacle, and if such be the case, a smoothing strategy must be employed to render $\rho$ artificially-continuous, such as is done in [65].

Consequently, we prefer to employ the standard version of the DI method from now on.

4.6.4 DI method: an iterative scheme

With the so-called Born approximation initialization:

[TABLE]

the iterative scheme becomes

[TABLE]

Owing to the fact that this scheme is not very informative concerning the constant-density assumption issue, it will be replaced by another, more appropriate, iterative scheme in a subsequent section.

4.6.5 DI method: the 2D case of cylindrical symmetry

Let the origin $O$ be located within the obstacle, whence the cylindrical coordinates $r,\theta,z$ are appropriate, with $z$ the ignorable coordinate. Then, $\mathbf{x}=(r,\theta)$, $\mathbf{x^{\prime}}=(r^{\prime},\theta^{\prime})$, the closed boundary $\Gamma$ becomes the closed curve $\gamma$ and $d\Omega$ becomes $d\varpi$. $\gamma$ is assumed to be describable by the single-valued function

[TABLE]

and

[TABLE]

so that

[TABLE]

The $2\pi$-periodicity (in terms of $\theta$) of the boundary makes it legitimate to assume

[TABLE]

whence

[TABLE]

Henceforth, we assume that the constitutive properties of the cylindrical obstacle, depend, at worst, only on $r$ so that

[TABLE]

[TABLE]

Furthermore, the Green’s function admits the representations

[TABLE]

wherein

[TABLE]

so that the DI equation becomes, after projection

[TABLE]

4.6.6 The DI method: case in which the boundary of the ”obstacle” is a circle, $k^{[1]}=k^{[1]}(r)$ and $\rho^{[1]}$ is a constant as a function of $r$ and $\theta$

We now assume that either $\gamma$ is a circle of radius $a$ or the obstacle is enclosed within a circular cylinder (now termed ”obstacle”) whose boundary (in the sagittal plane) is the circle $r=a=\max_{\theta\in[0,2\pi[}r_{\gamma}(\theta)$ and that the medium within this enclosure cylinder is now termed $M^{[1]}(\rho^{[1]},c^{[1]})$ (actually a mixture of the media of the host and that of the original obstacle ; see sect. 4.7 for more details on this issue) whence the fact that the medium outside this ”obstacle” is $M^{[0]}(\rho^{[0]},c^{[0]})$ whose constitutive parameters do not, as before, depend on position. Finally, we make the assumption that $\rho^{[1]}$ (but not $c^{[1]})$) does not depend on position. It ensues that;

[TABLE]

[TABLE]

wherein $H(\zeta)$ is the Heaviside distribution. The consequence of the last expression is

[TABLE]

in which $\delta(\zeta)$ is the Dirac delta distribution. It follows that the DI equation takes the form:

[TABLE]

or, on account of the sifting property of the Dirac distribution, and the continuity of $\frac{p_{l,r^{\prime}}(r^{\prime})}{\rho(r^{\prime})}$ for $r=a$,

[TABLE]

4.6.7 The DI method: case in which the boundary of the ”obstacle” is a circle, and $k^{[1]}$ and $\rho^{[1]}$ are constants as a function of $r$ and $\theta$

We had

[TABLE]

With $k^{[1]}(r)=k^{[1]}=$cte., $D_{l}(r)$ becomes

[TABLE]

First consider the internal field. The boundedness condition on the field within the obstacle and the wave equation within the latter imply that it is reasonable to make the ansatzes

[TABLE]

so that

[TABLE]

whence, on account of (288),

[TABLE]

Employing (227)-(228) gives

[TABLE]

wherein

[TABLE]

It follows from (297) that

[TABLE]

so that

[TABLE]

wherein

[TABLE]

We found previously (294) that

[TABLE]

whence

[TABLE]

wherefrom

[TABLE]

or

[TABLE]

which agrees, once again, with the exact, reference solution (111)and (237).

Next, consider the external field. Again, on account of (297):

[TABLE]

which yields, after the employment of (227)-(228)

[TABLE]

Eq. (303) entails

[TABLE]

Moreover,

[TABLE]

whence

[TABLE]

or

[TABLE]

which agrees, once again, with the exact, reference DD-SOV solution (110) and (243).

4.6.8 DI method: the constant density solution for the ”obstacle”

Everything is the same as in the previous section except that we now assume that we are in the constant-density situation, i.e.,

[TABLE]

so that, assuming

[TABLE]

we find, in the same manner as previously:

[TABLE]

Recall by (306) that

[TABLE]

so that (317) is, of course, a consequence of (318) because $e_{l}^{<}=0$ when $\rho^{[1]}=\rho^{[0]}$.

The key feature of (318) is the separation, in its denominator, of the term ($e_{l}^{<}$) depending on the contrast of mass density density from the other term ($d_{l}^{<}$) which does not depend on the mass density. This crucial feature emerges clearly and naturally in the DI formulation, but not in the DD-SOV, BI-BI, DI-BI and BI-DI formulations.

4.6.9 DI method: an iterative solution for the ”obstacle” useful for small differences of density

Eqs (318) and (317) entail

[TABLE]

We have

[TABLE]

or

[TABLE]

with

[TABLE]

Note that $\epsilon$ is a dimensionless (potentially-) small parameter.

If the product $\|\epsilon f_{l}\|<1$ then it is appropriate to employ the geometric series

[TABLE]

to represent $B_{l}^{[1]}$, whereby the iteration formula, initialized by

[TABLE]

is obtained:

[TABLE]

Thus, the smaller is $\epsilon$ (i.e., the closer the two densities are to each other), the closer is $B_{l}^{[1](j)}$ to the constant-density solution $\mathcal{B}_{l}^{[1]}$, and the lesser the amount of iterates that are necessary to obtain a decent approximation to $B_{l}^{[1]}$. This will be verified numerically later on.

Let us now turn to $A_{l}^{[1]}$. We found previously from (310)-(311)that:

[TABLE]

[TABLE]

so that (312) becomes

[TABLE]

with

[TABLE]

and, in the constant-density case

[TABLE]

Eq. (319) entails

[TABLE]

so that, on account of (330),

[TABLE]

or, employing once again the geometric series representation

[TABLE]

whence the iterative scheme, with the intialization

[TABLE]

for the obtention of successive approximations of $A_{l}^{[1]}$

[TABLE]

from which we notice that the $j$th-order correction to $A_{l}^{[1]}$ is of order $\epsilon^{j}$, which, for small $\epsilon$, is the same as the $j$th-order correction to $B_{l}^{[1]}$, since the latter is also of order $\epsilon^{j}$.

4.6.10 Reconciling various publications dealing with the DI method

The publications [55, 51, 57], which apparently deal with the same problem of scattering by an obstacle having space variable-dependent constitutive properties, seem to differ from each other, and especially from our own formulation of the DI method. We show hereafter that they are all equivalent under certain prescribed conditions (actually the same as the ones assumed in our DI-BI method):

[TABLE]

wherein $\rho^{[1]}$ is generally different from $\rho^{[0]}$ and

[TABLE]

The starting point is the fundamental relation (249) of the BI-DI method, i.e.,

[TABLE]

From Green’s first identity we obtain

[TABLE]

or, after recalling the partial differential equation (177) satisfied by $p^{[1]}$, and on account of (336)-(337),

[TABLE]

whence

[TABLE]

wherein

[TABLE]

Eq. (341) agrees with (3.5)-(3.7) of [57] if we recognize the fact that Martin’s Green’s function is the negative of our Green’s function. The 2D version of (341) also agrees with (2) of [51] if we recall the correspondences of the constitutive parameters of the 2D acoustic problem with those of the 2D SH elastic wave problem (see sect. 2.10). For the same reason, the 2D version of (341) agrees with the 2D version of (12) in [55]. However, it should be stressed that the authors of [55] assume that their obstacle is homogeneous which means, in particular, that their equivalent of $\alpha$ is a constant with respect to the spatial coordinates (this being also our starting hypothesis), whereas the author of [57] seems to conclude that his result holds also for $\alpha$ dependent on $\mathbf{x^{\prime}}$.

We should mention that (341) agrees with (2) of [59] provided that it is understood that the authors of [59] assume piecewise constant (in $\mathbb{R}^{2}$) compressibility $\kappa$ instead of our assumption of piecewise constant (in $\mathbb{R}^{2}$) wavespeed $c$.

Finally, it should be pointed out that the expression of the DI for 3D elastic wave scattering given in [67] appears to derive from the Mal-Knopoff integral equation (attributed by the authors of [67] to a study of Kupradze) for the particular case of variable $\rho$ assuming constant $\lambda$ and $\mu$ in $\mathbb{R}^{3}$.

4.7 The reasons, in general 2D acoustic scattering problems, for choosing the configuration of two homogeneous, isotropic media separated by a boundary of canonical shape

A general acoustic scattering problem can take several elementary forms (this list is non-exhaustive):

(1) the obstacle is composed of several smaller, homogeneous objects of arbitrary shape separated one from the other,

(2) the obstacle is composed of only one object that is inhomogeneous and of arbitrary shape,

(3) the obstacle is composed of only one object that is homogeneous and of arbitrary shape,

the obstacle being, by definition, a subdomain of $\mathbb{R}^{n}~{};~{}n=1,2,3$, within (i.e., its interior) which the constitutive properties are different from those of the region (i.e., its exterior of infinite extent) that prevail (i.e., the constitutive properties) in the absence of the obstacle.

In 3D acoustic scattering problems, it is convenient to situate the origin somewhere near the geometric midpoint of the group of objects or within the boundary of the single object so that either the group, or single object, can be thought to be contained within a sphere (it is not necessary for this sphere to be the smallest one able to contain the group or single object. We then define $\Omega_{ext}$ and $\Omega_{int}$ as the regions exterior and interior to this sphere (henceforth termed the ’obstacle’) and $\Gamma_{ext-int}$ as the boundary of the latter. By definition, the constitutive properties of the medium in $\Omega_{ext}$ are those of the homogeneous medium filling the host region that prevails in the absence of the obstacle. Thus, the constitutive parameters of the medium in $\Omega_{ext}$ are $\rho^{[ext]}$, $c^{[ext]}$, which are both constants with respect to position in $\Omega_{ext}$. The question is: what mass density and wavespeed should we assign to the necessarily-inhomogeneous (except in case (3) when the boundary of the object is a sphere coinciding with the minimal sphere) medium occupying $\Omega_{int}$?

In 2D acoustic scattering problems (of more interest in this study), it is convenient to situate the origin somewhere near the geometric midpoint of the group of cylindrical objects (i.e., infinitely long in the $z$ direction) or within the boundary of the single object so that either the group, or single object, can be thought to be contained within a circular cylinder (it is not necessary for this circular cylinder to be the smallest one able to contain the group or single object. We then define $\varpi_{ext}$ and $\varpi_{int}$ as the regions exterior and interior to this circular cylinder (henceforth termed the ’obstacle’) in its cross section plane and $\gamma_{ext-int}$ as the boundary of the circular disk $\varpi_{int}$. By definition, the constitutive properties of the medium in $\varpi_{ext}$ are those of the homogeneous medium filling the host region that prevails in the absence of the obstacle. Thus, the constitutive parameters of the medium in $\varpi_{ext}$ are $\rho^{[ext]}$, $c^{[ext]}$, which are both constants with respect to position in $\varpi_{ext}$. The question is again: what mass density and wavespeed should we assign to the necessarily-inhomogeneous (except in case (3) when the boundary of the object is a circular cylinder coinciding with the minimal circular cylinder) medium occupying $\varpi_{int}$?

Consider the circular cylinder in its sagittal plane, i.e., the disk $\varpi_{int}$, within which there are two regions, one ($\varpi_{int}-\varpi_{1}$) in which the constitutive parameters $\mathfrak{c}_{j}~{};~{}j=1,2$ are $\mathfrak{c}_{j}^{[0]}$ (constants with respect to $\mathbf{x}$) and the other ($\varpi_{1}$ in which the constitutive parameters are $\mathfrak{c}_{j}^{[1]}(\mathbf{x})$ (i.e., generally position-dependent). A general mixing formula for the parameters $\mathfrak{c}_{j}$ is

[TABLE]

To simplify matters, we shall consider $\mathfrak{c}_{j}^{[1]}~{};~{}j=1,2$ to also be constants with respect to position which amounts to the supposition that the composite medium within the disk is binary, i.e., composed of areas in which the physical properties $\mathfrak{c}^{[0]}$ are constants and other areas in which the physical properties $\mathfrak{c}^{[1]}$ are also constants, different from $\mathfrak{c}^{[0]}$. Consequently,

[TABLE]

wherein $\phi=\frac{\varpi_{1}}{\varpi_{int}}$ is the filling fraction, or the ratio of the area $\varpi_{1}$ occupied by the material having parameters $\mathfrak{c}^{[1]}~{};~{}j=1,2$ to the total area in the disk $\varpi_{int}$.

Until now, we have been rather vague as to the physical nature of the parameters $\mathfrak{c}^{[1]}~{};~{}j=1,2$. The most natural choice is:

[TABLE]

whence

[TABLE]

Another choice is:

[TABLE]

whence

[TABLE]

Eqs. (345)-(346) and (347)-(348) are the simplest (there exist more complex) versions of the so-called series and parallel (respectively) mixing formulae [82] for binary obstacles.

Now let us summarize what we have done. We began with a group of objects or single object of arbitrary shape, each object having properties $\mathfrak{c}_{j}^{[1]}$ different from those, $\mathfrak{c}_{j}^{[0]}$ (constants as a function of position), of the host. For the obvious reason, made evident in sect. 4, that it is much simpler to treat a scattering problem involving a single, homogeneous obstacle of canonical shape than one involving several (or one) heterogeneous or homogeneous objects of arbitrary shape, we decided, aided by this possibility in the DI method, to enclose the original group of objects or single object of arbitrary shape by a virtual boundary of canonical shape (spherical in 3D, circular cylindrical in 2D). We termed the domain interior to this virtual boundary to be that of ’the obstacle’. This obstacle is generally-heterogeneous because it comprises regions having the physical properties of the previous objects and another region having the physical properties of the host. We saw in sect. 4, that the integral methods turn out to be much more easily exploitable if this obstacle, enclosed by the virtual boundary of canonical shape, is homogeneous. Finally, we proposed a method for homogenizing the obstacle which appeals to a simple mixture formula. Thus, the original problem reduces to that of a (single) homogeneous obstacle enclosed within a virtual boundary of canonical shape, this problem having been treated in depth in various subsections of sect. 4.

It can be objected that the problem we finally address is not the one we proposed to solve at the outset. This objection is founded, but we found no other way to treat the original problem if the goal is to do this by other-than-purely-numerical means. In fact, this was precisely our goal as regards our quest for answers regarding the constant mass-density assumption question.

A second remark has to do with the ’closeness’ of the solution of the scattering problem of the ’new’ obstacle (i.e., the homogeneous one with canonical-shaped virtual boundary) to that of the problem dealing with the ’old’ obstacle (i.e., the one with several objects or one object with arbitrary shape). Homogenization, and, in particular, its mixture formula version, has been shown to give rise to very useful (in both the inverse and forward problem contexts [68, 99, 100]), and often precise, approximations of the solution, especially at low frequencies. This fact constitutes another, important, justification of the procedure adopted in our study.

5 Numerical results

5.1 The numerical procedure and choice of parameters

The aim of the computations was to compare numerically the acoustic scattering amplitudes $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$, based on the constant-density assumption (i.e., $\epsilon=0$ in sect. 4.6.9): (a) to the exact (i.e., DD-SOV) $A_{l}^{[0]},~{}B_{l}^{[1]}$ expressions (110)-(111))respectively of the amplitudes, as well as (b) to the approximations $A_{l}^{[0](j)},~{}B_{l}^{[1](j)}~{};~{}j=1,2,3$ of the amplitudes, all of which do not rely on the constant-density assumption. $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ are simply (110)-(111), in which $\rho^{[1]}$ is taken to be equal to $\rho^{[0]}$. $A_{l}^{[0](j)},~{}B_{l}^{[1](j)}~{};~{}j=1,2$ are given in (324)-(325) and (334)-(335) and $\epsilon$ is the mass density contrast parameter defined in (322).

The graphs are relative to the external amplitudes $\mathcal{A}_{l}^{[0]},~{}A_{l}^{[0]},~{}A_{l}^{[0](j)}$ and internal amplitudes $\mathcal{B}_{l}^{[1]},~{}B_{l}^{[1]},~{}B_{l}^{[1](j)}$ as a function of: frequency $f$, $\rho^{[1]}$ and $c^{[1]}$, the other parameters being fixed at the following values: $\rho^{[0]}=1000~{}Kgm^{-3}$, $c^{[0]}=1500~{}ms^{-1}$, $a=0.1~{}m$, $B_{l}^{[0]}=1$. Usually, the scattering amplitude order $l$ is taken to range from 0 to 3, but at high frequencies, we also give results for much larger $l$.

5.2 Variation of $\rho^{[1]}$ for various orders $l$: case $\rho^{[1]}>\rho^{[0]}$ and low $f$

We observe substantial qualitative differences, that increase with the mass density contrast, between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$ respectively except for very small differences of $\rho^{[1]}$ from $\rho^{[0]}$. Better agreement is attained, both quantitatively and qualitatively, by means of $A_{l}^{[0](1)},~{}B_{l}^{[1](1)}$, and near-coincidence with $A_{l}^{[0]},~{}B_{l}^{[1]}$ is obtained by means of $A_{l}^{[0](2)},~{}B_{l}^{[1](2)}$ respectively.

5.3 Variation of $\rho^{[1]}$ for various orders $l$: case $\rho^{[1]}<\rho^{[0]}$ and low $f$

The same remarks apply to this series of figures as in the series of sect. 5.2.

5.4 Variation of $\rho^{[1]}$ for various orders $l$: case $\rho^{[1]}>\rho^{[0]}$ and higher $f$

The same remarks apply to this series of figures as in the series of sect. 5.2.

5.5 Variation of $\rho^{[1]}$ for various orders $l$: case $\rho^{[1]}<\rho^{[0]}$ and higher $f$

The same remarks apply to this series of figures as in the series of sect. 5.2.

5.6 Variation of (low-to-medium) frequency $f$ for various orders $l$: case of fairly-large ($\epsilon=0.231$) mass density contrast and large \big{|}\Im(c^{[1]})\big{|}

We observe, between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$ respectively, substantial qualitative, and even quantitative, differences, that are particularly pronounced near the extrema of the amplitudes. Better agreement is attained, both quantitatively and qualitatively, by means of $A_{l}^{[0](1)},~{}B_{l}^{[1](1)}$, and near-coincidence with $A_{l}^{[0]},~{}B_{l}^{[1]}$ is obtained by means of $A_{l}^{[0](2)},~{}B_{l}^{[1](2)}$ respectively.

5.7 Variation of frequency $f$ for various orders $l$: case of relatively-small ($\epsilon=0.091$) density contrast and large \big{|}\Im(c^{[1]})\big{|}

As expected the differences between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$ respectively are less-pronounced than in sect. 5.6. Otherwise, the same remarks apply here than in the previous section.

5.8 Variation of (low-to-medium) frequency $f$ for various orders $l$: case of relatively-small ($\epsilon=0.091$) small density contrast and small \big{|}\Im(c^{[1]})\big{|}

Due to the lesser attenuation in the obstacle, the differences between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$ respectively are observed to be more-pronounced than in sect. 5.7. Otherwise, the same remarks apply here than in the previous section.

5.9 Variation of low-to-medium frequency $f$ for various orders $l$: case of fairly-large ($\varepsilon=0.231$) mass density contrast and large wavespeed contrast for $\Re(c^{[1]})>c^{[0]}$

We observe, between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$ respectively, substantial qualitative, and even quantitative, differences, that are particularly pronounced near the extrema of the amplitudes. Better agreement is attained, both quantitatively and qualitatively, by means of $A_{l}^{[0](1)},~{}B_{l}^{[1](1)}$, and near-coincidence with $A_{l}^{[0]},~{}B_{l}^{[1]}$ is obtained by means of $A_{l}^{[0](2)},~{}B_{l}^{[1](2)}$ respectively.

5.10 Variation of low-to-medium frequency $f$ for various orders $l$: case of fairly-large ($\varepsilon=0.231$) mass density contrast and small wavespeed contrast for $\Re(c^{[1]})>c^{[0]}$

We observe, between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$ the same qualitative, and somewhat aggravated quantitative, differences as in sect. 5.9. Otherwise, the same comments as in the previous section apply here.

5.11 Variation of low-to-medium frequency $f$ for various orders $l$: case of fairly-large ($\varepsilon=0.231$) mass density contrast and small wavespeed velocity contrast for $\Re(c^{[1]})<c^{[0]}$

We observe, between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$ the aggravated qualitative, and similar quantitative, differences as in sect. 5.10. Otherwise, the same comments as in the previous section apply here.

5.12 Variation of high frequency $f$ for various orders $l$: case of fairly-large ($\varepsilon=0.231$) mass density contrast, small velocity contrast and large \big{|}\Im(c^{[1]})\big{|}

We observe, between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$, rather large qualitative, and quantitative, differences. These differences are reduced between the exact amplitudes and $A_{l}^{[0](1)},~{}B_{l}^{[1](1)}$, and even more so, between the exact amplitudes and $A_{l}^{[0](2)},~{}B_{l}^{[1](2)}$.

A remarkable aspect of these figures is the evidence in the internal amplitudes of the occurrence of a resonance, followed by an asymptotic regime in which these amplitudes gently tend to zero. The resonances set in at $f=\sim 20,\sim 30,\sim 50~{}KHz$ for $l=5,11,20$ respectively. These resonances are associated with the excitation of normal modes (as is testified by the fact that the denominator of $B_{l}^{[1]}$ vanishes or is nearly nil at the resonance frequencies), as explained in [96]. The asymptotic regime seems to set in near $80~{}KHz$ for all $l$. Beyond this frequency, the external field is a quasi-periodic function of $f$, with slowly-decreasing envelope.

5.13 Variation of frequency near the resonance frequency of the $l=20$ normal mode, for various $\Im(c^{[1]})$

We again observe, between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$, rather large qualitative, and quantitative, differences. These differences are reduced between the exact amplitudes and $A_{l}^{[0](1)},~{}B_{l}^{[1](1)}$, and even more so, between the exact amplitudes and $A_{l}^{[0](2)},~{}B_{l}^{[1](2)}$.

To give further support to the normal mode resonance hypothesis, we progressively diminished \big{|}\Im(c^{[1]})\big{|}. This resulted, as expected, in the increase of the height of the resonance peak. This peak is not infinitely-high for $\Im(c^{[1]})=0$ because, even in the absence of material losses, there are pseudo-losses due to radiation damping (a part of the incident energy is always scattered to outer space (i.e., lost)) even when the obstacle acts like a closed cavity.

5.14 Variation of frequency $f$ near the former $l=20$ normal mode resonance frequency for various $\Re(c^{[1]})$

We observe qualitative and quantitative differences between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$. These differences are reduced between the exact amplitudes and $A_{l}^{[0](1)},~{}B_{l}^{[1](1)}$, and even more so, between the exact amplitudes and $A_{l}^{[0](2)},~{}B_{l}^{[1](2)}$.

The increase of the velocity contrast is seen to progressively eliminate what was previously (when $c^{[1]}=1400~{}ms^{-1}$) identified as being the $l=20$ normal mode resonance near $f=50~{}KHz$. This does not necessarily mean that the resonances cease to exist althogehter, since they might show up for other $l$ at other frequencies when the velocity contrast is large (but we have not verified this).

5.15 Variation of $\Re(c^{[1]})$ for various $l$: case of $\epsilon=0.231$ density contrast at $f=20~{}KHz$

We observe fairly-substantial qualitative and quantitative differences between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$. These differences seem to increase somewhat with the velocity contrast (but less than with mass density contrast observed previously) and, as previously, are reduced between the exact amplitudes and $A_{l}^{[0](1)},~{}B_{l}^{[1](1)}$, and even more so, between the exact amplitudes and $A_{l}^{[0](2)},~{}B_{l}^{[1](2)}$.

5.16 Variation of $\Re(c^{[1]})$ for various $l$: case of $\epsilon=-0.429$ density contrast at $f=20~{}KHz$

We observe very-substantial qualitative and quantitative differences between $\mathcal{A}_{l}^{[0]},~{}\mathcal{B}_{l}^{[1]}$ and $A_{l}^{[0]},~{}B_{l}^{[1]}$. These differences seem to increase somewhat with the velocity contrast (but less than with mass density contrast observed previously) and, as previously, are reduced between the exact amplitudes and $A_{l}^{[0](1)},~{}B_{l}^{[1](1)}$, and even more so, between the exact amplitudes and $A_{l}^{[0](2)},~{}B_{l}^{[1](2)}$.

6 Conclusions

The approach, based on brute-force resolution of the partial differential equations (as they stand) of wave scattering, is certainly useful in contemporary studies and applications dealing with both forward and inverse scattering problems. The principal drawback of this approach is that it is purely-numerical, and thus gives little insight as to what parameter (or group of parameters), contribute the most to the different phenomena that intervene in wave scattering. A convincing example of this is the phenomenon of normal mode resonances which may be hard to detect in the numerical results if the frequency is not varied with very small steps, and certainly harder to explain from purely-numerical evidence. Another, practical, drawback of the numerical PDE approach is that it is computationally-intensive, even in 2D, scalar versions employed to get a first grasp on what are fundamentally 3D and/or vectorial problems, the latter requiring enormous computer resources out of reach of most research groups. This negative aspect of the brute-force PDE approach is particularly acute in inverse scattering problems because the tendency is more and more to solve them by the nonlinear full-wave inversion (NLFWI) method which appeals to optimization algorithms that require the resolution of (usually) many forward problems.

An alternative to the brute-force numerical resolution of the PDE is the integral equation (IE) approach. Actually, as we have shown herein, there exist many variants of this approach, the most appealing of which appears to be the one relying on domain integral (DI) representation of the scattered field because it is the easiest to employ when the obstacle is composed of several objects (thus is discontinuously-heterogeneous) or is continuously-heterogeneous. The other IE methods are interesting too because they involve less computations when the numerical stage is attained, but they require, at the outset, that the objects composing the obstacle be homogeneous with respect to at least one of the constitutive parameters (usually the mass density in acoustic and elastodynamic problems, permeability in electromagnetic problems). In fact, when this assumption of homogeneity is made, it also greatly simplifies the DI formulation, a fact that was recognized very early by those concerned with biophysical, marine acoustical, and geophysical applications. The principal reason for having adopted the DI approach in these applications was, and still is, that it leads, via the Born approximation and asymptotic forms of the Green’s function, to what is called diffraction tomography (DT), which constitutes an explicit, functional rather than numerical (although, somewhat distorted) form of solution of the inverse problem. Of course, the IE solutions (not necessarily of the DI variety) can also be incorporated in NLFWI schemes to speed up the numerical resolution of inverse problems.

Other than affording the advantages of simplicity and computational economy, one suspects that the constant mass density assumption (which is not usually invoked in the computational PDE approach) must have some drawbacks when employed in conjuction with an IE technique. The most obvious one is that, in the inverse problem context, it is not possible, by this means, to obtain a quantitative, realistic retrieval of the mass density of an obstacle such as a female breast [35, 101] simply because this density is taken to be that of water. or of the surrounding tissue. The second drawback of the constant-density assumption is that it is based on blind faith, especially in the geophysical context.

Consequently, the question we addressed herein was: how reliable is the constant-density assumption, particularly in the 2D acoustical and geophysical contexts in which the governing PDE is the Bergman PDE (equivalent to the 2D SH Mal-Knopoff PDE) that one might want to solve by an IE method?. As shown herein, the IE method can take a variety of forms: BI-BI, DI-BI, BI-DI and DI-DI, the latter of which exists in several versions which we found to be equivalent. All these IE methods, which have been employed to numerically solve a variety of both inverse and forward scattering problems, do not enable to answer the aforementioned question because they do not yield tractable mathematical (rather than purely-numerical) solutions for 2D obstacles of general shape and composition. Moreover, all the IE techniques except the DI-DI one, rely on the initial assumption that the density of the obstacle be spatially-invariant, although different from that of the host medium (this being less-radical than assuming the density of both the obstacle and host to be spatially-invariant and equal, i.e., the constant density assumption). Thus, to compare the answers furnished by the four forms of the IE technique, we were obliged, at the outset, to make the assumption of a constant-density obstacle. But this was not sufficient, because none of the four IE techniques enable the obtention of a mathematically-tractable solution for a macroscopically-homogeneous (this term is used either for a medium containing no macroscopic objects or for a porous medium such as a solid network of connected microscopically-thin rods within a fluid for which a homogenization scheme enables its constitutive properties to be approximated [64, 100] by those of an ’effective’ homogeneous fluid or solid) obstacle of general shape. Consequently, we imposed the additional assumption that the shape of the obstacle be canonical, this meaning that the governing PDE is separable in a given coordinate system and that the boundary of the obstacle can be represented by an equation $\mathfrak{k}=$const., where $\mathfrak{k}$ is one of the coordinates in the chosen coordinate system.

In fact, we chose the cylindrical coordinate system $r,\theta,z$ so that the obstacle is a cylinder [64] whose boundary is circular, i.e, $r=a=$const for all $\theta\in[0,2\pi$. For this type of macroscopically-homogeneous obstacle, the Bergman PDE can be solved [64, 97] in closed form via domain decomposition and separation of variables (DD-SOV) so that, at worst, one thus disposes of a reference solution to which can be compared the solutions furnished by the various IE methods. We say ’at worst’ because it is not evident that these IE methods can yield closed-form solutions even for such a homogeneous canonical-shaped obstacle so that the comparisons might have been no more than numerical in nature. It turns out (this was the result of a great deal of what is offered in this study) that all the IE schemes enable to re-generate the exact DD-SOV sollution which gives some legitimacy to these schemes, but still does not answer the question of the reliability of the constant-density supposition.

It could be thought that since we dispose of an exact solution, for the two constant-density media separated by a circular boundary, in convenient mathematical form, why not just mathematically compare this solution to its counterpart obtained by taking the two densities to be equal, or what is more interesting, nearly equal? Actually, there did not appear to exist an easy way to do this (particularly for the nearly-equal density case), but it was found that the DI-DI method, contrary to the other three IE methods, enables one to generate the solution in a form in which the density contrast $\epsilon$ emerges in a clearly-identifiable manner. From this type of solution, we were able to show that the amplitudes of the external and internal cylindrical waves comprising the scattered field each take the form of a series in powers of $\epsilon$ and that the first term of each such series is none other than the corresponding (external or internal wave) amplitude resulting from the constant-density assumption.

This series form of the solution for the non-constant density problem thus constitutes the mathematical tool, supplemented by numerical examples extracted therefrom, by which we were able to answer the question of the reliability of the constant-density assumption, at least as concerns the model of a macroscopically-homogeneous obstacle having a circular boundary. Finally, we gave arguments that justify the employment of this model to yield useful information concerning the constant-density assumption even when the obstacle is neither macroscopically-homogeneous (but can be homogenized) nor enclosed within a physical circular boundary (but can be enclosed within a virtual circular boundary, a procedure that gives rise to no difficulties in the DI-DI method).