Jetting mechanisms in bubble-pair interactions

The understanding of bubble dynamics is important for a variety of technological applications ranging from ship propellers¹ to modern medicine.^2–5 Therefore, it has been extensively studied by means of experimental and numerical research.^6–12 Hawker et al.¹³ investigated the free-field shock-induced collapse and jetting mechanism of a single bubble. They identified three characteristic phases of the collapse: the incident shock, water-hammer shock, and vortex ring. Johnsen et al.¹⁴ studied a similar setup but with a rigid wall positioned close to the bubble. They described a re-entrant jet, triggering a strong water-hammer pressure shock, which may lead to very high wall pressures.

Using a viscous fluid instead of the wall, tissue-like materials can be modeled.^15–19 The viscosity of the tissue-like material damps the reflection of the incoming incident shock.¹⁷ Hence, medical procedures, such as sonoporation, can be investigated.²⁰ 

Chew et al.²¹ experimentally investigated the collapse behavior of bubbles created by spark discharge. For a wide range of bubble distances and ratios of radii, they found that bubbles oscillating in phase generate jets toward each other. Inversely, out-of-phase bubbles generate jets away from each other. Similarly, Han et al.²² studied the collapse of a bubble pair with an axisymmetric description using a boundary element method based on a potential flow theory. They identified the three different collapse mechanisms of “coalescence,” “jet toward,” and “jet away.” Matching the results of Chew et al.,²¹ they found in-phase bubbles jetting toward and out-of-phase bubbles away from each other.

Recently, Hopfes et al.²³ proposed a novel experimental method for the investigation of pure, non-condensable gas bubbles. In their approach, bubbles are placed into a Gelrite^™ Gellan Gum mixture. The so-called gelatin acts as a proxy for water while allowing to create resting gas bubbles of desired size at atmospheric pressure. Within this setup, they study the free-field shock-induced collapse of two differently sized bubbles in four distinct configurations: shooting through, weak interaction, reversing collapse, and collapse toward. They find good agreement with the theory of Chew et al.²¹ for most but not all of their configurations. Our motivation for the numerical study of bubble-pair collapse is twofold. At first, we want to gain insight into mechanisms leading to the distinct collapse behaviors observed by Hopfes et al. Second, we want to further analyze postulates of Chew et al. about a critical phase difference steering the jetting behavior in the case of initially quiescent bubble pairs. For this purpose, we conduct six three-dimensional (3D) simulations with $17×109$ cells of effective resolution. Next to the four cases by Hopfes et al., two new cases are considered. The simulations are carried out using ALPACA.^24–26 Herein, we use a finite volume approach with fifth-order low-dissipation shock-capturing weighted essentially non-oscillatory (WENO) reconstruction.²⁷ The convective interface fluxes are calculated via the Harten, Lax, and van Leer approximate Riemann solver with contact restoration (HLLC).²⁸ A third order total-variation-diminishing (TVD) Runge–Kutta scheme is used for explicit time integration.^29,30 The gas–liquid interface is described via a sharp-interface level-set method³¹ following re-initialization.³² The two-phase Riemann problem is solved using a linearized Riemann solver.³³ A block-structured multiresolution scheme³⁴ with adaptive local time-stepping (ALTS)³⁰ enables full local and temporal adaptivity.

The remainder of this work is organized as follows: Sec. II presents the governing equations, the employed numerical methods, and the compute environment. In Sec. III, we present the chosen problem setup in which we also examine a single-bubble collapse. We investigate and discuss the bubble-pair collapse in Sec. IV. Concluding remarks and an outlook on future work are given in Sec. V.

The bubble collapse is described by the compressible Navier–Stokes equations (NSEs). Due to the small bubble sizes in our study, however, we neglect gravitational forces to obtain

Therein, $U=(ρ, ρu, ρE)T$ is the vector of conservative states, in which ρ denotes the density, u is the velocity vector, and $E=e+|u|2/2$ is the specific total energy with the specific internal energy e.

The right-hand side X accounts for inter-phase exchange of momentum and energy. The convective and viscous flux densities are

with the pressure p and the identity matrix I. The heat flux is denoted by q. Due to the short timescales, however, we neglect heat transfer allowing to drop this term. The viscous stress tensor T is given by

with the dynamic viscosity μ.

Accounting for the two-phase interaction, we rewrite Eq. (1) in the flux-based integral form as

The two immisicible phases fill a material volume $V$⁠. Subscripts l and g represent the gaseous and liquid phases, respectively. The interface surface $S$ separates the phases. We define the interface normal vector $nS$ pointing toward $Vl$⁠. The convective and viscous interface flux densities $FSc$ and $FSμ$⁠, respectively, are identical to the convective and viscous flux densities in Eq. (2), except that they are evaluated at the phase interface. Capillary forces are taken into account by the stress balance in the interface normal direction, described in more detail in Sec. II B.

We close the set of equations with the stiffened-gas³⁵ equation-of-state (EOS) as it describes both liquid and gaseous phases. Using the polytropic coefficient κ and by defining a background pressure $p∞$⁠, it reads

Setting $p∞=0$⁠, the equation degenerates to the ideal-gas equation.

We conduct our simulations with the sharp-interface multiresolution framework ALPACA.^25,26 Here, we briefly recite the most important numerical concepts used within this solver. For a detailed description as well as thorough validation of the framework, the interested reader is referred to the given literature. Note, ALPACA provides a wide variety of numerical methods to choose from. Hence, we shortly comment on the selected configuration of numerical building-blocks as well.

Section II A is discretized via a Cartesian finite-volume approach. Cubic cells with cell volume $(Δx)3$ are used. We compute the convective flux vector $Fc$ using the HLLC approximate Riemann solver²⁸ and employ a fifth-order WENO stencil²⁷ for the inherent cell-face reconstruction. The cell-face viscous flux densities are obtained by computation of velocity gradients in the cell center with a fourth-order central stencil and a subsequent reconstruction to cell-faces using a fifth-order WENO stencil.

We compute the time-evolving fluid interface $S$ using a level-set method with conservative convective³¹ and viscous interface interaction.³⁶ The level-set function,

indicates the distance to the interface. From this function, the interface is geometrically approximated by flat surfaces in each cell. We update the interface shape and position by

where the interface normal is given by

The interface normal velocity $uS$ is determined from a linearized two-material Riemann problem.³³ The stress balance in the interface normal direction with surface tension effects reads as

where $Π=−pI+T$ denotes the stress tensor and σ is the coefficient of surface tension.

After the advection step, the re-initialization equation

is solved iteratively in pseudo time τ to recover the signed-distance property $|∇·ϕ|=1$⁠.³² 

We advance the simulation in time using an explicit third-order TVD Runge–Kutta scheme.²⁹ The time step size $Δt$ is restricted to ensure stability accounting for acoustic, viscous, and capillary effects as stated in Ref. 25 with a Courant–Friedrichs–Lewy (CFL) number of CFL = 0.6.

We rely on ALPACA's multi-resolution (MR) compression^26,37 and its ALTS scheme¹⁹ to reduce the computational load. Deitering et al.³⁸ showed that adaptive mesh refinement (AMR) and MR schemes reach comparable compression rates. The latter, however, refines smooth small-scale flow phenomena such as rarefaction fans more reliably. Still, our highly resolved simulations require the usage of high performance computing (HPC) computer hardware. We use a designated partition of the of the CoolMUC-2 Linux Cluster at the Leibniz Supercomputing Center.³⁹ The system consists of five 40-way Intel Xeon Gold 6230 compute nodes at 2.1 GHz.

The bubble collapse is triggered by a weak planar shock wave with $p1=1×105 Pa$ and $p2=3×105 Pa$⁠, where the subscripts 1 and 2 denote quantities up- and downstream of the shock wave, respectively. This shock wave is consistent with the experiments of Hopfes et al.²³ and yields a Mach number of $MS=1.000 14$⁠. According to Ref. 23, the surrounding fluid is modeled as water with the background pressure and polytropic coefficient taken from Ref. 14. All other parameters are chosen following Ref. 40. The material parameters are listed in Table I. Note that the post-shock conditions are defined by the Rankine–Hugoniot relations⁴¹ for the chosen Mach number of $MS=1.000 14$⁠. The Reynolds and Weber numbers quantify the influences of viscous and capillary forces, respectively, in relation to inertial forces. Following Ref. 15, we use the pressure jump $PS=p2−p1$ to define both non-dimensional numbers. The diameter of the smaller bubble of the weak interaction case is used as length scale $D0,min$⁠. The Reynolds number,

and the Weber number,

indicate considerable contributions of both viscous and capillary forces.

We simulate the collapse of an air bubble in water to assess the numerical setup and the material parameters and obtain an estimate of the required resolution and computational load. We conduct the simulation in three resolutions $R0/Δx=25, 50, and 100$⁠, wherein R₀ indicates the initial bubble radius. Exploiting spherical symmetry suffices to simulate a quarter of the bubble. Apart from the symmetry planes, zero-gradient boundary conditions for all required states are prescribed at all external boundaries. Figure 1 shows a two-dimensional (2D) projection of the computational domain.

Note that we slightly adapt the setup for this reference simulation to match with the original work of Johnsen et al.¹⁴ Hence, we set $p2/p1=353$⁠. We use the material parameters introduced in Sec. III A. Only the density of water is set to $ρ1,l=998 kg/m3$ to be consistent with Johnsen et al. Post-shock quantities follow from the Rankine–Hugoniot relations.⁴¹ For comparison with Ref. 14, we non-dimensionalize the velocity by $cl=1647 m/s$⁠, pressure by $ρlcl2$⁠, and time by $R0/cl$⁠.

Figure 2 depicts the velocity field and the corresponding pressure field over time for the collapse using the finest grid resolution. With this figures, we introduce the labeling used throughout this manuscript: capital letters refer to the velocity field and small letters to the pressure field. Additionally, we will use numbered arrows to highlight certain features of the flow fields. Regarding the velocity field, the color-scale shows its magnitude while the glyphs indicate its direction.

In subfigure (a), we see the shock-wave (1) having passed the bubble and the transmitted shock (2) as described by Hawker et al.¹³ The pressure field exhibits a transition from post-shock states to the lower pressure inside the bubble (3). In the following, we call this region of decreasing pressure the pressure-gradient region (PGR) of a bubble. In subfigures (B) and (b), the collapse becomes non-spherical due to the baroclinic vorticity generation. After that, a jet forms in the direction of the shock-wave advection (C, 4). Due to its collapse, the pressure inside the bubble in (b) has risen above the pressure level caused by the transmitted shock. The bow shock (C) and (c) and the water-hammer shock (D) and (d) as described by Hawker et al.¹³ are also clearly visible.

A time history of the jet velocities is given in the upper part of Fig. 3. We observe convergence for the jet velocity with increasing numerical resolution. During the first phase of the collapse until $tcl/R0≈7$⁠, all curves match very well. After that, the reference collapse of Ref. 14 evolves more gradual and produces a lower peak jet velocity compared to our simulations. Note, Johnsen et al. used a diffusive interface model rather than the sharp level-set description employed here. Although we consider viscous and capillary forces, which are expected to slow down the jet, it reaches a higher velocity in our simulation. Hence, we attribute the slight disagreement to the different interface treatment. The lower part of Fig. 3 shows the non-dimensionalized bubble volume over time. We find consistent temporal evolutions for all resolutions in good agreement with Ref. 14.

The collapse behavior of bubble pairs is characterized by three non-dimensional parameters:

• size ratio:

  $S=RLRS,$

  (13)
• relative distance:

  $γ=dRL+RS≥1,$

  (14)
• phase difference:

  $Δθ=|(t1tosc,1−t2tosc,2)|−Δttosc,1,$

  (15)

defined according to Ref. 21. Herein, d is the distance between the bubble centers, and R_L and R_S are the radii of the large and small bubbles, respectively. The phase difference $Δθ$ is defined for experimental configurations in which bubbles have to be created in water first, e.g., by spark discharge. Hence, $Δt$ is the time difference between nucleation of both bubbles. t₁ and t₂ denote the elapsed time from nucleation of the respective bubble to the first completed collapse of one of the bubbles. $tosc,1$ and $tosc,2$ are the oscillation periods of the respective bubble. As in the initial configuration of Ref. 23, the bubbles are quiescent, $Δt=0$ and Eq. (15) degenerates to

Figure 6 depicts the initial configuration of a bubble pair. Depending on these parameters, Hopfes et al.²³ have identified the following four distinct cases of bubble-pair interactions:

• Shooting through can be observed, if two bubbles with a large size ratio are located close to each other, cf. Fig. 4(a). This leads to a spherical collapse of the large bubble with the concurrent collapse of the small bubble. The bubbles merge and form a new spherical bubble. During the rebound phase, puncturing is observed as a result of the induced jet from the collapse of the small bubble. This creates a new small bubble splitting off from the large one at the opposite side of the initial small bubble.

• Weak interaction is characterized by a large size ratio between the two bubbles and a large relative distance, cf. Fig. 4(b). As the name suggests, the interaction is limited, i.e., a jetting behavior cannot be observed. The general evolution is dominated by the large bubble. Its radius changes periodically due to collapse and subsequent rebound. The behavior of the small bubble is determined by two superimposed movements. First, the radius of the small bubble also changes periodically but with higher frequency. At the same time, its center follows the motion of the large one, i.e., during a collapse of the large bubble, the small one moves toward it and vice versa during a rebound. This leads to a splitting of the small bubble, whereas the large bubble is little affected by the presence of the small bubble.

• When two bubbles with size ratio at around S = 1.7 are located next to each other, a reversing collapse occurs, cf. Fig. 4(c). Initially, both bubbles contract spherically. The faster collapse of the smaller bubble leads to a non-spherical deformation. This produces a jet toward the other bubble, which quickly disintegrates the smaller bubble. In turn, it expands away from the larger one. The larger bubble forms a jet directed away from the smaller bubble. The jet punctures the larger bubble and produces a new small bubble, which moves away from the large bubble.

• Collapse toward can be observed, if two bubbles of similar size are located in a large relative distance, cf. Fig. 4(d). During the first collapse, both bubbles contract spherically. When rebounding, they become non-spherical and develop a liquid jet toward each other. This case is similar to non-spherical shock-induced bubble collapse and bubble collapse near a wall.

Figure 5 shows a map for these bubble-interaction mechanisms in the γ-$Δθ$-plane. The two lines represent the critical phase differences found by Chew et al.²¹ The solid line represents cases with similar sized bubbles. Cases with differently sized bubbles are categorized by the dashed line. For cases above and below respective the line, Chew et al.²¹ have observed jets away from and toward the bubbles, respectively. Additionally, two intermediate cases are denoted in Fig. 5, which are located among the original four cases. We have chosen the cases in an attempt to provoke a superposition of the previously observed collapse mechanisms. Also, the intermediate case two is located in between Chew's critical lines for similarly and differently sized bubbles. This makes it a curious case to test the validity of their taxonomy. Each mechanism is demonstrated with an exemplary parameter set, see Table II. For the four collapse mechanisms of Ref. 23, the bubble sizes are chosen to be consistent with the experiments. Regarding the intermediate cases, the parameters are selected such that the resulting non-dimensionalized parameters are located in the γ-$Δθ$-plane as described above. Furthermore, the bubble sizes in the intermediate cases are chosen to be similar to the original four cases.

Figure 6 sketches the simulation setup for the bubble-pair simulations in terms of the bubble radii and d. The domain has a size of $13 RL,max×26 RL,max×13 RL,max$⁠, where $RL,max$ is the radius of the large bubble of the shooting through case, as this is the largest simulated bubble. Note that the fluid properties for all bubble pair interactions are presented in Sec. III A. Due to the low Mach number, the bubbles respond to the elevated post-shock pressure rather than the pressure gradient and, hence, can be positioned either way relative to the shock-wave. Furthermore, exploiting the spherical symmetry of the problem, we simulate a quarter of the domain and, hence, the bubbles. Regarding asymmetries in the experiments by Hopfes et al.,²³ we believe imperfect alignment of the bubbles with the normal of the incident shock wave to be a driver of such asymmetries. The boundary conditions are set as in Sec. III B.

In the following, the term near side refers to the side of the bubble pointing toward the other one. The term far side refers to the side of the bubble, which points away from the other bubble. Additionally, the bubble, which is initially closer to the shock wave, is called an upstream bubble. The other bubble is called a downstream bubble. Furthermore, in each case, a purely descriptive distinction between the small and the large bubble is made. The small bubble of one case may be of similar size as the large bubble in another case. The terms do not imply an absolute measure.

In the shooting through case, the resolution is $R0S/Δx=37$ and $R0L/Δx=154$⁠. Figure 7 shows the velocity field and the corresponding pressure field. We observe both bubbles merging during the first collapse (A)–(D). In this process, gas from the small bubble is injected into the large one while the connecting funnel between the bubbles acts like a nozzle (B) and (C). The gas is accelerated to velocities of up to $150 m/s$ within a gas jet shooting through the large bubble during its first collapse (C)–(F). During the first rebound of the remaining bubble, the gas jet hits the far side and forms a funnel-shaped protrusion [(F) and (G), 1]. This protrusion disappears during the second collapse (I)–(K). The pressure field exhibits a cyclic behavior. In the collapse phases, the pressure inside the bubbles increases, which also leads to a pressure rise in their vicinity (a)–(f). The subsequent rebound causes a decrease in pressure inside the bubbles leading to a pressure decrease in their surroundings (f)–(h).

For this scenario, we see a distinctive difference to the experiments of Hopfes et al.²³ Based on their experimental imaging, they report a liquid jet leading to the funnel-shaped protrusion, which is visible both in subfigures (10–12) of Fig. 4(a) and subfigures (G) and (H) of Fig. 7. In our simulations, however, the evolution is rather caused by a gas jet. Unlike the experiments, we do not observe a bubble being shot out of the protrusion. We assume our simulation to be under-resolved A reason may be the numerical resolution that may be insufficient to capture this phenomenon, as the radius of a new bubble would be resolved by approximately five cells only. The bulge on the downstream side of the remaining bubble in (K) and (k) can also be observed in the experiments.

We explain the observed behavior as follows. Initially, the bubbles are close to each other, so that the small bubble is positioned almost fully in the PGR of the large bubble. Hence, the high post-shock pressure cannot be fully established between them (a). Due to its size, the small bubble collapses faster. Thereby, it is exposed to a pressure gradient due to the still lower pressure of the large bubble on its near side and the elevated post-shock pressure on its remaining surface (a). This leads to the collapse of the small bubble toward the large one, whereby both bubbles merge (B) and (b). Due to the pressure difference, the gas from the small bubble flows into the large bubble (B) and (b) and creates the gas jet (C)–(F). The jet pulls the interface of the bubble inward at the near side [(F) and (G), 2]. As soon as the jet vanishes, the interface moves back outward during the next collapse [(H) and (I), 3]. This leads to the bulge visible both in (K) and (k) and the experiments.

Following the categorization of Chew et al.,²¹ this case lies in the jet-away region. For the respective combination of phase difference $Δθ=0.76$ and relative distance $γ=0.99$⁠, they have found jet-away cases with liquid jets.²¹ Although we do not observe a liquid jet, the gas jet points away from the small initial bubble, matching Chew's terminology.

We observe some anomalies in the velocity field which we attribute to numeric grid effects (I) which can be mitigated with higher numerical resolution⁴² and do not invalidate the overall evolution.

Now, we investigate the weak interaction case. In this case, the bubbles are resolved with $R0S/Δx=26$ and $R0L/Δx=102$⁠. Figure 8 shows the velocity field and the corresponding pressure field.

We observe the small bubble to follow the motion of the large one: during the first collapse of the large bubble, the small one moves toward it [(A)–(E), 1]. After that, the small bubbles moves away from the large bubble while the latter rebounds [(F)–(H), 2]. The same behavior is observed during the second collapse and rebound of the large bubble [(I)–(M), 3 and 4]. It is noticeable that with the beginning of its first rebound (F) and (f), the large bubble completely dominates the pressure field. Until this point, the radius of small bubble changes periodically at its own frequency (A)–(E) and, hence, determines the pressure field in its direct vicinity (a)–(e). After that, the oscillation frequency of the small bubble is in phase with one of the large bubbles. Additionally, the small bubble cannot establish an own pressure field in its vicinity anymore.

In contrast to the results of Hopfes et al.,²³ the reported splitting of the small bubble cannot be observed in our simulations. In accordance with the experiments, we see the small bubble following the oscillations of the large one without affecting it. Neither the experiments nor the simulation show jetting behavior of either bubble.

The collapse behavior is determined by the large relative distance. It leads to only the edges of the PGR of each bubble overlapping. Hence, the pressure difference between the near and far sides of both bubbles after the shock wave has passed is small (a). This leads to a spherical first collapse of both bubbles. Due to its size, the smaller bubble contracts faster than the large one. Therefore, the pressure inside the small bubble rises faster than in the large bubble. Hence, the small bubble is exposed to a higher pressure on its far side than on its near side (b). This leads to an asymmetric velocity field inside the small bubble (B) and (C) and its motion toward the large one. Meanwhile, the pressure inside the large bubble has risen as well (e) and (f), leading to an inversion of the velocity field inside the small bubble (E)–(G). This generates a movement away from the large bubble (F)–(H). During the second collapse and rebound of the large bubble, both bubbles oscillate in phase.

In the taxonomy of Chew et al., also no jetting is observed for $γ=1.74$ and $Δθ=0.74$ as in this case. Hence, the bubbles may be too far from each other to develop a jet after all.

Next, we investigate the reversing collapse case with a resolution of $R0S/Δx=64$ and $R0L/Δx=111$⁠. Figure 9 shows the velocity field and the corresponding pressure field.

During its first collapse, the smaller bubble moves and deforms in the direction of the large one (B, 1). We observe a gas jet shooting through the bubble and toward the large one (C, 2). Due to the reversing motion, the near-side interface of the smaller bubble moves inward and collides with the jet [(D) and (E), 3]. After that, a weak gas jet shoots from the near side of the large bubble toward its far side [(F)–(H), 4]. Similarly, a weak liquid jet punctures the smaller bubble [(G) and (H), 5].

For this case, Hopfes et al.²³ reported the name-giving reversing motion and two liquid jets. The first jet shoots from the smaller bubble toward the large one; the second jet penetrates the large bubble and causes a small bubble to detach from it. In our simulations, we also observe the reversing motion. Similar to the previous shooting through case, both jets appear as gas jets here. The second jet is too weak to reach the far side of the large bubble and, therefore, does not cause bubble detachment. The weak liquid jet puncturing the smaller bubble in (H) has not been reported by Hopfes et al. However, their experimental images hint at its existence, cf. the marked areas in frames 13–15 of Fig. 4(c). Apart from that, Hopfes et al. described parts of the smaller bubble splitting off during its first collapse, which is not visible here.

Examining the velocity and pressure fields, we see the PGR of both bubbles overlapping so that the high post-shock pressure cannot establish between them (a). Due to the faster collapse, the smaller bubble sees the high post-shock pressure on its far side and the still lower pressure of the large bubble on its near side (b). This pressure difference causes an acceleration and, hence, the movement and deformation in the direction of the larger bubble (B) and (C) and (b) and (c). The far side of the smaller bubble moves further inward, which leads to the gas jet shooting through the bubble toward the large one (C). Meanwhile, the large bubble has collapsed leading to a high pressure in its vicinity (c), which in turn leads to an inversion of the motion of the smaller bubble (C) and (D). Therefore, the near side of the smaller bubble moves inward and collides with the gas jet (E) and (F). Apart from that, the reversing motion creates the weak liquid jet, which punctures its far side (G) and (H). During its first collapse and rebound, the large bubble is exposed to a pressure gradient, as on its near side the small bubble generates a pressure, which is above the pressure on the far side of the large bubble (b) and (d) and (e). Due to this pressure difference, the large bubble shows a non-spherical behavior (D)–(F) and (d)–(f). This causes the weak gas jet-away from the smaller bubble visible in (F)–(H).

Our simulations show similarities between the weak interaction case, cf. Sec. IV B, and the current one. In both cases, the pressure difference between the near and far sides of the smaller bubble induces an internal momentum. Hence, the small bubble moves toward the large one. The pressure increase in the surrounding of the large bubble inverses this motion, cf. Figs. 9(c) and 9(d) and 8(f) and 8(g). Due to the smaller γ here, the inversion is stronger. This is the root cause for the differing jetting and collapse behavior of the two cases.

Regarding the critical phase differences found by Chew et al.,²¹ this case is in the jet-away region. In the simulations, we first see a gas jet-toward (C)–(E), then a liquid and a gas jet-away (F)–(H). This mixed behavior can be explained with its position close to the critical phase difference in Fig. 5.

The velocity field and the corresponding pressure field of the collapse toward case are depicted in Fig. 10. The resolution is $R0S/Δx=86$ and $R0L/Δx=93$⁠.

During the first two oscillations, the bubbles move closer to each other. We observe a spherical collapse and rebound of both bubbles (A)–(C). The second collapse, however, becomes non-spherical (D)–(H). After that, diagonal liquid jets penetrate the far side of both bubbles [(H) and (I), 1 and 2]. They are reversed in the subsequent rebound only to reappear in the following collapse [(L), 3 and 4]. Additionally, the upstream bubble forms an axial liquid jet [(K)–(M), 5]. Subsequently, a part of the downstream bubble is split off (N). Afterward it is punctured by a liquid jet [(N)–(O), 6].

Comparing with Ref. 23, we also observe a transition from spherical to non-spherical collapse. However, in our simulations, the first collapse and rebound are spherical, whereas in the experiments, the first rebound is non-spherical already. We can confirm the liquid jets toward each other reported by Hopfes et al. However, instead of the diagonal behavior they describe axial jets. Nevertheless, the experimental images hint at the existence of diagonal jets, see the marked areas in frames 14 and 15 of Sec. III C.

Our simulation shows how the large γ imposes the collapse behavior, as only the edges of the PGRs overlap. Only a small pressure difference between the near and far sides of the bubbles establishes (a). Hence, the first collapse and rebound of both bubbles are spherical (B) and (C), (b) and (c). However, due to the small pressure difference, the velocity field inside the bubbles is not symmetric (B). For this reason, the bubbles move toward each other (C) and (c). The reduced distance increases the pressure difference between the near and far sides of each bubble (c). Therefore, the asymmetry in the velocity field inside the bubbles is amplified (D). This leads to a non-spherical second collapse and subsequent non-spherical rebound (D)–(F) and (d)–(f). The again reduced distance increases the pressure difference between the near and far sides of each bubble further (F) and (f). Together with the non-spherical shape, this yields the diagonal liquid jets during the next collapse (H) and (I) and the gas flowing radially outward during the subsequent rebound, which amplifies the asymmetry (I). The bubbles still move toward each other. This leads to stronger jets with every collapse (L) and (M) and (N) and (O) and eventually to the splitting and puncturing of the smaller bubble.

Following the classification of Ref. 21, this case lies in the jet-toward area, see Fig. 5. The appearing diagonal jets in [(H), (L) and (M)] as well as the axial jets in (K)–(M) and (N)–(O) confirm this postulate.

The velocity and pressure fields of this case are depicted in Fig. 11. The resolution is $R0S/Δx=39$ and $R0L/Δx=133$⁠.

We observe the small bubble following the motion of the large one during the first collapse [(A)–(C), 1] and subsequent rebound [(D) and (E), 2] of both bubbles. The rebound of the large bubble induces a liquid jet penetrating the small one [(D)–(F), 2]. After that, a liquid jet forms at the far side of the small bubble and punctures it [(G)–(I), 3]. Apart from a deformation in the direction of the smaller bubble during the second collapse and rebound (G)–(I), the large bubble remains unaffected by the interaction.

Examining the velocity and pressure fields, we see the PGRs of both bubbles overlap so that the post-shock pressure cannot establish between them (a). Hence, the smaller bubble moves and deforms in the direction of the large one (B), as observed in the reversing collapse case, see Sec. IV C. As the far sides of both bubbles are exposed to a larger pressure gradient than the near sides, they experience a stronger acceleration during the first collapse (b) and (C). The subsequent rebound of the large bubble reverses the motion of the small one (D). This leads to a liquid jet penetrating the small bubble and in turn to a non-spherical shape of it (D) and (E). During the second collapse, the far sides of both bubbles are exposed to a higher pressure gradient again (f). In the case of the large bubble, this leads to an asymmetric collapse and subsequent rebound (H) and (I). The small bubble, however, develops a liquid jet leading to its puncturing (G)–(I).

Comparing this case with the previous ones, we mainly observe features of weak interaction and reversing collapse, cf. Secs. IV B and IV C, respectively. As it remains largely unaffected by the interaction, the behavior of the large bubble corresponds to the weak interaction case. The initial behavior (B), the reversing motion (D), and the first liquid jet (D) and (E) of the small bubble are characteristics for the reversing collapse case. However, the second liquid jet does not appear in reversing collapse. This jet forms due to the acceleration of the far side of the bubble caused by the pressure gradient (G)–(I). The same mechanism leads to the non-spherical second collapse and subsequent rebound of the large bubble (F)–(K). This behavior was observed during the collapse toward case. Notably, the combination of reversing collapse- and collapse toward-like behavior leads to a bidirectional liquid jetting of the smaller bubble, which has not been observed in the previous cases.

One could expect shooting through-like behavior in this case as it lies between this and the weak interaction case in Fig. 5. However, the bubbles seem to be located too far apart from each other to show this behavior. This matches the results of Hopfes et al.²³ as they report a mixed case of reversing collapse and weak interaction for this combination of γ and $Δθ$⁠.

Following the classification of Chew et al.,²¹ this case lies in the jet-away area. Hence, the first jet follows their proposal. The second jet, however, is a jet-toward, as it is a collapse toward-like mechanism.

The resolution of this case is $R0S/Δx=71$ and $R0L/Δx=102$⁠. Figure 12 shows the velocity field and the corresponding pressure field.

We observe the bubbles moving closer to each other and becoming non-spherical during the first two collapses (A)–(D). During the third collapse, a diagonal liquid jet forms at the far side of the small bubble [(G), 1]. The subsequent rebound of the large bubble reverses the motion of the small one and induces a liquid jet penetrating the near side of the small bubble [(G)–(I), 2]. During the following collapse, both bubbles develop axial liquid jets on their far sides [(K)–(M), 3 and 4]. In the case of the small bubble, the jet punctures the bubble (M).

Our simulations show that only the edges of the PGRs of both bubbles overlap (a). This leads to a spherical first collapse and subsequent rebound (A)–(C). However, the pressure gradient on the far sides of both bubbles is larger than on the near sides. Hence, the velocity field becomes asymmetric, and the bubbles move closer to each other. During the second collapse and rebound, both bubbles become non-spherical as the gas flows radially outward (C) and (D). Additionally, the bubbles have moved even closer to each other. As in the collapse toward case, the interaction of the bubbles becomes stronger with decreasing distance. Hence, the smaller bubble develops a diagonal liquid jet on its far side (G). Due to the size difference, the bubble oscillations become more and more out of phase with the ongoing interaction. This leads to the large bubble reversing the motion of the small one (G) and (H), which in turn causes a liquid jet penetrating the small bubble (G)–(I). As the distance between the bubbles has further decreased, both of them form a liquid jet during the next collapse (K)–(M).

In this case, we see strong similarities with collapse toward, cf. Sec. IV D. The movement toward each other while becoming non-spherical during the first oscillations and the liquid jets toward each other are characteristic for this case. Additionally, we also observe a diagonal jet here. However, as the bubbles have a larger size ratio than in collapse toward, they also have a larger $Δθ$⁠. This causes the reversing motion and the liquid jet-away characteristic for the reversing collapse case. As in the intermediate case one, the combination of reversing collapse- and collapse toward-like behavior leads to a bidirectional liquid jetting of the smaller bubble.

Regarding the critical phase difference found by Chew et al.,²¹ this case lies in the jet-toward area. Hence, the collapse toward-like jets follow their proposal. However, the jet-away originating from the reversing motion does not fit into this classification.

We have investigated six collapsing bubble pair configurations with respect to the induced jetting mechanisms. The collapse toward case, see Sec. IV D, is characterized by large γ and small $Δθ$ values. During the first two oscillations, the bubbles move toward each other. The closer the bubbles get, the higher the pressure difference between the near and far side of each bubble. Apart from that, the bubbles become more asymmetric with each completed collapse and subsequent rebound. As soon as they are about $γ=1.34$ away from one another, they develop liquid jets directed toward each other. Before axial jets form, we observe diagonal jets, which might be caused by the asymmetry of the bubbles. Such jets have not been stated explicitly by Hopfes et al.;²³ nevertheless, their experimental images hint at their existence. As the collapse toward case lies in the jet-toward area, the direction of the jets fits the classification of Chew et al.²¹ 

Different jetting mechanisms are observed for the reversing collapse case, see Sec. IV C. It is characterized by small γ and moderate $Δθ$ values. At first, the phase difference imposes a pressure difference between the bubbles and, subsequently, causes jetting inside the small bubble toward the large one. Then, the oscillation of the larger bubble reverses the motion of the smaller one. This induces a liquid jet-away which leads to the puncturing of the smaller bubble. Finally, the higher pressure inside the smaller bubble causes a weak gas jet-away inside the larger bubble. Hence, two of the three jets fit the taxonomy of Chew et al.,²¹ as this case lies in the jet-away area. However, the investigated configurations differ from Chew's experiments in three ways. First, Chew et al. only described the jets, which finally puncture the bubbles and second, distinguished only between pure jet-toward and jet-away cases. Third, they exclusively described liquid jets. Regarding the reversing collapse case, this third difference excludes the observed gas jet-toward from the experiment's scope.

The intermediate case two shows characteristics of both previously discussed cases, cf. Sec. IV F. This matches the expected behavior from the problem setup, see Fig. 5. On the one hand, as for the reversing collapse case, we observe liquid jet-away behavior penetrating the smaller bubble. On the other hand, we see diagonal and axial liquid jets-toward as in the collapse toward case. Following the classification of Chew et al.,²¹ this case lies in the jet-toward area. Hence, the bidirectional liquid jetting of the smaller bubble does not comply with their theory. As mentioned above, however, Chew et al. described solely the jets, which puncture the bubbles. For this reason, this case could also fit their taxonomy when only the final set of jets-toward is taken into account.

Finally, we have investigated three cases with large $Δθ$ values. In combination with large γ values, the weak interaction case occurs, see Sec. IV B. In this case, the bubbles are too far from each other to develop a jet. The large size difference prevents an approaching of the bubbles as observed in the collapse toward case. This fits both the results of Hopfes et al.²³ and Chew et al.,²¹ as the latter do not report jetting cases for this combination of $Δθ$ and γ values.

Putting the bubbles closer together, we obtain the intermediate case one, cf. Sec. IV E. This case shows mixed behavior of the weak interaction and the reversing collapse case. Nevertheless, we also observed a collapse toward-like mechanism. It generates a liquid jet-toward, which ultimately punctures the small bubble. This bidirectional liquid jetting of the smaller bubble is similar to the one observed in the intermediate case two and does also not fit the taxonomy of Chew et al.,²¹ as this case lies in the jet-away region. Hence, the non-compliant jet is the one finally puncturing the smaller bubble. For this reason, we cannot make the same argument as in the intermediate case two and report a difference between our simulations and the classification of Chew et al.

For small γ values, the shooting through case is observed, see Sec. IV A. Here, the bubbles are so close together that the small one collapses into the large one. This leads to gas advecting from the small bubble into the large one. Therein, the merged section of the interface acts as a nozzle and infers gas velocities $>150 m/s$⁠; three times higher than the velocities observed in the other cases. This is the name-giving gas jet shooting through the large bubble away from the small bubble. Since this case lies in the jet-away region, it fits the taxonomy of Chew et al.²¹ 

Altogether, we observed two major differences to the experiments of Hopfes et al.²³ First, they described liquid jets in the shooting through and reversing collapse case, which appear as gas jets here. Second, we do not see satellite bubbles in shooting through and reversing collapse. Additionally, Hopfes et al. reported a splitting of the smaller bubble during weak interaction and reversing collapse.

Regarding the first difference, Hopfes et al. stated that liquid jets cannot be observed directly in their results but are deduced from indications for them.²³ So it seems possible that they have, in fact, observed gas jets. On the other hand is the gas jet in the shooting through case caused by the merged section of the interface acting as a nozzle. Hence, if the bubbles do not merge at the beginning of the collapse, e.g., because of a slightly larger distance between them, another collapse mechanism leading to a liquid jet might be triggered. Note that Hopfes et al. calculated the initial distance of the bubbles from a two-dimensional image.²³ Therefore, the value is only accurate if the centers of the two bubbles are in the same plane.

The absence of bubble splitting and formation of satellite bubbles might be explained by an insufficient numerical resolution. Another explanation could be the used material parameters, as the gas jet is too weak to rip-off a part of the bubble into a satellite one. These material parameters are approximations obtained from the literature and might differ more strongly from the experiments than intended.

We have simulated three-dimensional shock-induced collapse dynamics of bubble pairs following the experimental setup of Hopfes et al.²³ using high-order numerical methods. We found that the main driver of their motion is a combination of the pressure difference between the near and far sides of the bubbles and their phase difference.

Depending on the combination of these parameters, cf. Fig. 5, we observed different jetting mechanisms. Detailed descriptions of the observed phenomena are given for each case, and overarching trends are summarized. In the following, we reflect the main insight provided by this study.

Contrary to the interpretation of experimental imaging by Hopfes et al.,²³ we observed gas jets rather than liquid jets in the shooting through and reversing collapse case. Furthermore, we introduced two new cases of bubble-pair interaction—intermediate case one and intermediate case two—to demonstrate the superimposition of two or more of the collapse mechanisms found by Hopfes et al. This superimposition leads to a new phenomenon which we call bidirectional jetting as both jets toward and jets away occur. In addition to the axial liquid jets in the collapse toward case, our simulations revealed the existence of diagonal liquid jets. These diagonal jets could arise due to asymmetry in the collapse of the bubbles.

For the majority of the observed jets the theory of Chew et al.²¹ held. The non-jetting weak interaction case also fit their description. However, we observed some jets showing a different behavior than predicted by their taxonomy. These jets appeared during reversing collapse, intermediate case one and intermediate case two.

Surprisingly, all these cases are located well within a distinct region in the γ-$Δθ$-plane according to their classification. Apparently, Chew's taxonomy fits well at the boundaries of the proposed region. In the core region, however, the collapse mechanisms cannot be separated and rather appear to be superimposed. The observed jets differ from Chew's classification. A more thorough parameter study analysis seems justified. Such investigations could also reveal new insight, e.g., finding a criterion that decides between the emergence of a gas jet and a liquid jet in the shooting through case.

Apart from that, we saw further differences in comparison to the experiments of Hopfes et al.,²³ which could be mitigated with a higher numerical resolution. These include missing satellite bubbles in the shooting through and reversing collapse cases and the small bubble of the weak interaction case not splitting. An increased numerical resolution could also reveal additional small-scale details. Furthermore, we observed only minor grid-resolution effects, which do not affect the physical conclusions drawn from our results.⁴² 

Aside from the mentioned tasks, further work will investigate the influences of the material parameters. Applying parameters to model materials, such as hard gelatin, would be of interest for medical applications.^15,43

The authors have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 667483). The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. for funding this project by providing computing time on the CoolMUC-2 Linux-Cluster at Leibniz Supercomputing Centre.

The authors have no conflicts to disclose.

Benedikt Biller: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal). Nils Hoppe: Conceptualization (equal); Software (equal); Supervision (equal); Validation (equal); Writing – review and editing (equal). Stefan Adami: Project administration (equal); Supervision (equal); Writing – review and editing (equal). Nikolaus A Adams: Funding acquisition (equal); Resources (equal); Writing – review and editing (equal).

The data that support the findings of this study are openly available in mediaTUM at https://mediatum.ub.tum.de/1459250?show_id=1663469&id=1459250&change_language=en, Ref. 44.