Theoretical study on bubble dynamics under hybrid-boundary and multi-bubble conditions using the unified equation

This paper aims to use the unified bubble dynamics equation to investigate bubble behavior in complex scenarios involving hybrid free surface/wall boundaries and interactions between multiple bubbles. The effect of singularity movement on the unified equation’s form is analyzed after deriving the bubble pulsation equation using a moving point source and a dipole, followed by discussions on the effect of migration compressibility-related terms on the bubble dynamics. In addition, the present study accounts for the impact of hybrid boundaries, including crossed and parallel boundaries, by introducing a finite number of mirror bubbles for the former and an infinite number of mirror bubbles for the latter. Spark bubble experiments and numerical simulation are conducted to validate the present theory. The application of the unified equation in multi-bubble interactions is exemplified by computing a spherical bubble array containing more than 100 uniformly distributed cavitation bubbles under different boundary conditions. The bubble cluster-induced pressure peak can reach nearly two times or even higher than that of an individual bubble, highlighting the damage potential caused by cavitation bubble clusters.

As early as 1917 and 1949, Rayleigh [31] and Plesset [32] established the bubble pulsation equation in a incompressible flow field.In 1970, Hicks [33] modified the Rayleigh-Plesset equation by considering incompressible migration.Keller et al. [34] derived the bubble pulsation equation in compressible flow fields based on the wave equation in 1980.Chahine [35] included the bubble migration in Keller equation by adding an incompressible migration term in 2009.These theories have played an important role in the bubble dynamics field.However, because bubbles are typically not isolated, their dynamic behavior is obviously affected by boundary effects, multiple bubbles, flow field environment, gravity, and many other conditions.These factors can cause bubbles to have pulsation and compressible migration.However, the above-mentioned theories did not fully consider bubble pulsation and migration in a compressible flow field.
To more accurately predict bubble dynamics in theory, Zhang et al. [36] constructed the basic solution of the moving singularity of the wave equation from the most basic mathematical principles while considering the compressible migration effect of bubbles in 2023.Through a series of mathematical and physical calculations and experimental verification, they derived a unified equation for bubble dynamics with an elegant mathematical form that can simultaneously consider boundary effects, multiple bubbles, flow field environment, gravity, migration, compressibility, viscous force, and surface tension.The unified equation enables a theoretically accurate prediction of bubble dynamic behavior, guiding experiments and predictions of new physical phenomena.
Extensive studies on the dynamics of single bubbles near a single boundary, including the rigid wall [37][38][39], free surface [40,41,41], and elastic boundary [42][43][44], have uncovered and clarified many interesting physical phenomena: re-entrant jet, annular jet, and water spike on the free surface.In recent years, bubble dynamics under the action of multiple bubbles or multiple boundaries have also received much attention.Two bubbles generated simultaneously produce two liquid jets in opposite directions [45,46], while the collapse patterns of bubbles generated unsynchronously are highly correlated with the phase difference [47,48].Some special phenomena, such as the mutual penetration between bubbles and fine thin jets, are found by Liang et al. [47] and Tomita et al. [48].When the bubble buoyancy effect is considered, the bubble jet direction is deflected upward [49,50], altering the bubble migration direction.When the third bubble is included, it enhances or weakens the strength of the bubble jets generated by the other two bubbles depending on the phase and size difference of the three bubbles [51].In addition, for bubble clusters formed by more bubbles, the bubbles inside the cluster are accelerated in collapsing by the surrounding bubbles and generate a much stronger pressure pulse than that of a single bubble [52][53][54].Different ambient flow fields could significantly alter the dynamic characteristics of bubble clusters [55,56].For bubble dynamics near multiple boundaries, when two rigid walls affect the bubble dynamics in different directions, the bubble jet is directed toward the closer wall [57][58][59].Furthermore, when the bubble buoyancy effect is considered, the jet may turn from downward to upward upon improving the buoyancy effect [60].When both a free surface and a sidewall are present at similar distances to the bubble, the bubble is more affected by the free surface, producing a downward but slightly inclined jet [61][62][63][64].[65] investi-gated three boundaries, including the free surface, sidewall, and bottom wall, and linked the various bubble collapse patterns (e.g., jetting and quasi-spherical collapsing) to distance parameters and the buoyancy effect.
The above studies used numerical simulations and experiments to investigate bubble dynamics in different environments in most cases.Numerical methods such as boundary element [66][67][68][69], finite element [70][71][72], finite volume [73][74][75], and smooth particle hydrodynamics methods [76,77] have been greatly developed in the past decades for predicting bubble behavior and revealing its underlying mechanism.Experimental methods such as laser pulse [78][79][80][81], electric discharge [42,[82][83][84] and small charge underwater explosions [19,85,86] are mature.However, when handling farfield problems or a large number of bubbles, both simulations and experiments can be costly.In such cases, theoretical models prove superior for investigating the volume oscillation, migration, and field pressure of a large number of bubbles.
In this study, we derive a unified equation using a moving point source and dipole, and the effects of singularity motion and compressible bubble migration are discussed.We then further apply the theory to bubble dynamics near a hybrid boundary, including crossed and parallel boundaries, and compare the results with experiments and numerical simulations.The dynamics of three-dimensional bubble arrays are also discussed by computing a spherical bubble cluster containing over 100 cavitation bubbles under different boundary conditions, using the derived equation.Computed results show that the bubbles in the cluster could have greater damage potential to the surrounding structure than a single bubble.
2 Derivation of the unified equation for bubble dynamics

Point source and dipole for moving bubble
In the compressible flow, the conservation equations for mass and momentum are and where ρ, u, p, and τ represent the pressure, density, velocity, and deviatoric stress tensor of the fluid, respectively; F is the body force.
Following the previous works, we assumed that the fluid around the bubble is weakly compressible [36], and the wave equation is obtained: in which C is sound speed in water and φ is the velocity potential.
In this study, the origin of a spherical coordinate system o − rθϕ is set to the initial bubble center.Assume that the bubble remains spherical during the whole process.θ = 0 denotes the bubble migration direction.In the previous study [36], we obtained a unified bubble equation by superimposing the velocity potential induced by a point source and a diple.In this study, we further derive the unified bubble equation by considering the movement of the two singularities.First, a source is set at the origin with the strength of f 0 (t), and the basic solution of the wave equation can be written as follows: As the point source and dipole move with an arbitrary velocity of v(t), the basic solution can be obtained by the superposition of infinite ϕ f 0 as follows: where r t denotes the vector pointing from the source at t − |r t |/C to r.However, providing the analytic solution of r t is challenging because a nonlinear equation is involved.With reasonable simplification, the migration velocity v can be treated as a constant as the singularity propagates and exerts influence from the bubble center to the bubble surface.During this period, the displacement of the singularities Rv/C is a small quantity relative to the bubble radius R. Figure 1 shows a sketch of the singularities and the propagation of their influence.R denotes the bubble radius at the present time t.The position of the singularities at the time t − R/C is the same as that of the origin at the time t (blue dot).Thus, the singularity-emitted influences at t − R/C reach the bubble surface (r = R) at t. Therefore, the distance from the present location of the singularities (red dot) to the origin O(t) is Rv/C.In addition, the velocity potential at r is determined by the singularity-emitted influences at time t − |r t |/C.The singularity position at that time is indicated by the orange dot.The distance between the location of the singularities at the present time t and that at t − |r t |/C can be easily obtained as |r t |v/C.Thus, we have the equation with a single unknown, |r t |, as follows: where ε = v/C; t 0 is a reference time to offset the origin to the bubble center such that t − t 0 = 0.However, when calculating the derivative of d and D, t − t 0 should be retained.Then, |r t | can be obtained by solving the above equation: Substituting equation (7) into equation ( 5), we can obtain the expression of the moving source: where D is expressed as follows: For a moving dipole with a strength of q(t), the wave equation can be solved by taking the limit as follows: where e denotes a unit vector along the migration direction of the bubble.Thus, the solution of the wave equation for a source with the strength f and a diple with the strength q both moving with a constant velocity v can be expressed as follows: where Q is written as Verifying that equation (11) satisfies the wave equation is a straightforward task.Next, the velocity components in r and θ directions at the bubble wall can be expressed as follows: and respectively, where f ′ and q ′ denote the derivative of f and q with respect to their arguments, respectively.

Single bubble dynamics equation
The kinetic boundary condition [36] for oscillation and migration on the bubble surface can be expressed as follows: respectively, where V denotes the bubble volume and ∂V denotes the bubble surface.Ṙ is the time derivative of the bubble radius.Denoting and substituting equation ( 13) into above two equations, we have respectively.Using a perturbation method on equation (18) with perturbation parameter of 1/C, we obtain the expression of with only the terms up to the order C −1 reserved.Differentiating both sides of equation ( 17) with respect to time, we have The above equation is in exactly the same form as that in the previous work [36].dF/dt is an operator that can be used to consider the effects of multiple bubbles, boundaries, and the flow field environment.Let Ma = Ṙ/C, then equation ( 20) can be expressed in the form of Mach number To obtain the specific expression of dF/dt in equation ( 20), the Bernoulli equation is employed: where H is the enthalpy difference; H = P b P a ρ −1 dp, P b is the fluid pressure on the outside of the bubble surface; P a is the ambient pressure at the bubble position (P a = P B + P E ; P B is induced by boundaries or other bubbles; P E denotes extra pressures including the hydrostatic pressure); ρ is the liquid density.In this study, we reserve only the first-order term of H with respect to (P b − P a )/(ρC 2 ) and then H = (P b − P a )/ρ.
On the bubble surface, the normal stress balance can be expressed as follows: where P g , σ, and τ n are the inner bubble pressure, surface tension, and normal viscous stress.P g is calculated with the equation of state of the internal gas-the adiabatic ideal gas equation was adopted in this study for simplicity, as suggested by Lee [87] and Wang [88].According to the Stokes assumption on the viscous stress of Newtonian fluid, the normal viscous stress can be expressed as follows: where µ is the fluid viscosity.Retaining only the effect of the point source, we have the following relation evaluated at the bubble surface As for the velocity divergence term in equation ( 25), employing the linearized equation of the state of the external fluid, we have Thus, plug the above two equations into equation ( 25), we have When the high-order terms are dropped, equation ( 24) agrees with that used in Zhang et al. [36].Similar to the work of Zhang et al. [36], dF/dt can be derived as follows: Substituting equation ( 29) into equation ( 20), the bubble oscillation equation is obtained: Expanding the above equation, we could get the following differential equation: which includes the effect of the compressible bubble migration on bubble pulsation that is ignored in the Keller-Miksis equation.Equation ( 31) is a typical second-order ordinary differential equation, which is often solved using the fourthorder Runge-Kutta method to obtain the temporal variation of the bubble radius.Another variable v in equation ( 30) is obtained using the bubble migration equation [36]: where C a and C d are the added mass coefficient and drag coefficient, respectively; the values of C a and C d are taken as 1.0 and 0.5, respectively; To include the propagation of pressure waves [36], the pressure of a point r in the flow field induced by a bubble is calculated as follows: 3 Modeling bubble dynamics near hybrid boundaries

Bubble dynamics at a corner crossed by two boundaries
A type of common hybrid boundary is the corner.Near a corner formed by two crossed boundaries, bubble dynamics are influenced by the boundaries in two distinct directions.Predicting bubble behavior near corners at arbitrary angles is challenging due to the lack of physical mapping between corners at different angles.However, when dealing with two semi-infinite boundaries at an angle of α = π/k (where k is a positive integer), the boundary effect can be accurately modeled by incorporating a finite number of mirror bubbles with a specific spatial distribution based on image theory.
To account for the effects of two boundaries using image theory, the positions of mirror bubbles on both sides of any boundary need to be completely symmetrical, which requires the position of the mirror bubbles to be arranged reasonably.The spatial distributions of mirror bubbles for two typical corners of two walls (α = π/4 and π/2) are presented in Fig. 2. The coordinate system is the same as that in Fig. 1, in which the origin is set at the initial bubble center and θ = 0 denotes the direction of bubble migration.The black circle represents the real bubble while others denote mirror bubbles.Define o N as the position vector of the real bubble N. Let us describe the plane of boundary i with n i • r + b i = 0, where n i denotes the outer normal vectors of boundary i and b i is a constant.The position of each mirror bubble is obtained by the symmetry of the other bubbles about the boundary.For example, mirror bubbles 1 and 2 are obtained from the real bubble symmetric about boundaries 1 and 2, respectively; mirror bubbles 3 and 4 are obtained from mirror bubbles 1 and 2 symmetric about boundaries 2 and 1, respectively.Here, we denote the position vector of the mirror bubble i of bubble N as o Ni : The black circle denotes the real bubble, and the others denote mirror bubbles.The solid lines represent the boundaries.The origin of the coordinate system is at the initial bubble center and the direction of bubble migration is consistent with θ = 0, which is the same as that in Fig. 1.
Let the reflection coefficients of boundaries 1 and 2 be ω 1 and ω 2 , respectively.The reflection coefficient is 1.0 for the rigid boundary and -1.0 for the ideal free surface.Thus, the coefficient of the mirror bubble ξ Ni can be expressed as follows: in which ξ N0 = 1.
As the bubble dynamics are affected by all the image bubbles, equations of multi-bubble interaction are employed [36] to consider the influence of crossed boundaries.The velocity u and derivative of velocity potential φ induced by the bubble N at r are as follows: and Considering the effect of mirror bubbles and other real bubbles, the velocity u B and pressure P B of the background flow field of real bubble M can be written as follows: (38) and respectively.L is the number of the real bubbles.By substituting the above two equations into equations ( 30) and (32) for the real bubble, we have the equations for the bubble dynamics considering the crossed boundary effect in the compressible flow.
Then, the pressure in the flow field is the sum of that induced by the real bubble and all the mirror bubbles: in which N all is the number of all bubbles including real and mirror bubbles.
As C tends to infinity, the above model reduces to the bubble dynamics equations for bubble M in the incompressible fluids: and Equation ( 42) reduces to the incompressible multi-bubble equation, as presented in the previous work [36] when C a = 0.5 on the left side.

Bubble dynamics between two parallel boundaries
As the angle of the crossed boundary is zero, the two boundaries become parallel and n 1 + n 2 = 0, as shown by boundaries 1 and 2 in Fig. 3. Using image theory to simulate the effect of two parallel boundaries on a bubble in-between presents a challenge, as an infinite number of mirror bubbles is theoretically required to represent the boundary condition (the maximum value of i in equations ( 34) and ( 35) tends to infinity).To overcome this challenge, we gradually increased the number of mirror bubbles and calculated the dynamic behavior of the real bubble.We continued this process until the calculation results reached a steady state, where further increases in the number of mirror bubbles no longer significantly affected the bubble dynamics results.The number of mirror bubbles is defined as m.At this point, we deemed the reproduction of the two parallel boundaries' effect to be satisfactory.We define the upper boundary as boundary 1 and the lower boundary as boundary 2.
To illustrate the convergence of the result as the number of mirror bubbles is sufficiently large, a typical case is considered with ω 1 and ω 2 set to 1.0.The results are nondimensionalized with the maximum radius of the bubble in the free field, fluid density (ρ = 1000 kg/m 3 ), and pressure at infinity as reference values.The initial conditions for bubble pulsation are P * 0 = 20, R * 0 = 0.06 and Ṙ * 0 = 5.98, in which the superscript '*' represents the dimensionless physical quantity.The initial dimensionless distance between the bubble and boundary 1 is 3.0, while the distance from boundary 2 is 2.0.The number of mirror bubbles m ranges from 1 to the case where the peak pressure of bubble pulsation converges, as shown in Fig. 4. As m increases, the asymmetry effect of the mirror bubbles on both sides of the boundaries decreases, and the difference between the calculated results reduces.As m is small, the error of the calculation results is relatively large since the excess mirror bubbles significantly impact the real bubble.If we focus only on the radius and period of the bubble, the calculation results converge within 0.1% after m reaches 14.After m reaches 30, the peak pressure achieves convergence within a tolerance of 0.1%.Thus, we take m as 30 in the subsequent calculations.

Bubble dynamics between two parallel boundaries and near a sidewall
Further, when the boundary contains both crossed and parallel features, more mirror bubbles are involved and the pressure and velocity of the background flow field of bubbles are a linear superposition of mirror bubbles and other real bubbles.For example, when boundary 3 is introduced as illustrated in 3, the positions of the bubbles at both sides of the three boundaries should be perfectly symmetrical.However, when the angle between two of the boundaries is obtuse, the image method cannot equate the boundary effect.Thus, only the situation where the sidewall is perpendicular to the free surface and bottom could be modeled, that is, a semi-truncated rectangular channel.Under this condition we need to symmetrically move all the bubbles under the influence of the two parallel boundaries to the other side of the third boundary, as shown in Fig. 3. Thus, the position vector and coefficient of the mirror bubble on the other side of boundary 3 are as follows: and in which n 3 is the outward normal vector and ω 3 is its reflection coefficient; b 3 is a constant; i ranges from 1 to m in equations 43 and 44.As the bubble 2m + 1 is symmetry with the real bubble about boundary 3, we have o The velocity and pressure of the background flow field of the real bubble M can be written as follows: (45) and  4 Results and discussion

Bubble dynamics near hybrid boundaries
In the work of Chahine et al. [35,89], the Keller equation was modified by adding v 2 /4 to account for the effect of bubble migration on bubble pulsation, similar to the bubble dynamics equation in the incompressible fluid [33,90].However, this approach does not consider the effect of compressible bubble migration.In this study, we investigated the effect of bubble migration on bubble oscillation in a compressible flow by removing the terms Ṙv 2 /4C and Rvv/2C in equation ( 31) and compared the results with the original equation.Figure 5 presents a comparison between theory and experiment, where a bubble was generated by the explosion of 1.125 g PBX at a distance of 0.4 m from the free surface.Further details on the experimental system and the initial calculating conditions can be found in the previous work [36].The red solid line shows the result of the original equation, while the blue line represents the result when only the term v 2 /4 is retained among the three migration-related terms.The radii and period of the red and blue lines differ significantly, especially in the second cycle of the bubble, indicating that the compressible terms Ṙv 2 /4C and Rvv/2C affect the energy loss of the bubble caused by the fluid compressibility.The deviation between the theoretical migrations and the experiment occurs mainly Figure 5 Comparison between the theoretical and experimental results of (a) bubble radius and (b) vertical bubble migration.The bubble is generated by an underwater explosion of 1.125 g PBX at a distance of 0.4 m from the free surface.The experimental come from the work of Zhang et al. [36].
in the second and later cycles and could be due to complex factors, including the non-spherical bubble behavior induced by gravity or boundaries, such as jetting, splitting, and phase transition.In addition, we also calculated the time histories of the radius and migration for a large-charge underwater explosion bubble, and compared them with the results from the Euler finite element method (EFEM) [71,91], as shown in Fig. 6.EFEM is a numerical method used to solve the compressible Euler equation system.It proves to be effective in simulating deformable bubbles.The bubble is generated by 500 Kg TNT explosives at a depth of 40 m.Compared to the result obtained by directly adding v 2 /4 and the result from Keller equation, the bubble radius and migration in the present theory agree better with the numerical simulation results.The two terms related to compressible migration obviously affect the dynamic characteristics of the bubbles.The result in the present theory agrees well with that of EFEM, especially in the second cycle.Therefore, the present theory, with all three migration-related terms, provides a more accurate approach for accounting for the effect of bubble migration on bubble pulsation in compressible fluids simultaneously.
To validate the present theory in predicting bubble dynamics near hybrid boundaries, we conducted two experiments involving spark-induced bubbles and compared the bubble radii and vertical migration between theory and experiment, as illustrated in Fig. 7 and Fig. 8.In the first case, we calculated the bubble dynamics near two perpendicular walls.In the experiment, to account for the effect of the vertical sidewall, we simultaneously generated two bubbles of the same size to equate a vertical rigid wall in the mid-pendant of the two bubbles, according to image theory.The initial bubble is located 33 mm from the sidewall surface and 34 mm from the bottom.The initial conditions for spark-induced bubbles in the theory are determined using the method proposed by Wang [92]: the maximum radius R max and the minimum pressure P min are required at the moment of maximum bubble volume.P min can be determined by the ratio R max2 /R max where R max2 is the maximum bubble radius in the second cycle.We performed a backward integration starting from the moment of maximum bubble volume with zero wall velocity until the initial condition is obtained.The initial conditions of the bubble in this case are as follows: P 0 = 1.2 MPa, R 0 = 3.35 mm and Ṙ0 = 95 m/s.We also provided time-history curves of the bubble radius and migration when the sidewall and bottom are present individually.The combined effect of the two walls significantly prolongs the bubble pulsation period, leading to a significant deviation of Keller equation in predicting the bubble period.Our theory is in excellent agreement with the experimental results.
In the second case, we compute the bubble dynamics near the free surface and vertical sidewall at a distance of 49.5 mm from the free surface and 35 mm from the sidewall.Similarly, the sidewall is equivalently represented by two horizontally arranged spark-induced bubbles of the same size in the experiment.The initial conditions of the bubble are P 0 = 1.2 MPa, R 0 = 2.98 mm and Ṙ0 = 90 m/s.Considering only the effect of the free surface underestimates the pulsation period of the bubble, while considering only the sidewall overestimates it.The vertical migration of the bubble is mainly caused by the Figure 6 Comparison between the theoretical and numerical results of (a) bubble radius and (b) vertical bubble migration.The bubble is generated by an underwater explosion of 500 Kg TNT explosive at a depth of 40 m.The effect of the free surface is ignored in both theoretical and numerical results as the distance between the bubble and free surface exceeds 5 times the maximum bubble radius.
free surface.Due to the neglect of boundary effects and bubble migration, Keller equation cannot accurately calculate the period and energy loss of the bubble.Our theoretical reproduce the evolution of the bubble radius and migration in the experiment well.
To validate the influence of parallel boundaries in our theory, we compared them with numerical simulations using EFEM.The results are shown in Figure 9, which includes the bubble radius, migration, and pulsation pressure.The bubble is initially located between two parallel boundaries, namely the free surface and the bottom wall, with a distance of 1.5 m from both boundaries.The initial pulsation conditions of the bubble are: P 0 = 60 MPa, R 0 = 32 mm, and Ṙ0 = 0 m/s.The maximum radius of the bubble is 0.365 m, resulting in a dimensionless distance of about 4.1 between the bubble and the two boundaries.The present theory exhibits good agreement with the EFEM results regarding the radius, period, and migration of the bubble at this distance parameter.While the EFEM-calculated pulsation pressure is slightly smaller than the theoretical prediction, this discrepancy can be attributed to numerical dissipation.Overall, our theory reproduces the time-dependent process of bubble pulsation pressure well compared with the traditional Keller equation.
Finally, we calculate the dynamics of three bubbles under the combined action of three boundaries, including the sidewall, bottom, and free surface, and compare the experimental results with the theoretical values, as shown in Fig. 10.In the experiment, to equate the effect of the sidewall, six sparkinduced bubbles are generated simultaneously by an identical voltage, and they are horizontally spaced at a distance of 40 mm.The three bubbles on the left are marked as bubble 1, bubble 2, and bubble 3, as illustrated in Fig. 10(a).These configurations correspond to a sidewall on the right located 20 mm to bubble 3.All bubbles are initially positioned 49 mm from the free surface and 45 mm from the bottom.Due to the relatively free flow of fluids on the left side of bubble 1, bubble 1 is the first to collapse to its minimum volume.The pulsation of bubble 3 has the maximum pulsation period due to the influence of the right sidewall and two bubbles on the left.Throughout the entire pulsation process, all bubbles migrate downward in the vertical direction because the buoyancy effect of the bubbles is too weak to counteract the influence of the free surface and bottom.The initial radius for the bubbles in theory are: R 01 = 2.36 mm, R 02 = 2.25 mm and R 03 = 2.41 mm.The initial inner bubble pressure and pulsation velocity of the three bubbles are 1.2 MPa and 100 m/s, respectively.The evolutions of the radius and migration of bubbles in the experiment agree well with the computed results in the present theory, as illustrated in Fig. 10(b) and (c).

Multiple bubble dynamics near boundary
We first conducted an experiment with nine spark-induced bubbles distributed in the three-dimensional space and compared the obtained experimental values with theoretical results, as shown in Figure 11.Fig. 11(a) illustrates the relative positions of nine bubbles on the plane captured using a high-speed camera, where bubbles 1 and 2 are marked as representative bubbles of the outermost and second circles, respectively.Bubble 3 is identified as the central bubble.The nine bubbles are distributed in three planes.The outermost cycle of bubbles is situated in the closest plane to the field of view, whereas bubbles 2 and 3 are generated 48 mm and 22 mm away from the plane of bubble 1, respectively.Both the outermost and second circles of bubble array comprise four uniformly distributed bubbles on a 2 × 2 grid, as illustrated in Fig. 11(a).The first circle's grid interval is set at 116 mm, while that of the second circle is 49 mm.As the positions of the nine bubbles are nearly symmetrical in the photographed plane, all the bubbles migrate toward the plane center.Note that some black dots are presented in the fourth frame, which are cavitation bubbles induced by pressure variation due to the rapid collapses and rebounds of bubbles.The outer bubble has a smaller pulsation period than the inner bubble due to the relatively free flow of fluid on the side away from the inner bubble.All bubbles have the same initial inner pressure and pulsation velocity in theory: P 0 = 1.2 MPa and Ṙ0 = 120 m/s.The initial radius of bubbles are: R 01 = R 03 = 1.94 mm and R 02 = 1.88 mm.Considering the slight differences in bubble dynamics during most of the first pulsation cycle, the other bubbles in the outermost and second cycles have the same initial conditions as bubble 1 and bubble 2, respectively.Figure 11(b) and (c) presents a comparison of the time-history of the bubble radius and vertical migration between theory and experiment.The present theory reproduces the dynamic features of bubbles observed in the experiment well, including the pulsation period of bubbles, accelerated contraction of bubble 3 at the end of collapse and bubble migration direction.
Next, as an application of the unified equation to the multiple bubbles, we discuss the dynamics of a spherical bubble cluster with 105 cavitation bubbles below a rigid wall.The bubble cluster comprises four layers of bubbles.The innermost layer has only one bubble located at the cluster center, while the outermost layer has 58 bubbles distributed evenly on a sphere with a radius of 195 mm.The radii of the middle two layers are 81 mm and 135 mm, respectively.The layer with a radius of 81 mm contains 14 bubbles uniformly distributed on it, and the layer with a radius of 135 mm contains 32 bubbles.The bubble cluster center is located 300 mm away from the wall.All bubbles have the same initial conditions: P 0 = 1.2 MPa, R 0 = 1.0 mm and Ṙ0 = 180 m/s.We designate three representative bubbles as bubbles 1, 2, and 3, as shown in Fig. 12(a).The maximum radius that a pulsating bubble can reach in a free field is defined as the reference length (7.6 mm).Fig. 12(b) presents the spatial distribution of bubbles at two typical moments in the dimensionless system, where the color denotes the amplitude of the vertical mi-    gration.At t * = 1.666, all bubbles are in their first pulsation cycle and show no significant migration tendency.As the distance from the wall increases, the migration tendency of the bubbles becomes weaker.However, at t * = 4.598, most of the bubbles, except for those still in the first pulsation cycle, migrate significantly toward the wall, resulting in the entire bubble cluster moving closer to the wall.This is reflected by the color bar values in Fig. 12(b).We also plot the radius of bubbles 1, 2, and 3 over time in Fig. 12(c) and the pressure evolution of the measuring point in Fig. 12(d).Bubble 3 has the shortest pulsation period, while bubble 1 has the longest, which is attributed to the difference in the number of surrounding bubbles, including real ones and mirror ones on the other side of the wall.During the collapse phase, the oscillation of bubble 1 accelerates because its surrounding bubbles rebound.These observations are similar to those of a planar bubble array observed in the previous study [36].The bubble energy formula [93] is employed to calculate the dimensionless bubble energy, and we found that the energy of bubble 1 in the first cycle increased from 4.19 when oscillating individually to 9.21 when in the bubble cluster.This indicates that some bubbles absorb the energy of other bubbles during bubble interaction.This also causes the wall pressure induced by some bubbles in the cluster to be more than two times higher than the pressure induced by bubble 1 alone.
Finally, we compute the dynamics of the bubble cluster under the hybrid-boundary condition, comprising a free surface and a vertical rigid wall, as shown in Fig. 13.The initial conditions of the bubbles are the same as those in Fig. 12. Bubbles 1, 2, and 3 are marked in Fig. 13(a).The bubble cluster center is 52 from the free surface and 79 from the rigid wall in the dimensionless system.Fig. 13(b) illustrates the relative bubble size and vertical bubble migration at two typical moments.The vertical migration of bubbles gradually increases as the distance from the free surface increases.The bubbles at the top of the bubble cluster exhibit downward migration due to the influence of the free surface, while the bubbles at the bottom migrate upward due to the weaker effect of the free surface.At t * = 0.950, the bubbles are all in the initial expansion stage with relatively uniform sizes.However, at t * = 4.598, the bubbles inside the cluster are relatively larger, which can also be reflected in the radius-time curve of the bubbles in Fig. 13(c).The maximum radius and period of bubble 3 significantly exceed those of bubbles 1 and 2, as the pulsation of the central bubble is slowed down by the surrounding contracting bubbles at the late stage of the first cycle.In Fig. 13(d), we present the wall pressure induced by the bubble cluster and bubble 1 alone.The peak pressure on the wall does not increase as much as that in Fig. 12 due to the effect of the free surface.Nevertheless, the wall pressure induced by the bubble cluster significantly exceeds that of the single bubble 1, with its pressure peak reaching close to twice that of single bubble 1.

Summarization and Conclusion
In this study, the unified equation for bubble dynamics [36] was rederived using the moving point source and dipole.We found that the form of the bubble pulsation equation was the same as that reported by Zhang et al. [36] considering the motion of singularities.The theory is applied to situations involving multiple bubbles and hybrid boundaries.The study's findings are summarized as follows: 1.The compressibility of bubble migration has an important effect on bubble dynamics behavior.The present theory fully considers the coupling between compressible bubble pulsation and migration in a compressible flow field.
2. A finite number of mirror bubbles are included to model the bubble dynamics near crossed boundaries.An infinite number of mirror bubbles must be introduced to account for the effect of parallel boundaries.Through calculations with various numbers of mirror bubbles for the parallel boundaries, all the results converge to within 0.1% when the number of mirror bubbles on the side of a boundary reaches 30.A good achievement is obtained between the theoretical result and the experiment and simulation, which demonstrate the capability of the present theory to predict bubble dynamics near hybrid boundaries.
3. A spherical bubble cluster containing more than one hundred cavitation bubbles is computed using the present theory under two boundary conditions: one is below a rigid wall and the other is near a vertical rigid wall below the free surface.The theoretical results indicated that some bubbles in the cluster could produce pressure peaks of more than double those of the individual bubble on the surrounding structures.
This study is based on the assumption of spherical bubbles, which could be used to predict the main characteristics of bubbles, including bubble size, period, migration, and pulsation pressure.For strong asymmetrical deformations of bubbles, such as splitting, breaking into pieces, and dissolving into the liquid, advanced numerical and experimental methods are required to conduct in-depth research.

Figure 1
Figure1Geometric relations between an arbitrary location r (marked by the green dot), the moving singularities, and the spherical bubble.The blue dot denotes the origin of the coordinate system o − rθϕ at time t.The singularities at t − R/C coincides with the blue dot.The red and orange dots represent the locations of the singularities at t and t − |r t |/C, respectively.

Figure 2
Figure2Schematic diagram depicting the image theory and mirror bubbles for the crossed boundaries with two angles (a) α = π/4 (b) α = π/2.The black circle denotes the real bubble, and the others denote mirror bubbles.The solid lines represent the boundaries.The origin of the coordinate system is at the initial bubble center and the direction of bubble migration is consistent with θ = 0, which is the same as that in Fig.1.

Figure 3
Figure 3 Schematic diagram depicting the mirror bubbles set to simulate the boundary effect containing crossed and parallel boundaries.Three solid lines represent three boundaries with boundary 3 perpendicular to the other two boundaries.The black circles on the left side of boundary 3 denote the real bubble, while the rest are mirror bubbles distributed alternately about boundaries 1 and 2. Circles on the right side of boundary 3 denote the mirror bubbles of all bubbles on the left side.

Figure 4
Figure 4 Convergence of (a) the bubble radius and (b) the flow field pressure with the number of mirror bubbles.The bubble is between two infinite parallel rigid walls (boundaries 1 and 2).The dimensionless distance between the two boundaries is 5.0 with the initial bubble 3.0 away from boundary 1, and the dimensionless initial conditions for the bubble are P * 0 = 20, R * 0 = 0.06 and Ṙ * 0 = 5.98.

Figure 7
Figure 7 Spark-induced bubble experiment and its comparisons of theoretical and experimental results near the bottom and vertical sidewall.The vertical sidewall on the right side is equivalently represented by two horizontally arranged spark-induced bubbles of the same size in the experiment.The distances of the initial bubble from the bottom and sidewall are 34 mm and 33 mm, respectively.(a) Experimental images at different moments.Frame width: 64 mm.(b) Bubble radius (in theoretical calculations: P 0 = 1.2 MPa, R 0 = 3.35 mm and Ṙ0 = 95 m/s.).(c) Vertical bubble migration.

Figure 8
Figure 8 Spark-induced bubble experiment and comparisons of theoretical and experimental results near the free surface and vertical sidewall.Two bubbles of the same size are generated in the experiment to equate a vertical sidewall in the mid-pendant of the two bubbles.The distances of the initial bubble from the free surface and sidewall are 49.5 mm and 35 mm, respectively.(a) Experimental images at different moments.Frame width: 65 mm.(b) Bubble radius (in theoretical calculations: P 0 = 1.2 MPa, R 0 = 2.98 mm and Ṙ0 = 90 m/s.).(c) Vertical bubble migration.

Figure 9
Figure 9  Comparison between the theoretical and EFEM simulation results for a bubble between the free surface and bottom.The distance between the two boundaries is 3 m and the initial bubble is equally distant from the free surface and bottom.In both theoretical and numerical calculations: P 0 = 60 MPa, R 0 = 32 mm and Ṙ0 = 0 m/s.(a) Bubble radius.(b) Vertical bubble migration.(c) Pulsation pressure (the measurement point is set at 1 m away from the initial bubble in the horizontal direction).

Figure 10 Figure 11
Figure10Comparison of the interaction of the three spark-induced bubbles with the theoretical results near the vertical sidewall, bottom and free surface.The effect of the sidewall is included in six horizontally equidistantly distributed bubbles induced by an identical voltage in the experiment.(a) Selected experimental images at typical moments.The sidewall is 20 mm from bubble 3. The three bubbles are initially located 49 mm from the free surface and 45 mm from the bottom.The distance between bubble 2 and bubble 1, as well as bubble 3, is 40 mm.Frame width: 119 mm.(b) Bubble radius (in theoretical calculations: R 01 = 2.36 mm, R 02 = 2.25 mm, R 03 = 2.41 mm, P 01 = P 02 = P 03 = 1.2 MPa and Ṙ01 = Ṙ02 = Ṙ03 =100 m/s).(c) Vertical bubble migration.
This work is funded by the National Natural Science Foundation of China (52088102, 51925904), the National Key R&D Program of China (2022YFC2803500, 2018YFC0308900), Finance Science and Technology Project of Hainan Province (ZDKJ2021020), and the Xplore Prize.

Figure 12
Figure 12 Computed dynamics of a spherical bubble cluster with 105 cavitation bubbles below a rigid wall.(a) Schematic diagram of the bubble cluster in the xz plane.The black bold solid line denotes the rigid wall on which a pressure measurement point is placed directly above the bubble cluster.(b) The relative size and vertical migration of bubbles at two typical moments (P * 0 = 11.55,R * 0 = 0.13, Ṙ * 0 = 17.65).(c) Bubble radius for three typical bubbles.(d) Comparison of pressure at the measuring point between the bubble cluster and single bubble 1.

Figure 13
Figure 13 Computed dynamics of a spherical bubble cluster with 105 cavitation bubbles near a vertical rigid wall and a free surface.(a) Schematic diagram of the bubble cluster and boundary conditions in the yz plane.The rigid wall is located at y * = 79 and the free surface is located at z * = 0.The bubble cluster center is at x * = y * = 0 and z * = −52.The pressure measurement point on the rigid wall is horizontal to the bubble cluster center.(b) Relative size and vertical migration of bubbles at two typical moments.The initial conditions for bubble pulsation are the same as those in Fig. 12. (c) Bubble radius for three typical bubbles.(d) Comparison of pressure at the measuring point between the bubble cluster and single bubble 1.