Generalization and stabilization of exact scattering solutions for spherical symmetric scatterers

In a previous work the authors described a fast high-fidelity computer model for acoustic scattering from multi-layered elastic spheres. This work is now extended with a scaling strategy significantly mitigating the problem of overflow and thus expanding the useful frequency range of the model. Moreover, new boundary conditions and loads are implemented. Most important are the fluid-fluid and solid-solid couplings, which allow a completely general layering of the scattering object. Sound hard and sound soft boundary conditions are implemented for solids and fluids respectively. In addition to the existing acoustic excitation, mechanical excitation in the form of point-excitation and surface excitation are implemented. Attenuation in the form of hysteresis damping as well as viscous fluid layers are also included. Several numerical examples are included, with the purpose of validating the code against existing reference solutions. The examples include air bubbles, a coated steel shell and a point-exited steel shell.


Introduction
Computer models of acoustic scattering from three-dimensional elastic objects including internal structure are important tools in underwater acoustics, as well as in many other fields such as non-destructive testing, seismology, noise control and medical imaging.
In underwater acoustics such scattering models allow the study of realistic scattering scenarios, both from marine life [1,2] and man-made objects [3,4,5,6,7]. Present work studies scattering by objects in the presence of interfaces, which enable modeling of partially and completely buried objects. Being able to predict the acoustic response of objects aids in the development of methods for detection and identification, and reduces the need for expensive experiments [8,9,10,11,12]. Another application is noise and vibration reduction, where sandwich structures with sound absorbent cores are used to improve the acoustic insulation [13,14,15].
Many computational techniques have been developed over the years for studying scattering from fluidloaded objects, ranging from fast approximate methods to computer intensive high-fidelity models. In the development of such methods, reference models are an invaluable tool for validating the accuracy and valid ranges of the different methods [16,17].
Exact solutions are only known for a few simple shapes: the infinitely long cylinder and the sphere. There is therefore a need for reference models for more complex objects of simple shapes included inside layered acoustic media; for this purpose, we have developed an open sourced 1 toolbox for scattering by multi-layered elastic spheres [18]. The model computes scattering by an elastic object submerged in an infinite fluid : An outer spherical aluminum shell with a coating layer is immersed in water (unbounded domain), the interior being flooded with water and containing another steel shell also being flooded with water but has also an air bubble in its center. This illustrates some of the generalization being made; both fluid-fluid and solid-solid couplings are made for the multi layered model. medium, waveguide effects are not considered. The model is exact in the sense that it is based on analytical expressions but requires numerical calculations to compute a truncated series with accuracy depending only on the precision used.
The present work solves some limitations of the model described in [18] and extends its capabilities. New boundary conditions are implemented, including fluid-fluid and solid-solid coupling (e.g. as in Figure 1); this removes the former restriction on the layering to be alternating fluid and solid, such that completely general layering is now allowed. Sound hard and sound soft boundary conditions are also implemented, allowing perfectly rigid and pressure release domains to be modelled. In order to simulate realistic materials attenuation needs to be included, and attenuation in the form of viscous fluid layers and hysteresis damping are implemented for this purpose. Source types now include both acoustical (plane wave and point source) and mechanical (point force and surface) excitation.
Another issue that is addressed is the overflow problem that occurs at high ka-numbers due to the exponential behaviour of Bessel functions for large order. The resulting instabilities may hide significant physics and render the solution useless. A scaling strategy is applied that mitigates this problem and expands the useful frequency range significantly.
The novelty of this work is thus the stabilization and generalization of the work in [18] both in terms of domain composition but also type of medium and boundary conditions. Besides some small alternation to the notations and the inclusion of scaling functions (for stabilization) the new scheme is similar to that of [18]. The main achievement in this work is arguably the stabilization as this scaling is completely absent in the literature. The main motivation for this work (i.e. new boundary conditions) is to enable verification of numerical codes although the physical problem in and of itself is of interest.
The article is organized as follows: Section 2 introduces full Navier-Navier coupling conditions allowing completely general layering and implements viscous fluid layers through the linearized Navier-Stokes equation. The problem of overflow is treated in Section 3, where a scaling strategy is implemented to mitigate the exponential behavior. The resulting stable solutions are given in Section 4. Several numerical examples are given in Section 5: Scattering by gas bubbles in a liquid; empty spherical steel shell with soft coating; spherical steel shell with point excitation; and viscoelastic sphere. The examples serve to validate the implemented code by a comparison with existing benchmark solutions.

General solution
The restriction in [18], that the scattering object needs to consist of alternating fluid and solid domains is removed. This forces us to alter the notation slightly. Denote by M the number of media (solids and fluids) present such that we have M − 1 spherical interfaces at radii R m , m = 1, . . . , M , separating the solid and fluid layers (where R m > R m+1 ). If the innermost media has support at the origin, R M = 0. The domain number m is denoted by with the convention R 0 → ∞ for unbounded domains. We will also generalize the fluid domain to be viscous using linearized equations of continuity and the Navier-Stokes equation for a barotropic, viscous non-heat-conducting compressible fluid [19,20,21] ∂ρ m ∂t wherep m represents the deviation of pressure from its mean value,v m is the velocity,ρ m is the time-varying part (as opposed to the constant equilibrium density ρ m ) of the total mass density,ρ tot,m = ρ m +ρ m , and µ m and µ b,m are the shear and the bulk coefficient of viscosity, respectively. Finally, c m is the ideal speed of sound evaluated at ambient conditions. Equations (2) to (4) can be combined into a single equation for the velocity field The time-dependent velocity field,v m , can be written in terms of the time-dependent displacement field,ȗ m , byv The displacement field,ȗ m , in the solid domain is governed by (as in [18]) where the bulk modulus, K m , and the shear modulus, G m , are defined by the Young's modulus, E m , and Poisson's ratio, ν m , as Equations (5) and (7) in the frequency domain are obtained through Fourier transform 2 as 2 For a time dependent functionΨ (x, t), its Fourier transform is [18] Ψ where ω = 2πf is the angular frequency, f the frequency and i = √ −1 is the imaginary unit. That is, we use the breve notation, ·, for the time-dependent Fourier transforms of functions in the frequency domain. and respectively. Equation (9) can be written in terms of the fluid displacement using the Fourier transformation of Equation (6) v m = −iωu m (11) such that For fluid domains we define which makes Equation (12) equivalent to Equation (10). This enables us to use same expressions for fluid and solid domain when expressing the general solution.
The solution for the displacement field in domain m, may be decomposed as 3 where the displacement potentials φ m and ψ m solves respectively. The angular wave numbers a m and b m are given by For non-viscous fluids we use the classical notation for the wave number, k m = ω cm = a m , as c m = c s,1,m in this case.

Remark 1. The angular wave numbers may for fluid domain be written in terms of the fluid parameters as
The following relations are obtained from Equations (11), (14) and (15) 3 We use displacement-based potentials also for the fluid as was done in [22] (in contrast with the velocity based potentials used in [21,20]). 4 Note that [19,20] use the (linear Taylor expansion) approximation am ≈ ω cm 1 + i ωµm The displacement is given in spherical coordinates as u m = u r,m e r + u ϑ,m e ϑ with components As in [18] Equations (15) and (16) can be written in spherical coordinates as Using separation of variables, each of these equations can be reduced to a couple of spherical Bessel and Legendre equations, with the associate Legendre polynomials of zero and first order and spherical Bessel functions as solutions (as described in [18]). The solution can be written as 5 , where m,n are coefficients to be found and 6 The derived functions (e.g. the displacement field u r,m ) can be found in Appendix A. Hysteresis damping [21, p. 146] can be included by modifying Youngs modulus with an addition of an imaginary loss factor for solid domains. However, [21] assumes the loss factor to be the same in both the longitudinal and transverse wave velocities. We thus generalize by whereη 1,m andη 2,m are the loss factor in longitudinal and transverse wave velocities, respectively. For non-viscous fluids this modification is equivalent with but many authors [21,23] do not use the square root for attenuation in non-viscous fluids (i.e., c m → For unbounded domains Ω 1 we impose the Sommerfeld radiation condition [24] Note that the Einstein's summation convention will be used throughout this work. Moreover, as in [18] the spherical coordinate system (r, ϑ, ϕ) is used throughout this work. 6 Here, jn(x) and yn(x) are the spherical Bessel functions of first and second kind, respectively, h n (x) is the Hankel function of first kind, and Pn(x) are the Legendre polynomials.
as r → ∞ uniformly inx = x r . We now introduce prescribed (e.g. Appendix B) incident potential fields φ inc,m and ψ ϕ,inc,m that solves the same equations and has the same derived formulas for the displacement fields, u r,inc,m , and stress fields, σ rr,inc,m , as φ m and ψ m , respectively (except, possibly, for the Sommerfeld radiation condition in Equation (33)). For the derived quantities, such as u r,m , we simply add "inc" to the subscript (i.e. u r,inc,m ). For all applications investigated herein, ψ ϕ,inc,m = 0. The coupling conditions (Neumann-to-Neumann, NNBC) for spherical symmetric objects are given by σ rr,m + σ rr,inc,m − (σ rr,m+1 + σ rr,inc,m+1 ) = 0, pressure boundary condition, If R M = 0, the following boundary conditions may be implemented on the innermost interface with z M being the impedance [25,26]. The impedance condition is used only for non-viscous fluids.
As the Hankel functions of the first kind satisfies the Sommerfeld condition in Equation (33), these will be the radial basis functions for the outermost domain (that is, A arising from evaluating the boundary conditions Equations (35) to (41) at all interfaces (cf. [18]). All solutions are given in terms of spherical Bessel functions which have exponential behavior for large order (as illustrated in [18]). This is a problem as one computationally run into "0 · ∞" evaluations (overflow). As we will later see, we can scale the coefficients to include this exponential behavior to mitigate this problem. To find proper scales, a discussion of asymptotic behavior of the Bessel functions is in order.
m Domain Media Free parameters Non-zero coefficients 1 and To eliminate exponential behavior in this case, the scaling e −| Im z| can be used for the Bessel functions of first and second kind, and e −iz for the Hankel function of the first kind.
We consider now the uniform asymptotic expansions for Bessel functions of large order [28, 9.3.38, p. 368]. Define z(ζ) through 2 3 Then the uniform asymptotic expansion of the Bessel function of the first kind is given by [28, 9.3.35-9.3.37, and finally, starting with u 0 = v 0 = 1, [29] should use separate expansions for when the order is roughly the same size as the magnitude of the argument. Proper treatment of the transition regions of the asymptotic expansion of the Bessel functions requires usage of another expansion as outlined in [29].

Remark 2. "Transition regions"
The implementation of the Airy functions in [30] provides the option of using the scales exp( 2 3 y 3/2 ) and exp(−| 2 3 Re y 3/2 |) for the Airy functions Ai(y) and Bi(y), respectively, in order to eliminate the exponential behavior (same scales are used for the first derivatives). This motivated scales s (i) n (x), i = 1, 2, 3, for the Bessel functions (j n (y) and y n (y)) and the Hankel function of first kind (h (1) n (y)), respectively, where with ν = n + 1/2.
The plots are visually indistinguishable validating the correct scales s (i) n (z) for the Bessel functions Z (i) n (z). The usage of Equations (46) to (48) will be limited to the inside of the "cup" 7 in these plots as n must be larger than the magnitude of z for these formula to apply. However, the scales motivated by the same formulas works remarkably well in other parts of the domain as well. Further evidence for this can be obtained using high precision toolboxes like the multiprecision computing toolbox Advanpix. A test was performed (using 100 digits precision) in the same domain used in Figure 3 which resulted in the bounds 2.718 · 10 −6 |Z (i) n (z)/s (i) n (z)| 6.657 · 10 −2 , i = 1, 2, 3, for (n, z) ∈ D using 41 × 41 × 41 uniformly sampled points (for reference: max i,n,z |Z (i) n (z)| ≈ 9.341 · 10 29 375 , i = 1, 2, 3, (n, z) ∈ D). The "0/0" evaluations are here discarded.

Computationally stable solutions
In the following we transform the expressions in [18] to mitigate computational instability based on the discussion in previous section. The scales are always applied on the final formulas in [18] (that is, after differentiation and manipulations). Define the intermediate radii bỹ we can write the solution as where 7 Defined by

Remark 4.
For large n and |x − y| even the combined evaluations (Remark 3) will yield overflow, and the strategy presented herein only extends the range of parameters at which the solution may be computationally evaluated (for a given floating point precision).

Remark 5.
For a fixed x one can show that g (1) n (x) ∼ x 2n , g (2) n (x) ∼ 2n |x| and g (3) n (x) ∼ 2n x , as n → ∞, such that g (i) n should not be considered to have exponential behavior in this context.

The linear system of equation in Equation (42) is now modified to bẽ
whereH n andC n are modifications of H n and C n taking into account the scaling functions w n (x, y) and s (i) n (x), respectively.
In Figure 4 the setup used in figure 6 in [18] (with the S135 model with SHBC on the inner sphere and parameters in Table 1) is used to illustrate the improvement of this scaling. The cut-off at which machine epsilon precision is no longer obtained 8 is increased from k 1 R 1 ≈ 107 to k 1 R 1 ≈ 614. This limit is extended by roughly a factor 10 if quadruple-precision is used (which increases REALMAX). At the end however, exponential functions must still be evaluated which still limits computational evaluations for high enough frequencies. To extend the computational frequency range one needs a representation with more bits for the exponent than the standard 8-bit representation for the exponent for double precision. In Figure 4d quadruple precision is used and covers the frequency range for most applications. 8 By this we mean when the ratio between term Nε and its sum in absolute value is less than machine epsilon precision.      Fig. 1 and Fig. 4 in [31].
Parameter Description

Numerical examples
To give further evidence for the correctness of the implemented code, comparison to existing benchmark solutions by Sage [31], Skelton [21], Hetmaniuk [32], and Ayres [25], will be presented.
It is customary to present results in the far-field. For the scattered pressure p 1 , it is defined by withx = x/|x|. From the far-field pattern, the target strength, TS, can be computed. It is defined by where P inc is the amplitude of the incident wave at the geometric center of the scatterer (i.e., the origin). Unless stated otherwise, the direction of the incident wave will be given by the vector d s = [0, 0, 1] .

Sage benchmark problem
Sage [31] considers an air bubble and uses the parameters in Table 2. Sage defines the cross-section area by  Moreover, the spectrum has been sampled closely, revealing small (less significant) eigenmodes not shown by Sage. It seems as if the y-axis in [31] is scaled by 1/π, and so this is also done here. A good agreement is here observed.

Skelton benchmark problem
Skelton [21] considers a steel spherical shell with a hysteretic loss factor, scattering an incident plane wave directed along the positive x 3 -axis, and uses the parameters in Table 3. In Figure 6 several cases are compared with the data from Figure 10.4 to Figure 10.7 in [21] 10 . Good agreements are observed in all cases considered. The case of an empty steel shell without coating is considered for higher frequencies in Figure 7a and higher still in Figure 7b. The latter plots illustrate the importance of the scaling strategy presented herein as the solution without scaling (computed from code based on [22]) does not converge to machine epsilon precision due to overflow computations. In [21] the equally spaced dips in the range f ∈ [2 kHz, 10 kHz] and the "hump" round f = 14 kHz is explained to stem from first symmetric Lamb wave and first antisymmetric Lamb wave, respectively. However, both of these are actually symmetric Lamb waves, and the phenomena round f = 140 kHz (not covered in [21]) represents antisymmetric Lamb waves.

Hetmaniuk benchmark problem
Hetmaniuk [32] considers a steel spherical shell with the parameters in Table 4

Ayres benchmark problem
Ayres [25] considers a rubber spherical shell with the parameters in Table 5. The loss factors for the   Figure 12 and Figure 17 in [32].   Figure 12 and Figure 17 in [25].

Parameter Description
Parameter Description       rubber sphere is here defined byη (65) Three pairs of α and β was considered and the resulting scaled far-field pattern is plotted in Figure 9a and Figure 9b. Both plots are compared with the data from Figure 1 in [25] 12 but is sampled a lot more closely to reveal more eigenmodes, and good agreement is yet again observed.

Conclusions
In this work the solution in [18] has been generalized and extended with new boundary conditions and loads, and generalized to include materials with damping. Moreover, a scaling strategy is presented which extracts exponential behaviour from the Bessel functions in order to deal with overflow issues. The problem is significantly mitigated and extends the domain of computation by an order of magnitude for the frequency. Using quadruple precision, the problem is eliminated for practical applications at the expense of the computational efficiency. This significantly extends the range for the parameters at which the solution can be computationally evaluated.
Several numerical examples has been investigated and verified against reference solutions in the literature. The example from Skelton illustrates the importance of the scaling strategy presented herein as previously hidden important characteristics has now been reviled.

B. The incident wave
At the interfaces it is assumed that the incident wave may be decomposed as inc,m,n Z (i) n (ξ m (r)), where the coefficients are found from inc,m,n Z (i) n (ξ m (r)) = 2n + 1 2 π 0 φ inc,m (r, ϑ)Q (0) n (ϑ) sin ϑ dϑ.
A plane wave has such an expansion (it is most sensible to imposed this only in the outermost unbounded domain, that is φ inc,m = 0 for m > 1) φ inc,1 = Φ inc,1 e ia1r cos ϑ = Φ inc,1 ∞ n=0 Q (0) n (ϑ)(2n + 1)i n Z (1) n (ξ 1 (r)), inc,1,n = Φ inc,1 (2n + 1)i n and A (i) inc,1,n = 0 for i = 2, 3. The incident wave due to a point source located at x s = r s e 3 (in domain m such that R m < r s < R m−1 ) can be more efficiently evaluated then what was presented in [18]. (2n + 1)P n (cos ϑ) j n (a m r)h (1) n (a m r s ) r < r s j n (a m r s )h (1) n (a m r) r > r s .

(B.4)
It is also possible to implement a point force mechanical excitation located at x s = r s e 3 where r s = R m for a given interface m [21] φ inc,m | r=Rm = Φ inc,m δ(ϑ) where H is the Heaviside function.
14 The relation between Φ inc,m and the amplitude of the incident wave in a non-viscous fluid used in [18], P inc,m , is given by P inc,m = ρmω 2 Φ inc,m