Time-domain immersed-boundary simulation of acoustic propagation between two spherical gas bubbles^a)

Acoustic scattering from bubbles is a fundamental research topic that continues to have a broad range of applications. The pulsation resonance frequency, damping, scattering, and extinction cross sections of bubbles are important for various acoustic applications¹ with the linear acoustic regime for bubbles encompassing ocean and underwater applications. Role of bubble geometry and multiple bubbles interactions remain topics of ongoing analytical,^2–10 numerical,^11,12 and experimental^13,14 research.

Linear resonance frequencies of a single bubble have been predicted with an empirical shape factor method^2–4 and subject to comparisons to immersed-boundary methods.^11,12 Extensions for the resonance of multiple bubbles have received less attention. Feuillade and Werby⁵ and Feuillade^6,7 provided analytical solutions of two and three same-sized bubbles with the same- or opposite-phase oscillation patterns. Gaunaurd et al.⁸ studied two different-sized bubbles subject to an incident plane wave at arbitrary angles. Valier-Brasier and Conoir⁹ extended the two-bubble scattering for viscous and thermal damping. However, in these previous studies, the corresponding interior gas motion was not investigated. Moreover, these analytical solutions are limited with respect to bubbles in close contact or touching.

Extensions of single spherical bubble formulas¹⁰ to multi-bubble systems have been attempted within limits,^3,4 spurning the development of numerical solutions based on acoustic wave propagation.^15,16 The acoustic fields at both the exterior and interior of the bubbles for two in-contact bubbles are still less studied due to the difficulty of dealing with the geometric contact in the numerical methodologies. Experimental^13,14 comparisons with the analytical and numerical results have been primarily for same-sized bubbles.

We investigate a two-bubble system, with different sizes and distances between bubbles, such that the acoustic pressure and velocity profile—both outside and inside the bubble—and the resonance frequency of the scattering field, are solved using a combination of the finite-difference time-domain (FDTD) methods and the immersed-boundary (IB) method.^11,12 Specific cases of two in-contact bubbles illustrate the advantages of the method and applications related to fish swim bladders. This method can incorporate the entire domain of multiple geometries simultaneously, and facilitate interface contact between bubbles.

In the above expressions, $P 0$ is the ambient pressure in air, or ambient internal pressure of the bubble, $γ=1.4$ is the heat capacity ratio, $α a= c a 2 γ P 0$⁠, and $c a$ is the speed of sound in the air. Viscous and thermal effects are neglected at the current stage. A full description of the computational method and validation have been reported previously in Refs. 11 and 17. In this study, $ρ w=1050 kg/ m 3$⁠, $c w=1531 m/s$⁠, $c a=340 m/s$⁠, and $ρ a$ is calculated through the ideal gas relation, $ρ a= γ P 0 c a 2$⁠, where $P 0= 10 6 Pa$⁠.

The computational domain incorporates a Cartesian coordinate system (⁠ $x$⁠, $y$⁠, $z$⁠). The geometry for the bubbles is two spheres, where a sphere has the radius $a′$ with the center of the sphere at the negative side of the $y$-axis as front bubble and another sphere of radius $a$ with the center of the sphere at the positive side of $y$-axis as back bubble. An artificially infinite inviscid fluid medium surrounds the bubbles by applying perfectly matched layers (PMLs). An incident acoustic plane wave with frequency $f i$ and pressure amplitude $P i$ is traveling in the positive $y$-direction and impinging upon the bubbles.

Figure 1 shows a side view of how the simulations of the computational domain of the $y$- $z$ plane at $x=0$ are set up. The plane wave is specified as a monochromatic sinusoidal wave with a fixed wavenumber $k=2π f i/ c w$ from the negative $y$-direction of the computational domain located at $10a′$ away from the center of the front bubble. Table 1 shows the cases presented in the paper. Other cases were simulated during the study but are not presented in the paper for the brevity. This value of $ka$ is selected to be larger than 1 to study the higher modal oscillations.^18,19 The distance between the centers of the two bubbles is $d$ while the centers of the bubbles are located on the $y$-axis at $0 , ± d / 2 , 0 m$⁠. The receiver arrays at the interior of the bubble are parallel to $x$⁠, $y$⁠, and $z$-axes, intersecting at the center of the bubble and for another arrangement, the two bubbles are along the $z$-axis instead of the $y$-axis. The results from an alternative geometry related to fish swim bladders are shown in a larger-scale simulation such that the resonance frequencies can be focused on the monopole and the relative effect of size difference. This large-scale simulation will be discussed in a separate section. The computation grid mesh is a Cartesian grid mesh with $Δx$⁠, $Δy,$ and $Δz$ set as 1/10 of the radius of the bubble in each direction, $Δx$ =  $Δy$ =  $Δz$ = 0.00001 m = 0.1  $a$⁠. The time step, $Δt$⁠, is set to be $2× 10 − 10$ s. The small-time step corresponds to a small Courant–Friedrichs–Lewy (CFL) number equal to $0.03$ and provides sufficient resolution for the waves at the higher frequency parameters.^11,17 The order of accuracy of the scheme is second order in both space and time.^11,17 An incident sinusoidal plane wave of the amplitude magnitude equal to 0.000 01 Pa is investigated for the linear acoustic regime and generated from the negative y-boundary. All the other boundaries are set as perfectly matched layer (PML) boundaries, as Fig. 1 shows, which represents the free space on these boundaries. The computational domain size has $50$ grids of PML in the $x$- and $z$-direction and $250$ grids of PML in the positive $y$-direction, where the size varies based on the distance between bubbles, $d$⁠. The interior background pressure of the bubbles is $P 0= 10 6 Pa,$ not only for the best visualization of the resonance responses,^11,12 but also for an assumed condition of a few hundred feet under the surface of the sea.^20,21 A study of the numerical capability for the interior background pressure was reported earlier using the same numerical package.¹² 

We investigate linear resonance response from two spherical bubbles with respect to variations in size and distance extending previous work.^11,12 The acoustic velocity is shown over the receiver arrays with respect to the resonance frequency. The bubble boundary is assumed fixed, an approximation suitable for linear acoustics.^11,12 For bubbles with radii $a=100 μm$ on the $y$-axis with a variable distance $d$ between the centers of two bubbles, the resonance responses are collected at the centers of the bubbles^2–7,11,12 with six receiver arrays at the interior along the three axes.

Linear resonance frequency responses are investigated for the underwater target applications. In Sec. 3.1, frequencies related to monopole, dipole, and quadrupole oscillation are determined. Multiple resonance peaks correspond to incident waves, symmetric/anti-symmetric modes, and higher modal frequencies.⁷ The symmetric resonance is generated by the in-phase oscillation and the anti-symmetric resonance follows the opposite-phase oscillation. The lowest frequency peak is the symmetric resonance frequency of the bubble system, namely, $f r L$⁠, with another peak corresponding to the anti-symmetric resonance frequency, $f r S$⁠. We highlight the related resonance frequency responses since the $f r L$ and $f r S$ are close to larger and smaller bubbles' natural frequency.¹² Higher order modes are the dipole mode frequency, $f d$⁠, and quadrupole mode frequency, $f q$⁠, respectively. The incident wave frequency $f i$ is apparent in power spectral density (PSD). Table 1 shows all the cases simulated and analyzed for this study with the Minneart frequency² of the smaller bubble $f r S , M$ provided as a reference.

In cases 1–4, bubbles of the same size with different distances are simulated following the results reported in Valier-Brasier and Conoir.⁹ Larger-scale cases for two different-sized bubbles and a lower frequency incident wave were also computed. These simulations were constructed such that the relative radius of the smaller bubble $ka=0.1$ and the monopole mode oscillation is obtained. Only vertical cases are presented due to the similarity to tandem ones. In addition, case 8 replaces the larger bubble at the positive $z$-axis with a prolate spheroidal bubble as a preliminary model of a fish bladder. The $a′$ represents both the major axis and the minor axis of the prolate bubble in this case.

The interior acoustic field for each individual bubble was computed as in our previous studies¹¹ for a three-dimensional distribution of the interior pressure and velocity. Two different-sized bubbles¹² show the non-uniform interiors with off-center concentrations of pressure and velocity.

Resonance frequency shifts reported⁶ from multiple scattering theory are reproduced for two same-sized bubbles for the configuration in Ref. 9, located on the $y$-axis along the wave propagation direction. The geometry is selected such that the resonance frequencies of closely located up to touching bubbles are investigated. Even with viscous and thermal effects neglected, a simulation comparison has been made appropriately with the study of Valier-Brasier and Conoirt⁹ using the same configuration in Fig. 1. The numerical setting follows Sec. II. The bubble at (0, −d/2, 0) m is the front bubble and the bubble located at (0, d/2, 0) m is the back bubble. The incident wave propagates from the negative $y$-direction with incident frequency $f i$ =  $3062 kHz$ and interior background pressure $P 0$ =  $10 6 Pa$⁠. Under conditions with $ka≫1$⁠, higher order oscillation modes may occur. The situation where a distance between the two bubbles corresponds to the wavelength when $d=5a$ is case 4 in Table 1.

Figures 2(a) and 2(b) show the PSDs obtained from the centers of each bubble and a zoomed in view of the monopole resonance peak in Fig. 2. The symmetric resonance peaks, $f r L$⁠, are comparable to those reported by Feuillade;^6,7 however, the anti-symmetric resonance frequency, $f r S$⁠, has shifted to a larger value compared with the analytical solution with a decrease in the distance b. This shift is close to the resonance shift reported in Ref. 9, indicating the limits of the multiple scattering approximations in truncation.⁶ In case 1, the two bubbles are touching; analytical solutions fail due to the singularity in the multiple scattering theory. The $f L S$ is still close to the predicted symmetric resonance frequency but the $f r S$ is much higher than previously predicted. In case 4, due to the larger distance, only the front bubble experiences the scattered wave of the back bubble while the back bubble did not receive a signal reflected from the front bubble. These results at this separation distance have been previously studied^6,7 and verified numerically.⁹ Extremely long-distance cases are examined in which the bubbles behave as uncoupled oscillators. Higher order oscillation responses are found, notably, $f d$ and $f q$⁠, which are dipole and quadrupole related modes.

For the initial study, the symmetric monopole-type oscillations are studied on the frequency domain as $f r$ and the anti-symmetric is neglected due to its absence in large distance cases and the low magnitude in PSD. The magnitudes of the PSD of sound velocity component at each direction are $|U|$⁠, $|V|$⁠, and $|W|$⁠. These will be referred to as velocity levels. The pressure level (not shown) at both $f r S$ and $f r L$ shows a notable similarity to the interior pressure distribution at the monopole mode in a single bubble (Fig. 7 in Ref. 11). The sound velocity level captured on the $y$-axis is shown in Fig. 3 for each case, with the magnitude of sound velocity at the resonance frequency normalized by the incident amplitude $U i= P i / ρ w c w=6.221× 10 − 12 m/s$⁠. The array location is non-dimensionalized by $a$⁠. A conceptual diagram of the geometry is shown on the top right with the red lines indicating the array locations.

Figures 3(a) and 3(b) show the pressure with colored contours and velocity vectors arrows for case 2, $d=2.2a$⁠, on the $y$- $z$ plane. The velocity magnitude is scaled for improved visualization with a reference velocity vector shown at the top right. Areas of small velocity magnitude indicate negligible deformations of the bubble gas–water interface. In both Figs. 2(a) and 2(b), the front bubble has a dipole-like pressure distribution while the back bubble exhibits dipole and quadrupole responses. The velocity vectors reveal that there are multiple oscillation centers in coupled, spherical bubbles due to the contribution of higher order modal oscillations with a time lag from radiative damping. The dipole and quadrupole oscillations correspond with source or sink centers while the monopole can be mainly discerning as spatially symmetric displacements equidistant from the interface boundary. Larger velocity magnitudes appear on the boundary of the bubble due to the combination of strongly coupled bubble oscillations and interior wave reflection from a symmetric boundary. The interior velocity along the wave propagation direction is shown in Figs. 3(c)–3(f).

In Figs. 3(c)–3(f), multipoles are captured from monopole (solid lines), dipole (dashed lines), and quadrupole (dash-dotted lines). All velocity profiles at the monopole are shifted in the front bubble. In the back bubble, the velocity profiles are closely symmetrical to the front bubble except for case 4, which also indicates that the back bubble does not experience the reflection wave from the front bubble for the time duration of the simulation. Notably, the dipole velocity profiles only shift when the two bubbles are close, namely, $d=2a$ and $d=2.2a$⁠, demarcating the regime higher order scattering interactions. Greater distances tend to have less influence on dipole mode oscillations compared to the quadrupole modes. However, quadrupole oscillations seem less stable for increasing separation, as dipoles dominate the coupling. These shifts can be seen in $d=2a$ but become less obvious at greater separations.

In this section, an alternative geometry is studied with the bubbles located on the $z$-axis with a fixed distance between the centers of the bubbles which is representative of typical anterior and posterior fish bladder chambers.²² A smaller bubble of radius $a=0.01 m$ with its center located at the negative axis of the $z$-direction at $0 , 0 , − 0.015 m$ along with a larger bubble of radius $a′$ with its center located at the positive axis of the $z$-direction at $0 , 0 , 0.015 m$⁠. The distance between the centers of the two bubbles is fixed for $d=0.03 m$⁠. The incident acoustic frequency $f i$ is $2437 Hz$⁠. For parameters, the dimensionless radius of the smaller bubble is fixed for $ka=0.1$⁠. This value of $ka$ is selected for a monopole response. For the specially designed case 8, the larger spherical bubble has been replaced by a prolate spheroidal bubble with its minor axes are $0.5a$ on $x$- and $y$-axes and its major axis is $4a$ on the $z$-axis. The center of the prolate bubble has been changed to $0 , 0 , 0.035 m$ such that two bubbles are still touching together. The Cartesian grid sizes are $Δx$ =  $Δy$ =  $Δz$= 0.001 m = 0.1  $a$ and the time step is set to be $2× 10 − 8$ s which maintains the same CFL number as in Sec. 3.1. The grid size is selected as $1/10$ of the smaller bubble radius where the larger bubble may have finer spatial resolution. The magnitudes of the interior acoustic velocity level at both $f r S$ and $f r L$ are provided for each of the bubbles for cases 5–7 in Table 1. All the physical variables are determined at the resonance frequencies, $f r S$ and $f r L$⁠, in these figures with these values normalized as in Sec. 3.1.

Figure 4(a) illustrates the magnitude of resonance the $z$-direction velocity component, $|W|$⁠, collected from the interior of both larger and smaller bubbles with the receiver location normalized by $a$⁠. A diagram on the top shows the receiver array locations. Three cases are selected for $a ′=a$⁠, $a ′=1.5a$⁠, and $a ′=2a$⁠. Note that in case 7, two bubbles in contact are computed with the immersed-boundary method. Receiver arrays are shown in the top right geometry diagram as red lines along the $z$-direction across the bubbles' centers. The receiver location is nondimensionalized with respect to the radius of the bubbles, $a$ or $a′$⁠, and the relative location to the center of the larger bubble. The velocity level results on the $x$- and $y$-axes are similar to those of a single bubble with a small shift in the larger bubble at $f r S$⁠. All velocity levels at $f r L$ (and $f r$ for case 5) rise from the center of the bubble to the top end demonstrating that the larger bubble tends to oscillate at a “semi-monopole” type. For case 7, where the two bubbles are contacting, both $f r S$ and $f r L$ have significantly high levels at $z/a=−0.5$⁠. This phenomenon is likely due to the relatively large oscillation amplitude of the smaller bubble influenced by the resonance oscillation of the larger bubble. At $f r S$⁠, the velocity level drops when the receiver location reaches the top of the bubble. A similar phenomenon also occurs in case 6 at $f r S$ but the velocity level undergoes limited changes near the top of the bubble. For two identically sized bubbles, the oscillating center moves closer to the other bubble instead of remaining at the center of the respective bubbles. For two different-sized bubbles, velocity level distributions at $f r L$ decrease for receiver locations away from the larger bubble. The dips for the velocity levels at $f r S$ shift away from the center of the bubble towards the bottom of the smaller bubble and with larger shifts for corresponding larger relative radius ratios. The highest velocity level is at the contact or the closest point on the bubbles for all cases, corresponding with the larger velocity vector at the contact point.¹² The smaller bubbles still maintain the “semi-monopole” type oscillation, but the oscillation centers are shifted due to the interaction with the larger bubble.

Figures 4(b) and 4(c) show two touching bubbles (as in case 7) with the pressure given as colored contours and the vector arrows depicting acoustic velocity on the $y$- $z$ plane. Note that the arrows are scaled for enhanced visualization with a reference velocity vector at the bottom right. The small velocity illustrates the small deformation of the bubbles. Likewise, the interior pressure distribution at any time instance is spatially non-uniform even though the overall pressure might be considered approximately uniform. The oscillating center in the larger bubble has shifted along $y$-direction for symmetric smaller bubble responses. Larger velocity occurs near the two bubbles' boundary for increased contact, differing from close proximity.¹² Inclusion of viscous and thermal damping as well as boundary deformation will be among the subjects for future work. The opposite phase occurs, but it is not as obvious, such that only in-phase contours are reported.

Figure 4(d) illustrates a prolate spheroidal bubble touching with a spherical bubble to model a fish swim bladder captured at a time instance where both bubbles are expanding. Note that the volumes of both bubbles are the same. The oscillation center and the pressure difference in the prolate bubble is clearly shown. It is shown that the velocities are of larger magnitude at the two tips of the prolate, similar to a single prolate bubble¹¹ while the smaller spherical bubble oscillates at an approximate monopole type. Similar to case 7, at and near the touching point, the velocities increase dramatically and are pushing towards the smaller bubble, showing a forced oscillation from the prolate bubble.

The resonance response of two spherical bubbles is investigated with a FDTD method combined with the immersed-boundary method. With the immersed-boundary method, the pressure and velocity of the interior gas inside the bubble, which has been neglected in previous numerical studies, are simulated together with the outside scattering field. Agreement is found with approximate solutions for the symmetric and anti-symmetric modes for two identical-sized spherical bubbles;¹² more comparisons with experiments will be done in future studies. The oscillation center at the dominant resonance frequency shifts closer to the second bubble while the oscillation center at the minor resonance frequency tends to separate. Alternative arrangements encompassing large and small scales reveal the direction of the incident wave has a limited effect. A non-spherical touching case has been shown to further demonstrate the applicability of the simulation. Acoustic scattering from multiple bubbles with various shapes, touching, or in contact bubbles, has applications for fish swim bladders (dual chamber) and bubbles entrapped in sediment.

J. S. A. acknowledges support from U. H. Manoa (sabbatical leave) and Office of Naval Research #N0014-20-1–2651. The authors have no conflicts of interest to declare.