Gas-encapsulated droplets have recently been promoted as an effective technique for fluid transport. Shock waves are herein proposed as an instant release mechanism for the encapsulated fluid, which subsequently discharges into the surroundings. This release process relies on the intricate bubble dynamics and droplet response to the shock driving, which are discovered through numerical and theoretical investigations. The key factors involved in the process, such as the complex shock pattern, pressure amplification, and the generation of a sheet jet cascade, are characterized. These observations are further supported by analytical models derived to predict the water hammer pressure, sheet jet velocity, and droplet drift.

The use of encapsulated fluid structures^{1} has been of particular interest in the fields of targeted drug delivery,^{2} fluid transport,^{3} chemical microreaction,^{4} and mixing^{5} processes. The interaction of such layered fluid formations with acoustic waves is being intensively studied due to their numerous biomedical applications such as medical ultrasound imaging^{6} and therapy.^{7} A multi-layer encapsulation approach has been broadly adopted in the field of controlled drug delivery where a payload, entrapped by an outer shell layer, is released using ultrasonic excitation. The outer layer of the delivery vehicle typically serves to ensure stability and long circulation times, while the core medium contains the active agent.^{8} Similar to ultrasound, shock waves induced externally by a shock wave generator^{9} or originating from nearby cavitating bubbles,^{10} which can also trigger the unloading of the payload. Gas-encapsulated liquids are three-fluid systems consisting of a droplet-embedded gas bubble immersed in an ambient liquid. When the ratio of the gaseous shell thickness *h* to the bubble diameter $db,0$ is lower than $10\u22123$, the encapsulation structure is termed an antibubble. In the absence of coating, antibubbles possess very short lifetimes^{11} and poor stability against micellar complexes or external pressure fluctuations.^{12} However, recently reported liquids encapsulated within thick-shelled gas bubbles ($h/db,0>10\u22121$)^{13} show promising controllable characteristics in terms of reproducibility, scalability, high stability, and resistance to external shear stresses, which make these structures a strong attractive candidate for liquid transport in various industrial applications.

This Letter investigates the interaction of a planar shock wave with a thick-shelled gas-encapsulated drop, promoting shock waves as an instantaneous, controllable, and efficient alternative mechanism for the droplet release. Here, shock wave lithotripsy is chosen to trigger the droplet release due to the relative simplicity of its modeling. To illuminate this uncharted area, a detailed numerical and theoretical study is conducted aiming at a comprehensive physical understanding of the bubble-drop-shock interplay.

The interaction of a planar shock wave with an air-encapsulated water droplet is simulated using a seven-equation conservative multicomponent compressible flow model that accounts for both viscous and capillary effects. The fluid components are assumed immiscible. An interface-capturing scheme is used, combining the flow model with a shock-capturing finite-volume method. This is implemented in the open-source code ECOGEN,^{14} which has been validated for several compressible multiphase flow configurations,^{15–17} including gas bubble dynamics problems^{17–19} and droplet breakup due to high-speed flows.^{15,20–22} More details on the governing equations and numerical methods are given in the supplementary material. The problem is treated using a two-dimensional, axisymmetric formulation. The liquid phase is modeled by the stiffened-gas equation of state $p\u2113=(\gamma \u2113\u22121)\rho \u2113el\u2212\gamma \u2113\pi \u221e$, where $el$ is the internal energy, $\gamma \u2113=6.12$ is an empirically determined constant, and $\pi \u221e=3.43\xd7108\u2009Pa$ is a constant representing the attraction between water molecules.^{23–25} The air follows the ideal gas law. Thermal effects are neglected. The planar shock wave is initialized inside the domain and travels from left to right. For a given incident shock (IS) Mach number *M*, the initial flow field is determined from the Rankine–Hugoniot jump relations for the stiffened-gas equation of state using a downstream density of $998\u2009\u2009kg/m3$ and a pressure of 1 atm. Throughout the Letter, subscripts 1, 2, b, and d refer to the pre-shocked state, post-shocked state, bubble, and droplet variables, respectively. Pressure, density, sound speed, and velocity parameters are, respectively, denoted by $pi,\u2009\rho i,\u2009ci$, and $ui$, where *i *=* *1 or 2. The gas-encapsulated droplet is initially at rest and located at the center of the bubble. The initial droplet to bubble size ratio is denoted $\epsilon =rd,0/rb,0$, where *r* is the radius. The bubble radius of $rb,0=250\u2009\u2009\mu m$ is kept constant for all simulations. This corresponds to the current order of magnitude reached experimentally.^{13} The excess pressure inside the bubble and the droplet over the ambient pressure are computed according to the Young–Laplace equation. The surface tension, *γ*, between air and water is 72 mN/m. A symmetric boundary condition is applied at the lower boundary of the computational domain, which represents the droplet axis. Non-reflective boundary conditions are enforced on the remaining boundaries. The high spatial and temporal accuracy at a reasonable computational cost are ensured by a three-level adaptive mesh refinement, which is adapted to track all interfaces and flow discontinuities. The simulations used a Cartesian grid and are carried out in a $12rb,0\xd730rb,0$ rectangular computational domain. Following a grid resolution study, 112 points per droplet diameter are used to solve the gas-encapsulated droplet system (see the supplementary material for more details).

At the very early stages, the incident shock (IS) interacts with the proximal bubble wall resulting in a three-wave system consisting of the incident shock, the reflected wave (ISr), and the transmitted shock (TS), as shown in Fig. 1(a). As a result of the water-to-air acoustic impedance ratio being greater than one, the ISr is a rarefaction wave. The proximal bubble wall deforms as the shock propagates across its interface, resulting in a non-spherical bubble collapse. The droplet is not affected by the early shock-bubble interaction. When the transmitted shock reaches the droplet, it is reflected (TSr) and transmitted through the interface. The TSr wavefront then reflects at the proximal bubble wall, thus generating a reflection of the TSr shock (TSrr). A close examination of the TS–TSr wave system over time shows the formation of an irregular Mach reflection consisting of the TS and TSr wavefronts, and a Mach stem (MSt) all interconnected by a triple point as seen in Fig. 1(b). The encapsulated droplet starts to deform as it interacts with the incident post-shock flow, that is, when the proximal bubble wall reaches the front wall of the droplet. As time proceeds, the continuous collapse of the bubble increases the exposure of the droplet surface to the incident flow, and, consequently, the hydrodynamic forces exerted on it. Thus, the droplet continuously flattens and drifts in the streamwise direction. At the bulk–droplet–bubble boundary, the bubble's edge forms a cusp where a supersonic sheet jet starts to develop as a consequence of the strong water-hammer shock (WS) generated by the impact of the upstream bubble wall on the droplet's front face [Fig. 1(b)]. Note that the WS propagating in the surrounding liquid (WSa) is distinguished from the one propagating in the droplet (WSb). The WSb also affects the droplet as it leads to the formation of low-pressure regions [Fig. 1(c)]—likely to trigger cavitation—by experiencing internal reflection at the droplet boundary (WSbr). The supersonic sheet jet induces a detached bow shock (BS) ahead of its tip that adopts a mushroom-like shape. When interacting with the bubble wall, the detached bow shock is both transmitted (TBS) and reflected (BSr).

Out of all these shock waves, the WS is by far the most violent event which can yield extremely high pressures. Classically, water-hammer experiments rely on the impact of a liquid droplet on a solid substrate.^{26} In the present configuration, however, the surrounding fluid in the post-shock state plays the role of the solid substrate, and the impact velocity, herein denoted $r\u0307b|\alpha =0$, corresponds to the proximal bubble wall velocity at impact time, $ti$, for the polar angle displayed in Fig. 1(a), *α* = 0. Combining the post-shock flow conditions derived from the Rankine–Hugoniot jump relations with the analytical expression for the water-hammer pressure inside the unshocked water droplet $pwh=\rho dcdr\u0307b|\alpha =0\rho 2c2/(\rho 2c2+\rho dcd)$ yields

Assuming a spherical bubble collapse, the bubble speed can be derived from the Rayleigh collapse equation^{27} as $r\u0307b|\alpha =0=[2(p2\u2212pb,0)(\epsilon \u22123\u22121)/(3\rho 2)]1/2$, where $pb,0$ is the initial pressure inside the bubble. Combining it with Eq. (1) and expressing *ρ*_{2} and *p*_{2} for a given *M* through the Rankine–Hugoniot jump relations result in an analytical expression for the water-hammer pressure inside the droplet. Figure 2(a) compares the analytically computed water-hammer pressure with that extracted from the numerical simulation for *ε* of 0.3 and 0.8. While agreement on the order of magnitude is found, the results reveal significant discrepancies that increase for decreasing *ε*. This can be explained by comparing numerical measurements of the bubble-to-drop wall impact velocities with predictions from the Rayleigh collapse equation that exhibits deviations, as shown in Fig. 2(b). Indeed, the upstream part of the bubble starts to deform upon impact, resulting in a decrease in sphericity with $rd,0$ and with increasing *p*_{2}. This confirms that the spherical collapse assumption overestimates the bubble wall velocity because of the faster dynamics of such an event.

To account for the non-spherical nature of the bubble collapse, $r\u0307b|\alpha =0$ can be approximated by an axisymmetric segmented formulation of the Keller–Miksis model (KMEs),^{28,29} where the initially spherical bubble is divided into *n* segments of radius $r|\alpha $ and angle $\alpha \u2208[0,\pi /(n\u22121),\u2026,\pi ]$. The radial evolution for each of them is computed by numerically solving the Keller–Miksis equation (KME). Although the KME assumes spherical symmetry of the bubble, such segmentation can conveniently approximate the spatial pressure differences along the bubble interface caused by the shock propagation while assuming purely radial flow.^{29} The segments are coupled by the total gas volume of the bubble $Vb=Vr+(2\pi \delta \alpha /3)r3|\alpha \u2009sin\u2009(\alpha )$, where $Vr$ is the volume of the remaining segments and $\delta \alpha $ is the angle discretization step. Here, the canonical KME is modified to account for the encapsulated droplet by modeling the polytropic process of the gas volume as $\xi =(Vb,0\kappa \u2212Vd,0\kappa )/(Vb\kappa \u2212Vd,0\kappa )$, with the polytropic exponent $\kappa =1.33$. The initial volumes of the bubble and the droplet, respectively, $Vb,0$ and $Vd,0$, are taken for perfect spheres from their respective radii. The problem is initiated by setting the radius and radial velocity of each segment to $r0|\alpha =rb,0$ and $r\u03070|\alpha =0$, respectively. The speed of sound and the density of the liquid are set to the pre-shocked conditions, $c1$ and $\rho 1$, and the over pressure, $Ps$, is set to 0. The shock wave front propagation is modeled by applying the post-shock conditions to a segment located at *α* once the elapsed time from the contact of the shock with the bubble is greater or equal to $rb,0[1\u2212cos\u2009(\alpha )]/(Mc1)$. Beyond this instant, the speed of sound and density, for a specific segment, are set to post-shocked conditions, *c*_{2} and *ρ*_{2}. The shock wave pressure, $Ps$, is set to the shock pressure, and the initial velocity of this segment, describing the local bubble interfacial speed, is set as the velocity of the rarefaction wave ISr^{30} given by $r\u0307|\alpha =2u2Z2/(Z2+Zb)$,^{31} where $Zi=\rho ici$ is the characteristic specific acoustic impedance. Note that $Z2\u226bZb$, so that the initial velocity reduces to $r\u0307|\alpha \u22482u2$, until the wall speed exceeds it.

At the end, the KMEs equation reads

where $r\xa8|\alpha $ is the radial acceleration and *H* is the specific enthalpy, given by

The vapor pressure $pv$ and viscosity *μ* are taken for water at 20°. The ambient pressure *p*_{1} is set to 1 atm. Figure 3 displays a direct comparison between the KMEs model and a simulation computed with ECOGEN, which reports very good agreement on the bubble collapse prediction. Note that the bubble simulated by the spherical Rayleigh–Plesset equation (RPE) is already in the rebound phase and expanding. This is due to the spherical assumption, which overestimates the bubble wall velocity, as non-spherical collapses are weaker.^{30,32}

Determining the bubble interface speed $r\u0307b|\alpha =0$ in Eq. (1) using the KMEs model yields a semi-analytical expression for the water-hammer pressure inside the droplet. Numerical artifacts, linked to the segmented formulation, arise for segments located between ($7\pi /4,\pi /4$) and are slightly visible in Fig. 3(right), resulting in an over-estimation of the speed. Note that these artifacts are likely to increase for smaller droplets and weaker shock waves. To minimize this spurious numerically induced velocity fluctuation in the analysis, the impact velocity at *α* = 0 is instead approximated as the average of the speed of every segment in the above-mentioned range [solid lines in Fig. 2(b)]. The shaded regions display $\xb12\sigma $, where *σ* is the standard deviation associated with the above-mentioned averages. Despite remaining discrepancies as *ε* increases, both water-hammer pressure and impact velocity computed with KMEs show a much better agreement with the values extracted from the numerical simulations compared to the analytical model, as shown in Figs. 2(a) and 2(b). As expected from their dependency on bubble sphericity, differences between the analytical and the semi-analytical water-hammer pressure become consistently larger with *M* and decreasing *ε*. In the range of *ε* and *M* herein investigated, the droplet influence on the bubble dynamics, when exposed to a sustained shock wave, is negligible.

The high pressure resulting from the water-hammer shock acts on the bubble cusp, driving the surrounding liquid into the gas bubble and forming a sheet jet [Fig. 1(b)]. This process is similar to sheet jets observed in classical shock-bubble interactions, where the driving water hammer shock is produced by two opposite sides of the bubble interface coming into contact.^{33} At the very early stages of the sheet jet development, $ujet\u2248ucusp\u2248uWSa$, where $ujet$ is the jet speed, $ucusp$ is the bubble cusp velocity, and $uWSa$ is the WSa wavefront speed. Noting that $Zb/Z2\u226a1$ and $Zd/Z2\u22481,\u2009ujet$ can be evaluated using the approximate linear dependency $uWSa=c2+2up,WSa$,^{34} where $up,WSa$ is the WSa particle velocity approximately equal to $2u2$.^{31} Figure 2(c) compares the theoretically predicted jet velocity against the average values extracted from the numerical simulations. Their good agreement confirms the independence of the jet speed from the drop-to-bubble size ratio *ε* and implies that the jet dynamics is mainly driven by the liquid flow parameters behind the WSa. When the sheet jet reaches the bubble wall and disrupts it, an additional water-hammer shock is generated. Similar to the first, this second water-hammer shock produces yet another sheet jet that propagates inside the bubble in the direction of WS_{2} before deviating because of the pressure gradient induced by the initial shock wave (IS) [Fig. 1(c), fourth frame]. This results in a curved jet. Once again, when the second jet reaches the bubble wall, a third water-hammer shock is initiated, resulting in a third sheet jet [Fig. 1(c), sixth frame]. This recurrent jet formation process is referred to as the jet cascade. As seen in Fig. 1(c), the jet cascade overcomes the bubble collapse process, leading to bubble fragmentation into several toroidal bubbles. As the jet cascade propagates, the jets converge toward the bubble axis where they eventually interact and generate a very high pressure region at the rear of the bubble. For $\epsilon =0.4$, such pressure amplification yields values of 6.6 and 26.7 GPa for *M *=* *1.1 and *M *=* *2.1, respectively.

After the shock wave has impacted the proximal bubble wall, the droplet becomes exposed to the surrounding liquid. Consequently, the droplet's center-of-mass drifts along the shock's direction of propagation. Considering the droplet as a non-deformable spherical particle of constant mass and assuming that both pressure forces from the incident post-shock flow and bubble are applied to the cross-sectional droplet area, the acceleration of the droplet is trivially derived from Newton's second law and reads $ad=3(p2\u2212pb)/(4\rho drd,0)$. Now, assuming a constant acceleration of the droplet, the displacement of its center-of-mass, normalized to its initial radius, reads $\Delta x/rd,0=(t/2tr)2$, where $tr$ is the characteristic transport time of a shocked gas-encapsulated droplet equal to $rd,0\rho d/(p2\u2212pb)$, which can be derived from the Rankine–Hugoniot relations as

Figures 4(a) and 4(b), respectively, display the dimensional and dimensionless displacement of the droplet's center-of-mass in its axial direction over the normalized time $t/tr$ and compare the analytical model with numerical results from simulations. The numerical data for every shock wave Mach number are found to collapse on the same curve and to be in excellent agreement with the analytical model.

In addition to internal forces, the droplet is affected by the propagating water-hammer shock. Because of the large water-to-air acoustic impedance ratio, the droplet interface behaves as a perfect mirror, virtually trapping the WSb energy inside the droplet.^{22} The confined WSb, thus, experiences nearly total reflection (WSbr) at the droplet interface, yielding the focusing of expansion waves. The droplet is, therefore, locally exposed to a tensile force likely to rupture the liquid and generate cavitation. According to vapor nucleation theory, cavitation in pure water occurs for pressures below −134 MPa at 300 K.^{35} For heterogeneous nucleation thresholds in impure water, the literature is less specific but reports much lower threshold tensions, ranging from −0.1 to −1 MPa.^{36–39} For *M *=* *1.1 and $\epsilon =0.4$, the pressure inside the droplet locally drops down to −3.04 MPa, that is, by more than twice the pressure threshold for homogeneous cavitation, which reasonably suggests the occurrence of shock-induced cavitation inside the droplet through either homogeneous or heterogeneous nucleation. The regions within the droplet falling below −1 and −134 MPa are illustrated at two successive times [Fig. 1(c), fifth and sixth frames]. Because ECOGEN does not account for phase change, the cavitation process is not further discussed.

While the liquid core deforms, it mixes with the surrounding liquid. In the unlikely absence of cavitation, the droplet dissemination is briefly quantified by measuring the droplet's mean density of distribution in the bulk liquid *D.*^{40} The $\u27e8D0\u27e9/\u27e8D\u27e9$ ratio converges asymptotically to zero [Figs. 4(c) and 4(d)], indicating complete mixing. Increasing *M* for a constant *ε* results in a more efficient mixing as a result of the more intense driving flow field. When the droplet deformation is dominated by the jet cascade-induced high-pressure region located at the rear of the bubble [last two frames of Fig. 1(c)], the exponential shape of the dissemination shows an oscillation, seen at $t/tr\u22483$ in Fig. 4(d). This dampens when the inertial forces overcome the high-pressure region effects on the droplet deformation.

To summarize, these results lay the ground work for the shock-induced release of a gas-encapsulated droplet. The underlying dynamics is found to be dominated by the complex shock patterns within the bubble-drop system and the production of an extremely high water-hammer pressure. This results in an intricate converging jet cascade as well as cavitation in the droplet. Further research should focus on (i) effects of various geometrical configurations such as non-concentric bubble–droplet systems, (ii) identifying the contribution of cavitation in the droplet release dynamics, (iii) quantifying scale effects in shock-induced drug unloading, and (iv) ultrasound as an alternative release mechanism for the encapsulated drops. For such gas-encapsulated droplets to become viable vehicles for drug delivery applications, a technique to produce stable gas-encapsulated drops at the low microscale should also be demonstrated.

See the supplementary material for additional details on the numerical model and methods as well as the animation corresponding to Fig. 1(c).

The authors acknowledge the financial support from ETH Zurich and the ETH Zurich Postdoctoral Fellowship programme.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Luc Biasiori-Poulanges:** Conceptualization (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (lead); Methodology (lead); Software (lead); Validation (equal); Writing - original draft (equal); Writing - review and editing (equal). **Guillaume T. Bokman:** Investigation (equal); Validation (equal); Writing - original draft (equal); Writing - review and editing (equal). **Enea Baumann:** Investigation (equal); Writing - original draft (equal); Writing - review and editing (equal). **Outi Supponen:** Funding acquisition (equal); Investigation (equal); Supervision (lead); Validation (supporting); Writing - original draft (supporting); Writing - review and editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.