Space-time domain solutions of the wave equation by a non-singular boundary integral method and Fourier transform

TL;DR

This paper presents a high-precision, stable boundary integral method combined with Fourier transform to efficiently solve the space-time wave equation for complex scattering problems, accurately capturing multiple scattering effects.

Contribution

It introduces a non-singular boundary integral approach for Helmholtz problems and applies Fourier transform for efficient space-time wave solutions, improving stability and accuracy.

Findings

Exact enforcement of radiation boundary conditions

High numerical stability for small degrees of freedom

Efficient response analysis to various incident pulses

Abstract

The general space-time evolution of the scattering of an incident acoustic plane wave pulse by an arbitrary configuration of targets is treated by employing a recently developed non-singular boundary integral method to solve the Helmholtz equation in the frequency domain from which the fast Fourier transform is used to obtain the full space-time solution of the wave equation. The non-singular boundary integral solution can enforce the radiation boundary condition at infinity exactly and can account for multiple scattering effects at all spacings between scatterers without adverse effects on the numerical precision. More generally, the absence of singular kernels in the non-singular integral equation confers high numerical stability and precision for smaller numbers of degrees of freedom. The use of fast Fourier transform to obtain the time dependence is not constrained to discrete time steps and is particularly efficient for studying the response to different incident pulses by the same configuration of scatterers. The precision that can be attained using a smaller number of Fourier components is also quantified.

Figures

Tables

Table 1. Table 1: The percentage root mean squared (RMS) relative errors at all positions on the scattered wave between the scattered wave calculated by using the indicated number of Fourier amplitudes and by using all Fourier amplitudes. The results are for the observation points at ( 0 , 0 , − 4 ​ a ) 0 0 4 𝑎 (0,0,-4a) and ( 0 , 0 , − 1.1 ​ a ) 0 0 1.1 𝑎 (0,0,-1.1a) . The ratio ϵ ≡ | P ~ ​ ( ω m ) | / | P ~ m ​ a ​ x | italic-ϵ ~ 𝑃 subscript 𝜔 𝑚 subscript ~ 𝑃 𝑚 𝑎 𝑥 \epsilon\equiv|\tilde{P}(\omega_{m})|/|\tilde{P}_{max}| denotes the magnitude of the Fourier amplitudes relative to the maximum amplitude.

Number of		% RMS error at	% RMS error at
Fourier amplitudes	$ϵ$	$(0,0,-1.1a)$	$(0,0,-4a)$
5	$>20\%$	4.4	4.9
10	$>3.5\%$	1.6	1.7
20	$>0.9\%$	0.5	0.6

Equations

$$\rho \frac{\partial \mathbf{u}}{\partial t}+\rho \mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p,$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0,$$

$${\rho }_0\frac{\partial \mathbf{u}}{\partial t}=-\nabla p,$$

$$\frac{\partial {\rho }_1}{\partial t}+{\rho }_0\nabla \cdot \mathbf{u}=0.$$

$$p={\left(\frac{\partial p}{\partial \rho }\right)}_{{\rho }_0}{\rho }_1\equiv c^2{\rho }_1.$$

$$\left[{\nabla }^2-\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}\right]p(\mathbf{x},t)=0.$$

$$\frac{\partial f(\mathbf{x},t)}{\partial x}\to \frac{f^{(n)}(i+1,j,k)-f^{(n)}(i-1,j,k)}{2\Delta x}$$

$$\frac{\partial f(\mathbf{x},t)}{\partial t}\to \frac{f^{(n+1)}(i,j,k)-f^{(n-1)}(i,j,k)}{2\Delta t}.$$

$$c\Delta t\le {\left[\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}+\frac{1}{\Delta z^2}\right]}^{-1/2}.$$

$$c_0({\mathbf{x}}_0)p({\mathbf{x}}_0,t_0)={\int }_{-\infty }^{\infty }dt{\int }_{S}dS(\mathbf{x})\left[\frac{\partial p(\mathbf{x},t-t_0)}{\partial n}\mathcal{G}-p(\mathbf{x},t-t_0)\frac{\partial \mathcal{G}}{\partial n}\right].$$

$$\mathcal{G}(\mathbf{x}-{\mathbf{x}}_0,t-t_0)=\frac{\delta (t-t_{ret})}{|\mathbf{x}-{\mathbf{x}}_0|},t_{ret}\equiv t_0-|\mathbf{x}-{\mathbf{x}}_0|/c$$

$$\left[{\nabla }^2-\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}\right]\mathcal{G}(\mathbf{x}-{\mathbf{x}}_0,t-t_0)=-4\pi \delta (\mathbf{x}-{\mathbf{x}}_0)\delta (t-t_0).$$

$$h(\mathbf{x},t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }H(\mathbf{x},\omega )\exp (-i\omega t)d\omega .$$

$${\nabla }^{2}P(\mathbf{x},\omega )+k^{2}P(\mathbf{x},\omega )=0,k^2\equiv {\omega }^2/c^2.$$

$${\int }_{S+S_{\infty }}\left[P(\mathbf{x})-P({\mathbf{x}}_0)g(\mathbf{x})-{\left(\frac{\partial P}{\partial n}\right)}_{0}f(\mathbf{x})\right]\frac{\partial G({\mathbf{x}}_0,\mathbf{x})}{\partial n}dS(\mathbf{x})$$

$$={\int }_{S+S_{\infty }}\left[\frac{\partial P(\mathbf{x})}{\partial n}-P({\mathbf{x}}_0)\nabla g(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})-{\left(\frac{\partial P}{\partial n}\right)}_0\nabla f(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})\right]G({\mathbf{x}}_0,\mathbf{x})dS(\mathbf{x}).$$

$$G(\mathbf{x},{\mathbf{x}}_0)=\frac{\exp (ik|\mathbf{x}-{\mathbf{x}}_0|)}{|\mathbf{x}-{\mathbf{x}}_0|}$$

$${\nabla }^{2}f(\mathbf{x})+k^{2}f(\mathbf{x})=0,f({\mathbf{x}}_0)=0,\nabla f({\mathbf{x}}_0)\cdot {\mathbf{n}}_0({\mathbf{x}}_0)=1$$

$${\nabla }^{2}g(\mathbf{x})+k^{2}g(\mathbf{x})=0,g({\mathbf{x}}_0)=1,\nabla g({\mathbf{x}}_0)\cdot {\mathbf{n}}_0({\mathbf{x}}_0)=0.$$

$$4$$

$$-{\int }_{S+S_{\infty }}\left[P(\mathbf{x})-P({\mathbf{x}}_0)g(\mathbf{x})-{\left(\frac{\partial P}{\partial n}\right)}_{0}f(\mathbf{x})\right]\left[\frac{\partial G({\mathbf{x}}_p,\mathbf{x})}{\partial n}-\frac{\partial G({\mathbf{x}}_0,\mathbf{x})}{\partial n}\right]dS(\mathbf{x})$$

$$+{\int }_{S+S_{\infty }}\left[\frac{\partial P(\mathbf{x})}{\partial n}-P({\mathbf{x}}_0)\nabla g(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})-{\left(\frac{\partial P}{\partial n}\right)}_0\nabla f(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})\right]\left[G({\mathbf{x}}_p,\mathbf{x})-G({\mathbf{x}}_0,\mathbf{x})\right]dS(\mathbf{x})$$

$$\displaystyle p^{inc}({\bf x},t)\equiv\tilde{p}(\tau)=\left\{\begin{array}[]{ccc}\qquad 0,\qquad\qquad\qquad\quad-4\pi N_{c}<\tau<-2\pi N_{c}\\ \sin(\tau)\exp(-\alpha|\tau|),\qquad-2\pi N_{c}<\tau<2\pi N_{c}\\ \qquad 0,\qquad\qquad\qquad\qquad 2\pi N_{c}<\tau<4\pi N_{c}\\ \end{array}\right.$$

$$\xi /a$$

$$\eta /a$$

$$f(\mathbf{x},\omega )$$

$$g(\mathbf{x},\omega )$$

$$I_1\equiv {\int }_{S_{\infty }}\left[\frac{\partial g(\mathbf{x})}{\partial n}G({\mathbf{x}}_0,\mathbf{x})-\frac{\partial G({\mathbf{x}}_0,\mathbf{x})}{\partial n}g(\mathbf{x})\right]dS=-4\pi ,$$

$$I_2\equiv {\int }_{S_{\infty }}\left[\frac{\partial f(\mathbf{x})}{\partial n}G({\mathbf{x}}_0,\mathbf{x})-\frac{\partial G({\mathbf{x}}_0,\mathbf{x})}{\partial n}f(\mathbf{x})\right]dS=0.$$

$$I_{11}=\frac{i}{2}{\int }_{S_{\infty }}\left[e^{ikr\cos \theta }-e^{-ikr\cos \theta }\right]\frac{e^{ikr}}{r}k\cos \theta dS,$$

$$I_{12}=-\frac{1}{2}{\int }_{S_{\infty }}\left[e^{ikr\cos \theta }+e^{-ikr\cos \theta }\right]\frac{e^{ikr}}{r^2}(ikr-1)dS.$$

$$I_1=2\pi \frac{ik}{2}{\int }_0^{\pi }\left[e^{ikr\cos \theta }(\cos \theta -1)+e^{-ikr\cos \theta }(-\cos \theta -1)\right]\frac{e^{ikr}}{r}{r}^2\sin \theta d\theta .$$

$$I_1=2\pi (ik){\int }_{-1}^1{e}^{ikru}(u-1)e^{ikr}rdu,$$

Full text

††preprint: JASA/123

I Introduction

The space-time solution of the scalar wave equation underpins the prediction of the scattering of acoustic waves by discrete targets and is relevant to applications that range from noise suppression to seismic exploration to ultrasonic therapy. There are a number of complementary approaches to finding general numerical solutions to the problem. One of these is based on a direct solution of the wave equation by replacing derivatives in the spatial and time variables by finite differences Wang (1966); Yee (1966), referred to as the finite difference time domain approach. This method is also used extensively in the study of electromagnetic scattering in the time domain Taflove (1988). Another approach to solving the wave equation is to extend the conventional boundary integral method to the time domain using the time-dependent Green’s function to represent the spatial solution in terms of values of the wave function on the boundaries of scatterers and the time evolution is treated by time marching Groenenboom (1983). Recently there is renewed theoretical interest in the stability of the time dependent solutions of the wave equation at large times particularly for the canonical problem of scattering by a sphere in an infinite spatial domain for which space and time variations can be represented analytically in terms of infinite series of spherical harmonics and Bessel functions with time-dependent coefficients Greengard et al. (2014); Martin 2016a ; Martin 2016b .

In this paper, a recently developed non-singular boundary integral formulation of the solution of the wave equation in the frequency domain is used as the basis of constructing the solution in the time domain by Fourier transform. In this non-singular formulation, the usual singularities of the surface integrals have been eliminated analytically Klaseboer et al. (2012); Sun et al. (2015). This confers a number of advantages in that high accuracy in the evaluation of the surface integrals can be achieved with simple quadrature and with fewer surface nodes. In addition, the algorithm remains stable even when the wave number is very close to the resonant values Sun et al. (2015). The non-singular nature of the integrals means that field values near boundaries can be evaluated directly without the need for further steps to avoid numerical divergences that are characteristic of the traditional boundary integral method.

Since the solution at any time value can be found directly by Fourier transform, it is not subjected to error accumulation effects of time marching methods. Furthermore, this approach is particularly suited for exploring the effects of different incident pulsed waves on a fixed configuration of scatterers. This is because once the boundary integral solution that depends on the scatterers is found, the scattering by different types of incident pulses can be found directly by a process of linear computational complexity.

Before giving details of our approach, it is instructive to review the characteristics of the existing finite difference time domain method and time marching solution of the conventional boundary integral method of finding the space-time solution of the wave equation, as well as to touch on recent theoretical studies of the time dependence of series solutions of the wave equation – this is done in Section II. Our non-singular boundary integral formulation combined with the Fourier transform method is introduced in Section III. Results for the scattering of a plane wave pulse by different targets are given in Section IV. The paper closes with a discussion of possible future directions of our approach to obtaining space-time domain solutions of the wave equation in other applications.

II Overview of existing methods

The equations that govern acoustic wave propagation are obtained by combining the Euler momentum equation

[TABLE]

with the continuity equation

[TABLE]

that relate the position, ${\bf x}$ and time, $t$ dependent density, $\rho$, velocity, ${\bf u}$ and pressure, $p$ Wang (1966).

For small amplitude oscillations, the density is written as $\rho\equiv\rho_{0}+\rho_{1}$ and only terms linear in $p$, ${\bf u}$ and the small density deviation, $|\rho_{1}|\ll\rho_{0}$, from the constant mean density, $\rho_{0}$, are retained in Eqs. 1 and 2 to give a pair of first order partial differential equations:

[TABLE]

[TABLE]

These equations can be closed by introducing the material constitutive equation characterized by the speed of sound, $c$

[TABLE]

Eliminating ${\bf u}$ and $\rho_{1}$ from Eqs. 3, 4 and 5, gives the wave equation for the pressure, $p$

[TABLE]

II.1 Finite difference time domain

The finite difference time domain solution of the acoustic wave equation follows the method that was developed for solving the propagation of Maxwell’s electromagnetic equations Yee (1966). Instead of working directly with the wave equation, Eq. 6, the equivalent system of first order equations, Eq. 3 to 5, is solved at discrete values: $f({\bf x},t)\rightarrow f^{(n)}(i,j,k)$ with indices for the Cartesian spatial nodes ${\bf x}=(x,y,z)\rightarrow(i,j,k)$ and time steps $t\rightarrow n$. Partial derivatives are approximated by first order central differences Wang (1966):

[TABLE]

with similar approximations for partial derivatives with respect to the other spatial coordinates $y$, $z$ and to time, $t$

[TABLE]

So given boundary and initial conditions, the space-time solution is obtained by time marching using Eqs. 7 and 8 to solve Eqs. 3 to 5.

Although the finite difference time domain algorithm is conceptually straightforward, there are a number of technical issues that require careful implementation, see Taflove (1988) for details. For instance, to ensure convergence, the step sizes in time, $\Delta t$ and space, $\Delta x,\Delta y,\Delta z$ are constrained by the condition:

[TABLE]

Numerical dispersion effects associated with the relative orientation of the spatial grid and the direction of propagation can arise. If there are changes in the spatial grid density within the solution domain, necessitated for example by differences in characteristic length scales of the problem, care needs to be exercised to avoid unphysical reflections at the boundary between regions of the different grid densities. Physical boundaries between different media are assumed to conform to the stepwise nature of the grid. If the problem domain is infinite, then an ‘outer’ boundary needs to be constructed with boundary conditions that will satisfy the Sommerfeld radiation boundary condition at infinity Sommerfeld (1912) so as to avoid unphysical reflections back into the solution domain.

II.2 Time marching with conventional boundary integral methods

The boundary integral equation formulation of the solution of the wave equation avoids the task of solving the wave equation in a 3D spatial domain due to the elliptic nature of the problem. Instead, it is only necessary to determine values of the function and its normal derivative on the boundary surfaces, $S$ that enclose the problem domain. This reduces the dimension of the problem by one. The solution of the wave equation, Eq. 6 at the space time point $({\bf x}_{0},t_{0})$ can be expressed as a surface integral over the surface, $S$, involving the function, $p$ and its normal derivative $\partial p/\partial n\equiv\nabla p\cdot{\bf n}$, where the surface normal ${\bf n}({\bf x})$ points out of the solution domain Groenenboom (1983)

[TABLE]

The constant $c_{0}({\bf x}_{0})$ is the solid angle subtended at ${\bf x}_{0}$: $c_{0}({\bf x}_{0})=4\pi$, if ${\bf x}_{0}$ lies within the solution domain and $c_{0}({\bf x}_{0})=0$, if ${\bf x}_{0}$ lies outside the domain. If ${\bf x}_{0}$ is on the boundary surface $S$, the solid angle subtended at the boundary, $c_{0}({\bf x}_{0})$, will depend on the details of the local geometry of the boundary that is relevant, as for instance, in numerical implementations where the boundary is represented by piecewise continuous surface elements.

The kernel of this integral equation is the time-retarded Green’s function

[TABLE]

that satisfies

[TABLE]

Given initial conditions and boundary data for $p$ (Dirichlet problems) or $\partial p/\partial n$ (Neumann problems), Eq. 10 can be solved by a time marching method. The divergence of the Green’s function ${\cal{G}}$ at ${\bf x}={\bf x}_{0}$ means that Eq. 10 contains singularities that require careful treatment in numerical evaluations of the surface integrals. The most common approach is to simply represent the surfaces $S$ by a number of planar elements and assume the function $p$ is constant within each element. The use of higher order surface elements to represent the surface with more accuracy introduces further complexities because of the presence of the singularities in the kernel.

Even though the physical acoustic problem is well-behaved on domain boundaries that are the surfaces of scatterers, the singular behavior of the Green’s function, ${\cal{G}}$ that originates from the mathematical formulation of Eq. 10, means that the precision to which function values can be computed near boundaries or in problems in which surfaces are close together may be compromised, and will require additional effort to resolve.

II.3 Stability of series solutions

Both the finite difference time domain method discussed in Section II.1 and the conventional boundary integral formulation summarized in Section II.2 use time marching to track the time evolution. This process accumulates global numerical error as time progresses. Recently, there has been detailed theoretical analysis of the stability of the solution at large times for the solution of the wave equation outside a single sphere in an infinite domain. In this simpler problem, the spatial variation can be represented analytically as infinite series comprised of spherical harmonics and Bessel functions. The time dependent coefficients of such expansions turn out to grow exponentially with increasing order. As a consequence, a loss of significant figures will result from cancellations between terms of growing magnitude at large times unless new formulations are used to calculate the coefficients Greengard et al. (2014); Martin 2016a ; Martin 2016b .

It is therefore attractive to be able to retain the lower spatial dimensionality of the boundary integral approach in combination with a different way to treat the time evolution that does not accumulate global error associated with time marching. In the next section, a recently developed non-singular boundary integral formulation will be used that eliminates all the singular behavior that arises from the Green’s function. A fast Fourier transform of the frequency domain solution is used to circumvent the error accumulation characteristics of time marching solutions. Further efficiencies can be gained by focusing only on frequency components in the dominant part of the power spectrum of the incident wave.

III Non-singular boundary integral Fourier transform method

The twin objectives of the present non-singular boundary integral Fourier transform method to obtain space-time solutions of the wave equation are to retain the lower spatial dimensionality feature of the boundary integral method while avoiding the use of time marching to track the time evolution. The boundary integral approach automatically satisfies the Sommerfeld radiation at infinity exactly. Since there is no need to represent the 3D domain by a grid, it is not necessary to be concerned with numerical dispersion issues associated with grid based methods. In the absence of the need to have a fixed spatial grid, there is flexibility to accommodate special characteristics in the shapes of the scatterers.

A recently developed non-singular boundary integral method Klaseboer et al. (2012); Sun et al. (2015) to solve the Helmholtz equation in the frequency domain is employed. With this method, the usual singularities associated with the conventional boundary formulation are eliminated analytically. This makes it easy to use quadratic surface elements to represent the geometric features of boundaries more faithfully and to do so with a smaller number of degrees of freedom to minimize the problem size. Rather than simply assuming a constant function value for each surface element, quadratic interpolants are used to represent functional variations within each element since there are no singular integrals to complicate such an approach. Consequently, surface integrals can be evaluated efficiently with simple quadrature. Also, the absence of singularities in the kernel means that function values on or near the boundaries can be calculated without restrictions or possible loss of precision.

The following convention is used to define the Fourier representation of a function, $h({\bf x},t)$ in space and time in terms of its Fourier transform $H({\bf x},\omega)$:

[TABLE]

In the frequency domain, the wave equation for the pressure in Eq. 6 becomes

[TABLE]

The non-singular boundary integral equation for $P({\bf x},\omega)\equiv P({\bf x})$ is Klaseboer et al. (2012); Sun et al. (2015)

[TABLE]

where dependence on $\omega=kc$ in all functions is suppressed to ease the notation. The Green’s function, $G$ is given by

[TABLE]

and the functions $f({\bf x})$ and $g({\bf x})$ can be any convenient solution of the equations

[TABLE]

[TABLE]

with ${\bf n}({\bf x}_{0})$ being the outward normal at ${\bf x}_{0}$.

The integrals in Eq. III are taken over the surfaces, $S$ of scatterers and over the surface at infinity, $S_{\infty}$ that together enclose the 3D solution domain. If $f({\bf x})$ and $g({\bf x})$ obey Eqs. 17 and 18, the terms containing $f({\bf x})$ and $g({\bf x})$ and their gradients in Eq. III will cancel the singular behavior of $G$ and $\partial G/\partial n$ at ${\bf x}={\bf x}_{0}$. The choices adopted here for $f({\bf x})$ and $g({\bf x})$ are given in Appendix A. In addition, the integrals over the surface at infinity, $S_{\infty}$ in Eq. III can also be evaluated analytically, see Appendix B. As the surface integrals in Eq. III are not singular, they can be evaluated efficiently and accurately using quadrature, without the need to interpret them as principal value integrals.

It is also worthy to note that the solid angle, $c({\bf x}_{0})$ that appears in the conventional boundary integral formulation, see Eq. 10, has been eliminated in our non-singular formulation, Eq. III. This is advantageous in practical numerical implementations because there is no longer the need to be concerned with calculating the solid angle at ${\bf x}_{0}$ that depends on the details of the local surface geometry of the surface elements.

A numerically robust way to evaluate the pressure $P({\bf x}_{p},\omega)$ at a point ${\bf x}_{p}$ in the solution domain that may be arbitrarily close to a boundary is to use the following expression (again for brevity, the dependence on $\omega=kc$ is suppressed in all functions)

[TABLE]

where ${\bf x}_{0}$ is a point on the surface that is closest to ${\bf x}_{p}$. Again all integrals in Eq. III are free of singularities and well-behaved even as ${\bf x}_{p}\rightarrow{\bf x}_{0}$, and therefore, can be evaluated using the standard Gauss quadrature. Proofs of these results are given in Sun et al. (2015).

Having found the pressure, $P({\bf x}_{p},\omega)$ in the frequency domain, the space-time solution, $p({\bf x}_{p},t)$ can be found by taking the inverse Fourier transform, Eq. 13 using the discrete fast Fourier transform method Cooley et al. (1965).

IV Results - scattering of a plane wave pulse

We consider the scattering of an incident wave that comprises an infinite periodic train of plane wave pulses that replicates a fundamental waveform. Consider an example of the fundamental wave pulse with the form, where $\tau\equiv k_{0}(z-ct)\equiv(2\pi/\lambda_{0})(z-ct)$

[TABLE]

and travels in the $z$ direction. This pulse, $\tilde{p}(\tau)$ has $2N_{c}$ oscillatory cycles modulated by the constant $0<\alpha<1$ and pre- and post-padded by zero amplitudes to make a total non-dimensional width, $\tilde{w}=8\pi N_{c}$, as shown in Fig. 1. This fundamental wave pulse is then replicated to create an infinite periodic train as illustrated in the inset of Fig. 1b. By sampling the incident wave at $N_{f}$ evenly spaced points in the interval $-4\pi N_{c}<\tau<4\pi N_{c}$ then gives a set of discrete Fourier components, $P^{inc}({\bf x},\omega)\equiv\tilde{P}(\omega)$. In Fig. 1a, such an incident wave, $\tilde{p}(\tau)$, is shown with $2N_{c}=4$ oscillatory cycles, each of wavelength $\lambda_{0}$ sampled at $N_{f}=128$ values. Since $\tilde{p}(\tau)$ is an odd function, its Fourier transform $\tilde{P}(\omega)$ is an even function so there are only $N_{f}/2=64$ unique values of the amplitude $|\tilde{P}(\omega_{m})|$ as shown in Fig. 1b. The aliasing properties of the discrete Fourier representation will produce the infinite wave train shown in the inset of Fig. 1b.

With the sampling rate shown in Fig. 1, one would, in general, solve the wave equation, Eq. 14, in the frequency domain at each of the 64 values of $k_{m}=\omega_{m}/c$ using the non-singular boundary integral method given by Eq. III. Then the inverse discrete fast Fourier transform of these solutions will give the total pressure in the space-time domain that is a sum of the incident and scattered components: $p({\bf x},t)=p^{inc}({\bf x},t)+p^{scat}({\bf x},t)$. As the Fourier spectrum $|\tilde{P}(\omega_{m})|$ shown in Fig. 1b is dominated by only a few values of $\omega_{m}$, then as shall be seen, accurate results can be constructed from the boundary integral solutions obtained only at the $k_{m}=\omega_{m}/c$ values that correspond to say just the 10 largest $|\tilde{P}(\omega_{m})|$ amplitudes.

In the non-singular boundary integral equation, Eq. III, for the pressure in the frequency domain, we need to specify the functions $f({\bf x})$ and $g({\bf x})$ that satisfy Eqs. 17 and 18 that will also determine the value of the integral over the surface at infinity, $S_{\infty}$ in Eq. III. These technical details are given in the Appendices.

IV.1 Spherical scatterers

Numerical results are given for the space-time domain solution of the scattering of the plane wave pulse given in Fig. 1 by a sphere of radius, $a$. The width of the fundamental pulse, $w$ in Fig. 1a is taken to be $w=20.1a$ so that the parameters in Eq. 23 are $k_{0}a=2\pi a/\lambda_{0}=16\pi/20.1$ and the modulating constant is taken to be $\alpha=0.1$. The choice $w=20.1a$ is to avoid having the wave number being too close to the resonant value for the sphere. Two types of boundary conditions are considered: a ‘soft’ sphere that corresponds to the boundary conditions: $p=0$ and a ‘hard’ sphere that is specified by $\partial p/\partial n=0$ on the surface.

The results are presented with $t=0$ (and the first time step) defined as the moment the leading edge of the oscillatory part of the pressure pulse first come into contact with the sphere as illustrated in Fig. 2. By the 78th time step, the oscillatory part of the incident wave has completely passed over the sphere and so no further scattering will take place thereafter. However, until the 128th time step, that is also the last time step in the fundamental period, the scattered wave will continue to travel away from the sphere. In all video files in the supplementary material, the animation starts at time zero, as described above, and continues for 128 time steps.

The fast Fourier transform of the incident pulse was taken with 128 sampling points and since the pulse is an odd function, this gives 64 distinct frequencies: $k_{m}a=(2\pi/20.1)m$, $m=0,1,2...63$, at which the non-singular boundary integral equation, Eq. III, has to be solved. To account for the above time convention, the incident wave must be adjusted by a multiplicative phase factor $e^{i\beta}$, with $\beta=-k_{m}a[w/(4a)-1]$ in the boundary integral equation.

As in previous work Sun et al. (2015), the sphere surface was represented by 500 quadratic elements and 1002 nodes. Quadratic interpolation was used to represent the variation of the function within each element to construct a linear system from Eq. III that was solved using a direct method. Although this is less efficient, it does have the advantage that if the form of the incident pulse is changed, there is no need to solve the boundary integrals again, provided that the wave width $w$ remains the same.

For ‘soft’ sphere: $p=0$ or ‘hard’ sphere: $\partial p/\partial n=0$ boundary conditions on the surface, the non-singular boundary integral equation was solved for the scattered field, $p^{scat}$ that obeyed the Sommerfeld radiation condition at infinity, and on the sphere surface, $p^{scat}$ was given in terms of the incident field, $p^{inc}$: $p^{scat}=-p^{inc}$ for the soft sphere or $(\partial p^{scat}/\partial n)=-(\partial p^{inc}/\partial n)$ for the hard sphere. The time sequence of wave amplitudes in Figs. 3 and 4 shows the space-time variation of the scattered and total wave as the incident pulse within a $20a\times 20a$ square in the $xy$-plane as the incident wave traverses the sphere that is located at the origin. In view of the symmetry of the problem, both the scattered pressure wave and the total pressure wave in the $yz$-plane can be displayed within the same figure. Videos of animations of these results are available in the electronic supplement (Video 1 and 2). Although the boundary condition for the scattered wave is not spherically symmetric, the scattered wavefronts become spherical as they travel away from the sphere.

Our results obtained from quadratic elements were compared with those obtained using 774 linear elements and 389 nodes with linear interpolation for the variation of the function on each element. Results from using quadratic or linear elements agree to better than 2 significant figures. Since an analytical solution of the wave equation, Eq. 14 as a series expansion in terms of spherical harmonics and Bessel functions is available, see for example Doinikov (1994), it can also be ascertained that the present non-singular boundary integral solution with the stated numbers of elements and nodes are correct to better than 2 significant figures.

In Fig. 5, the scattered wave is shown as a function of time at a position 4 radii from the sphere center: ${\bf x}=(0,0,-4a)$ or at just one-tenth of a radius from the sphere surface at ${\bf x}=(0,0,-1.1a)$ obtained by using just 5, 10 or 20 terms of the largest Fourier amplitudes $|\tilde{P}(\omega_{m})|$ out of the 64 amplitudes to construct the time behavior. These results are consistent with the familiar notions that with Fourier representations, quite acceptable answers can be obtained using just the dominant frequencies of the incident pulse to construct the space-time solution of the wave equation.

The results in Fig. 6 provide a simple illustration of the space-time variation of scattering and interference in the presence of two spherical scatterers. The identical spheres of radius $a$, with ‘hard’ boundary conditions $\partial p/\partial n=0$ are placed at a distance $3a$ between centers along the direction of the incoming pulse. The pulse has $2N_{c}=10$ cycles and a total width $w=20.1a$ that corresponds to $k_{0}a=2\pi a/\lambda_{0}=40\pi/20.1$ for the incident pulse. The snap shots of the scattered and total waves shown in Fig. 6 illustrate the appearance of temporal and positional dependence of constructive and destructive localized interference that can occur in between the spheres, Fig. 6b, or on the far surface of the downstream sphere, Fig. 6c. The amplification effect can be up to 80% larger than the maximum amplitude of the incoming pulse. The entire animated sequence of the scattering is available in the online supplementary material (Video 3).

IV.2 Wave focusing

To illustrate the transient wave focusing effect of a ‘hard’ axisymmetric acoustic bowl whose surface with coordinates $(\xi,\eta)$ are constructed by taking the closed curve defined by the parametric equation over $0\leq\theta<2\pi$:

[TABLE]

and rotating the curve about the $\eta$-axis. This axis of symmetry is then oriented at an angle of $0.15\pi$ radian relative to the direction of propagation of the incident pulse. The pulse with width $w=20.1a$ is similar to that in Fig. 1 but has $\alpha=0.001$ and $2N_{c}=6$ oscillatory cycles so that $k_{0}a=(2\pi/\lambda_{0})a=24\pi/20.1$ in Eq. 23.

The space-time variation of the focusing and directional effects are illustrated in Fig. 7 with an oblique view in Fig. 8. Animations of these results are available in the online supplementary material (Video 4 and 5). The focusing effect of this ‘hard’ reflecting bowl creates time-varying high amplitude pressure hotspots close to the bowl surface that can amplifiy the maximum amplitude of the incident pulse by about 4.5 times. Relative to spherical scatterers in Figs. 3, 4 and 6, the higher directional and intensity effects in the far field due to focusing is evident well after the incident pulse has passed the scatterer.

V Conclusion

This paper demonstrates an approach to finding the space-time dependent solution of the scalar wave equation in the context of acoustic scattering and propagation that does not involve time-marching. The scattering of incident pulses of finite spatial extent and time duration were used as illustrative examples. The method builds on the recent development of a non-singular boundary integral formulation of the solution of the wave equation in the frequency domain, the Helmholtz equation Sun et al. (2015) and then uses fast Fourier transform Cooley et al. (1965) of the result to obtain answers in the time domain. This is fundamentally different to the time-marching method based on a fully finite difference representation of the wave equation in the space and time domains or the time-marching solution based on the conventional boundary integral formulation using time-retarded Green’s functions.

Although a boundary integral based approach may appear to require many solutions of the Helmholtz equation in the frequency domain, it was shown that quite accurate results could be obtained by only using those frequencies that have relatively large amplitudes in the power spectrum of the incident pulse. By exercising judicious choice in selecting the physically important frequencies, considerable savings in computational effort can be achieved. Furthermore, the use of the non-singular boundary integral also means that higher precision can be obtained with fewer degrees of freedom. Also the concern with numerical handling of singularities that invariably accompany traditional boundary integral formulations is eliminated.

This Fourier transform approach to obtain the time domain solution also affords the flexibility to obtaining solutions at any required time point and therefore is not subjected to the space-time stability constraints, such as that in Eq. 9, on the time step size in time-marching methods.

For a fixed configuration of scatterers, the solution matrix of the boundary integral equation can be stored so that exploring the effects of varying the incident wave can be carried efficiently without having to solve the boundary integral equation again. In contrast, with time-marching methods, changing the indident wave will require an ab initio solution.

The present approach therefore provides a viable alternative to the established time-marching methods of solving acoustic problems in the space-time domain given its different characteristics. With recent reformulation of the Maxwell’s equations for electromagnetic scattering in terms of coupled scalar wave equations for the Cartesian components of the electric field, ${\bf E}$ and the scalar function, $({\bf x}\cdot{\bf E})$ Klaseboer et al. (2017); Sun et al. (2017), this approach has broader applications for finding time dependent solutions beyond the field of acoustics.

Acknowledgements.

This research was supported in part by an Australian Research Council Discovery Project Grant to DYCC.

Appendix A The functions $f({\bf x},\omega)$ and $g({\bf x},\omega)$

The functions $f({\bf x},\omega)$ and $g({\bf x},\omega)$ needed to ensure that the boundary integral equation, Eq. III is non-singular can be any functions that satisfy Eq. 17 and 18. The choice, with ${\bf k}=k\;{\bf n}({\bf x}_{0})=(\omega/c)\;{\bf n}({\bf x}_{0})$,

[TABLE]

is one possibility. It is easy to verify that this will ensure the absence of singularities in the integrands in Eq. III.

Appendix B The integral over the surface at infinity, $S_{\infty}$

The integral over the surface at infinity, $S_{\infty}$ in Eq. III has two separate contributions from the functions $f({\bf x})$ and $g({\bf x})$. This Appendix gives the derivation of the key results:

[TABLE]

[TABLE]

The integral over $S_{\infty}$ gives rise to a term $4\pi P(\mathbf{x}_{0})$ on the left hand side of Eq. III. Although this resembles the term with the solid angle in the conventional boundary integral method (Eq. 10), it is of a totally different origin. Note also that if different $f({\bf x})$ and $g({\bf x})$ functions are chosen, the $S_{\infty}$ integrals will be different, see for example Sun et al. (2015).

The results in Eqs. 27 and 28 can be derived as follows. Without loss of generality, $S_{\infty}$ can be taken as the surface of a sphere with radius, $r=|\mathbf{r}|=|\mathbf{x}-\mathbf{x}_{0}|$ and centered at $\mathbf{x}_{0}$, in the limit as $r\rightarrow\infty$. The outward surface normal is $\mathbf{n}=\mathbf{r}/r$, so that $\mathbf{k}\cdot\mathbf{n}=k\cos\theta$ and $\mathbf{k}\cdot(\mathbf{x}-\mathbf{x}_{0})=kr\cos\theta$. By using Eqs. 16, 25 and 26 and their derivatives with respect to the unit normal vector $\mathbf{n}$, one can write $I_{1}=I_{11}+I_{12}$ where

[TABLE]

[TABLE]

In the limit $r\rightarrow\infty$, $(ikr-1)\rightarrow ikr$, so using the explicit form $dS=2\pi r^{2}\sin\theta d\theta$ in $I_{11}$ and $I_{12}$ give

[TABLE]

By using $u=\cos\theta$ and $du=-\sin\theta d\theta$, one can write

[TABLE]

that can be simplified by an integration by parts as

[TABLE]

where the second term vanishes as $r\rightarrow\infty$, giving the analytical result

[TABLE]

Following similar steps, the $I_{2}$ integral is

[TABLE]

By using the same change of variable, $u=\cos\theta$, one finds

[TABLE]