3.1. Potential flow solution

The potential flow solution is obtained by solving the inviscid part of the linearised momentum equation

(3.2)\begin{equation} \frac{\partial {\boldsymbol u}_{\phi}}{\partial t} ={-}\frac{1}{\rho}\,\boldsymbol{\nabla} p. \end{equation}

The magnitude of the pressure gradient field (or acceleration field) at short infinitesimal times governs the dynamics of bubble collapse at a finite time. To obtain the acceleration, we introduce the dimensionless pressure $\tilde {p}=({p-p_{0,L}})/({p_\infty - p_{0,L}})$. Note that in the linear problem, the interface does not move, therefore surface tension effects at short times only introduce a correction on the scaling prefactor $(p_\infty -p_{0,L})$ upon which the solution depends. The divergence of (3.2) leads to a Laplace equation ($\Delta \tilde{p} = 0$) that can be solved with appropriate boundary conditions for $\tilde {p}$. In particular, we impose $\tilde {p}=1$ far away from the bubble, $\tilde {p}=0$ at the bubble interface, and ${{\rm d}\tilde {p}}/{{\rm d}n} = 0$ at the solid boundary. The solution of this equation depends only on the geometry of the bubble, which is fully described as a function of the contact angle $\alpha$. This simplified representation of the fluid fields at short times is referred as the ‘free surface model’ in this paper.

In figure 2(a), we start presenting the isocontours of the magnitude of the non-dimensional acceleration ($|\boldsymbol {\nabla } \tilde {p}|$) field obtained for a spherical cap bubble with contact angle $\alpha = {3}{\rm \pi} /{4}$, while figure 2(b) shows the variation of the same along the wall in the liquid phase with respect to the distance from the contact line obtained from the free surface model. The isocontours in the left half are obtained numerically from the free surface model, while those in the right half are obtained from the DNS of the Navier–Stokes equations (using the set-up shown in figure 4a), accounting for the presence of a gas with non-zero density ($\rho _l/\rho _g = 10^{-3}$) and a finite value of viscosity ($\mu _l/\mu _g = 10^{-2}$) at ${t U_c}/{R_c} = 2.07 \times 10^{-6}$. Good agreement between the two models supports the argument that (2.2) is a good representation of the fluid dynamics at very short times. We distinguish clearly two separate regions: the far field where the pressure gradient contours are hemispherical caps centred at the axis of symmetry, and the near field where the contours for pressure gradient are intricate and diverge on approaching the triple contact point. Next, we characterise each of these regions.

Figure 2.

(a) Isocontours of the magnitude of acceleration for a bubble with $\alpha = 3{\rm \pi} /4$. The isocontours in the left half are obtained from the free surface model, and those in the right half are obtained from DNS. (b) Non-dimensional acceleration magnitude along the wall obtained for a bubble with $\alpha =3{\rm \pi} /4$ using the free surface model.

3.1.1. Near field

The numerical solution of the interface acceleration at the initial time is shown in figure 3(a), where we use the angle with respect to the solid wall $\theta _w$ to parametrise the interface position. The origin $\theta _w = 0$ corresponds to the point of contact between the interface and the solid wall, and $\theta _w={\rm \pi} /2$ is the intersection of the interface with the axis of symmetry. Figure 3(a) depicts the local non-dimensional acceleration magnitude at bubble interface $\vert {\boldsymbol a} \vert /a_0$ for four representative cases with $\alpha ={\rm \pi} /3, {\rm \pi}/2, 2{\rm \pi} /3,3{\rm \pi} /4$. The value of the reference acceleration is chosen as $a_0 = ({p_\infty - p_{0,L}})/{\rho R_c}$, which is the acceleration magnitude for a spherical bubble with radius $R_c$ and a known pressure difference. As expected, in the case $\alpha = {\rm \pi}/2$, the numerical solution recovers the Rayleigh–Plesset solution, and the non-dimensional acceleration is uniform and equal to 1 all over the interface. When $\alpha < {\rm \pi}/2$, the interface acceleration (thus velocity) tends to zero at the contact point, and therefore it ceases to move even in the potential flow problem with a slip wall. For $\alpha > {\rm \pi}/2$, the appearance of the singularity is evident as the acceleration diverges as we approach the $\theta _w\to 0$ limit. To interpret this result, one must keep in mind that the general solution of the Laplace equation in spherical coordinates can be obtained by separation of variables that leads to an infinite series. This series sometimes diverges at the edge where Dirichlet and Neumann boundary conditions meet due to the appearance of a singularity (Dauge Reference Dauge2006). This singularity was reviewed extensively in the past (Landau & Lifshitz Reference Landau and Lifshitz1959; Blum & Dobrowolski Reference Blum and Dobrowolski1982; Steger Reference Steger1990; Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten1997; Li & Lu Reference Li and Lu2000), but studies where this singularity appears in bubble dynamics problems are rare. The asymptotic solutions close to the point where homogeneous Dirichlet and Neumann boundary conditions meet can be expressed alternatively by using the general solution of the Laplace equation (Li & Lu Reference Li and Lu2000; Dauge Reference Dauge2006; Yosibash et al. Reference Yosibash, Shannon, Dauge and Costabel2011), and take the form

(3.3)\begin{equation} \tilde{p}_s= \sum_k^\infty C_k \tilde{r}_s^{b_k} \cos\left(b_k \theta_s \right), \end{equation}

where $\tilde {r}_s=r_s/R_c$ is the non-dimensional distance from the contact line, and $b_k=({{\rm \pi} }/{\alpha }) (k+1/2)$. Taking the derivative with respect to the normal of the interface, we readily find the interface acceleration magnitude near the contact line as

(3.4)\begin{equation} \vert {\boldsymbol a}_I \vert = \frac{1}{\rho} \left| \frac{\partial p}{\partial n} \right| \approx a_0 C_0 b_0 \tilde{r}_s^{b_0 - 1}. \end{equation}

This expression exhibits a singularity at the contact point ($r_s \to 0$) when $b_0 = {{\rm \pi} }/{2\alpha } < 1$ (or $\alpha > {\rm \pi}/2$), implying that the acceleration at the triple contact point is infinite. In these conditions, the first term in the expansion is the leading-order term that eventually dominates the solution in the region $r_s \le R_c$. This can also be seen clearly in figure 2(b), where the fitting curve shown with a solid line is obtained using the first term of the series solution (i.e. parameters $C_0$ and $b_0$ in (3.4)). We repeat the fitting procedure for various $\alpha$ to verify that the numerical values of $b_0$ match well with theoretical predictions, and find $C_0$ that is a constant of order 1 that increases slightly with $\alpha$ (figure 3b).

Figure 3.

(a) Non-dimensional acceleration magnitude ${\vert {\boldsymbol a}_I\vert }/a_0$ along the interface parametrised using the angle $\theta _w$ (measured in the counter-clockwise direction from the point of contact of wall and the axis of symmetry) for $\alpha = {{\rm \pi} }/{3},{{\rm \pi} }/{2},{2}{\rm \pi} /{3},{3}{\rm \pi} /{4}$. Dots represent the numerical solution from the free surface model, and the thick lines are predictions using (3.4) evaluated at the interface. (b) Exponent $b_0$ and coefficient $C_0$ obtained from fitting the pressure gradient obtained from the free surface model along the wall for different $\alpha$. (c) Non-dimensional averaged interface acceleration magnitude as a function of the contact angle using different methods: analytical expression given by (3.6) (blue line), the solution from the free surface model (yellow line), the DNS solution of the Euler equations (black crosses), and the DNS solution of the Navier–Stokes (NS) solver for ${Re}=100$ (red crosses).

Focusing on the characterisation of the regimes where the singularity appears, figure 3(a) shows clearly that the numerical solution (dots) is well described by (3.4) evaluated at the interface (solid lines), showing excellent agreement close to the contact point, i.e. small values of $\theta _w$. Near the axis of symmetry (e.g. $\theta _w\to {\rm \pi}/2$ and $\tilde {r}_s \approx 1$), the errors are visible and the first term in the series does not suffice to describe accurately the acceleration field.

3.1.2. Far field solution

The far field flow created by the bubble corresponds to a punctual sink sitting at the intersection between the wall and the axis of symmetry. The integration of the momentum equation in the radial coordinate provides the magnitude of the pressure gradient generated by a punctual sink at an arbitrary distance from the sink as

(3.5)\begin{equation} \left| \frac{{\rm d}\tilde{p}_{{far}}}{{\rm d}\tilde{r}_w} \right| = \frac{\vert \bar{\boldsymbol a}_I \vert}{a_0}\,\frac{1 + \cos(\alpha)}{\tilde{r}^2_w}, \end{equation}

where $\vert \bar {\boldsymbol a}_I \vert$ is the magnitude of averaged acceleration along the interface, which determines the strength of the punctual sink. Figure 2(b) shows clearly that this equation captures well the decay of the pressure gradient at distance $\tilde {r}_w$ far from the interface. Because the first term of the series, i.e. (3.4), predicts the interface acceleration reasonably well, we obtain an estimation of the averaged acceleration magnitude for bubbles with $\alpha > {\rm \pi}/2$ as

(3.6)\begin{equation} \vert \bar{\boldsymbol a}_I \vert = \left| \frac{1}{S_b} \int \frac{-1}{\rho_l}\,\frac{\partial p}{\partial n}\,{\rm d}S_b \right| = a_0 C_0 b_0\,G(\alpha), \end{equation}

where $S_b$ stands for the bubble interface surface, and $G(\alpha )$ is a geometrical factor that can be computed numerically, assuming that the first term in the series is indeed the leading-order term along the entire bubble interface. Figure 3(c) shows the averaged non-dimensional interface acceleration magnitude obtained using this model, which compares well with the numerical results of the free surface model and the solution of both the viscous and inviscid Navier–Stokes equations. The results of the simplified model are obtained using the value of $C_0$ computed numerically and reported in figure 3(b). For $\alpha > {\rm \pi}/2$, $\vert \bar {\boldsymbol a} \vert$ increases with $\alpha$. The three proposed models agree well for the angles tested, implying that the free surface model as well as the simplified expression given in (3.6) capture well the averaged bubble response at short times for sufficiently large Reynolds numbers.

3.2. Viscous correction

The potential flow solution is formally exact only at the initial time $t = 0$ when the interface is at rest. As soon as the interface is set in motion, a thin boundary layer develops immediately, regularising the flow close to the contact point. We use the results from the DNS of the Navier–Stokes equations at short times to investigate the influence of viscosity on the dynamic response of the bubbles.

We solve for the full set of (2.1a) where, for the sake of simplicity, we neglect surface tension effects. The equations are solved using the All-Mach implementation (Popinet Reference Popinet2015; Fuster & Popinet Reference Fuster and Popinet2018), where the one-fluid approach is used to solve the modelling equations. This solver has been used successfully already to investigate some aspects of the interaction of an oscillating bubble and a rigid tube by Fan, Li & Fuster (Reference Fan, Li and Fuster2020), and interactions of collapsing bubbles with a free surface by Saade et al. (Reference Saade, Jalaal, Prosperetti and Lohse2021). The set-up for numerical simulations is shown in figure 4(a), where a spherical-cap-shaped bubble is initialised in a square domain of size 100 times the initial radius of curvature of this bubble. We use such large domains so that the pressure waves reflected from the boundaries do not affect the solution; moreover, the grid is coarsened away from the bubble to diffuse these waves. The simulations are axisymmetric about the left-hand boundary, and the bubble is initially in contact with the bottom boundary where a non-slip boundary condition is applied. The far-away boundaries are given constant far-field pressure conditions and a homogeneous Neumann condition for velocity. The grid is refined progressively near the interface to a minimum grid size $\Delta x_{min}/R_c = 0.00038$. The maximum time step in the simulation is governed by an acoustic Courant–Friedrichs–Lewy restriction, which in our case is 0.4. The initial pressure field inside the domain is given by (2.2), the velocity field is given as zero, and the total energy is given by the equation of state.

Figure 4.

(a) The numerical setup for DNS. (b) Results from the DNS for $\alpha = 2 {\rm \pi}/3$ and ${Re} = 100$: time-averaged interface acceleration in the direction parallel to the wall, $a_{r,I} = u_{r,I}/t$, as a function of the distance from the wall at five different times (in colour where $t U_c/R_{c,0}\in [7.27 \times 10^{-5},3.6 \times 10^{-4}]$). For reference, we include the potential flow solution given by (3.4) evaluated at the interface (solid black line). The inset represents a zoom in to the viscous boundary layer generated close to the wall.

In figure 4(b), we characterise the bubble motion at short non-dimensional time $t U_c/R_{c}$ by showing the evolution of the time-averaged interface acceleration magnitude parallel to the wall, defined as $a_{r}(t) = u_r(t)/t$, as a function of the normal distance from the wall. The results shown with coloured dots are obtained from DNS for $\alpha = 2{\rm \pi} /3$ and ${Re} = 100$, where the colour scale represents the non-dimensional time varying between $7.27 \times 10^{-5}$ and $3.6\times 10^{-4}$. The inset shows the clear development of a boundary layer very close to the wall; the thickness of this layer is defined by the height for which the velocity is maximum. Outside this region, the free surface potential flow solution (3.4), shown with a solid black line, predicts the interface velocity obtained from DNS relatively well, therefore the potential flow solution is a good approximation of the collapse process outside the viscous boundary layer. Inside the viscous boundary layer, the interface acceleration is sensitive to the slip length imposed (see Appendix A) and also to the movement of the contact line. In this region, the solution of the free surface model at $t=0$ cannot be extrapolated to predict the flow field near the contact line at short times. Remarkably, although the contact line motion is grid-dependent, the average acceleration and the maximum jet velocity are shown to be independent of the slip length imposed (see Appendix A). This fact, together with the slight dependence of the average acceleration on ${Re}$ (figure 3c), confirms that for large enough Reynolds numbers, the bubble dynamic response is governed mainly by the potential flow solution and not by the particular contact line motion model used. In these conditions, the free surface potential flow model at $t=0$ is a useful tool to predict the average dynamics of the interface at short times.

In order to characterise the jet formation process, we extract the maximum interface velocity magnitude from the DNS data, which we name the jet velocity ${u}_{jet}$. This velocity is fitted using a power-law function of the form (figure 5a)

(3.7)\begin{equation} \frac{u_{jet}}{U_c} = A \left( \frac{t U_c}{R_c} \right)^q, \end{equation}

where $q$ is a coefficient that becomes smaller than 1 when the acceleration is singular at $t=0$. A theoretical estimation of coefficient $q$ accounting for viscous effects can be obtained from the singular solution described in the previous subsection as follows. The jet velocity after a short time $t_s$ can be approximated from the initial acceleration evaluated at the interface at height equal to the thickness of the boundary layer $z=\delta (t_s)$. The $\delta$ is given approximately by the solution of the Stokes problem for the flow near the flat plate that is started impulsively from rest, therefore we can impose that $\delta (t_s) \approx C_\delta \sqrt {\nu t_s}$, where $C_\delta$ is a constant of order unity. In this case, the jet velocity velocity is estimated as

(3.8)\begin{equation} u_{jet} (t_s) = a_{r,I} \left(t=0, z_I=\delta(t_s)\right) t_s, \end{equation}

where we assume that the interface decelerates quickly as soon as it enters inside the viscous boundary layer. Under these assumptions, the coefficient $q$ is obtained readily as

(3.9)\begin{equation} q=\frac{1}{2} + \frac{\rm \pi}{4\alpha}. \end{equation}

Figure 5.

(a) Evolution of jet velocity for different contact angles from DNS (dots) and fitting using (3.7) (solid lines). (b) Exponent $q$ from DNS fitting and predictions of (3.9).

As shown in figure 5(b), this model is capable of capturing the evolution of the jet velocity at short times found from DNS simulations accounting for viscous effects when $\alpha > {\rm \pi}/2$. The exponent $q$ (figure 5b) decreases as $\alpha$ increases because the acceleration gradient in the normal direction increases with $\alpha$, as seen in figure 3(a). Note that the DNS results confirm that $q < 1$ for $\alpha > {\rm \pi}/2$, and therefore the jet acceleration is singular at $t=0$ even in the presence of viscous effects.

Another interesting quantity to describe the overall interface dynamics is the vorticity generated during the collapse. At short times, vorticity is concentrated in a sheet vortex at the wall and at the interface. By assuming a zero stress condition at the interface for $t = 0$, the only non-zero component of the vorticity can be written as

(3.10)\begin{equation} \omega_\phi ={-} 2\,\frac{\partial {u}_r}{\partial \theta}. \end{equation}

It follows from the potential flow solution (figure 3a) that the sign of the acceleration gradient changes depending upon initial contact angle, i.e. if it is greater or less than ${\rm \pi} /2$, reversing the sign of the vortex sheet generated at the interface. This behaviour is indeed observed from the numerical simulations as shown in figures 6(a,b), where the strength of the vortex sheet is plotted with saturated colour maps. The change in dominant colour represents the opposite sign of vorticity at the interface between the two cases.

Figure 6.

The vorticity is shown in colour map and velocity vectors obtained from DNS after a short time (${t U_c}/{R_c} = 0.1$) for ${Re} = 100$, for (a) $\alpha = 2 {\rm \pi}/3$, and (b) $\alpha = 5 {\rm \pi}/12$.

4.1. High-Reynolds-number regime

Having described the short-time dynamics, in this section we investigate the long-time consequences of this singular response using the experiments and DNS. The results for two representative cases $\alpha > {\rm \pi}/2$ and $\alpha < {\rm \pi}/2$ shown in figures 7(a,b) are discussed now.

Figure 7.

Snapshots of bubble shapes are shown for the two representative cases; the numerical bubble shapes, shown with red curves, are overlaid at same times and scaled to the same length as in the experiments. (a) Case where $\alpha > {\rm \pi}/2$; snapshots are taken every $0.1$ ms. (b) Case where $\alpha < {\rm \pi}/2$; snapshots are taken every $0.125$ ms.

4.1.1. Dynamics of bubbles with an effective contact angle larger than ${\rm \pi} /2$

First, we discuss the case where the singularity is present (e.g. $\alpha > {\rm \pi}/2$). The generation of bubbles with $\alpha > {\rm \pi}/2$ can be investigated using the set-up described briefly in Appendix B, and Bourguille et al. (Reference Bourguille, Bergamasco, Tahan, Fuster and Arrigoni2017) and Tahan et al. (Reference Tahan, Arrigoni, Bidaud, Videau and Thévenet2020). The laser is focused at the bottom of a liquid tank where a water drop is attached. The cavitation in the water drop induces a shock wave in the aluminium plate that leads to the appearance of multiple bubbles in the water tank that interact to form a flat bubble with shape similar to a spherical cap at maximum volume (figure 7a). The snapshots are taken every $0.1$ ms; the interface from the second snapshot is extracted and fitted with an approximate spherical cap which gives $R_c = 7.56\times 10^{-3}$ m and $\alpha = 0.727 {\rm \pi}$.

We reproduce the experimental conditions numerically using the numerical set-up described in § 3.2. The minimum mesh size is $\Delta x = 60\,\mathrm {\mu }{\rm m}$, the far-field pressure $p_\infty$ is 1 atm, and the pressure inside the bubble is set to the low value 0.1 atm, which is chosen by trial and error to match the numerical collapse time with the collapse time in experiments. We consider both surface tension and viscous effects in the numerical simulations, with ${Re} \approx 75\,600$ and ${We} \approx 10\,500$. The numerically obtained bubble shapes plotted with red contours, and scaled to bubble size in snapshot 2, are subsequently overlaid on the experimental snapshots after the same times. Very good agreement is seen between the numerical and experimental bubble shapes; small differences persist because of the simplification of the bubble shape during the expansion, and the influence of gravity, mass transfer and other effects that are not considered in the numerical simulations.

At the beginning of the collapse phase, the highest interface velocity is developed at the edge of the viscous boundary layer (see figure 4b), which leads to the appearance of an annular jet that is generated parallel to the wall, and a mushroom-like shape of the interface contour (see figure 7(a) and supplementary movie 1). These results are consistent with previous numerical works of Lauer et al. (Reference Lauer, Hu, Hickel and Adams2012) and Koukouvinis et al. (Reference Koukouvinis, Gavaises, Georgoulas and Marengo2016). Remarkably, when the collapse is strong enough and the jet reaches the axis of symmetry, a stagnation point appears there, and a secondary upward jet normal to the wall is generated. The re-entrant jet observed for $\alpha > {\rm \pi}/2$ is not conventional in cavitation and generates vortex ring structures similar to those observed by Reuter et al. (Reference Reuter, Gonzalez-Avila, Mettin and Ohl2017) that are persistent in nature and can travel large distances in comparison to bubble size. The generation of this vortex ring is illustrated numerically in figure 8(a), where the colour maps show the vorticity field. A clear re-entrant annular jet is observed, followed by a mushroom like structure, and the bubble detaches from the wall; consequently a jet in the upwards direction leads eventually to formation of a vortex ring. In the experiments, we visualise this by adding dye in the bottom of the tank (figure 8b) during the collapse of a flat bubble. This vortex ring can travel to a long distance and induce unusual long-range effects like free surface waves and jetting (see supplementary movie 3).

Figure 8.

(a) Bubble interface in black curve and the vorticity field in the colour map as obtained from DNS, shown at (consecutively from left to right) ${t U_c}/{R_c}= 0,0.19,0.3,0.81,1.00,1.45,1.95,4.90$. The annular jet parallel to the wall reaches the axis, resulting in the jet being directed away from the wall, which later generates the vortex ring. The results are plotted for $\alpha = 2{\rm \pi} /3$, ${Re} = \infty$ and ${We} = \infty$. (b) Visualisation of the liquid flow field obtained by adding dye at the bottom of the tank, resulting from collapse of the bubble, shown at every 2.5 ms.

4.2. Finite-Reynolds-number effects

Popinet & Zaleski (Reference Popinet and Zaleski2002) showed that for small bubbles collapsing in the vicinity of a rigid wall, the jets can be suppressed due to the viscous effects. In this subsection, we use the DNS to predict quantitatively the range of size for which the observations reported in this paper apply. We take an ideal case of an air bubble with $p_{b,0} = 0.1$ atm and constant contact angle equal to $120^{\circ }$ collapsing in water at an ambient far-field pressure. In this situation, the characteristic velocity of the process stays constant and equal to $10\,{\rm m}\,{\rm s}^{-1}$, and the initial radius of curvature is then left as the unique parameter controlling the values of the Reynolds and Weber numbers.

In figure 9(a), we show the maximum non-dimensional velocity reached during the collapse as a function of the Reynolds number. For Reynolds numbers above a given critical value, ${Re} > {Re}_c= 100$, the influence of the Reynolds number on the peak velocities is only marginal. However, when the Reynolds number is below this critical value, the jet velocity drops dramatically further, decaying with decreasing ${Re}$. This sudden change in the peak velocities is controlled by the appearance of jetting, which is not visible for ${Re} \leq 100$ (see figure 9b). For a low-pressure air bubble collapsing in water at atmospheric pressure, these results reveal that the observations reported (reversed re-entrant jet) are applicable for bubbles larger than $R_c > 10\,\mathrm {\mu }{\rm m}$. We must still emphasise that the results reported here are limited to spherical cap bubbles. In reality, other kinds of behaviour can be observed depending on the mechanism by which the bubbles are generated and the bubble shape reached at the instant of maximum expansion. If the bubble shape is not a spherical cap, or the microlayer is not able to drain during the bubble collapse, or the bubble size is comparable to pits, then it can lead to a variety of other collapse and jetting dynamics discussed by Lechner et al. (Reference Lechner, Lauterborn, Koch and Mettin2020), Reuter & Ohl (Reference Reuter and Ohl2021) and Trummler et al. (Reference Trummler, Bryngelson, Schmidmayer, Schmidt, Colonius and Adams2020), respectively.

Figure 9.

(a) Peak of the non-dimensional velocity as a function of Reynolds number for an air bubble at $p_0=0.1$ atm collapsing in water at atmospheric pressure. (b) Interface contours as a function of the Reynolds number at the instant of minimum radius.