2. Experimental and numerical methods

Our experimental set-up consists of an amplified titanium–sapphire femtosecond laser emitting at 800 nm with a pulse duration of 65 fs, pulse energy of 3 mJ and a repetition rate of 1 kHz. A Pockels cell, synchronized with the laser, selects a smaller repetition rate going from a fraction of a Hz to several Hz. The laser beam, of typical radius 2 mm, is then focused with a long working distance microscope objective (50$\times$ N.A. 0.5 (Olympus SLMPLN series)) in a volume of liquid (ultrapure water from a Millipore/Milli Q system) contained in a rectangular fused quartz cell (Hellma 101.150$\phi$s) with a square transparent base of 1 cm$^2$. The pulse energy can be varied using a half-wave plate combined with a polarizing beamsplitter from 2.5 to 30 $\mathrm {\mu }$J. The objective was mounted on a stage for a precise positioning of the focal point with respect to the liquid–air interface. The laser beam waist was estimated to be roughly 1 $\mathrm {\mu }$m. The cell walls were rendered hydrophobic using OTS (octadecyltrisiloxane) in Toluene solution to minimize the effect of the meniscus on the visualization. A high speed camera (Phantom V 640) equipped with a magnifying lens (Navitar) was then used to visualize the dynamics of both bubble and jet formation. Typical frame rates were 110 000 frames s$^{-1}$ with an exposure time of 1 $\mathrm {\mu }$s. The spatial resolution was 3.5 $\mathrm {\mu }$m pix$^{-1}$. The acquisition of the images from the camera was triggered using a train of pulses delivered by an Agilent generator, itself synchronized to a recurring starting pulse powered by the Pockels cell driver.

For the numerical simulations, we use the open-source package OpenFOAM (Jasak, Jemcov & Tukovic Reference Jasak, Jemcov and Tukovic2007), a finite-volume solver for partial differential equations, such as the Navier–Stokes equations. For the time integration we use an Euler scheme with the time step dynamically set to ensure a maximum Courant number less than 0.25. We use compressibleInterFOAM, a solver for two compressible, non-isothermal immiscible fluids, which uses a volume of fluid method to solve for the volume fraction of the liquid phase $\gamma$ with the transport equation (Miller et al. Reference Miller, Jasak, Boger, Paterson and Nedungadi2013; Koch et al. Reference Koch, Lechner, Reuter, Köhler, Mettin and Lauterborn2016; Zeng, Gonzalez-Avila & Ohl Reference Zeng, Gonzalez-Avila and Ohl2020)

(2.1)\begin{align} \frac{\partial \gamma}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\gamma \boldsymbol{u}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\gamma (1 - \gamma) \boldsymbol{U}_{\boldsymbol{r}} &= \gamma (1-\gamma) \left( \frac{1}{\rho_g}\frac{{\rm D} \rho_g}{{\rm D} p_g} - \frac{1}{\rho_l}\frac{{\rm D} \rho_l}{{\rm D} p_l} \right) \frac{{\rm D} p}{{\rm D}t}\nonumber\\ &\quad + \gamma \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}, \end{align}

where $\boldsymbol {u}$ is the velocity field, $\boldsymbol{U}_{\boldsymbol{r}}$ is the relative velocity between the two phases, $t$ is time, $p$ is the pressure, $\rho$ is the density and the subscripts $g$ and $l$ represent the gas and liquid phases, respectively. The final term on the left-hand side is an artificial term that counteracts the numerical diffusion of the interface, keeping it sharp (Rusche Reference Rusche2003; Koch et al. Reference Koch, Lechner, Reuter, Köhler, Mettin and Lauterborn2016). The compressibilities ${\rm D} \rho / {\rm D} p$ are determined by the equation of state for the gas and liquid phases, which for simplicity we choose to be the ideal gas and fluid laws, neglecting higher-order effects of temperature and pressure (Miller et al. Reference Miller, Jasak, Boger, Paterson and Nedungadi2013; Koch et al. Reference Koch, Lechner, Reuter, Köhler, Mettin and Lauterborn2016). Numerical simulations using compressibleInterFOAM have recently been performed to study cavitation and jetting in thin gaps (Zeng et al. Reference Zeng, Gonzalez-Avila and Ohl2020).

Heat and mass transfer are neglected, and the governing equations for the motion of the fluid are then the continuity equation

(2.2)\begin{equation} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} ) = 0, \end{equation}

and the compressible Navier–Stokes equations

(2.3)\begin{equation} \frac{\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \boldsymbol{u}) ={-} \boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \tau + \boldsymbol{F}, \end{equation}

where $\tau$ is the viscous stress tensor $\tau = \mu [ \boldsymbol {\nabla } \boldsymbol {u} + \boldsymbol {\nabla } \boldsymbol {u}^T - \frac {2}{3} (\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}) \boldsymbol{\mathsf{I}}]$, with $\mu$ the dynamic viscosity and $\boldsymbol{\mathsf{I}}$ the identity tensor; $F$ is the force due to surface tension which is modelled using a continuous-surface-force method which preserves the sharpness of the interface in order to accurately estimate the curvature (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992; Koch et al. Reference Koch, Lechner, Reuter, Köhler, Mettin and Lauterborn2016).

Equations (2.1)–(2.3) are solved together in a rectilinear domain. The extent of the three-dimensional (3-D) simulation domain was 400 $\mathrm {\mu }$m in both horizontal directions ($x$, $y$) and ranged from 800 to 1200 $\mathrm {\mu }$m in the vertical direction ($z$), with the liquid–air interface ($z=0$) fixed at 300 $\mathrm {\mu }$m above the domain bottom. We used a uniform resolution of ${{\rm d} x} = 2$ $\mathrm {\mu }$m. The boundary conditions are zero gradient at all domain boundaries. More information on the sensitivity to domain size and resolution can be found in Appendix A.

Simulating the exact experimental situation, including the laser–fluid interaction, is complex. To simplify our approach, we assume that this interaction simply nucleates a small bubble of initial high pressure $P$ beneath the surface at $z = -H$ ($x = y = 0$). We use a spherical bubble with initial radius $R_0 = 10$ $\mathrm {\mu }$m, which remains roughly spherical during the initial dynamics. This is significantly different from experiments where the bubbles observed experimentally are elongated in the direction of the beam (the focal area is a vertically elongated ellipsoid due to the necessarily weak numerical aperture of the long working distance objective). As we show, despite these simplifications, our numerical simulations capture the essential features of the experimental observations.

4. Analysis and discussion

While there exists detailed analysis and theory for the initial jet formation due to the expanding bubble (Blake & Gibson Reference Blake and Gibson1987; Robinson et al. Reference Robinson, Blake, Kodama, Shima and Tomita2001), there appears to be little understanding of the second, thicker jet. Previous work has argued that the mechanism for the second jet differs little from that of the first jet (Patrascioiu et al. Reference Patrascioiu, Fernández-Pradas, Palla-Papavlu, Morenza and Serra2014b): the secondary bubble expansion pushes the fluid away from the free surface into the air. However, in both experiments and simulations, the second jet emerges when the bubble is imploding and recoiling away from the interface, before it has begun its second expansion (as in figure 1 for $t \lesssim 25$ $\mathrm {\mu }$s and figure 5a). If the second jet were the product of the second bubble expansion, then one would expect a delay even longer than the several microseconds between the first bubble expansion and first jet (see figure 1c) since the bubble is even further from the surface when the second jet emerges. The absence of any delay points to the failure of the argument that the second jet arises from the second expansion.

A qualitative test of the hypothesis of the importance of the secondary expansion for the second jet can be made by artificially removing the bubble when it reaches its minimum size in a simulation (see figure 3). When the bubble reaches its minimum size at $t \simeq$ 14 $\mathrm {\mu }$s, after its implosion, we set the phase in the submerged bubble to water, leaving the velocity, pressure and temperature fields as they are. In this way the second expansion is effectively, albeit artificially, removed. A comparison between a simulation with (figure 3a) and without (figure 3b) the second expansion shows only minute differences. The second jet still emerges and engenders a large droplet. If the second expansion is not responsible for the second jet, its source must be sought elsewhere.

Figure 3.

Vertical slices of the simulations at different times for $P = 200$ bars and $H = 40$ $\mathrm {\mu }$m. The liquid (water) phase is dark and the gas (air) phase is light. Until time $t \simeq 14$ $\mathrm {\mu }$s, when the bubble size reaches its first minimum, the two simulations (a) and (b) are the same with dynamics as in figure 1 for $t \lesssim 25$ $\mathrm {\mu }$s. Then for (b), the bubble is artificially removed by setting the air phase below the interface to water and the simulation allowed to evolve. In this way the second expansion of the bubble is removed. A comparison of the second jet with (a) and without (b) the second expansion shows that the second jet is much the same in its absence.

Here, the additional information provided by the simulations becomes essential. In figure 4 we show the velocity and pressure fields of a representative simulation where $P = 200$ bar and $H = 40$ $\mathrm {\mu }$m. These fields are consistent with the explanation that the first jet is created by the initial bubble expansion pushing the fluid upward, the primary flow. The pressure gradient is directed radially away from the bubble, and the velocity is likewise directed radially outward. During the initial formation of the second jet, however, the velocity and pressure fields are significantly different. The pressure is much larger below the surface than at the surface, and the largest gradient occurs when the bubble is at its minimum size ($t=14$ $\mathrm {\mu }$s). The collapse of the bubble, which becomes toroidal, generates a rapidly retracting counter-jet and a substantial inflow into the region between the bubble and the surface ($t = 13$ $\mathrm {\mu }$s). This creates the high-pressure region below the surface as well as a saddle flow ($t = 14$ $\mathrm {\mu }$s). The lateral flow into the region between the bubble and the surface, which we term the secondary flow and which is established before the bubble expansion, appears to be the source of the second jet. Less than 1 $\mathrm {\mu }$s later, the second jet has already emerged, well before the bubble re-expands. In figure 5(a) we show the time of the second jet emergence $T_{jet}$ vs the time that the bubble collapses to its minimum size $T_{min}$. The near perfect correlation between the two times for all pressures, energies and depths reinforces our argument that the second jet is related to the bubble contraction, and the secondary flow it engenders, and not the bubble re-expansion. Identifying the secondary flow as the mechanism behind the second jet formation, as done qualitatively in Chen et al. (Reference Chen, Yu, Su, Chen and Chen2013), permits a simple analysis of the second jet velocity scaling.

Figure 4.

Vertical slices of the simulated flow fields at $t = 1$, 13, 14 and 15 $\mathrm {\mu }$s for $P = 200$ bars and $H = 40$ $\mathrm {\mu }$m, focusing on the region near the interface. (a) Plots of the liquid (dark) and air (light) phases with velocity vectors super-imposed. The velocity vector lengths at $t = 1$ $\mathrm {\mu }$s correspond to the distance a fluid particle will travel in 1 $\mathrm {\mu }$s, while the velocity vector lengths for the remaining times correspond to the distance a fluid particle will travel in 2 $\mathrm {\mu }$s. (b) Colour plots of the pressure field $p$. The phase boundaries are indicated by thick black lines. (c) Colour plots of the horizontal velocity $u_y$. (d) Colour plots of the vertical velocity $u_z$. The first column ($t = 1$ $\mathrm {\mu }$s) shows the bubble several $\mathrm {\mu }$s prior to the emergence of the first jet. The remaining columns are before ($t = 13$ $\mathrm {\mu }$s), during ($t = 14$ $\mathrm {\mu }$s) and after ($t = 15$ $\mathrm {\mu }$s), the emergence of the second jet, which is accompanied by a substantial secondary flow. A substantial secondary flow accompanies the emergence of the second jet, with the flow drawn in by the bubble collapse and counter-jet feeding the second jet and creating the large pressure gradient that drives it.

Figure 5.

(a) Plot of the second jet emergence time $T_{jet}$ vs the time when the submerged bubble reaches its minimum size $T_{min}$. The uncertainties on the experimental and numerical results correspond to the time between recorded times, 9 and 1 $\mathrm {\mu }$s, respectively. (b) Plot of the first $U_1$ and second $U_2$ jet velocity vs $H$ for fixed $P = 400$ bar from simulations. Both vary approximately linearly near a threshold value of $H_c \simeq 90$ $\mathrm {\mu }$m.

Mass conservation determines a relationship between the secondary inflow velocity and the second jet velocity, $U_2$. Although the degree of contribution to the secondary flow from the bubble collapse and counter-jet is unclear, we attribute the large downward flow in figure 4(a) ($t \geq 13$ $\mathrm {\mu }$s) to the counter-jet, which we assume has a velocity comparable to the first jet $U_1$. By considering a cylindrical volume of liquid of diameter $D_c$ and length $l_c$, we can write $2 {\rm \pi}(D_c/2) l_c U \sim {\rm \pi}(D_c/2)^2 U_1$, which yields

(4.1)\begin{equation} U \sim U_1 \frac{D_c}{4 l_c}. \end{equation}

We anticipate that $D_c$ corresponds to the maximum nucleated bubble size $D_{max}$ while the vertical length scale $l_c$ is as yet unknown. (Additional evidence for a single horizontal length scale proportional to $D_{max}$ is given in figure 6(a), where we show that the droplet diameter $D \propto D_{max}$.) In figure 5(b) we plot both $U_1$ and $U_2$ vs $H$ for simulations at $P = 400$ bar. Both velocities decrease approximately linearly with $H$ and approach zero at $H \sim 90$ $\mathrm {\mu }$m, beyond which no jetting is observed. The ratio between $U_2$ and $U_1$ is thus constant and given by $D_{max}/4 l_c$ (4.1), which after measuring $D_{max}$ yields $l_c$. In figure 6(b) we show that $l_c \simeq D_{max}$ for several values of $P$. The bubble dimensions thus set the relation between $U_2$ and $U_1$.

Figure 6.

(a) Plot of the droplet size $D$ vs horizontal maximum bubble diameter $D_{max}^h$ from simulations. Combining the results from several depths indicates that $D \propto D_{max}^h$. (b) Plot of first $U_1$ and second $U$ jet velocity vs initial depth $H$ from simulations for three different initial pressures $P = 200$, 300 and 400 bar. For all pressures, both $U_1$ and $U$ decrease nearly linearly to zero at a critical depth $H^*$. The ratio of the slopes of these two lines is $D_{max}/4l_c$, as argued in the main text. Inset: $D_{max}/l_c$ vs $P$. The ratio is close to unity for all $P$, suggesting that there is only a single dominant length scale.

In figure 7(a,b) we show that the linearity of $U_2$ with $H$ up to a threshold value $H_c$, as in figure 5(b), is generic for both experiments with fixed $E$ and for additional simulations with fixed $P$. Unsurprisingly, the critical depth $H_c$ increases while the slope $\partial U / \partial H$ decreases with an increase of $E$ or $P$. The similarities between these graphs suggest that the experimental and numerical results can be brought into correspondence with a suitable normalization. Since we have argued that the second jet is associated with bubble collapse, we estimate a characteristic bubble collapse velocity as $U_{\sigma } = \sqrt {\sigma /\rho D_{max}}$, where $\sigma$ is the surface tension between water and air, and normalize $U_2$ with $U_\sigma$. To account for the dependence on a critical depth $H_c$, we plot $U_2/U_\sigma$ vs $(H_c - H)/D_{max}$, which is similar to but distinct from the standoff distance $2H/D_{max}$ used in previous work (Pearson et al. Reference Pearson, Cox, Blake and Otto2004) by virtue of the inclusion of $H_c$ necessary for collapse. As figure 7(c) shows, the data collapse satisfactorily. Although the bubbles in the simulations are roughly spherical, the bubbles in experiments are strongly elliptical, so to account for this asymmetry we used the maximum bubble horizontal width $D^h_{max}$ and vertical length $D^v_{max}$ as the relevant length scales. In the simulations these two length scales are nearly the same. If different, as in experiments, we expect the vertical length scale $D^v_{max}$ to affect the critical depth $H_c$, while we expect the horizontal length $D^h_{max}$ to control the bubble collapse speed. Thus $D^v_{max}$ is used to normalize $(H_c - H)$ and $D^h_{max}$ is used to normalize the velocity. With these choices the experiments and simulations collapse, but separately. To bring them into accord with each other, we introduce a parameter $\beta \simeq 0.48$ which is used only for the simulations to artificially modify the length scales $D^v_{max}$ and $D^h_{max}$. Since the value of $\beta$ is close to the ratio of the $D^h_{max}$/$D^v_{max}$ for experiments (see Appendix B), this suggests that we have fully captured the behaviour of both simulations and experiments with the normalization in figure 7(c), excepting a dependence on the bubble asymmetry at its maximum size.

Figure 7.

The second jet velocity $U_2$ vs initial depth $H$ for (a) experiments with several fixed $E$ and (b) simulations with several fixed $P$. All curves reveal the same linearity of $U_2$ with $H$ and a threshold $H_c$. (c) Second jet velocity $U_2$ from experiments and simulations normalized by $U_\sigma = (\sigma /\rho \beta D_{max}^h)^{1/2}$ vs a normalized distance from the interface $(H_c-H)/\beta D_{max}^v$. Here, $\beta = 0.48$ (close to the experimental ratio $D_{max}^h/D_{max}^v \simeq 0.41$) brings the experiments into accord with the simulations, including initially simulations where the initial bubble shape was elliptic (see Appendix D). (d) Phase diagram for jetting using the two normalized control parameters $H/R_0$ and $PR_0/\sigma$. The interface between jetting and no-jetting scales as $H/R_0 \propto (PR_0/\sigma )^{1/2}$. Inset: plot of $\sigma H_c^2$ vs $P R_0^3$ with the same symbols as the main figure. A linear calibration fit to the simulations determines $PR_0^3$ for experiments, with $P$ from (3.1) then yielding $R_0$.

We add further evidence to the correspondence between experiments and simulations with a combined phase diagram in figure 7(d) indicating the parameter ranges for jetting and no jetting. For the simulations we systematically varied $H$ and $P$ to determine when jetting was observed, yielding a jetting to no-jetting boundary empirically found to be $H/R_0 \sim (PR_0/\sigma )^{1/2}$. The scaling of the boundary can be understood by recasting the dependence as $\sigma H_c^{2} \sim PR_0^3$, a relationship born out in the inset of figure 7(d). This simply states that, for jetting to occur, the bubble energy must overcome a surface energy (Pearson et al. Reference Pearson, Cox, Blake and Otto2004; Chen et al. Reference Chen, Yu, Su, Chen and Chen2013). We note that the length scale of the surface energy $H_c$ is itself simply related to the nucleated bubble dimensions at maximum expansion (see figure 10d). The surface in question is thus that which will protrude from the initially flat surface due to a bubble of size $D_{max} \propto H_c$, with a curvature also set by these length scales. To include the experimental data in the phase diagram requires knowledge of $P$ and $R_0$, which we are not able to measure directly. Using a line calibrated to the simulation data in the inset of figure 7(d), we determine the value of $P R_0^3$ corresponding to the value of $\sigma H_c^2$, which we can determine experimentally. Using $P$ estimated with (3.1), the value of $R_0$ turns out to be similar to the beam waist, $\sim$1 $\mathrm {\mu }$m, as found in Schaffer et al. (Reference Schaffer, Nishimura, Glezer, Kim and Mazur2002). The experimental data can then be brought into correspondence with the numerical data in a satisfactory manner as shown in figure 7(d). Taken together, figure 7(c,d) provides a comprehensive understanding of the jet formation using LIFT as well as the correspondence between laser pulse energy and the engendered pressure giving rise to the cavitation bubble.