3.1. Evolution of a circular bubble

To establish a theoretical understanding of the bubble evolution, we will also assume a circular bubble, an assumption that may be partly justified by the a posteriori observation that the amplitude of the branching structures is significantly smaller than the scale over which the pressure field varies. Before the bubble is formed, mass conservation of the liquid implies that the added volume under the ball as it is rising must be balanced by a corresponding influx. We may write the volume change within a radius $r$ as

(3.1)\begin{equation} \dot{V} (r)= {\rm \pi}r^2 \dot{h} = 2{\rm \pi} r (h+\Delta h )\,u(r). \end{equation}

The increase in separation $\Delta h = r^2/(2R)$ follows from the geometric observation that the angle $\alpha$ in figure 1(b) occurs in several places. This implies that the gap-averaged fluid velocity takes the form

(3.2)\begin{equation} u(r) ={-}\frac{\dot{h} r/2}{h + r^2/(2R)} . \end{equation}

This velocity decays to 0 at $r=0$ for symmetry reasons, and has a maximum at $r=\sqrt {2hR}$. Using the corresponding maximum velocity $u_m= \dot {h}\sqrt {R/(8h)}$ and separation $h= 1\ \mathrm {\mu } {\rm m}$ (see figure 2b), we may compute the upper bound of the Reynolds number $Re=u_mh/\nu <0.5$, where $\nu$ is the kinematic viscosity of water, and $\dot {h}$ and $h$ are estimated from figure 2(b).

Since $Re$ is small, we may take the Stokes equation to govern the fluid pressure $P(r)$. Moreover, since the radius of curvature of the streamlines ${\sim }R$ is much larger than the gap, we may apply a lubrication limit of this equation. This equation will rely on the velocity field of (3.2), and may be written as

(3.3)\begin{equation} \frac{\partial P}{\partial r} = 6\mu \dot{h}\,\frac{r}{(h+r^2/(2R))^3}, \end{equation}

where $\mu$ is the dynamic viscosity. Integrating this from $r$ to $\infty$, where the fluid has the atmospheric background pressure $P_0$, yields

(3.4)\begin{equation} P(r) = P_0 - \frac{3\mu R \dot{h}}{h^2(1 +r^2/(2Rh))^2} \end{equation}

for the pressure field that exists before the bubble is formed.

Once a bubble of radius $r_c(t)$ is formed, the conservation of liquid volume inside a radius $r$ may be described as

(3.5)\begin{equation} \dot{V} (r) = {\rm \pi}(r^2-r_c^2) \dot{h} = 2{\rm \pi} r_c (h+ \Delta h(r_c) ) \dot{r}_c - 2{\rm \pi} r (h+ \Delta h (r)) u, \end{equation}

where we have neglected the mass loss due to the vapour that is formed inside the bubble. This gives the velocity

(3.6)\begin{equation} u(r) = \frac{h+r_c^2/(2R)}{h+r^2/(2R)}\,\frac{r_c}{r}\,\dot{r}_c - \frac{(r^2 -r_c^2 )\dot{h}/(2r)}{h+r^2/(2R)} . \end{equation}

With this velocity field the lubrication approximation of the Stokes equation becomes

(3.7)\begin{equation} \frac{\partial P}{\partial r}= \frac{6\mu \left( (r^2-r_c^2)\dot{h} - 2r_c \dot{r}_c (h+r_c^2/(2R)) \right) }{(h + r^2/(2R))^3r} . \end{equation}

Since the vapour pressure at room temperature is ${\sim }0.03\ \mathrm {atm}$, which will turn out to be three orders of magnitude smaller than the typical liquid pressures, we take it to be zero and integrate (3.7) from $r=r_c$, where $P=0$, to $r= \infty$, where $P=P_0$.

Following ch. 11.2 in Israelachvili (Reference Israelachvili2011), the order of magnitude of the van der Waals forces between a sphere of radius $R$ and a surface may be estimated as ${\sim }1\ \mathrm {\mu } \text {N}$ if the separation is taken as $h=10$ nm. The vertical force on the plate due to surface tension is $2{\rm \pi} r_c \sim 4\ \mathrm {\mu } \text {N}$ if $r_c=10\ \mathrm {\mu } \text {m}$. On the other hand, the local forces due to the negative water pressures $P \sim -10\ \mathrm {atm}$ is $F_P\approx {\rm \pi}r_c^2 P \sim 100\ \mathrm {\mu } \text {N}$, so that the hydrodynamic force dominates the van der Waals and capillary forces by two orders of magnitude as the bubble is first recorded.

However, the capillary pressure drop $P_c$ across the water vapour interface – which is due to the surface tension and the curvature in the vertical direction, and depends on the dynamic wetting angle $\alpha _w \approx 50^\circ$ of the water – is not negligible at the very start of the bubble growth. The pressure $P_c = 2 \sigma \cos \alpha _w/h \approx 15\ \mathrm {atm}$ when $h=40\ \mathrm {nm}$, corresponding to frame (2) of figure 2(c), quickly drops to $P_c \approx$1 atm somewhere between the $h$ values of frames (2) and (3) of figure 2(c). For simplicity, we will neglect $P_c$ in the calculation of the $r_c(t)$ evolution. As we will see, this limits the model to describe the later part of the bubble evolution (shown in frames (2)–(4) of figure 2c) when the capillary pressure becomes small compared to the viscous pressure drop.

Integration then gives

(3.8)\begin{equation} \frac{P_0R}{6\mu \dot{h}} = K(r_c) - I(r_c)\left( \frac{r_c^2}{R^2} + \frac{2r_c \dot{r}_c h}{R^2\dot{h}} \left( 1 + \frac{1}{x_c} \right)\right), \end{equation}

where we have introduced the integrals

(3.9)\begin{equation} \left.\begin{gathered} I(r_c)=R^3\int_{r_c}^{\infty} \frac{{\rm d}r}{r(h+r^2/(2R))^3} = \frac{1}{2}\left( \frac{R}{h} \right)^3 g(x_c),\\ K(r_c) = R \int_{r_c}^{\infty} \frac{{\rm d}r\,r}{(h+r^2/(2R))^3} = 2\left( \frac{R}{r_c} \right)^4 \frac{1}{(x_c+1)^2}, \end{gathered}\right\} \end{equation}

where $x_c=2Rh/r_c^2= 2{\rm \pi} Rh/A$, $A(t)={\rm \pi} \,r_c^2(t)$, and we have introduced the function

(3.10)\begin{equation} g(x_c) = \ln (x_c+1) - x_c\,\frac{1+3x_c/2}{(x_c+1)^2}. \end{equation}

An equation similar to (3.8) has been derived for the lubrication flow between parallel plates by Tirumkudulu & Russel (Reference Tirumkudulu and Russel2003). Now, solving (3.8) for $\dot {r}_c$ brings it to the form $\dot {r}_c = F(h,\dot {h},r_c)$, which allows us to solve the numerical scheme

(3.11)\begin{equation} r_c \rightarrow r_c + F(h,\dot{h}, r_c) \,{\rm d}t, \end{equation}

where $h=h(t)$ and $\dot {h} = \dot {h} (t)$ are the measured values. The result of such a calculation for a single experiment is shown in figure 4(a), while figure 4(b) shows the early-time result for all experiments.

Figure 4.

(a) The bubble radius as a function of time when $h(t)$ and $\dot {h}$ are taken as the measured values from a typical experiment (dots) along with the calculation based on (3.11) (line). Cavitation time $t_c$ is estimated as the lowest root of a parabolic fit of the data. (b) Comparison between the measured and calculated values of the first recorded bubble radii from all the experiments; the first value is shown by blue dots, and the second by green triangles. The red line has slope 1.

3.2. Pressure field

With the values of $r_c(t)$ at hand, (3.7) may be used to calculate the pressure field $P(r)$ for $r>r_c$. Integrating this equation from $r$ to $\infty$ gives

(3.12)\begin{equation} P(r) = P_0 - 6\mu \left( \left(K(r)-I(r)\,\frac{r_c^2}{R^2}\right) \frac{\dot{h}}{R} - \frac{2r_c \dot{r}_c h}{R^3}\,I(r) \left( 1+\frac{1}{x_c}\right) \right), \end{equation}

where the argument $r$ in $I$ and $K$ means that $x= 2Rh/r^2$ must be used in place of $x_c$ in (3.9).

3.3. Scaling behaviour from constant $\dot {h}$ theory

The rather weak agreement between theory and measurement that is displayed in figure 4(a) is greatly improved by the introduction of non-dimensional variables and the consideration of their scaling behaviour.

The extension of a piezoelectric element has a linear relationship with the applied voltage, and in the experiments, this voltage is increased linearly with time during the initial ramp. Therefore, the height $h$ would increase linearly in time as well, save for the effects of elastic deformations of the support system. Such an elastic response is clearly present, as is evident from the oscillations in $h(t)$, shown in figure 2(a). However, figure 2(b) shows that during the existence of the bubble, a linear increase of $h(t)$ may be a reasonable approximation, which we will apply in the following. This is also justified by the possibility of constructing future, more rigid experimental set-ups.

Introducing the characteristic time scale

(3.13)\begin{equation} \Delta t = \sqrt{\frac{12\mu R}{P_0 \dot{h}}} \end{equation}

and area

(3.14)\begin{equation} A_0 = \frac{24 \mu}{P_0\,\Delta t}\,{\rm \pi} R^2 , \end{equation}

we may rewrite (3.8) in terms of the dimensionless variables $t'=t/\Delta t$ and $A'=A/A_0$. The point of doing this is (i) to establish a dynamic equation with relevant dimensionless groups, (ii) to obtain the scaling behaviour relating the control parameters to the bubble area and lifetime, and (iii) to identify the limits to this scaling behaviour.

By assumption, $h(t)=\dot {h} t$, so the theoretical time of cavitation is

(3.15)\begin{equation} t_c = \frac{h_c}{\dot{h}}. \end{equation}

In terms of the dimensionless variables, the different terms of (3.8) then take the form

(3.16)\begin{equation} \left.\begin{gathered} \frac{P_0R}{6\mu \dot{h}}= 8 \left(\frac{{\rm \pi} R^2}{A_0} \right)^2,\\ \frac{2r_c \dot{r}_c h}{R^2\dot{h}}=\frac{A_0}{{\rm \pi} R^2}\,t'\,\frac{{\rm d}A'}{{\rm d}t'},\\ \frac{r_c^2}{R^2}= \frac{A_0}{{\rm \pi} R^2}\,A',\\ x_c= \frac{ t'}{A'},\\ K(x_c,A') = 2 \left(\frac{{\rm \pi} R^2}{A_0} \right)^2\frac{1}{A'^2(x_c+1)^2},\\ I(x_c) = 4 \left(\frac{{\rm \pi} R^2}{A_0} \right)^3 \frac{g(x_c)}{t'^3}. \end{gathered}\right\} \end{equation}

Combining the above results, (3.8) may be written as

(3.17)\begin{equation} \frac{{\rm d}A'}{{\rm d}t'} = \frac{t'}{t'+A'} \left( \frac{t'^2}{g(x_c)} \left( \frac{1}{2(t'+A')^2} - 2\right) - \frac{A'}{t'}\right) , \end{equation}

where all dimensional quantities, such as $P_0$, $R$ and $\mu$, have dropped out.

This equation may be solved numerically subject only to the initial condition $A'(t'_c) = A'_c$. The dependence on $A'_c$ may be neglected as long as $A_c \sim h_c^2$; replacing $A_c$ by its $A_c\rightarrow 0$ limit changes the subsequent $A(t)$ values by less than 1 %.

The non-dimensionalized initial times $t'_c =h_c/(\dot {h}\,\Delta t)$, which are obtained from the experiments, are shown in figure 5(b). Using the extreme values of this figure, $t'_c = 0.035$ and 0.15 (neglecting only the single $t'_c=0.25$ point), we obtain the non-dimensionalized areas $A'(t')$ shown by the red and green curves in figure 6. For the smaller $t'_c$ values, when $t_c' \le 0.001$, $A'$ has essentially converged to its $t'_c \rightarrow 0$ limit.

Figure 5.

(a) The pulling velocity $\dot {h} = \textrm {d}h/\textrm {d}t$ at the point of cavitation, along with the corresponding height $h$. The blue dots show the initial $h$ measurements, and the green dots show the values at the time immediately before the bubble is observed. The black dots show the values immediately after cavitation, and the red line is a linear fit to the latter data to guide the eye. (b) The effective non-dimensionalized initial time $t'_c =h_c/(\dot {h}\,\Delta t)$ as a function of $h_c$.

Figure 6.

Normalized bubble area versus normalized time for the extreme $t'_c$ values shown in figure 5(b), as well as the values $t'_c=0.01{-}0.0001$.

Figure 6 shows that $A'(t')$ varies significantly with $t'_c$ over the experimental range $t'_c=0.035{-}0.15$, and when $t'_c$ becomes larger than these values, the bubble will simply not grow. However, when the area (see figure 7) is normalized by its maximum value $A_m$, rather than $A_0$, and the time is normalized by the bubble lifetime $t_l$, rather than $\Delta t$, the dependence on $t'_c$ becomes much smaller. This is illustrated in figure 8(a) where the predicted variation of $A(t)/A_m$ with $t'_c$ is less than 5 % over the full experimental range. This variation is smaller than the experimental uncertainty, so the prediction is that the area may be written in the form

(3.18)\begin{equation} \frac{A(t)}{A_m} = f\left(\frac{t-t_c}{t_l}\right) \end{equation}

for a curve $f(x)$, which is the same for all the cavitation bubbles within the experimental uncertainty. The predicted data collapse of $A(t)/A_m$ is shown in figure 8(a), using the experimental data of figure 7.

Figure 7.

Bubble area versus the post nucleation time $t_c$. A parabolic fit to each of the data series from the 37 experiments yielded bubble nucleation times, maximum areas, and bubble lifetimes $t_l$. All the data were then rescaled to give the plot in figure 8(a).

Figure 8.

(a) Collapse of normalized bubble area versus normalized bubble lifetime for experimental data (dots) and theory (lines). The theory is shown for both the smallest and largest measured values of $t'_c =h_c/(\dot {h}\,\Delta t)$. (b) Maximum bubble area versus bubble lifetime from experimental data (dots), a linear fit to data of slope $-1.0$ (solid black line), the theory given by (3.19) (solid red line), and the theory using the asymptotic value $A'_m=0.0698$ (dashed black line). Note that the experimental areas and lifetimes have been computed by the maximum value and the difference of the roots of a quadratic fit of the experimental data.

The normalized bubble lifetime value $t'_l=0.770$ is obtained by inspection of the $t'_c<0.01$ graphs in figure 6. However, when $t'_c \gtrsim 0.01$, $t'_l$ and the maximum value $A'_m$ will depend on $t'_c$. The corresponding maximum area is

(3.19)\begin{equation} A_m = A_0\,A'_m(t'_c) = A'_m (t'_c)\,\frac{ 24 \mu {\rm \pi}R^2}{ P_0 }\,\frac{1}{\Delta t}, \end{equation}

where the replacement $\Delta t = t_l/t'_l$ gives

(3.20)\begin{equation} A_m = A'_m (t'_c)\,t'_l(t'_c)\,\frac{ 24 \mu {\rm \pi}R^2}{ P_0 }\,\frac{1}{ t_l} . \end{equation}

In the limit when $t'_c \ll 0.01$, we may replace $A'_m (t'_c)\,t'_l(t'_c) \rightarrow 0.0524$, and we are left with the asymptotic scaling $A_m \propto 1/t_l$.

On the other hand, we may also use (3.13) to express $A_m$ by $\dot {h}$ and $R$. This yields

(3.21)\begin{equation} A_m = 4A'_m{\rm \pi}\,\sqrt{\frac{3\mu \dot{h}}{P_0}}\,R^{3/2}, \end{equation}

which shows how the maximum bubble area depends on the radius of curvature of the lens. This is of some significance relative to the ‘lifted Hele-Shaw cell’ experiments for which $R\rightarrow \infty$. In this case, the prediction is that the bubble would always extend beyond the finite boundaries of the cell.

Figure 8(b) shows to what extent the predicted $A_m\propto 1/t_l$ scaling is fulfilled in the experiment. The theoretical prediction, which is shown in figure 8(b) (red line), is given by (3.19), which is obtained by setting $h_c=30$ nm and varying $\dot {h}$ to get a range of values for $\Delta t$ and $t'_c$ as given by (3.13) and (3.14), respectively. The corresponding asymptotic $t'_c\rightarrow 0$ result, which is shown by the dashed line, is given by (3.20). The solid black line which is a linear fit to the experimental data of slope $-1.0$, seems to indicate agreement with the asymptotic result. However, the uncertainty in the $t_l\gtrsim 1$ ms data may also yield agreement with the non-asymptotic theory (red line). The uncertainty is caused mainly by the fact that the experimental lifetimes are extracted from the curves of figure 7 by fitting these curves to parabolas and defining the lifetimes as the distance between their roots. This procedure may cause errors, in particular for longer lifetimes. However, for the $t'_c$ values where $A'_m$ and $t'_l$ have converged, i.e. when $t_l\lesssim 1$ ms, the predicted asymptotic scaling behaviour of $A_m$ is well satisfied by the experiments.

3.4. Pressure field and secondary nucleation

Figure 9 shows the pressure in the water outside the growing bubble at subsequent times. It is calculated from (3.12), using the measured values of $\dot {h}$ and the corresponding calculated values of $\dot {r}_c$. During the early expansion phase, figure 9 shows how the pressure outside the bubble occasionally falls below the initial cavitation pressure given by (3.4), and sometimes not. This variability in the water tension, and the fact that defects or impurities away from the centre may cause nucleation at pressures less negative than that required at the centre, may explain the occasional secondary nucleation events. However, the radius at which this outside pressure reaches its minimum is approximately $r_c \sim 20\ \mathrm {\mu } \textrm {m}$, a size that is reached within $\sim$0.05 ms. On the other hand, figure 2(b) shows that the largest initial $\dot {h}$ values exist for a time $\sim$0.5 ms. Hence the time available for secondary nucleation is an order of magnitude smaller than the time available for the primary nucleation, an observation that agrees well with the fact that only occasionally is such a secondary bubble observed in the experiments, as is illustrated in figure 9.

Figure 9.

(a,b) The water pressure outside the expanding bubble as a function of $r>r_c$ for the two different measured time series $h(t)$. Red curve: earliest $r_c$ value. Dashed curves: when the bubble has reached $r_c= 2.5\ \mathrm {\mu } \textrm {m}$, $5\ \mathrm {\mu } \textrm {m}$, $7.5\ \mathrm {\mu } \textrm {m}$, …, and $\dot {r}_c >0$. The blue dots represent the pressure given by (3.4) with $r=0$ at the time when the cavitation bubble is formed. (c) Image of a cavitation bubble with a secondary bubble nucleated outside the main bubble. Scale bar is $20\ \mathrm {\mu } \textrm {m}$.

3.5. Saffman–Taylor instability

As the bubble grows, a less viscous fluid (the water vapour) displaces the more viscous water, creating the conditions for the Saffman–Taylor instability to develop. Linear stability analysis (Feder, Flekkøy & Hansen Reference Feder, Flekkøy and Hansen2022) predicts that the fastest-growing perturbation has wavelength

(3.22)\begin{equation} \lambda_m = {\rm \pi}(h+ \Delta h) \sqrt{\frac{\sigma}{\mu \dot{r}_c}} , \end{equation}

where $h+\Delta h$ is the local gap width at $r=r_c$. Figure 10 shows $\lambda _m$ as a function of $h$ for all the experiments, using the $\dot {r}_c$ values obtained from (3.11) along with the measured $h$ values, and viscosity $\mu =10^{-3}\ \mathrm {Pa}\ \mathrm {s}$. Since the time step used in that calculation is smaller than the time between the frames at which $h$ is measured, the calculated $\dot {r}_c$ value that enters in (3.22) evolves numerically between each recorded $h$ value. This is why there are several $\lambda _m$ values for each $h$ value, effectively providing a measure of uncertainty in the predicted $\lambda _m$ values. The red dots, which are obtained by estimating the finger number from figure 3, fall below the approximately linear trend shown by the grey dots. The sign of this discrepancy is the same as that found by Lindner et al. (Reference Lindner, Derks and Shelley2005).

Figure 10.

The dominant Saffman–Taylor wavelength as a function of $h$. The blue and grey dots show the $\lambda _m$ values given by (3.22) using the $\dot {r}_c$ values computed from (3.11) and all the experimental $h$ values, with the blue dots corresponding to the experiment shown in figure 3. The red dots show the actual wavelength that is determined by counting the fingers in figure 3.

However, initially at small $h$ values, there is an approximate agreement between the measured and predicted $\lambda _m$ values, which strongly suggests that the finger structure is indeed caused by the Saffman–Taylor instability. At later times, the smaller $\lambda _m$ fingers will already have grown to a point where they will dominate the larger $\lambda _m$ perturbations, which come to have the larger growth rates too late to dominate the earlier growth.