A unified theory for bubble dynamics

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 insights into inter-bubble energy transfer and coupling of bubble-induced pressure waves.

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 researches.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.Lord 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][88][89][90][91][92][93][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 other researchers [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 co-workers [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  2 /4 (surfaceaveraged pressure, SAP, and  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 causation.
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 airgun-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][115][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 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 section II.The theory is extended to multiple-bubble interaction in section III and bubble dynamics near boundaries in section 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 are presented in section V.A discussion on the applicability is given in section 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 section VII.

II. DERIVATION OF THE UNIFIED THEORY FOR BUBBLE DYNAMICS
A. Basic equations of the theory Firstly, 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 , , u, and τ , respectively.The conservation equations for mass and momentum are and where F is the body force.Meanwhile, the sound speed  in the external fluid field is related to the pressure  and the fluid density  as Substituting equations ( 1) and (3) into the derivative of equation ( 2) with respect to time yields The sound speed  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 equation ( 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  that satisfies u = ∇.Substituting it into equation ( 4) and integrating from infinity to the current location while neglecting the variation of the body force and the terms containing  −2 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  − to the bubble center. = 0 points to the direction of the bubble migration velocity relative to the ambient flow denoted by v. v = e = v m − u a , where e is a unit vector along the direction of v, v m is the migration velocity, and u a is the ambient flow velocity.u a = u B + u E , where u B is the velocity induced by boundaries or other bubbles and u E represents extra velocities including the velocity of the incoming background flow u ∞ .
In this work, the bubble is assumed to oscillate and migrate while maintaining a spherical shape.We construct a solution of equation (5) with two singularities, i.e., a source and a dipole, the strength of which are denoted by  and , respectively.The velocity potential at an arbitrary location r and time  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 r is close to the bubble and considering that the migration velocity  is small compared with , 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 r.Thus, the velocity potential at r = (, ) and  is expressed as the summition of equations (A2) and (A3) after simplification, and reads in which  ′ denotes the derivative of  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  ′ denotes the derivative of  with respect to its argument.The kinetic boundary condition on the bubble surface for the oscillation and the migration can be expressed as and d d respectively, where  is the volume occupied by the bubble and  is its boundary, i.e., the bubble surface. denotes the time derivative of .Let's denote  () = −  ( − /) and () = ( − /).Substituting equation ( 7) into equations ( 9) and ( 10), we have and respectively.Employing the perturbation method on equation (12), we arrive at with the terms of the order  −2 ignored.Then, take the derivative of equation ( 11) with respect to , 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 /4, where  is the bubble volume acceleration, meanwhile, the left-hand side is equivalent to the corresponding driving force.d/d can be determined according to the physical problem considered, enabling equation (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 nonspherical 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 equation (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  g , the fluid pressure on the outer bubble surface  b , surface tension, viscous stress, and other additional terms.Here we express the equilibrium condition as where  is the surface tension and  the liquid viscosity. g 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,  g is composed of the non-condensable gas pressure and the saturated vapor pressure  v as where  0 is the initial pressure of the non-condensable gas,  the polytropic exponent and  0 the initial bubble radius.
To obtain d/d in equation (15), we introduce the dynamic boundary condition expressed by the Bernoulli equation in the moving coordinate system where v • u = (cos, −sin) • (  ,   ) and −1 d is the enthalpy difference at the bubble surface.Here,  a represents the ambient pressure at the bubble center. a =  B +  E where  B is the pressure induced by boundaries or other bubbles and  E represents extra pressures including the hydrostatic pressure  ∞ and acoustic pressures  A .With the motion of the fluid being an adiabatic process,  resolves into where  = ( b −  a )/( 2 ).
Note that the terms containing  in equation ( 18) will lead to non-spherical bubble deformation.To comply with the spherical assumption while taking into consideration their averaged effects, we integrate equation ( 18) over the bubble surface to eliminate  and obtains Finally, substitute the above into equation ( 15) and 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 equation (15) as where F =  2 (  2 /2 +  2 /4 + ) is to be determined according to specific conditions of the problem under consideration. 2 /2,  2 /4, and  indicate equivalent forces induced by bubble oscillation, migration, and the ambient flow field, respectively.The oscillation in a gravity field is described by equation (21) with the ambient flow ignored (u a = 0).In addition, neglecting both the ambient flow and bubble migration ( i.e. v = 0), equation ( 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  = 0,  a =  ∞ , and  is calculated with equation ( 19) retaining only the first term.On this basis, if the second term on the left-hand side is substituted with (1− /) /, equation (24) becomes the Gilmore equation.
For bubble dynamics in an incompressible flow, i.e.,  → ∞, equation ( 15) can be reduced to and equation ( 21) can be simplified as which further reduces to the Rayleigh-Plesset equation when the first term is neglected and  a =  ∞ .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, equation (21) involves two unknowns,  and .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  is the total mass of the gas inside the bubble, g is the gravity acceleration,  d is the drag coefficient,  p is the projected area of the bubble in the migration direction, and 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.X represents extra forces such as lift force and is not considered in this paper.Substituting  = ( b −  a ) / into equation ( 18) will give us  b .Then, the inviscid resultant force on the bubble can be expressed as where  a is the added mass coefficient and  = 4/3 3 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, equation ( 27) can be transformed into For a bubble oscillating and migrating in a still and unbounded flow, we have v = v m and ∇ a = g = (0, 0, −).Then equation (29) reduces to The bubble oscillation and migration can be obtained by solving equations ( 21) and ( 29) simultaneously, or by solving (21) and (30) if the ambient flow is not considered.

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  calculated by equation (19) retaining only the first term.Combining equations ( 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| −1 and the effect of migration are ignored.Note that  on the right-hand side of the above two equations including those in the subscripts should also be evaluated at  − (|r| − )/.Here, |r| −  denotes the minimum distance from the bubble surface to r. Substituting the above into equation (31) and ignoring the ambient velocity u a , one may obtain the formula for the pressure at r as Note that, the pressure disturbance in the fluid field with a distance  from the bubble center at time  is induced by the bubble at an earlier time,  − (|r| − )/, i.e., the bubble-induced disturbance is delayed in time due to fluid compressibility.In an incompressible flow where  approaches infinity, we have and Similarly, substituting them into equation (31), we have the degraded counterpart of equation (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 Section II.For an arbitrary point r in the fluid field, the velocity potential and the velocity induced by the migration of bubble  are of the orders of |o  − r| −2 and |o  − r| −4 , respectively, where o  is the center of bubble .Thus, only the effects of the oscillation of bubble  need to be considered.Hence, by taking the derivatives of equation ( 6) with respect to r and  with the solution of  substituted, the current velocity and the time derivative of the potential induced by bubble  at r can be expressed as and respectively, where   =  − (|r  | −   )/ is the initiation time of a disturbance induced by bubble  that later arrives at r at time .In the equations, the subscript  or  indicates the quantities of the corresponding bubble.
When considering the dynamics of bubble , 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 u B and the pressure  B of bubble  induced by other bubbles, considering the influences from all the other bubbles, as and where  is the total number of the interacting bubbles and o  is the center of bubble .Substituting the above two equations into equations ( 21) and ( 29) for bubble , and setting  as 1, 2, . . ., and , respectively, one may obtain a set of 2 equations that describe the oscillation and migration of the  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  approaches infinity, i.e., the external fluid is incompressible, we can obtain the reduced bubble dynamics equations as and where r   = o  − o  .The above equations with  = 1, 2, . . ., and  describe the dynamics of the multiple interacting bubbles in an incompressible flow.Particularly, for the dynamics of two interacting bubbles, i.e.,  = 2 and (, ) = (1, 2) or (2, 1), the above equations reduce to and

IV. MODELING BUBBLE DYNAMICS NEAR BOUNDARIES
Besides 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 e b • r +  = 0 where e b is the outward normal vector of the boundary and  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,  = −1.For bubble   , which is the image of bubble , it is required that respectively.Therefore, similar to equations ( 40) and ( 41), the flow velocity and the pressure at o  induced by other bubbles and all the images can be written as and Substitute the above two equations in to equations ( 21) and ( 29) for bubble , 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  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  b is the stand-off distance from the bubble center to the boundary.FIG.1: A comparison between the theoretical results and the experimental data [120] of the radius of an acoustic bubble in a free field.

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 Firstly, 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 m before being stimulated at  = 0 by a high-frequency acoustic wave, − amp sin(2  a ), with the pressure amplitude  amp = 125 kPa and the frequency  a = 29 kHz.In the theoretical calculations, we set the polytropic exponent as  = 1.4,and the saturated vapor pressure as  v = 2338 Pa (same for the subsequent all sections).We measure the effect of gravity with the buoyancy parameter  = √︁  max / ∞ in which  max 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 m before collapsing and continues to oscillate with a high frequency.The maximum radius after the first rebounding is merely 12 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 further details can be found in Zhang et al. [68].The bubble reached a maximum radius  max = 29.8mm 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  0 = 45 kPa,  0 = 2.07 mm, and  0 = 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  max and the minimum pressure  min when the bubble radius reaches  max . min can be determined by the ratio  max2 / max where  max2 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 Fig. 2(c) and (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 m × 4 m × 4 m cubic 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  max = 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  0 = 42.52 mm,  0 = 61.80m/s, and  0 = 304.13kPa.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 Fig. 3(c-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 3(c).Also, the migration calculated with our theory is more comparable to the experiment than that with the Geers-Hunter model, see 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 bubbles 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 mm 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 mm and 14.6 mm and the first oscillation periods were 3.05 ms and 2.83 ms, respectively.The initial radii used for theoretical calculation are  01 = 2.50 mm and  02 = 2.29 mm for bubble 1 and 2, respectively.The initial speeds of the bubble wall are  01 =  02 = 130 m/s and the initial pressure inside the bubbles are  01 =  02 = 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 Fig. 4(b-d ).The high-speed images of the bubble pair are shown in (b), where the two bubbles migrate towards 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 Fig. 4 (c) and (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 mm and 9.0 mm, respectively and the first oscillation periods were 2.92 ms 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  01 = 1.60 mm,  02 = 1.65 mm,  01 = 240 m/s,  02 = 110 m/s,  01 = 1.2 MPa, and  02 = 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 towards each other.The center of bubble 2 reciprocated during multiple oscillation cycles, first migrating away from bubble 1 and then towards.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 Fig. 5 (b) and (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  01 = 2.22 mm and  02 = 2.48 mm, and the initial oscillation speeds  01 = 180 m/s and  02 = 150 m/s, respectively.The initial internal pressure of both bubbles is 1.2 MPa.While oscillating in volume, the bubbles keep migrating towards 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 Fig. 6(b) and (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 4 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 towards 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:  01 = 1.41 mm,  02 = 1.36 mm, and  03 = 1.18 mm, respectively.The initial oscillation speeds are  01 =  02 = 180 m/s and  03 = 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 Fig. 7(b) and (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  b = 2.02 max , with the maximum bubble radius  max = 0.768 mm.The experimental setup is shown in Fig. 8(a) and high-speed images of the laser bubble oscillating while migrating towards 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  0 = 1.2 MPa,  0 = 0.121 mm, and  0 = 130 m/s in the same way as stated before.The results are plotted against the experimental data in Fig. 8(c) and (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 towards 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  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 [85,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  max = 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  = 0.126.The initial conditions of the bubble are  0 = 43.57mm,  0 = 61.83m/s, and  0 = 303.94kPa 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 Fig. 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 equation (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.0276s and 0.0283s) calculated by equation ( 56) and from the experiment, the bubble period calculated by the present theory, 0.0277s, 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 mm 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 downwards, affected by the free surface.The initial conditions for bubbles 1 and 2 in the theoretical calculation were  01 = 2.92 mm,  02 = 2.1 mm,  01 = 100 m/s,  02 = 180 m/s, and  01 =  02 = 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 Fig. 10(b) and (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  0 = 10 MPa,  0 = 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 Fig. 11(b) and (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.In addition, the above experiment was conducted in a reservoir where it was difficult to capture the bubble shape.Therefore, we present photos of an air-gun bubble produced in a water tank in Fig. 11(d) where it was better observed, though certain boundary effect may

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 are 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  0 = 0.173 max ,  0 = 100 ∞ ,  = 1.4,and  0 = 0. Fig. 12(a) depicts the theoretical and the numerical results of bubble radius variations with the standoff distance to the rigid boundary being  = 2 max .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  max ,  ∞ and , respectively.The scales for the other quantities are calculated accordingly, e.g., √︁  ∞ / for velocity and  max √︁ / ∞ 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.3 max , 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.4 max .

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  = √︁  ∞ // is taken as 0.013.The initial conditions of the two bubbles are given by the same  0 = 100 ∞ and  0 = 0.The initial radius of bubble 1 is  0 = 0.175 max 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   =  0,  0, /( − 1) where  = 1 or 2 indicates physical quantities of bubble 1 or 2, respectively.We then denote  * =  1 / 2 as the bubble energy ratio.The initiation of bubble 1 is marked as  = 0.The phase difference between the bubbles is defined as  = Δ * / * f where Δ * is the non-dimensional initiation time of bubble 2 and  * f 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  * 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 (5 times the  max 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  * −  parameter space.The migration is represented by the average relative speed of the two bubbles as in Fig. 14(a).The positive values represent the two bubbles eventually repelling each other and negative values represent attracting, with three demarcation lines  1 −  3 .The isoline  1 indicates that, at  * = 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  * and reaches approximately 0.55 at  * = 10, which means that a larger phase difference is required for the bubbles to repel each other when the bubble energy difference increases. 2 resides within 0.9 <  < 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  3 where  * 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  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  = 0.5 and  * = 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  * = 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 towards each other is prominent.The variation of the non-dimensional first period of bubble 1,  * , 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  1 and  2 ,  * 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. * reaches a minimum of 1.877 near  * = 1,  = 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  * leads to higher  * due to the weaker influence of bubble 2 with lower energy.Compared to  * , altering  results in more obvious changes in  * .As  approaches 1.0,  * grows until reaching  * f = 1.986; as  approaches 0, isoline  1 is crossed and  * exceeds  * f 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  2 ,  * varies in a very limited range and approximates  * f , since bubble 2 hardly affects the period of bubble 1 when  * is large.The scaled internal bubble pressure is reflected in Fig. 14(d ) by a nephogram of  * R =  * max  * , where  * max is the non-dimensional maximum pressure inside bubble 1 and  * the non-dimensional radius. * R 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  * R being equal to what bubble 1 induces in a free field,  * Rf , 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  * R is expected.By contrast, if bubble 1 is collapsing, the pressure wave will accelerate the collapse and  * R will increase.The energy transmission rate increases with  2 , 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  * R with  reduces with the increase in  * , since a higher energy ratio means a smaller size of bubble 2 and lower influence on bubble 1.For the region above  1 ,  * R remains greater than  * Rf .This indicates energy transmission from bubble 2 to bubble 1.The greatest  * R , i.e. the highest transmission, occurs in a region around  = 0.85,  * = 1, with  * R = 1.56 * Rf .This is because that the highest − 2  is found right in the above region.Therefore, the expansion of bubble 2 intensifies the collapse of bubble 1, leading to the high  * R .The Variation of bubble period, represented by the first oscillation period of bubble 1 subjected to interaction with bubble 2. On L1 and L2, the value equals the first period of bubble 1 in a free field without other bubbles.(d ) Variation of the bubble-induced pressure peak at a unit distance, indicated by the highest pressure peak reached inside bubble 1 multiplied by its radius.On L1, L2, and L3, the value equals that of bubble 1 in a free field without a second bubble.
required minimum  for  * R >  * Rf starts as about 0.5 at  * = 1 and increases almost linearly until reaching about 0.75 at  * = 10, as reflected by the demarcation line  1 .For the region below  1 , the pressure wave of bubble 2 arrives at bubble 1 at its expansion phase and does negative work to it.Thus,  * R falls below  * Rf , indicating energy transmission from bubble 1 to bubble 2. The minimum  * R is about 4.5 and appears around  = 0.1 and  * = 1 when  2  arrives at the maximum.Additionally, it's worth noticing that as  approaches 0,  * R may also exceed  * Rf for certain ranges of  * , corresponding to the region between  2 and  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  * R exceeds  * Rf .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., (, ) = (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 Fig. 15(a) and (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 is shown in Fig. 15(c) and (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  * = 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.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.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 1-dimensional spherically symmetric Euler equations: where  =  + 1 2  2 is the total energy and  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  s and , respectively.Thus, at  >  s , the fluid is undisturbed and at  <  is the gaseous product inside the bubble.These two parts provide the boundary conditions of the fluid within  <  <  s .In this work, we focus on the fluid between  and  s 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  is mapped linearly into the coordinate for the numerical solution  as which is also illustrated in Fig. 16, where  =  s −  is the length of the computational domain.Then we have where Note that an extra source term emerges in the scaled coordinate system.Then, the hydrodynamic pressure at  can be given by Here,   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 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 ( = 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  = , the fluid pressure of the water between  and  s equals the inner pressure of the bubble, i.e.,  g .The Mach number for the shock wave  =   / ∞ is calculated with the Huginiot condition Then, the initial conditions for the fluid density and velocity are given by and respectively.Here, the subscript ∞ indicates the physical quantities of the undisturbed water.The initial bubble radius  and the shock front radius  s at  =  are given by where  c is the radius of the explosive charge.In this paper, the value of  is determined by letting  =  c /100 at  = .
Within this framework, we can obtain the boundary conditions.For the boundary representing the material contact (the left boundary), we have û = 0 and  =  g .Thus, the flux is given by For the boundary representing the shock front (the right boundary), the flux is given by where  0 = ( ∞ +   )/( − 1) 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 <  <  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.

Initial condition for the unified bubble theory
In this work, we start the simulation with the unified bubble theory at the moment denoted by  e when  s reaches 10.The initial condition at  =  e is derived from the simulation of the initial phase prior to  e using the equations above.Let us denote  =  e 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  0 with the kinetic energy of the fluid field excluding the shock wave region: where  is a constant smaller than 1 to avoid the shock wave region in the integration.In this paper,  is taken as 0.65.

Shock wave far-field propagation
Besides 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.(B28) 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  > 10 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 G =  ( e ) for the left end, and G =  ( s ) for the right end. (B29)

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   elements and the ℓth element is bounded by  ℓ and  ℓ+1 .
Firstly, 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   .For element ℓ, the control equation (B7) of the near-field shock wave is then given in the Galerkin form where F is the flux at the boundary to be calculated according to specific boundary conditions.The matrix form is given by An explicit time marching scheme can be used to update K, e.g., the Runge-Kutta method.In this paper, 1 element with the  15 space was used first until  > 2.Then, 20 elements with the  3 space were used until  > 10.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, equation (B27) is solved with the same approach.Only 1 element and  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 3-Dimensional model should be used.

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 m and 30.0 m, respectively.The comparison of the experimental and the simulated results using the near-field and far-field methods are illustrated in Fig. 18 where a good agreement is achieved.

FIG. 2 :
FIG. 2: Underwater spark-generated bubble experiment and comparisons of theoretical and experimental results in a gravity field.(a) Schematic diagram of the experiment setup.(b) Selected sequential high-speed images of the spark-generated bubble.The capturing time is labeled below each image in milliseconds.Frame width, 51.62 mm.(c) Comparisons between the theoretical and the experimental results of bubble radius.(d ) Comparisons between the theoretical and the experimental results of the bubble migration.

FIG. 3 :
FIG. 3: Underwater explosion bubble experiment and comparisons of theoretical and experimental results in a gravity field.(a) Schematic diagram of the experiment setup.(b) Selected sequential high-speed images of the underwater explosion bubble.The capturing time is labeled below each image in milliseconds.Frame width, 0.490 m.(c) Comparisons between the theoretical and the experimental results of the bubble radius in the gravity field.(d ) Comparisons between the theoretical and the experimental results of the bubble migration in the gravity field.(e) Comparisons between the theoretical and the experimental results of the pressure pulse of the first collapse.The experimental pressure data is measured at the test point as in the schematic diagram.

FIG. 4 :
FIG. 4: The first experiment case of two interacting spark-induced bubbles and comparisons between theoretical and experimental results.(a) Schematic diagram of the experiment setup.(b) Selected sequential high-speed images of the spark-induced bubble pair.The capturing time is labeled below each image in milliseconds.Frame width, 149 mm.(c) Comparisons between theoretical and experimental results of bubble radius variation.(d ) Comparisons between theoretical and experimental results of bubble migration.

FIG. 5 :
FIG. 5: The second experiment case of two interacting spark-induced bubbles and comparisons between theoretical and experimental results.(a) Selected sequential high-speed images of the spark-induced bubble pair.The capturing time is labeled below each image in milliseconds.Frame width, 96.86 mm.(b) Comparisons between theoretical and experimental results of bubble radius variation.(c) Comparisons between theoretical and experimental results of bubble migration.

2 FIG. 6 :
FIG. 6: Comparison of the interaction of three spark-generated bubbles with the theoretical prediction.(a) Selected sequential high-speed images of the three spark-induced bubbles.The capturing time is labeled below each image in milliseconds.Frame width, 147 mm.(b) Comparison between theoretical and experimental results of bubble radius variation.(c) Comparison between theoretical and experimental results of the bubble migration in the vertical direction.

FIG. 7 :
FIG. 7: Comparison of the bubble dynamics in a cluster of sixteen spark-generated bubbles with the theoretical prediction.(a) Selected sequential high-speed images of the bubble cluster.The black dots in the seventh frame are cavitation bubbles triggered by the pressure fluctuation due to the collapses and rebounds of the spark bubbles.The capturing time is labeled below each image in milliseconds.Frame width, 151 mm.(b) Comparison between theoretical and experimental results of the bubble radius variation.(c) Comparison between theoretical and experimental results of the bubble migration in the vertical direction.

FIG. 8 :
FIG. 8: Laser-induced cavitation bubble experiment near a rigid wall and comparisons between theoretical and experimental results.(a) Schematic diagram of the experiment setup.(b) Selected sequential high-speed images of the laser-induced cavitation bubble.The capturing time is labeled below each image in microseconds.Frame width, 2.16 mm.(c) Comparisons between theoretical and experimental results of bubble radius variation.(d ) Comparisons between theoretical and experimental results of bubble migration.

FIG. 9 :
FIG. 9: Underwater explosion bubble experiment and comparisons of theoretical and experimental results near a free surface.(a) Schematic diagram of the experiment setup.(b) Selected sequential high-speed images of the underwater explosion bubble.The capturing time is labeled below each image in milliseconds.Frame width, 0.446 m.(c) Comparisons between the theoretical and the experimental results of the bubble radius considering gravity and boundary effect.(d ) Comparisons between the theoretical and the experimental results of the bubble migration considering gravity and boundary effect.(e) Comparisons between the theoretical and the experimental results of the pressure pulse of the first collapse.The experimental pressure data is measured at the test point as in the schematic diagram.

FIG. 10 :
FIG. 10: Comparison of the interaction of three spark-generated bubbles near the free surface with the theoretical results.(a) Selected sequential high-speed images of the three bubbles.The capturing time is labeled below each image in milliseconds.Frame width, 153 mm.(b) Comparison between theoretical and experimental results of bubble radius variation.(c) Comparison between theoretical and experimental results of the bubble migration in the vertical direction.

FIG. 11 :
FIG. 11: Air-gun bubble experiment and comparisons of theoretical and experimental results near a free surface.(a) Schematic diagram of the experiment setup.(b) Comparisons between the theoretical results of bubble radius variation versus time.(c) Comparisons between the theoretical and the experimental results of bubble-induced pressure variation versus time.(d ) Selected sequential high-speed images of the air-gun bubble generated in a water tank with the initial internal pressure and volume of the chamber being  0 = 5 MPa and  0 = 2458 cm 3 , respectively.The capturing time is labeled below each frame in milliseconds.Frame width, 113 cm.

FIG. 12 :
FIG. 12: Comparisons between the theoretical result and the numerical simulation.(a) Bubble radius variations near a rigid boundary and in a free field obtained using the present theory and different numerical methods.(b) The discrepancy between the theoretical and numerical results of the first bubble oscillation period as a function of the standoff distance to the boundary scaled by  max .

FIG. 13 :
FIG.13: Sketch for the configuration of the interacting bubble pair and the measuring point for fluid field pressure.

FIG. 14 :
FIG. 14: Dynamics of two interacting bubbles with varied initial energy ratio  * and phase difference .(a) Variation of bubble migration, represented by the relative speed between the two bubbles averaged over the first three and a half oscillation periods of bubble 1.The positive values indicate bubbles repelling each other and the negative attracting.L1, L2, and L3 represent where the relative speed is zero.(b) Typical time curves of the relative migration between the two bubble centers.(c)Variation of bubble period, represented by the first oscillation period of bubble 1 subjected to interaction with bubble 2. On L1 and L2, the value equals the first period of bubble 1 in a free field without other bubbles.(d ) Variation of the bubble-induced pressure peak at a unit distance, indicated by the highest pressure peak reached inside bubble 1 multiplied by its radius.On L1, L2, and L3, the value equals that of bubble 1 in a free field without a second bubble.

FIG. 15 :
FIG.15: Theoretical and numerical results of the dynamics of two interacting bubbles with phase and energy differences obtained using the present unified theory and the Eulerian Finite Element Method (EFEM), respectively.(a) Temporal variation of bubble radii for the bubble pair with  = 0 and  = 2 (the first scenario).(b) Temporal variation of the pressure induced by the bubble pair in the flow field with  = 0 and  = 2. (c) Temporal variation of bubble radii for the bubble pair with  = 0.5 and  = 1.6 (the second scenario).(d ) Temporal variation of the pressure induced by the bubble pair in the flow field with  = 0.5 and  = 1.6.

) 1 𝐿(
Denote û =  as the fluid velocity observed in the coordinate system of  .Thus,û =  −   − ).(B6)Denote ∂/ ∂ = / + (  + )/ as the derivative with respect to  when  is fixed.Then, the Euler equations can be rearranged and expressed in the following compact form:

FIG. 16 :
FIG.16: Sketch for the linear mapping between the physical domain () to the numerical solution domain () for the near-field solver.
B10) where  c is the radius of the explosive charge and , ,  1 ,  2 , and  are material related constants.If the surrounding fluid is liquid, we model it by the Tammann EoS  = ( − 1) −   .(B11) For far-field flow, we modeled the fluid with the Tait equation  = ( ∞ +   )   ∞  −   .(B20)Thus, the sound speed  satisfies 2 =   ∞ ( ∞ +   ) 1/ (  +   ) 1−1/ .(B21)Expanding  at  =  ∞ and denoting   =  −  ∞ , we have =  ∞ +     +  ( ) 2 , (B22)where  = ( − 1)/(2 ∞  ∞ ).The energy of the shock wave denoted by  can be evaluated by =   +   ≈ ( ∞  2 ∞ ) −1 ( ) 2in the far field.The shock wave energy is transported at the sound speed .Thus, the conservation of the shock wave energy can be expressed as that the change of wavelength of the shock wave is slow in the far field, we map  to  with the following relation =  −  e  s −  e , (B24)where  e is the left end of the computational domain and  s −  e is the length of the computational domain which is a constant and denoted by   .Approximating  by () −1   , we have replace  with ( ∞  2 ∞ ) −1  2 , and  with B22 neglecting  ( ) 2 , then equation (B25) can be expressed in a conservation-law form as ∂ ( ) 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 ∂ ( ) the flux  = ( ∞ −   )  + 3 4   ( ) 2 and the source  = − /.Then, the pressure history at  can be expressed as (, ) =  ∞ when  >  s ,  ∞ +   for else.

FIG. 17 :
FIG. 17: Experimental and theoretical results at the initial phase of an underwater explosion with the charge equivalent to 2.0 g TNT.(a) Initial bubble expansion and shock wave front propagation.(b) An image of the underwater explosion captured at 0.106 ms from detonation.The theoretical results are projected on the image with the green and the red lines representing the bubble wall and the shock wave front, respectively.

FIG. 18 :
FIG.18: Experimental and theoretical results of shock wave pressures for underwater explosions with the charges being equivalent to (a) 7.0 kg TNT measured at a distance of 4.0 m from the charge center and (b) 8.0 kg TNT measured at a distance of 30.0 m from the charge center.