In this work, we established a novel theory for the dynamics of oscillating bubbles such as cavitation bubbles, underwater explosion bubbles, and air bubbles. For the first time, we proposed bubble dynamics equations that can simultaneously take into consideration the effects of boundaries, bubble interaction, ambient flow field, gravity, bubble migration, fluid compressibility, viscosity, and surface tension while maintaining a unified and elegant mathematical form. The present theory unifies different classical bubble equations such as the Rayleigh–Plesset equation, the Gilmore equation, and the Keller–Miksis equation. Furthermore, we validated the theory with experimental data of bubbles with a variety in scales, sources, boundaries, and ambient conditions and showed the advantages of our theory over the classical theoretical models, followed by a discussion on the applicability of the present theory based on a comparison to simulation results with different numerical methods. Finally, as a demonstration of the potential of our theory, we modeled the complex multi-cycle bubble interaction with wide ranges of energy and phase differences and gained new physical insight into inter-bubble energy transfer and coupling of bubble-induced pressure waves.

## I. INTRODUCTION

Bubbles are ubiquitous in nature and of great significance in numerous fields of science and engineering such as marine science, ocean engineering, mechanical and material engineering, environmental and chemical engineering, and medicine and life science.^{1–6} Bubble dynamics is a fundamental scientific problem attracting widespread interest. The pulsating bubbles that undergo volumetric oscillations due to pressure imbalance have many important applications. For example, some acoustic and laser cavitation bubbles,^{7–14} which are typically millimeter or micrometer-sized oscillating bubbles, can facilitate ultrasonic cleaning,^{15,16} sonoluminescence,^{17–19} ultrasound therapy,^{20–27} drug/gene delivery,^{28–32} inkjet printing,^{33–35} and bubble propulsion.^{36,37} Underwater explosion/implosion bubbles and air-gun bubbles, which can be meter-sized oscillating bubbles, are central to underwater explosions/implosions^{38–44} and geophysical explorations,^{45–47} respectively. Hydraulic cavitation bubbles,^{48–55} entrained air bubbles/cavities,^{56–62} heat-generated vapor bubbles,^{63–65} and spark-generated bubbles,^{66–71} which also oscillate and are usually sized from millimeters to meters, can be crucial to the performance of turbines, propellers, waterborne vehicles, engines, reactors, and spark sound sources. There are also oscillating bubbles with ultra-small or large sizes, such as oscillating nanobubbles^{72,73} or the bubbles produced in an underwater volcano eruption.^{74} Today, new studies on bubble dynamics keep emerging and new applications^{75–84} of the oscillating bubbles are being discovered.

The dynamics of bubbles, of either natural or artificial sources, is complicated due to the existence of different scales, various boundaries, and multiple oscillation cycles that bring enormous challenges to theoretical, numerical, and experimental research. Theoretical research is crucial to understanding the physics of bubbles under different conditions, which is the foundation for the utilization of bubbles and the circumvention of their hazardous effects. The theoretical research on bubble dynamics dates back to the early 20th century. Rayleigh^{85} proposed a preliminary mathematical description of the collapse of a cavity in an infinite fluid field. Based on that, Plesset derived the classic Rayleigh–Plesset (RP) equation for bubble oscillation in an incompressible flow.^{86} The loss in bubble energy due to acoustic radiation, such as bubble collapse-induced pressure waves,^{87–94} is overlooked due to the incompressible assumption and the RP equation may not be suitable when the damping in bubble radius, period, or energy is important. Taking into consideration the weak compressibility of the fluid outside the bubble, Herring,^{95} Gilmore,^{96} Trilling,^{97} Keller and Kolodner,^{98} Hickling and Plesset,^{99} Flynn,^{100} Keller and Miksis,^{101} and Fujikawa and Akamatsu^{102} established weakly compressible bubble models. Among them, Keller's work in which the model was formulated using the wave equation and the incompressible Bernoulli equation enjoyed a most widespread application. Prosperetti and Lezzi^{103,104} developed a model for bubble dynamics considering the compressibility of the far-field fluid based on the perturbation theory. These models have played a major role in the theoretical research on bubble dynamics in a free field. However, bubbles, regardless of type or scale, are seldom isolated and always in coupling with different conditions giving rise to complex bubble dynamics behaviors. In this paper, we proposed a theory describing bubble oscillation and migration in a compressible fluid. The theory provides a framework allowing for the effect of different types of boundaries, bubble interaction, background flow, gravity, migration, fluid compressibility, viscosity, and surface tension and presents a unification of important bubble models such as the Rayleigh–Plesset equation, the Gilmore equation, and the Keller–Miksis equation.

A major concern in the study of theoretical bubble models is the coupling effect of bubble oscillation and migration. Hicks^{105} developed a bubble dynamics model considering the migration based on the incompressible assumption. Best^{106} developed an incompressible model, where a time-dependent dipole is inserted to allow for the migration of the spherical bubble. The migration-induced pressure variation was manifested by an additional term $v2/4$ (surface-averaged pressure, SAP, and *v* is the migration velocity). Similar methods were adopted by Brujan *et al.*^{107} and Seo *et al.*.^{58} Geers and Hunter^{108} proposed a doubly asymptotic approximation (DAA) model specialized for a single underwater explosion bubble in a free field considering fluid compressibility and gravity-induced migration. To the best knowledge of the authors, extensions of the model to incorporate bubble migration for interacting bubbles or bubbles near boundaries in a compressible flow are still not available. The theory proposed in the current study can consider the effect of compressibility on not only bubble oscillation but also migration of various causations.

Another important issue for the theoretical bubble models is the interaction between multiple bubbles. Different bubble models have been applied to predict the behavior of a bubble in a multi-bubble system. Harkin *et al.*^{109} established a model for the spherical oscillation and translational motion of a pair of interacting gas bubbles in an incompressible liquid using Weiss sphere theorem and the method of SAP. Bremond *et al.*^{110} extended the RP equation by considering the pressure modification induced by the surrounding bubbles and studied the weak interaction between two or more cavitation bubbles before the first bubble collapse. A similar scheme was also applied by Gonzalez-Avila *et al.*^{111} for cavitation bubbles and by Ziolkowski *et al.*^{112} and Li *et al.*^{113} for air-gun-bubble clusters in geophysical exploration. However, these previous works ignored the traveling time of pressure waves, which could be an oversimplification for the reproduction of the bubble-induced pressure field.^{42,114–116} To overcome this, we considered the time delay for the pressure wave propagation and our theory can reproduce some interesting phenomena that were not found with the previous models, such as the reflection of a pressure wave on bubble surfaces.

The effect of fluid boundaries on bubble behavior is another problem for the spherical bubble theories. To reproduce the bubble dynamics near a rigid wall, incompressible bubble models with the image method^{106,107} as well as the Gilmore model with the SAP method^{15,117} were used. The image method has also been applied to the free surface boundary.^{118} For more accurate prediction, we incorporated in our model the effect of the bubble-induced velocity fields, pressure waves, and their reflections at the boundaries on the bubble behavior. This is based on the consideration of fluid compressibility and traveling time of disturbances. It is worth noting that the negative pressure that is formed by the reflection of pressure pulses at the free surface can be captured by our theoretical model.

This paper is structured as follows. The present theory for bubble dynamics, including the unified equations for bubble oscillation and migration, is derived in Sec. II. The theory is extended to multiple-bubble interaction in Sec. III and bubble dynamics near boundaries in Sec. IV. A comparison with previous bubble models and experimental results with a variety in the scale and source of bubbles, boundaries, and other ambient conditions is presented in Sec. V. A discussion on the applicability is given in Sec. VI followed by a demonstration of the capability of the present theory via solving complex interaction between two bubbles with phase and energy differences. Finally, this work is summarized and conclusions are drawn in Sec. VII.

## II. DERIVATION OF THE UNIFIED THEORY FOR BUBBLE DYNAMICS

### A. Basic equations of the theory

First, we derive the fundamental equation of the unified bubble dynamics theory from the basic laws of conservation.

The density, pressure, velocity, and deviatoric stress tensor of the fluid outside the bubble are denoted as *ρ*, *p*, ** u**, and $\tau $, respectively. The conservation equations for mass and momentum are

and

where $F$ is the body force. Meanwhile, the sound speed *C* in the external fluid field is related to the pressure *p* and the fluid density *ρ* as

The sound speed *C* is usually much greater than $|u|$. Thus, we adopt the weakly compressible assumption and neglect the density gradient in the above equation. The oscillating bubble driven by inertia is usually of large Reynolds numbers and short oscillating cycles. When only a few cycles are considered, the effect of viscosity can be neglected^{101,103,119} here. Thus, the first item of the last term in Eq. (4) is dropped, and we will introduce the viscosity back to the system later from the dynamic boundary condition of the bubble surface. Assuming that the flow around the bubble is irrotational, then there exists the velocity potential $\phi $ that satisfies $u=\u2207\phi $. Substituting it into Eq. (4) and integrating from infinity to the current location while neglecting the variation of the body force and the terms containing $C\u22122$ except the first one, we then have the wave equation

For an oscillating bubble, migration may happen due to the anisotropy of the surrounding flow such as the pressure gradient caused by nearby boundaries or the gravity field. We tie the origin of a spherical coordinate system $o\u2212r\theta \varphi $ to the bubble center. *θ* = 0 points to the direction of the bubble migration velocity relative to the ambient flow denoted by ** v**. $v=ve=vm\u2212ua$, where $e$ is a unit vector along the direction of

**, $vm$ is the migration velocity, and $ua$ is the ambient flow velocity. $ua=uB+uE$, where $uB$ is the velocity induced by boundaries or other bubbles and $uE$ represents extra velocities including the velocity of the incoming background flow $u\u221e$.**

*v*In this work, the bubble is assumed to oscillate and migrate while maintaining a spherical shape. We construct a solution of Eq. (5) with two singularities, i.e., a source and a dipole, the strength of which are denoted by *f* and *q*, respectively. The velocity potential at an arbitrary location ** r** and time

*t*induced by the moving singularities can be found in Appendix A. However, it is hard to write down explicit expressions when they are moving at varying speeds. If we consider the flow field within a small range around the bubble, despite the influences from boundaries and other bubbles, the direction of bubble migration can be deemed constant and the velocity field axisymmetric to

*θ*= 0. When

**is close to the bubble and considering that the migration velocity**

*r**v*is small compared with

*C*, the variation of the location of the singularities can be neglected during the very short time when the influences are propagating from the singularities to an arbitrary location

**. Thus, the velocity potential at $r=(r,\theta )$ and**

*r**t*is expressed as the summation of Eqs. (A2) and (A3) after simplification and reads

in which $q\u2032$ denotes the derivative of *q* with respect to its argument.

Note that the wave equation governing the fluid field outside the bubble is not applicable to the fluid field inside the bubble and the propagation from the singularities to the bubble surface is a mathematical extension of the fluid field outside the bubble. Then, we have the two velocity components at the bubble surface, induced by the two singularities, as

and

respectively, in which $f\u2032$ denotes the derivative of *f* with respect to its argument. The kinetic boundary condition on the bubble surface for the oscillation and the migration can be expressed as

and

respectively, where *V* is the volume occupied by the bubble and $\u2202V$ is its boundary, i.e., the bubble surface. $R\u0307$ denotes the time derivative of *R*. Let us denote $F(t)=f(t\u2212R/C)$ and $Q(t)=q(t\u2212R/C)$. Substituting Eq. (7) into Eqs. (9) and (10), we have

and

respectively. Employing the perturbation method on Eq. (12), we arrive at

with the terms of the order $C\u22122$ ignored. Then, taking the derivative of Eq. (11) with respect to *t*, we have

which may be reorganized as

Equation (15) describes the oscillation of a migrating bubble and is the base of the present unified theory for bubble dynamics. Its right-hand side equals $V\xa8/4\pi $, where $V\xa8$ is the bubble volume acceleration; meanwhile, the left-hand side is equivalent to the corresponding driving force. $dF/dt$ can be determined according to the physical problem considered, enabling Eq. (15) to model complicated bubble dynamics under different conditions. Therefore, the present theory can be extended while retaining a unified, concise, and elegant mathematical form to incorporate different conditions. The above forms the basis for the subsequent work of this paper. Additionally, if we consider the non-spherical motion of a simply connected bubble, the velocity potential and the bubble surface shape may be expanded by the bubble deformation modes. If the motion of a toroidal bubble (bubble ring) is considered, we may need to move the singularities off the rotating axis and introduce an extra vortex to incorporate the velocity circulation around the toroidal section.

### B. Bubble oscillation equation

In this section, we derive the bubble oscillation equation by introducing dynamic boundary conditions into Eq. (15). The deformation of fluid particles during the bubble oscillation would introduce additional viscous stress components. Viscosity was neglected in bulk liquid in the previous derivation. However, within the proximity of the bubble boundary, the viscosity effect can be restored as a correction term to the equilibrium condition for normal stresses on the bubble surface, which involves the pressure of the internal gas on the inner bubble surface $Pg$, the fluid pressure on the outer bubble surface $Pb$, surface tension, viscous stress, and other additional terms. Here, we express the equilibrium condition as

where *σ* is the surface tension and *μ* the liquid viscosity. $Pg$ may be affected by many factors such as mass and heat transfer, wave effect, and chemical reactions, and can be determined according to the problem being considered. For simplicity, we assume that the oscillation process is adiabatic and that the internal gas pressure is uniform. Thus, $Pg$ is composed of the non-condensable gas pressure and the saturated vapor pressure $Pv$ as

where *P*_{0} is the initial pressure of the non-condensable gas, *γ* is the polytropic exponent, and *R*_{0} is the initial bubble radius.

To obtain $dF/dt$ in Eq. (15), we introduce the dynamic boundary condition expressed by the Bernoulli equation in the moving coordinate system

where $v\xb7u=(v\u2009cos\u2009\theta ,\u2212v\u2009sin\u2009\theta )\xb7(ur,u\theta )$ and $|u|2=ur2+u\theta 2$. $H=\u222bPaPb\rho \u22121dp$ is the enthalpy difference at the bubble surface. Here, $Pa$ represents the ambient pressure at the bubble center. $Pa=PB+PE$, where $PB$ is the pressure induced by boundaries or other bubbles and $PE$ represents extra pressures including the hydrostatic pressure $P\u221e$ and acoustic pressures $PA$. With the motion of the fluid being an adiabatic process, *H* resolves into

where $\varpi =(Pb\u2212Pa)/(\rho C2)$.

Note that the terms containing *θ* in Eq. (18) will lead to non-spherical bubble deformation. To comply with the spherical assumption while taking into consideration their averaged effects, we integrate Eq. (18) over the bubble surface to eliminate *θ* and obtains

Finally, substituting the above into Eq. (15), we have the bubble oscillation equation

Equation (21) is the core of the present theory that can be extended to various scenarios such as multiple-bubble interaction and bubble dynamics near boundaries. The above equation evolves into a unified and elegant mathematical form similar to Eq. (15) as

where $F\u0303=R2(R\u03072/2+v2/4+H)$ is to be determined according to specific conditions of the problem under consideration. $R\u03072/2,\u2009v2/4$, and *H* indicates equivalent forces induced by bubble oscillation, migration, and the ambient flow field, respectively.

The oscillation in a gravity field is described by Eq. (21) with the ambient flow ignored ($ua=0$). In addition, neglecting both the ambient flow and bubble migration (i.e., $v=0$), Eq. (21) is simplified as

which describes bubble oscillation in a free field. Equation (21) can be expanded as the following form:

It can be simplified to the Keller-Miksis equation when *v* = 0, $Pa=P\u221e$, and *H* is calculated with Eq. (19) retaining only the first term. On this basis, if the second term on the left-hand side is substituted with $(1\u2212R\u0307/C)RH\u0307/C$, Eq. (24) becomes the Gilmore equation.

For bubble dynamics in an incompressible flow, i.e., $C\u2192\u221e$, Eq. (15) can be reduced to

and Eq. (21) can be simplified as

which further reduces to the Rayleigh–Plesset equation when the first term is neglected and $Pa=P\u221e$. Therefore, our theory unifies the Rayleigh–Plesset equation, the Gilmore equation, and the Keller–Miksis equation.

### C. Bubble migration equation

Since the problem is migration-related, Eq. (21) involves two unknowns, *R* and *v*. Therefore, an additional equation is required to make the problem solvable. The momentum equation of the bubble migration in the gravity field can be written as

where *m* is the total mass of the gas inside the bubble, ** g** is the gravity acceleration, $Cd$ is the drag coefficient, $Ap$ is the projected area of the bubble in the migration direction, and $S(\xb7)=(\xb7)|\xb7|$ is a signed square operator. The second term on the right-hand side denotes the inviscid part of the hydrodynamic force exerted on the bubble by the surrounding fluid, and the third term approximates the viscous part.

**represents extra forces such as lift force and is not considered in this paper.**

*X*Substituting $H=(Pb\u2212Pa)/\rho $ into Eq. (18) will give us $Pb$. Then, the inviscid resultant force on the bubble can be expressed as

where $Ca$ is the added mass coefficient and $V=4/3\pi R3$ is the bubble volume and ** n** is the external normal vector.

Since the density of the gas inside the bubble is usually much smaller than that of the liquid outside the bubble, we may ignore the inertial force and the gravity force of the gas. Hence, Eq. (27) can be transformed into

For a bubble oscillating and migrating in a still and unbounded flow, we have $v=vm$ and $\u2207Pa=\rho g=(0,0,\u2212\rho g)$. Then, Eq. (29) reduces to

### D. Bubble-induced velocity and pressure fields

We can calculate the velocity and pressure fields induced by the bubble once the bubble motion is solved. The Bernoulli equation at an arbitrary location ** r** reads

where the third term is *H* calculated by Eq. (19) retaining only the first term. Combining Eqs. (6) and (20), one may obtain the relation between the physical quantities at ** r** and at the bubble surface as

and

where the higher order terms of $|r|\u22121$ and the effect of migration are ignored. Note that *R* on the right-hand side of the above two equations including those in the subscripts should also be evaluated at $t\u2212(|r|\u2212R)/C$. Here, $|r|\u2212R$ denotes the minimum distance from the bubble surface to ** r**. Substituting the above into Eq. (31) and ignoring the ambient velocity $ua$, one may obtain the formula for the pressure at

**as**

*r*Note that the pressure disturbance in the fluid field with a distance *r* from the bubble center at time *t* is induced by the bubble at an earlier time, $t\u2212(|r|\u2212R)/C$; i.e., the bubble-induced disturbance is delayed in time due to fluid compressibility. In an incompressible flow where *C* approaches infinity, we have

and

Similarly, substituting them into Eq. (31), we have the degraded counterpart of Eq. (34) in an incompressible fluid expressed as

## III. MODELING MULTIPLE-BUBBLE INTERACTION

Bubbles usually do not appear in isolation either in nature or in engineering applications, and the interaction between multiple bubbles may result in a significant alteration of the bubble behavior. Here, we derive the equations for the oscillation and migration of multiple bubbles based on the theory in Sec. II. For an arbitrary point ** r** in the fluid field, the velocity potential and the velocity induced by the migration of bubble

*N*are of the orders of $|oN\u2212r|\u22122$ and $|oN\u2212r|\u22124$, respectively, where $oN$ is the center of bubble

*N*. Thus, only the effects of the oscillation of bubble

*N*need to be considered. Hence, by taking the derivatives of Eq. (6) with respect to

**and**

*r**t*with the solution of

*f*substituted, the current velocity and the time derivative of the potential induced by bubble

*N*at

**can be expressed as**

*r*and

respectively, where $tN=t\u2212(|rN|\u2212RN)/C$ is the initiation time of a disturbance induced by bubble *N* that later arrives at $r$ at time *t*. In the equations, the subscript *M* or *N* indicates the quantities of the corresponding bubble. When considering the dynamics of bubble *M*, we interpret the effect of other bubbles on its ambient fluid field as disturbances in the velocity and the pressure fields. Thus, we may obtain the flow velocity $uB$ and the pressure $PB$ of bubble *M* induced by other bubbles, considering the influences from all the other bubbles, as

and

where *L* is the total number of the interacting bubbles and $oM$ is the center of bubble *M*. Substituting the above two equations into Eqs. (21) and (29) for bubble *M*, and setting *M* as 1, 2, $\u2026$, and *L*, respectively, one may obtain a set of 2*L* equations that describe the oscillation and migration of the *L* bubbles considering their mutual interaction. The pressure in the fluid field contains contributions from all the bubbles; therefore, the pressure at ** r** can be expressed as

When *C* approaches infinity, i.e., the external fluid is incompressible, we can obtain the reduced bubble dynamics equations as

and

where $rMN=oN\u2212oM$. The above equations with *M* = 1, 2, $\u2026$, and *L* describe the dynamics of the multiple interacting bubbles in an incompressible flow. Particularly, for the dynamics of two interacting bubbles, i.e., *L* = 2 and (*M*, *N*) = (1, 2) or (2, 1), the above equations reduce to

and

## IV. MODELING BUBBLE DYNAMICS NEAR BOUNDARIES

In addition to gravity and bubble interaction, nearby boundaries such as a rigid/solid wall or a free surface may also impose considerable influence on the bubble dynamics. In this paper, we incorporate the boundary effects by adopting the image theory in which the boundary condition can be satisfied by properly introducing image bubbles on the opposite side of the boundary. Then, the dynamics of a bubble is subjected to the influences of other real bubbles and all the image bubbles. Thus, we derive the equations of bubble oscillation and migration subject to the boundary effect by transferring the problem to a special multi-bubble interaction problem. Let us describe the boundary plane with $eb\xb7r+b=0$, where $eb$ is the outward normal vector of the boundary and *b* is a constant. Assume that the boundary can reflect the influence of a bubble with a reflection coefficient denoted by *α* depending on the property of the boundary. Particularly, for a rigid boundary, *α* = 1 and for an ideal free surface, $\alpha =\u22121$. For bubble *N _{i}*, which is the image of bubble

*N*, it is required that $oNi=oN\u2212(2oN\xb7eb+b)eb,\u2009RNi=RN,R\u0307Ni=\alpha R\u0307N,\u2009R\xa8Ni=\alpha R\xa8N$, and $vNi=\alpha [vN\u22122(vN\xb7eb)eb]$, respectively. Therefore, similar to Eqs. (40) and (41), the flow velocity and the pressure at $oM$ induced by other bubbles and all the images can be written as

and

Substituting the above two equations into Eqs. (21) and (29) for bubble *M*, we then have the equations for the dynamics of one or more bubbles considering boundary effect in a compressible flow. The above derivation can be further extended to various scenarios such as those with flexible, multiple, or hybrid boundaries. The pressure in the fluid field contains the contributions from all the bubbles and their images. Then, the pressure at ** r** can be written as

Note that when the reflection coefficient *α* of the boundary is negative, the pressure waves induced by the image bubbles become rarefaction waves that may fall within the cavitation limit of the liquid. Therefore, it may be necessary to apply a modification, e.g., the cutoff cavitation model, to the pressure to incorporate the cavitation effect. When *C* approaches infinity, i.e., the fluid outside the bubble becomes incompressible, the oscillation and migration equations can be rewritten as

and

and one may then obtain the bubble dynamics considering the boundary effect and bubble interaction in an incompressible flow. Particularly, for the dynamics of a single bubble interacting with a boundary, the above equations further reduce to

and

where $db$ is the stand-off distance from the bubble center to the boundary.

## V. COMPARISON BETWEEN THEORETICAL AND EXPERIMENTAL RESULTS

In this part, we validate our unified theory with experiments of acoustic bubbles, electric spark bubbles, underwater explosion bubbles, laser-induced cavitation bubbles, and air-gun bubbles that are varied in bubble scales, sources, boundaries, and other ambient conditions, as well as the results of classical bubble theories to demonstrate the advantages of the present theory.

### A. Bubble dynamics in a free field

First, we verify the correctness and validity of our unified theory for bubble dynamics using the experimental data of single sonoluminescence cavitation bubbles from Matula^{120} and results of the classical bubble theories from Plesset,^{86} Gilmore,^{96} Keller and Miksis,^{101} and Prosperetti and Lezzi^{103} concerning the dynamics of a bubble in a free field. In the experiment, the bubble was initially in an equilibrium state with a radius of 4.7 $\mu $m before being stimulated at *t* = 0 by a high-frequency acoustic wave, $\u2212Pamp\u2009sin\u2009(2\pi fat)$, with the pressure amplitude $Pamp=125$ kPa and the frequency *f*_{a} = 29 kHz. In the theoretical calculations, we set the surface tension as $\sigma =0.075$ N/m, the fluid viscosity as $\mu =0.001$ Pa s, the polytropic exponent as $\gamma =1.4$, and the saturated vapor pressure as $Pv=2338$ Pa (same for Secs. V B–V D). We measure the effect of gravity with the buoyancy parameter $\delta =\rho gRmax/P\u221e$ in which $Rmax$ is the maximum bubble radius. The effect of gravity is neglected since the buoyancy parameter *δ* is only 0.002. The theoretical and experimental results of the bubble radius variation are depicted in Fig. 1. Driven by the acoustic wave, the bubble expands to a maximum radius of 40 $\mu $m before collapsing and continues to oscillate with a high frequency. The maximum radius after the first rebounding is merely 12 $\mu $m, indicating a drastic loss in bubble energy. The above process is successfully predicted by our bubble theory and Gilmore's, Keller's, and Prosperetti's bubble equations, which can be proof of the validity of our theory for bubble dynamics in a free field neglecting gravity. Also, it is shown in Fig. 1 that the result of the Rayleigh–Plesset equation deviates from the experiment after the first bubble oscillation because the model neglects fluid compressibility and thus cannot correctly reproduce the damping in maximum radii of the acoustic bubble. Note that the result obtained from the Gilmore equation deviates from other compressible models, which coincides with the analysis presented by Prosperetti and Lezzi^{103} that the accuracy of the Keller model is better than the Gilmore model.

### B. Bubble dynamics in a gravity field

In this section, we continue to validate our theory for bubble dynamics in a gravity field using the experimental data of a spark bubble generated in a sealed container. The ambient pressure was reduced to 6.82 kPa to enhance the influence of gravity on the bubble behavior, resulting in significant bubble migration. The experimental setup is shown schematically in Fig. 2(a), and furthermore, details can be found in Zhang *et al.*^{68} The bubble reached a maximum radius $Rmax=29.8$ mm with the buoyancy parameter *δ* = 0.21. As shown in Fig. 2(b), the bubble oscillated with significant upward migration due to gravity. For theoretical calculation, the initial conditions of the bubble are set as $P0=45$ kPa, $R0=2.07$ mm, and $R\u03070=100$ m/s. In this paper, we determined the initial conditions for modeling spark and laser-induced bubbles with the present theory using the method introduced in Wang,^{121} which requires the maximum radius $Rmax$ and the minimum pressure $Pmin$, when the bubble radius reaches $Rmax$. $Pmin$ can be determined by the ratio $Rmax2/Rmax$, where $Rmax2$ is the maximum radius during the second oscillation cycle. A backward integration is performed starting from the bubble at the maximum expansion with zero wall velocity using the fourth-order Runge–Kutta method until the initial condition of the bubble is obtained. The added mass coefficient and the drag coefficient in the bubble migration equation are set as 1.0 and 0.5, respectively. The theoretical results are plotted against the experiment in Figs. 2(c) and 2(d). The present theory yields a better prediction of the bubble radius variation compared to Keller and Gilmore's equations for the second bubble oscillation cycle. It also correctly reproduces the strong bubble migration due to gravity, which is not obtained with the other two models. If we turn off the migration in our theory, the radius curve would be identical to that from the Keller equation. Thus, it is crucial to incorporate the effect of migration when buoyancy is significant. There is still a deviation between the present results and the experiment mainly after the first oscillation cycle that could be possibly due to multiple factors, including the gravity-induced non-spherical bubble behavior, such as jetting/splitting, spark-induced combustion/plasma, and liquid/vapor phase transition during the process. In addition, more accurate results may be achievable by changing the EOS for the present theory.

Next, we carry out further validation of the present theory with a small-charge underwater explosion bubble. The experiment setup is shown in Fig. 3(a). In the experiment, an explosive charge equivalent to 1.125 g PBX9501 (1.25 g TNT) was detonated at the center of a 4 × 4 × 4 m^{3} cubical tank filled with water. The water depth was adjusted so that the charge is 1.4 m below the free surface. The explosion generated a gas bubble with a maximum radius $Rmax=0.16$ m. The effect of the water surface on the bubble is unsubstantial due to the distance to the surface being significantly larger than the maximum bubble radius. The buoyancy parameter is *δ* = 0.12. Due to the buoyancy effect, the bubble migrated upward while oscillating, as shown in Fig. 3(b). The initial conditions of the bubble are obtained with the method given in Appendix B as $R0=42.52$ mm, $R\u03070=61.80$ m/s, and $P0=304.13$ kPa. The bubble dynamics are calculated using the present theory as well as the Keller equation and the Geers–Hunter model.^{108} The theoretical results are compared to the experimental data in Figs. 3(c)–3(e). The bubble radius variation predicted by our theory matches the experimental result with slightly higher accuracy than the other two models, as shown in Fig. 3(c). Also, the migration calculated with our theory is more comparable to the experiment than that with the Geers–Hunter model, see Fig. 3(d). The migration was not reproduced by the Keller equation. In addition, Fig. 3(e) shows the comparison between the theoretical and the experimental results of the pressure pulse induced by the first bubble collapse. The pressure test point in the experiment is at the same depth as the charge with a distance of 0.7 m. Compared with the results from the experiment and the present theory, the Keller equation significantly overestimates the bubble pulse due to the absence of the migration, while the Geers–Hunter model underestimates it because of the introduction of excessive dissipation.

### C. Dynamics of multiple interacting bubbles

In this section, we present experimental results of the dynamics of multiple interacting bubbles and further validated our theory by comparing the calculation to the experiment. In the first experiment, two bubbles were induced by electric sparks powered by a charged capacitor bank. We manually created two defects on a copper wire 0.21 mm in diameter that links the two poles of the capacitor bank. Upon discharge, two sparks were produced at the two defects due to increased local resistance, thus producing a bubble pair. Physical deviations in the two defects were intentionally created to give the bubble differences in size and phase. The distance is 250 mm from the bubbles to the water surface and 200–250 mm to nearby walls; thus, the bubbles are deemed far from boundaries considering the small radii (10–20 mm). A schematic diagram of the experiment setup is shown in Fig. 4(a). Two experiment cases were carried out where the two defects are horizontally placed at a distance of 93.6 and 42.7 mm, respectively, from each other, which are taken as the initial distances between the two bubbles in each pair. The bubbles are marked as bubbles 1 and 2 from left to right. In the first case, the maximum radii of the two bubbles were 16.0 and 14.6 mm and the first oscillation periods were 3.05 and 2.83 ms, respectively. The initial radii used for theoretical calculation are $R01=2.50$ and $R02=2.29$ mm for bubbles 1 and 2, respectively. The initial speeds of the bubble wall are $R\u030701=R\u030702=130$ m/s, and the initial pressure inside the bubbles is $P01=P02=1.2$ MPa. The calculation is conducted from bubble initiation to the third oscillation period after which the experimental data of migration and radius become hardly measurable. The experiment and calculations for this case are summarized in Figs. 4(b)–4(d). The high-speed images of the bubble pair are shown in (b), where the two bubbles migrate toward each other during oscillation due to the Bjerknes force, which is similar to the attraction acting on an oscillating bubble by a nearby rigid wall. Time curves of the radius and the bubble center position are plotted in Figs. 4(c) and 4(d), respectively. The results calculated with our unified theory are compared to that of the experiments as in (c) and (d), and a good agreement is obtained.

In the second experiment case, the energy difference between the two bubbles is increased. The maximum radii were 14.8 and 9.0 mm, respectively, and the first oscillation periods were 2.92 and 1.75 ms, respectively. The high-speed images of the bubble pair are shown in Fig. 5(a). The initial conditions for theoretical calculation are $R01=1.60$ mm, $R02=1.65$ mm, $R\u030701=240$ m/s, $R\u030702=110$ m/s, $P01=1.2$ MPa, and $P02=1.2$ MPa. Due to the size difference, the two bubbles changes from being in-phase to out-of-phase over time. The effect of one bubble on its in-phase partner resembles that of a rigid wall, but the effect becomes similar to a free surface when the bubbles are out-of-phase. Thus, contrary to the previous case, the bubble migration was no longer monotonically toward each other. The center of bubble 2 reciprocated during multiple oscillation cycles, first migrating away from bubble 1 and then toward. The present theory is capable of reproducing such a phenomenon. The prediction matches the experiment in bubble radius variation and captures the main features of the bubble migration, as shown in Figs. 5(b) and 5(c). Possible reasons for the deviation may include the non-spherical bubble behavior, such as jetting, and the heat and mass transfer over multiple oscillation cycles.

Next, we carried out an experiment on the interaction of three electric spark-generated bubbles and compared the dynamics to theoretical results, as shown in Fig. 6, which introduces higher complexity. The three bubbles are positioned at the vertices of a triangle with the top two bubbles being 118 mm apart and the lower bubble 115 mm from the other two. Since the behavior of the upper two bubbles is almost symmetrical to each other, we only compare the upper left bubble and the lower one, marked as Bubble 1 and Bubble 2, respectively. The initial radii of bubbles 1 and 2 for theoretical calculations are $R01=2.22$ and $R02=2.48$ mm, and the initial oscillation speeds $R\u030701=180$ and $R\u030702=150$ m/s, respectively. The initial internal pressure of both bubbles is 1.2 MPa. While oscillating in volume, the bubbles keep migrating toward each other, as shown in Fig. 6(a). The radius and the migration of the bubbles in the vertical direction in the experiment are well reproduced by the theoretical results, as shown in Figs. 6(b) and 6(c).

Further, the dynamics of a bubble cluster are modeled and compared with an experiment. The bubble cluster consisted of sixteen bubbles positioned uniformly on a 4 × 4 grid with an interval of 40 mm. All the bubbles were generated by the underwater discharge at the same instant. Considering the symmetry of the cluster, we only analyze the dynamics of three representative bubbles among them, marked as bubble 1, bubble 2, and bubble 3 as shown in Fig. 7(a). The bubble oscillation can be prolonged by neighbor bubbles with the same oscillation phase. Thus, bubble 3 surrounded by its neighbors from all four directions had the greatest oscillating period. On the contrary, bubble 1 located at the external corner was less affected by the other bubbles, which led to the shortest oscillation period. The bubbles away from the cluster center collapsed earlier and their re-expansion accelerated the final collapse of bubble 3. All the bubbles migrated toward the cluster center due to bubble interaction. In the theoretical calculation, the initial internal pressure of all the bubbles is 1.2 MPa. The initial radii of the three selected bubbles are as follows: $R01=1.41$, $R02=1.36$, and $R03=1.18$ mm, respectively. The initial oscillation speeds are $R\u030701=R\u030702=180$ and $R\u030703=220$ m/s, respectively. To simplify theoretical calculations, the initial radii and speeds of other bubbles are identical to the three bubbles in a symmetrical manner, which can be considered as reasonable given the slight differences in bubble behavior in most of the first cycle. The radius, period, and vertical migration from the simulation with the present theory are in good agreement with the experimental results, as shown in Figs. 7(b) and 7(c), demonstrating the theory's capability in dealing with complex bubble cluster dynamics.

### D. Bubble dynamics near a boundary

In this section, we validate our theory in the scenario of bubble oscillation near a rigid boundary. A cavitation bubble was generated at the atmospheric pressure by a Q-Switched Nd:YAG laser breakdown beneath a rigid boundary at a standoff distance of $db=2.02Rmax$, with the maximum bubble radius $Rmax$ = 0.768 mm. The experimental setup is shown in Fig. 8(a), and high-speed images of the laser bubble oscillating while migrating toward the rigid boundary at the top are shown in Fig. 8(b). For the theoretical calculation, we obtain the initial conditions of the bubble as $P0=1.2$ MPa, $R0=0.121$ mm, and $R\u03070=130$ m/s in the same way as stated before. The results are plotted against the experimental data in Figs. 8(c) and 8(d). When a bubble oscillates near a rigid wall, the flow is retarded by the presence of the wall. Thus, the bubble oscillation period increases with a decreasing standoff parameter. The theories of Keller and Gilmore do not allow for the boundary effect, which may lead to a discrepancy in the oscillation period and the maximum radius of the second cycle when compared to the experiment. The present theory has considered the boundary effect, and the calculated bubble radius and oscillation period are closer to the experimental data for multiple cycles, as shown in Fig. 8(c). During the bubble–wall interaction, the pressure between the bubble and the wall is smaller than that on the opposite side; thus, the bubble is pushed toward the wall, which is usually accompanied by a jet pointing to the wall. Our theory also takes into account the bubble migration induced by the boundary effect and the calculated migration curve is consistent with the experiment, as shown in Fig. 8(d). To further validate the present theory, we derive the Rayleigh time for a bubble oscillating near a rigid boundary by integrating the modified Rayleigh–Plesset equation

in time, which is similar to Refs. 85 and 86. Then, the Rayleigh-like bubble oscillation period near a rigid wall can be expressed as

in which *η* is the auxiliary variable. In this case, the Rayleigh-like period is 0.153 ms. Comparatively, the bubble period of the present theory is 0.156 ms, which is closer to the experimental result of 0.155 ms due to the consideration of multiple factors, e.g., the fluid compressibility, the internal pressure of the bubble, the surface tension, and the viscosity.

Additionally, we also compared the theoretical results of bubble dynamics near a free surface to the experimental data of a small-charge underwater explosion bubble. The experiment setup is shown schematically in Fig. 9(a). In the experiment, an explosive charge equivalent to 1.125 g PBX9501 (1.3 g TNT) was detonated at the center of the aforementioned tank with a water depth of 0.4 m to generate a gas bubble with a maximum radius $Rmax=0.169$ m. Due to the repelling effect of the free surface, the bubble migrated downward while oscillating, as shown in Fig. 9(c). The buoyancy parameter is $\delta =0.126$. The initial conditions of the bubble are $R0=43.57$ mm, $R\u03070=61.83$ m/s, and $P0=303.94$ kPa, which are obtained with the method given in Appendix B. The bubble dynamics are calculated using different models, and the results are compared to that of the experiments in Figs. 9(b) and 9(d). Compared to an infinite flow domain, the inertia of the liquid is reduced when a free surface is present; thus, the bubble can oscillate faster and have a smaller period. The pressure of the bubble gas is much smaller than the atmospheric pressure during most of the lifetime of the bubble; thus, the pressure gradient between the bubble and the free surface points upward and the bubble is repelled by the free surface, which is usually accompanied by a downward jet.

The Keller equation does not consider the boundary effect and bubble migration, while the Geers–Hunter model allows for migration but neglects the boundary effect. Consequently, there are deviations in the bubble radius and period to the experimental results, and a key feature that the bubble is repelled by the free surface is not reproduced. Compared to that, with the present model, the calculated bubble radius variation is in better agreement with the experiment and the migration of the bubble away from the free surface is also predicted. Similar to Eq. (55), it is easy to find that the Rayleigh-like bubble pulsation period near the free surface can be expressed as

Compared with the bubble periods (0.0276 and 0.0283 s) calculated by Eq. (56) and from the experiment, the bubble period calculated by the present theory, 0.0277 s, is also more accurate. In addition, Fig. 9(e) depicts the comparison between the theoretical and the experimental results of the pressure pulse induced by the first collapse of the bubble. The pressure is measured at the same depth as the charge center and 0.7 m away in the horizontal direction. The pressure peak and the pulse width calculated by the present theory match the experimental data. The negative pressure induced by the rarefaction wave that is formed by the reflection of the pressure pulse at the free surface is also successfully captured by the present theory.

We proceeded to compare the experimental data with the theoretical results obtained for multiple bubbles located near a free surface. In Fig. 10, we observe the progression of three bubbles in a linear arrangement near the free surface and compare the bubble radius and vertical migration between our theory and the experiment. At inception, all bubbles were positioned 47.6 mm away from the free surface. The left and right bubbles were initially located 52.8 and 55.8 mm away from the central bubble, respectively, which resulted in their almost symmetric configuration with respect to the central bubble. We focused our analysis on the dynamics of the left and central bubbles, which are identified as bubble 1 and bubble 2, respectively. Bubble 1 generated an oblique liquid jet and migrated in the same direction as bubble 2 under the combined action of the free surface and bubble 2. Bubble 2 experienced horizontal elongation during collapse, with both of its sides obstructed by the contracting bubbles on its left and right, while simultaneously migrating downward, affected by the free surface. The initial conditions for bubbles 1 and 2 in the theoretical calculation were $R01=2.92$ mm, $R02=2.1$ mm, $R\u030701=100$ m/s, $R\u030702=180$ m/s, and $P01=P02=1.2$ MPa. The present theory was found to accurately reproduce the time histories of bubble radius and migration in the experiment, as presented in Figs. 10(b) and 10(c).

In the last experimental verification, we applied our theory to model high-pressure air-gun bubble dynamics and compared the theoretical results against experimental data. When triggered, an air-gun rapidly released compressed air within a few milliseconds, forming a high-pressure bubble that expands and collapses. The experiments were conducted in a deep reservoir, and a schematic diagram of the experiment setup is given in Fig. 11(a). An air-gun was placed 1.8 m below the water surface. The pressure in the flow field was measured by a pressure sensor at the same depth with a distance of 2.0 m from the air-gun. For the theoretical analysis, we use an air-release model to assess the mass flow from the air-gun chamber to the bubble and the ideal gas law to update the internal pressure of the bubble. Details of the implementation can be found in previous works.^{122} Here, the initial internal pressure and volume of the chamber are $P0=10$ MPa, $V0=6227$ cm^{3}, respectively. The results of the present theory, as well as that of the Keller equation and the Gilmore equation, are compared to the experimental data in Figs. 11(b) and 11(c). The present theory is capable of considering bubble migration and the effect of the free surface and, therefore, produces a result in better agreement with the experimentally obtained pressure profile.

## VI. DISCUSSION

In this part, we discuss the applicability of our theory and then demonstrate its competence via implementation in the prediction of the complicated dynamic behavior of two interacting bubbles divergent in phase and energy, with new physics revealed.

### A. Applicability of the present theory

Here, we analyze and discuss the applicability of the present theory in an extended parameter space through a quantitative comparison to numerical simulations. A number of numerical methods have been successfully applied to bubble dynamics, such as boundary integral method (BIM)/boundary element method (BEM),^{119,122,123} finite difference method (FDM)/finite volume method (FVM)/Eulerian finite element method (EFEM),^{124,125} smooth partial hydrodynamics (SPH),^{126} and lattice-Boltzmann method (LBM).^{127,128} Here, we choose three representative methods, namely, BIM, EFEM, and SPH for the comparison. The effect of gravity, surface tension, and viscosity is ignored in both the theoretical and the numerical calculations to highlight the effect of boundaries on bubble dynamics. The initial conditions of the bubble are set as $R0=0.173Rmax,\u2009P0=100P\u221e,\u2009\gamma =1.4$, and $R\u03070=0$. Figure 12(a) depicts the theoretical and the numerical results of bubble radius variations with the standoff distance to the rigid boundary being $d=2Rmax$. Also, the results neglecting the boundary effect (i.e., in a free field) are presented. The boundary effect has a noticeable influence on the oscillation period of the bubble, and the results of our theory are in good consistency with the numerical simulation. Hereafter, symbols with “*” are non-dimensional physical quantities with the length, pressure, and density scales being $Rmax,\u2009P\u221e$ and *ρ*, respectively. The scales for the other quantities are calculated accordingly, e.g., $P\u221e/\rho $ for velocity and $Rmax\rho /P\u221e$ for time.

To further analyze the discrepancy between the two approaches quantitatively and evaluate the applicability of the present theory, we compare the results of the present theory to that of the BIM simulations within the first oscillation cycle of a bubble with varied standoff distances to a rigid boundary or a free surface, as depicted in Fig. 12(b). For scenarios with the rigid boundary, the discrepancy in period remains within 2%, except for when the standoff distance is smaller than $1.3Rmax$. For scenarios with the free surface, the discrepancy is higher than that for the rigid boundary at small standoff distances possibly due to nonlinear free surface motions (such as a water spike). In spite of that, the discrepancy drops to about 2% as soon as the standoff distance increases to $1.4Rmax$.

### B. Dynamics of interacting bubbles

Next, we further apply the present theory to complex scenarios of bubble dynamics to reveal new underlying physics. The interaction between two bubbles is fundamental to the study concerning the dynamics of multiple bubbles, which is a problem frequently encountered in a wide range of applications and troubles the researchers with sophisticated variations in bubble behavior influenced by many factors. Therefore, we have chosen the dynamics of a two-bubble system with different initial bubble energy ratios and initiation times as the objective. A sketch of the two-bubble system is shown in Fig. 13. Herein, we take the maximum radius of bubble 1 as the length scale and ignore the gravity. The Mach number $Ma=P\u221e/\rho /C$ is taken as 0.013. The initial conditions of the two bubbles are given by the same $P0=100P\u221e$ and $R\u03070=0$. The initial radius of bubble 1 is $R0=0.175Rmax$ and that of bubble 2 is adjusted to tune the energy ratio between the two bubbles. The parameters considered are defined as follows. The internal energy of a bubble at initiation is estimated as $Ei=P0,iV0,i/(\gamma \u22121)$, where *i* = 1 or 2 indicates physical quantities of bubble 1 or 2, respectively. We then denote $E*=E1/E2$ as the bubble energy ratio. The initiation of bubble 1 is marked as *t* = 0. The phase difference between the bubbles is defined as $\beta =\Delta t*/Tf*$, where $\Delta t*$ is the non-dimensional initiation time of bubble 2 and $Tf*$ is the non-dimensional first oscillation period of bubble 1 in a free field without bubble 2.

The parameter space is set as follows. The energy ratio $E*$ ranges from 1 to 10 and the phase difference *β* from 0 to 1. The initial distance between the two bubbles is set as a constant (five times the $Rmax$ of bubble 1 in a free field). We focus on three properties that reflect the main bubble dynamics, namely, the migration, the oscillation period, and the scaled internal pressure. The results are presented in Fig. 14 in an $E*\u2212\beta $ parameter space. The migration is represented by the average relative speed of the two bubbles as depicted in Fig. 14(a). The positive values represent the two bubbles eventually repelling each other and negative values represent attracting, with three demarcation lines $L1\u2212L3$. The isoline *L*_{1} indicates that, at $E*=1$ where the two bubbles are initiated with the same energy, the migration alters from attracting to repelling when *β* exceeds approximately 0.1. The threshold in *β* gradually grows with $E*$ and reaches approximately 0.55 at $E*=10$, which means that a larger phase difference is required for the bubbles to repel each other when the bubble energy difference increases. *L*_{2} resides within $0.9<\beta <1.0$ for different energy ratios, indicating a transition from repelling to attracting when the initiation of bubble 2 is delayed until the end of the first contraction stage of bubble 1. In the region indicated by *L*_{3} where $E*$ exceeds 6.6 and *β* is around 0.1–0.3, the migration turns from attracting to repelling again. In fact, in an extended area containing *L*_{3}, the bubble-center distance reciprocates around its initial value during bubble oscillation. An example of the corresponding complex relative migration–time curve is depicted in Fig. 14(b). An explanation is that with the large energy difference, the migration of bubble 2 is dominated by the surrounding fluid motion set into oscillation by bubble 1. In addition, the scenarios with the highest repelling amplitude occur around $\beta =0.5$ and $E*=1$. Here, the effect of bubble 1 on bubble 2 is similar to a free surface, and thus, the latter migrates away from the former. The scenarios with the highest attraction amplitude occur around *β* = 0 and $E*=1$; i.e., the two-bubble system is symmetrical. Here, the dynamics of a bubble mimics that near a rigid boundary, and hence, the migration toward each other is prominent.

The variation of the non-dimensional first period of bubble 1, $T*$, is shown in Fig. 14(c), which reflects the average energy of bubble 1 in the period in the ambient flow induced by bubble 2. On the isolines *L*_{1} and *L*_{2}, $T*$ equals the first period of bubble 1 in a free field without other bubbles. The pattern of Fig. 14(c) is similar to Fig. 14(a) but with opposite signs. The effects of bubble 2 on the pulsation period of bubble 1 can also be interpreted by analogy with the effects of a boundary. $T*$ reaches a minimum of 1.877 near $E*=1,\u2009\beta =0.45$, in which case bubble 2 initiates as bubble 1 approaches the maximum volume. In this scenario, the expansion of bubble 2 accelerates the contraction of bubble 1, similar to the effect of a free surface. From here, increasing $E*$ leads to higher $T*$ due to the weaker influence of bubble 2 with lower energy. Compared to $E*$, altering *β* results in more obvious changes in $T*$. As *β* approaches 1.0, $T*$ grows until reaching $Tf*=1.986$; as *β* approaches 0, isoline *L*_{1} is crossed and $T*$ exceeds $Tf*$ by a maximum of 5%. This is because the bubbles are in-phase and hinder each other, leading to a prolonged period similar to the bubble pulsation near a rigid boundary. In the region around *L*_{2}, $T*$ varies in a very limited range and approximates $Tf*$, since bubble 2 hardly affects the period of bubble 1 when $E*$ is large.

The scaled internal bubble pressure is reflected in Fig. 14(d) by a nephogram of $PR*=Pmax*R*$, where $Pmax*$ is the non-dimensional maximum pressure inside bubble 1 and $R*$ the non-dimensional radius. $PR*$ equals the pressure peak a bubble induces at a unit distance and implies the energy of bubble 1 at its first minimum volume. The three demarcation lines are the isolines for $PR*$ being equal to what bubble 1 induces in a free field, $PRf*$, when bubble 2 is absent. If bubble 1 is expanding when the pressure wave emitted by bubble 2 arrives, the pressure wave does negative work to bubble 1 and decelerates the expansion, in which case a smaller $PR*$ is expected. By contrast, if bubble 1 is collapsing, the pressure wave will accelerate the collapse and $PR*$ will increase. The energy transmission rate increases with $R2R\u0307$, which is measured on bubble 1 at the incidence of the pressure pulse induced by bubble 2. The negative value indicates the transmission from bubble 2 to bubble 1, while the positive value indicates the reversed transmission. The variation in $PR*$ with *β* reduces with the increase in $E*$, since a higher energy ratio means a smaller size of bubble 2 and lower influence on bubble 1. For the region above *L*_{1}, $PR*$ remains greater than $PRf*$. This indicates energy transmission from bubble 2 to bubble 1. The greatest $PR*$, i.e., the highest transmission, occurs in a region around $\beta =0.85,\u2009E*=1$, with $PR*=1.56PRf*$. This is because that the highest $\u2212R2R\u0307$ is found right in the above region. Therefore, the expansion of bubble 2 intensifies the collapse of bubble 1, leading to the high $PR*$. The required minimum *β* for $PR*>PRf*$ starts as about 0.5 at $E*=1$ and increases almost linearly until reaching about 0.75 at $E*=10$, as reflected by the demarcation line *L*_{1}. For the region below *L*_{1}, the pressure wave of bubble 2 arrives at bubble 1 at its expansion phase and does negative work to it. Thus, $PR*$ falls below $PRf*$, indicating energy transmission from bubble 1 to bubble 2. The minimum $PR*$ is about 4.5 and appears around $\beta =0.1$ and $E*=1$ when $R2R\u0307$ arrives at the maximum. Additionally, it is worth noticing that as *β* approaches 0, $PR*$ may also exceed $PRf*$ for certain ranges of $E*$, corresponding to the region between *L*_{2} and *L*_{3}. Here, the collapse of bubble 2 occurs at the late contraction stage of bubble 1. The pressure wave induced by the collapse of bubble 2 serves as an energy input to bubble 1 and that explains why $PR*$ exceeds $PRf*$.

To further discuss and demonstrate the capability of the present theory in the prediction of complex pressure waves, energy transition, and complex physical laws of bubble dynamics, we have chosen two typical scenarios from the above-mentioned parameter space, i.e., ($\beta ,E$) = (0, 2) and (0.5, 1.6), and modeled the corresponding bubble dynamics using the present theory. In the first scenario, the bubble inceptions are simultaneous, while in the second, bubble 2 is produced when bubble 1 is at the maximum expansion. We also compared the results to that of the numerical simulation using the EFEM. The radius variation and the bubble-induced pressure profiles are compared to the results of numerical simulations using the EFEM in 15. The test point where the pressure profiles are measured is indicated in Fig. 13, which is at the same depth as the bubble pair with a non-dimensional distance of 6.0 from the center of the configuration. The results of a single-bubble scenario where bubble 2 is removed are also added in the figure to highlight the influence of bubble 2.

The radius variation and pressure profiles for the first scenario with the simultaneously generated bubbles are given in Figs. 15(a) and 15(b), respectively. One may find that the maximum radius of bubble 1 in its second oscillation cycle has increased when bubble 2 is present. Using the bubble energy formulas,^{129} which is based on the maximum bubble radius, we calculated the dimensionless energy of the two bubbles and find that the maximum energy of bubble 1 increased significantly from 4.12 when oscillating alone to 4.68 when interacting with bubble 2. It indicates that bubble 1, in its first two oscillation cycles, absorbs the energy of bubble 2 during the bubble interaction. This is also reflected in the pressure profile, i.e., the first bubble period is prolonged, and the bubble collapse pressure is increased than that when bubble 2 is absent. Interestingly, the maximum radius of bubble 2 in different cycles also fluctuates. The maximum radius in the fourth cycle is noticeably higher than those in the second and the third. The dimensionless energy of bubble 2 in the fourth cycle reached 1.99, which is significantly higher than 1.54 and 1.55 of the two preceding cycles, indicating a gain in the energy of bubble 2, rather than simply damping. The energy gain is from the positive work done by the pressure waves induced by the collapse of bubble 1 occurring in the third and the fourth cycle of bubble 2. This is also reflected in the pressure profile; i.e., the pressure peak induced by the third collapse of bubble 2 is remarkably higher than that of the previous two collapses and comparable to that of the second collapse of bubble 1. The above indicates complex energy transmission between the two interacting bubbles.

The results of the second scenario where bubble 2 is generated when bubble 1 reaches its maximum expansion are shown in Figs. 15(c) and 15(d). The maximum radii for both bubbles are damping during oscillations. The peaks on the pressure profile fluctuate, but the peaks induced by the same bubble upon each collapse keep reducing. Intriguingly, a negative pressure peak is captured at about $t*$ = 1 in the theoretical result. It corresponds to the rarefaction wave generated after the pressure pulse induced by the initial expansion of bubble 2 reflecting on the surface of bubble 1. Such ghost reflections are also manifested in the subsequent process of the bubble interaction. In general, good agreement is found between the theoretical and the numerical (EFEM) results of the bubble radius variation and the complex bubble-induced pressure fluctuation in the flow field of both scenarios, which proves the capability of the present theory in explaining the sophisticated bubble dynamics.

## VII. CONCLUSION

Initiating from basic mathematical principles and physical equations, we established a novel and extendable theory for oscillating bubble dynamics. It can simultaneously consider the complex physical factors including boundaries, bubble interaction, ambient flow field, gravity, bubble migration, fluid compressibility, viscosity, and surface tension. At the same time, it retains a unified and elegant mathematical form when dealing with various bubble-related problems. The theory achieves the unification of classical bubble theories such as the Rayleigh–Plesset equation, the Keller–Miksis equation, and the Gilmore equation. We systematically compared the results of the present theory with those of the previous ones, the experiments, and the numerical simulations of bubble dynamics with a variety in bubble scales, sources, boundaries, and ambient conditions. It shows that the present theory has higher accuracy and applicability compared to the classical ones. It can accurately predict the key bubble dynamics features including bubble oscillation, migration, and collapse-induced pressure pulses under different conditions and shows potential in exploring more sophisticated physics and mechanisms behind bubble dynamics such as bubble interaction, coupling of bubble-induced pressure waves, and inter-bubble energy transfer. The theory presented in this work may provide new references for future explorations in bubble dynamics, cavitation, and other multiphase flow-related problems.

## ACKNOWLEDGMENTS

This work was funded by the National Natural Science Foundation of China (Nos. 51925904 and 52088102), the National Key R&D Program of China (Nos. 2022YFC2803500 and 2018YFC0308900), Finance Science and Technology Project of Hainan Province (No. ZDKJ2021020), the Heilongjiang Provincial Natural Science Foundation of China (No. YQ2022E017), and the Xplore Prize.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**A-Man Zhang:** Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (equal); Methodology (equal); Software (equal); Supervision (lead); Validation (equal); Visualization (equal); Writing – review & editing (equal). **Shi-Min Li:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Pu Cui:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Shuai Li:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – review & editing (equal). **Yunlong Liu:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request. The source codes used in the paper can be found at https://github.com/fslab-heu/unified-bubble-theory.git.

### APPENDIX A: SOLUTION OF A MOVING SINGULARITY FOR WAVE EQUATION

Assuming that there is a source at the origin with a strength of *f*_{0}(*t*), the solution of the wave equation can be written as

Making use of the linearity of the wave equation, the solution of the wave equation considering the movement of the source can be obtained by superposition of a set of $\phi f0$ as

where $rt$ is the vector pointing from the location of the source at $t\u2212|rt|/C$ to ** r**. It is not easy to obtain an explicit analytic expression for $rt$ unless the source is moving with a constant velocity. Once Eq. (A2) is obtained, the solution for a moving dipole with a strength of

*q*(

*t*) can be calculated subsequently with the following limit

### APPENDIX B: INITIAL CONDITIONS FOR UNDERWATER EXPLOSION BUBBLE

##### 1. Bubble formation and near-field shock wave propagation

In this part, we present an approach to model the initial phase of an underwater explosion bubble. Upon the formation of the explosion bubble, the fluid pressure around the bubble is very high and a strong shock wave emits from the bubble. The basic assumption of the wave equation is thus violated. Therefore, we solve the Euler equation in the conservation law form to treat this phase until the pressure and velocity are low enough for the initial condition of the unified theory for bubble dynamics. We assume that there is a gas bubble located in still water starting to expand due to high internal pressure. The thermal conduction and viscous effect are ignored. The fluid motion is described by the one-dimensional spherically symmetric Euler equations

where $E=\rho e+12\rho u2$ is the total energy and *e* is the specific internal energy per unit mass. The right-hand sides of the above equations represent the geometric source terms in a spherical symmetric coordinate system.

After the detonation, a shock wave forms at the explosive charge surface and propagates into the water. The gaseous product of the explosive forms a bubble that expands rapidly. The radii of the shock front and the bubble surface are denoted as $Rs$ and *R*, respectively. Thus, at $r>Rs$, the fluid is undisturbed and at *r* < *R* is the gaseous product inside the bubble. These two parts provide the boundary conditions of the fluid within $R<r<Rs$. In this work, we focus on the fluid between *R* and $Rs$ that is smooth and solvable with high-order methods, and avoid the difficulties in treating the multi-medium interface and the discontinuity at the shock front. The physical coordinate *r* is mapped linearly into the coordinate for the numerical solution *ζ* as

which is also illustrated in Fig. 16, where $L=Rs\u2212R$ is the length of the computational domain. Then, we have

We denote $u\u0302=d\zeta dt$ as the fluid velocity observed in the coordinate system of *ζ*. Thus,

We denote $\u2202\u0302/\u2202\u0302t=\u2202/\u2202t+(\zeta L\u0307+R\u0307)\u2202/\u2202r$ as the derivative with respect to *t* when *ζ* is fixed. Then, the Euler equations can be rearranged and expressed in the following compact form:

where

Note that an extra source term emerges in the scaled coordinate system. Then, the hydrodynamic pressure at *r* can be given by

Here, *p _{l}* is the pressure calculated with the equation of state of the surrounding fluid.

The Euler equation system is closed by the equation of states of the fluids. Typically, the internal gas produced by explosive is modeled with the JWL equation,^{130} which reads

where $Rc$ is the radius of the explosive charge and $A,B,R1,R2$, and *ω* are material-related constants.

If the surrounding fluid is liquid, we model it by the Tammann EoS

The initial conditions are calculated as follows. Initially, it is assumed that a high-pressure gas bubble is placed in a still fluid field. The internal gas is considered as still at the beginning for simplicity. Thus, a Riemann problem is formulated at the interface between the internal gas and the surrounding fluid. Here, we define the starting time (*t* = 0) of the explosion bubble dynamics as when the detonation wave reaches the explosive charge surface. We assume that the charge has fully converted into gaseous products and thus formed a bubble at this moment. We start the simulation from a short time *ε* after it. The fluid quantities are assumed to be constant between the bubble surface and the shock front. At $t=\epsilon $, the fluid pressure of the water between *R* and $Rs$ equals the inner pressure of the bubble, i.e., $Pg$. The Mach number for the shock wave $M=R\u0307s/C\u221e$ is calculated with the Hugoniot condition

Then, the initial conditions for the fluid density and velocity are given by

and

respectively. Here, the subscript $\u221e$ indicates the physical quantities of the undisturbed water. The initial bubble radius *R* and the shock front radius $Rs$ at $t=\u03f5$ are given by

where $Rc$ is the radius of the explosive charge. In this paper, the value of *ϵ* is determined by letting $L=Rc/100$ at $t=\u03f5$.

Within this framework, we can obtain the boundary conditions. For the boundary representing the material contact (the left boundary), we have $u\u0302=0$ and $p=Pg$. Thus, the flux is given by

For the boundary representing the shock front (the right boundary), the flux is given by

where $E0=(P\u221e+\gamma Pw)/(\gamma \u22121)$ is the total energy of the undisturbed fluid. Note that when the dynamic effect of the internal gas is considered, the Euler equation for $0<r<R$ should also be mapped to a scaled space and solved numerically. In such cases, its right boundary conditions can be connected to the conditions on the left boundary of the external fluid with a consistent flux. Conditions for the right boundary can be connected to the conditions for the left boundary of the external fluid with a consistent flux.

##### 2. Initial condition for the unified bubble theory

In this work, we start the simulation with the unified bubble theory at the moment denoted by $te$ when $Rs$ reaches 10*R*. The initial condition at $t=te$ is derived from the simulation of the initial phase prior to $te$ using the equations above. Let us denote $t=te$ as the moment when the near-field simulation ends. Then, it is straightforward that the initial condition for the unified bubble theory is chosen as

However, the energy of the oscillating bubble will be excessive, and the maximum radius overestimated if the above is adopted directly. The key reason is that the energy of the shock wave, which should have dissipated, has been included. In this paper, we prevent this by calculating $R\u03070$ with the kinetic energy of the fluid field excluding the shock wave region

where *a* is a constant smaller than 1 to avoid the shock wave region in the integration. In this paper, *a* is taken as 0.65.

##### 3. Shock wave far-field propagation

In addition to the bubble, the shock wave is also a major component of the physical process of an underwater explosion. The method to calculate near-field shock wave has been given in Appendix B 1. Nevertheless, when the far-field shock wave pressure is considered, this method can be very expensive because of the time step limitation related to sound speed. Thus, we propose a new method based on the weakly compressible assumption for quick estimation of far-field shock wave pressure. For far-field flow, we modeled the fluid with the Tait equation

Thus, the sound speed *C* satisfies

Expanding *C* at $p=P\u221e$ and denoting $\delta p=p\u2212P\u221e$, we have

where $Kc=(\gamma \u22121)/(2\rho \u221eC\u221e)$.

The energy of the shock wave denoted by *E* can be evaluated by $E=Ek+Ep\u2248(\rho \u221eC\u221e2)\u22121(\delta p)2$ in the far field. The shock wave energy is transported at the sound speed *C*. Thus, the conservation of the shock wave energy can be expressed as

Considering that the change of wavelength of the shock wave is slow in the far field, we map *r* to *ζ* with the following relation:

where $Re$ is the left end of the computational domain and $Rs\u2212Re$ is the length of the computational domain, which is a constant and denoted by *L _{p}*. Approximating

*u*by $(\rho C)\u22121\delta p$, we have

Let us replace *E* with $(\rho \u221eC\u221e2)\u22121p2$, and *C* with B22 neglecting $O(\delta p)2$, then Eq. (B25) can be expressed in a conservation law form as

The above equation describes the evolution of a pressure wave in a spherical symmetric problem. The second term on the left-hand side represents the linear transport of the wave energy. The third term is a nonlinear flux to reduce the slope of the pressure distribution curve after the shock front. The source term on the right-hand side incorporates the effect of the increase in the area of the shock front propagating outward. We can rewrite the above equation as

where the flux $G=(C\u221e\u2212R\u0307s)\delta p+34Kc(\delta p)2$ and the source $S=\u2212C\delta p/r$. Then, the pressure history at *r* can be expressed as

The initial conditions of the far-field simulation are determined by the final state of the near-field simulation. In this work, we take the far-field as $L>10R$, where we assume that the shock wave pressure is low enough for the linear assumption to become valid. Therefore, we start calculation for the far-field using the approach in this section instead of the one based on the Euler equation introduced in Appendix B 1. The L2 projection is used to calculate the pressure solution in the new polynomial space from the solution of the Euler equation. The boundary conditions are given by the following boundary flux:

##### 4. Discontinuous Galerkin method

The discontinuous Galerkin method^{131,132} is adopted to solve the near-field flow (the initial phase of bubble expansion) and the far-field shock wave propagation. The computational domain is discretized into *N _{dg}* elements, and the $\u2113$ th element is bounded by $\zeta \u2113$ and $\zeta \u2113+1$.

First, we take the control equation (B7) of the near-field flow in Appendix B 1 as an example. We approximate the solution by

where *β* is a complete set of the polynomial space *P ^{K}*. For element $\u2113$, the control equation (B7) of the near-field shock wave is then given in the Galerkin form

where $F\u0303$ is the flux at the boundary to be calculated according to specific boundary conditions. The matrix form is given by

where $Mij=\u222b\zeta \u2113\zeta \u2113+1\beta i\beta jd\zeta $ and $Bi=\u222b\zeta \u2113\zeta \u2113+1(F\u2202\beta i\u2202\zeta +S\beta i)d\zeta \u2212F\u0303\beta i|\zeta =\zeta \u2113\zeta =\zeta \u2113+1$. The numerical flux $F\u0303$ at $\zeta \u2113$ with $2\u2264\u2113\u2264Ndg$ is calculated with the HLLC^{133} Riemann solver.

An explicit time marching scheme can be used to update $K\u0307$, e.g., the Runge–Kutta method. In this paper, one element with the *P*^{15} space was used first until $L>2R$. Then, 20 elements with the *P*^{3} space were used until $L>10R$. If the air blast problem is considered, the Euler equation describing the explosion product should also be solved and the relay point between the near-field and far-field solvers should be later than the above for an underwater explosion problem because the compressibility of the air is much greater than that of the water.

Subsequently, one may calculate the initial conditions for the unified bubble theory based on Appendix B 2 to analyze the bubble motion or start the far-field simulation based on Appendix B 3 to predict the shock wave at a far distance. When solving the far-field shock wave propagation, Eq. (B27) is solved with the same approach. Only one element and *P*^{8} space are used since the profile of the shock wave is smooth enough. If the initial conditions of multi-bubble are going to be obtained at the same time and the interaction of their shock waves are considered, a full three-dimensional model should be used.

##### 5. Validation

Here, we validate the above methods with experimental measurements. Three underwater explosion experiments were carried out in a free field with explosives equivalent to 2.0 g, 7.0 kg, and 80.0 kg TNT, respectively. We captured the shock wave front and the initial bubble expansion in the first experiment using the high-speed camera, to which we compared our solution. A good match was obtained as depicted in Fig. 17. We also measured the detonation shock wave pressure histories in the second and the third experiments using a piezoelectric pressure sensor at distances of 4.0 and 30.0 m, respectively. The comparison between the experimental and the simulated results using the near-field and far-field methods is illustrated in Fig. 18, where a good agreement is achieved.