Bursting bubble in a viscoplastic medium

Abstract

When a rising bubble in a Newtonian liquid reaches the liquid–air interface, it can burst, leading to the formation of capillary waves and a jet on the surface. Here, we numerically study this phenomenon in a yield-stress fluid. We show how viscoplasticity controls the fate of these capillary waves and their interaction at the bottom of the cavity. Unlike Newtonian liquids, the free surface converges to a non-flat final equilibrium shape once the driving stresses inside the pool fall below the yield stress. Details of the dynamics, including flow energy budgets, are discussed. The work culminates in a regime map with four main regimes, all with different characteristic behaviours.

1. Introduction

Bubble bursting processes abound in nature and technology and have been studied for many years in fluid mechanics (Liger-Belair, Polidori & Jeandet 2008). For example, they play a vital role in transporting aromatics from champagne (Liger-Belair 2012; Vignes-Adler 2013; Ghabache et al. 2014, 2016), and pathogens from contaminated water (Poulain & Bourouiba 2018; Bourouiba 2021). The process is also responsible for forming sea spray through ejecting myriads of droplets (MacIntyre 1972; Singh & Das 2019). Bursting bubbles also play an important role in geophysical phenomena such as volcanic eruptions (Gonnermann & Manga 2007).

In Newtonian liquids, the bubble bursting mechanism is controlled by buoyancy, surface tension and viscosity. First, the air bubble (figure 1a) being lighter than the surrounding medium, rises and approaches the liquid–air interface (figure 1b). The thin film between the bubble and the free surface then gradually drains (Toba 1959; Princen 1963) and eventually ruptures, resulting in an open cavity (figure 1c, Mason 1954). The collapse of this cavity leads to a series of rich dynamical processes that involve capillary waves (Zeff et al. 2000; Duchemin et al. 2002) and may lead to the formation of a Worthington jet (Gordillo & Rodríguez-Rodríguez 2019). In some cases, the jet might break via a Rayleigh–Plateau instability, forming droplets (Ghabache et al. 2014; Ghabache & Séon 2016). The phenomenon is so robust that it even occurs in soft granular matter, when a rising bubble also bursts at the surface, leading to a granular jet (Lohse et al. 2004).

Figure 1. Schematics for the process of a bursting bubble: (a) a gas bubble in bulk. (b) The bubble approaches the free surface forming a liquid film (thickness $\delta$) between itself and the free surface. (c) A bubble cavity forms when the thin liquid film disappears.

Earlier work on bursting bubbles used boundary integral methods in an inviscid limit (Boulton-Stone & Blake 1993; Longuet-Higgins & Oguz 1995). However, the progress made using direct numerical simulation (DNS) tools for multiphase flows (Popinet 2003, 2009; Tryggvason, Scardovelli & Zaleski 2011) has resulted in models that take into account the effects of viscosity. In fact, some recent studies have revealed how the viscosity of a liquid affects the dynamics of bursting bubbles (Deike et al. 2018; Gordillo & Rodríguez-Rodríguez 2019).

For Newtonian liquids, Deike et al. (2018) have provided quantitative cross-validation of numerical and experimental studies. They have also given a complete quantitative description of the influence of viscosity, gravity and capillarity on the process, extending the earlier work of Duchemin et al. (2002). More recently, experiments and simulations, complemented by theoretical frameworks (Gañán-Calvo 2017; Gordillo & Rodríguez-Rodríguez 2019), have resulted in a profound understanding of the physics of bubble bursting in Newtonian fluids. Appendix B provides more details on the previous studies in the Newtonian limit and compares our results with those available in the literature.

Notably, despite many applications, such as in the food industry and geophysics, the influence of rheological properties on the collapse of bubble cavities is yet to be understood. Here, we study the dynamics of bursting bubbles in a viscoplastic medium using DNS. Viscoplastic or yield-stress fluids manifest a mix of solid and fluid behaviour. The materials behave more like an elastic solid below critical stress (yield stress); however, they flow like a viscous liquid above this critical stress. Readers can find detailed reviews on yield-stress fluids in Bird, Dai & Yarusso (1983), Coussot (2014), Balmforth, Frigaard & Ovarlez (2014) and Bonn et al. (2017).

Previous experiments and simulations have been reported for trapped bubbles in a viscoplastic medium (Dubash & Frigaard 2004; De Corato et al. 2019; Sun et al. 2020), rising bubbles in yield-stress fluids (Singh & Denn 2008; Tsamopoulos et al. 2008; Sikorski, Tabuteau & de Bruyn 2009; Mougin, Magnin & Piau 2012; Dimakopoulos, Pavlidis & Tsamopoulos 2013; Tripathi et al. 2015; Lopez, Naccache & de Souza Mendes 2018) and bubbles moving inside tubes filled with viscoplastic fluids (Jalaal & Balmforth 2016; Laborie et al. 2017; Zamankhan, Takayama & Grotberg 2018). We will show that the introduction of non-Newtonian properties can significantly influence the bursting behaviour of bubbles on a free surface. At moderate values of yield stress, the collapse of the cavity can still lead to the formation of a Worthington jet, but the droplet formation might be suppressed. At high yield-stress values, the unyielded region of the viscoplastic fluid can seize the collapse of this cavity, which leads to distinct final crater shapes.

The paper is organised as follows: § 2 describes the problem and the governing parameters. Section 3 provides a phenomenological analysis and § 4 presents the different modes of energy transfer during the viscoplastic bursting process. Section 5 presents the final equilibrium shapes. The work culminates in § 6 where we summarise the different regimes observed in the process of bursting in a phase diagram. The paper ends with conclusions in § 7.

2. Numerical framework and problem description

2.1. Governing equations

We consider the burst of a small axisymmetric bubble at a surface of an incompressible Bingham fluid. To non-dimensionalise the governing equations, we remove the length and velocity scales using the initial bubble radius $R_0$ and inertia–capillary velocity $V_\gamma$, respectively. Pressure and stresses are scaled with the characteristic capillary pressure $\tau _\gamma$ (see Appendix A). The dimensionless equations for mass and momentum conservation, for the liquid phase, then read

(2.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}

(2.2)\begin{gather}\frac{\partial\boldsymbol{u}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{uu}) ={-} \boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau} - \mathcal{B}o\,\hat{\boldsymbol{e}}_{\boldsymbol{\mathcal{Z}}}, \end{gather}

where $\boldsymbol {u}$ is the velocity vector, $t$ is time, $p$ is the pressure and $\boldsymbol {\tau }$ represents the deviatoric stress tensor. We use a regularised Bingham model with

(2.3)\begin{equation} \boldsymbol{\tau} = 2\min\left(\frac{\mathcal{J}}{2\|\boldsymbol{\mathcal{D}}\|} + Oh,Oh_{max}\right)\boldsymbol{\mathcal{D}},\end{equation}

where $\|\boldsymbol {\mathcal {D}}\|$ is the second invariant of the deformation rate tensor, $\boldsymbol {\mathcal {D}}$ and $Oh_{max}$ is the viscous regularisation parameter. The three dimensionless numbers controlling the equations above are the plastocapillary number $(\mathcal {J})$, which accounts for the competition between the capillary and yield stresses, the Ohnesorge number $(Oh)$ that compares the inertial–capillary to inertial–viscous timescales and the Bond number $(\mathcal {B}o)$, which compares gravity and surface tension forces:

(2.4a–c)\begin{equation} \mathcal{J} = \frac{\tau_yR_0}{\gamma},\quad Oh = \frac{\mu_l}{\sqrt{\rho_l\gamma R_0}},\quad \mathcal{B}o = \frac{\rho_l gR_o^2}{\gamma}. \end{equation}

Here, $\gamma$ is the liquid–gas surface tension coefficient and $\tau _y$ and $\rho _l$ are the liquid's yield stress and density, respectively; $\mu _l$ is the constant viscosity in the Bingham model. Note that in our simulations, we also solve fluid motion in the gas phase, using a similar set of equations (see Appendix A). Hence, the further relevant non-dimensional groups, in addition to those in (2.4a–c), are the ratios of density $(\rho _r = \rho _g/\rho _l)$ and viscosity $(\mu _r = \mu _g/\mu _l)$. In the present study these ratios are kept fixed at $10^{-3}$ and $2 \times 10^{-2}$, respectively.

2.2. Method

For our calculations, we use the free software program Basilisk C (Popinet & collaborators 2013–2021a; Popinet 2015). The code uses a volume of fluid (VoF) technique (Tryggvason et al. 2011) to track the interface, introducing a concentration field $c$, that satisfies the scalar advection equation. Hence, (2.1)–(2.2) and their counterparts for the gas phase are solved using a one-fluid approximation, where the surface tension acts as a body force, $\boldsymbol {f}_\gamma$, only at the gas–liquid interface (Brackbill, Kothe & Zemach 1992; Popinet 2009). In dimensionless form,

(2.5)\begin{equation} \boldsymbol{f}_\gamma = \kappa\delta_s\hat{\boldsymbol{n}} \approx \kappa\boldsymbol{\nabla}c, \end{equation}

where $\delta _s$ is the interface Dirac function, $\hat {\boldsymbol {n}}$ is a unit vector normal to the interface (Tryggvason et al. 2011) and $\kappa$ is the curvature of the interface, $z = S(r)$ (as discussed by Deserno 2004, pp. 14–16):

(2.6)\begin{equation} \kappa = \frac{\dfrac{\textrm{d}^2S}{\textrm{d}r^2}}{\left(1 + \left(\dfrac{\textrm{d}S}{\textrm{d}r}\right)^2\right)^{3/2}} + \frac{\dfrac{\textrm{d}S}{\textrm{d}r}}{r\left(1 + \left(\dfrac{\textrm{d}S}{\textrm{d}r}\right)^2\right)^{1/2}}. \end{equation}

In Basilisk C, the curvature in (2.6) is calculated using the height-function method. As the surface-tension scheme is explicit in time, the maximum time step is maintained at most at the oscillation period of the smallest wave-length capillary wave (Popinet 2009; Popinet & collaborators 2013–2021b). Note that the curvature defined above is, in fact, the dimensionless capillary pressure. Hence, in the text, the wave with the largest curvature is called the ‘strongest wave’.

Basilisk C also provides adaptive mesh refinement (AMR). We use this feature to minimise errors in the VoF tracer (tolerance threshold: $10^{-3}$) and interface curvature (tolerance threshold: $10^{-4}$). Additionally, we refine based on velocity (tolerance threshold: $10^{-2}$) and vorticity (tolerance threshold: $10^{-3}$) fields to accurately resolve the regions of low strain rates. For AMR, we use a grid resolution such that the minimum cell size is $\varDelta = R_0/512$, which dictates that to get similar results, $512$ cells are required across the bubble radius while using uniform grids. We have also carried out extensive grid independence studies to ensure that changing the grid size does not influence the results. Moreover, we employ free-slip and no-penetration boundary conditions for both liquid and gas at the domain boundaries. For pressure, a zero-gradient condition is employed at the boundaries. For cases where a Worthington jet breaks into droplets, an outflow boundary condition is employed at the top boundary to ensure that the drop does not bounce off the boundary. These boundaries are far away from the bubble (size of the domain is $8R_0$) such that they do not influence the process.

Note that our numerical method uses a regularised form of the Bingham constitutive equations (see (2.3) and Appendix A). Hence, we cannot resolve the exact position of the yield surface as $\|\boldsymbol {\mathcal {D}}\|$ is never precisely zero. However, we can safely assume that low values of $\|\boldsymbol {\mathcal {D}}\|$ will be associated with the plastic regions. In our simulations, $Oh_{max}=10^8$. We have ensured that our results are independent of this regularisation parameter (Appendix E.1). The regularisation of the constitutive model also forces us to choose a criterion for the stoppage time, $t_s$. In our simulations, we consider a significantly small cut-off kinetic energy to stop the simulations (see Appendix E.2 for details).

2.3. Initial condition

This initial shape of a bubble at a fluid–fluid interface (figure 1b) can be calculated by solving the Young–Laplace equations to find the quasi-static equilibrium state for an arbitrary Bond number, $\mathcal {B}o$ (see Lhuissier & Villermaux 2012; Walls, Henaux & Bird 2015; Deike et al. 2018; Magnaudet & Mercier 2020). As a starting point, for this study, we are only concerned with the limit of $Bo \rightarrow 0$, i.e. when capillary effects dominate the gravitational ones. We choose $\mathcal {B}o = 10^{-3}$ for all the simulations in this work. For this value, the initial bubble is nearly spherical in a surrounding Newtonian liquid. Note that the bubble sphericity is a crucial assumption (simplification) in our work. The actual initial shape of the bubble depends on its size ($\mathcal {B}o$), material properties ($\mathcal {J}$, $Oh$), the method of generation and its dynamics before approaching the interface (Dubash & Frigaard 2004; Dimakopoulos et al. 2013). Furthermore, for a bubble to rise in a viscoplastic medium, the buoyancy forces should be strong enough to yield the flow (Dubash & Frigaard 2004; Sikorski et al. 2009; Balmforth et al. 2014), i.e. $\mathcal {B}o \gg \mathcal {J}$. Hence, non-spherical and non-trivial shapes might be expected (Lopez et al. 2018). For such a limit, one should first solve the full dynamics of rising bubbles to achieve the correct initial condition for the bursting problem. Note that low $\mathcal {B}o$ bubbles could still form near a free surface in other situations. One example is the process of laser-induced forward transfer, in which a laser pulse generates a bubble near the free surface of a viscoplastic liquid (Jalaal et al. 2019).

For our given initial shape, the value of $\mathcal {J}$ varies between 0 and 64. This range is selected such that we will study a full range of yield-stress effects, from the Newtonian limit ($\mathcal {J}=0$) to a medium that barely deforms owing to a large yield stress ($\mathcal {J}=64$).

Following the common assumption in these types of problems (Deike et al. 2018; Gordillo & Rodríguez-Rodríguez 2019), we assume that the thin liquid cap of thickness $\delta$ (figure 1b) disappears at $t = 0$, resulting in the configuration shown in figure 1(c), i.e. the initial condition for our simulations. In figure 1(c), $(\mathcal {R}, \mathcal {Z})$ denotes the radial and axial coordinate system. Furthermore, $\mathcal {H}_i \approx 2$ is the initial bubble depth and $\theta _i$ is the initial location of the cavity-free surface intersection. Note that the curvature diverges at this intersection in such a configuration. We smooth the sharp edge using a small circular arc to circumvent this singularity, introducing a rim with a finite curvature $\kappa _0$ that connects the bubble to the free surface. We ensured that the curvature of the rim is high enough such that the subsequent dynamics are independent of its finite value (for details, see Appendix E.3).

3. Effects of yield stress on the bursting bubble

3.1. Phenomenology

This section describes the dynamics of bursting bubbles and the qualitative effects of the plastocapillary number $(\mathcal {J})$. Figure 2 illustrates four representative cases for this purpose (supplementary movies are available at https://doi.org/10.1017/jfm.2021.489). For a Newtonian liquid (figure 2a, $\mathcal {J} = 0$), the retraction of the rim leads to the formation of capillary waves.

Figure 2. Bursting bubble dynamics for different plastocapillary numbers. (a) $\mathcal {J} = 0.0$: a typical case with a Newtonian liquid medium, (b) $\mathcal {J} =0.1$: a weakly viscoplastic liquid medium in which the process still shows all the major characteristics of the Newtonian liquid, (c) $\mathcal {J} = 0.5$: a case of moderate yield stress whereby the jetting is suppressed, nonetheless the entire cavity still yields and (d) $\mathcal {J} = 1.0$: a highly viscoplastic liquid medium whereby a part of the cavity never yields. The left part of each panel shows the magnitude of the velocity field, and the right part shows the magnitude of the deformation tensor on a $\log _{10}$ scale. The transition to the black region (low strain rates) marks the yield-surface location in the present study. The time instances in this figure are chosen to show significant events throughout the process of bursting bubbles for different $\mathcal {J}$ numbers. For all the cases in this figure, $Oh = 10^{-2}$. Movies 1–3 are available in the supplementary material.

Part of these waves travels away from the cavity, forming regions of small strain rates (black dots in figure 2a: $t = 0.45$), which are advected with the train of capillary waves. Meanwhile the other part of the waves travel down the cavity (figure 2a: $t = 0.1$) and focuses on the bottom of the cavity (figure 2a: $t = 0.45$).

Consequently, a Worthington jet is formed as depicted in figure 2(a): $t = 0.65$. Furthermore, owing to the conservation of momentum, a high-velocity jet is also formed in an opposite sense to the Worthington jet, inside the liquid pool (figure 2a: $t = 0.65$). The Worthington jet can then break into multiple droplets owing to the Rayleigh–Plateau instability (Walls et al. 2015). In the Newtonian limit, the flow continues until the free surface is fully flat, that is, when the surface energy is minimised (figure 2a: $t = 4.00$).

The introduction of the yield stress, in general, slows down the flow owing to a larger apparent viscosity. Remarkably, even at large yield stresses, the early time dynamics near the retracting rim remain unchanged owing to the highly curved interface, as clearly shown in the first panels ($t = 0.1$) of figure 2(a–d). In contrast, the anatomy of the flow inside the pool is considerably affected owing to the yield stress. At low yield stresses ($\mathcal {J} = 0.1$ in figure 2b: $t = 0.1$), everywhere near the bubble cavity yields at early times. However, as the values of plastocapillary number increases, the size of the yielded region decreases ($\mathcal {J} = 0.5$ and $\mathcal {J} = 1.0$ in figures 2c and 2d, respectively).

Furthermore, at low values of $\mathcal {J}$, the flow focusing at the bottom of the cavity persists (figure 2b: $t = 0.50$), though, owing to the increased dissipation, it is less vigorous. As a result, the jet formed post-collapse is thicker, slower and less prominent (figure 2b: $t = 1.00$) as compared with the Newtonian case (figure 2a: $t = 0.65$). Notably, for small values of $\mathcal {J}$, a Worthington jet still forms and breaks up into droplets owing to the Rayleigh–Plateau instability.

Note that unlike the Newtonian case, where the final shape is always a flat free surface, a viscoplastic medium (i.e. finite $\mathcal {J}$) comes to a halt when stress inside the liquid drops below the yield stress. Hence, the final state can feature non-zero surface energy (figure 2b: $t \ge 4.00$).

At higher values of $\mathcal {J}$, the capillary waves are so damped that the flow focusing at the bottom of the cavity vanishes. At moderate $\mathcal {J}$ numbers (figure 2c where $\mathcal {J} = 0.5$), the capillary waves are still strong enough to travel over the entire cavity (figure 2c: $t = 0.25 - 0.80$). As a result, the entire cavity yields nonetheless, and the final shape still features a deep crater (figure 2c: $t = 1.60$). On further increasing $\mathcal {J}$ such that the yield stress equals the capillary stress ($\tau _y \sim \gamma /R_0$, i.e. $\mathcal {J} \sim O(1)$), the capillary waves do not yield the entire cavity (figure 2d: $t = 0.25 - 0.75$). Hence, the final shape is a deep crater that stores a large surface energy, contrary to the final shapes at small $\mathcal {J}$ values.

In this section, we mainly focus on the effect of yield stress (via $\mathcal {J}$) on the process of bursting bubbles. Appendix C contains the discussion on the effect of $Oh$. Furthermore, in the subsequent sections, we will discuss the features explained previously in more quantitative detail. Sections 3.2 and 3.3 then delineate the travelling capillary waves and the subsequent jet formation (or lack of it), respectively.

3.2. Capillary waves in the presence of yield stress

Capillary waves are critical in the bubble bursting process (Gordillo & Rodríguez-Rodríguez 2019). Initially, the breakage of the film and the retraction of the rim create a train of capillary waves of varying strengths (Gekle et al. 2009). However, sharper waves experience very high viscous damping. As a result, the wave focusing and jet formation are controlled by the strongest wave, which is not halted by viscous damping. We follow Gordillo & Rodríguez-Rodríguez (2019) and track the strongest wave by chasing the maximum curvature of the free-surface wave ($\|\kappa _c\|$). The location of this wave, on the cavity, is denoted by the angular position, $\theta _c$ (see inset in figure 3a).

Figure 3. Effects of viscoplasticity on the travelling capillary waves. (a) Variation of the location $(\theta _c)$ of strongest capillary with time. The grey dotted line denotes the Newtonian limit, $\theta _c -\theta _i \sim -V_\gamma t$ as described by Gordillo & Rodríguez-Rodríguez (2019). (b) Variation of the strength $(\|\kappa _c\|)$ of the strongest capillary wave with time. Snapshots of the deformation tensor modulus $\|\boldsymbol {\mathcal {D}}\|$ for (c) $\mathcal {J} = 0.2$ and (d) $\mathcal {J} = 1.0$. For all the cases in this figure, $Oh = 10^{-2}$. Movies 4–6 are available in the supplementary material.

For a Newtonian liquid $(\mathcal {J} = 0)$, at low Ohnesorge numbers (e.g. figure 3), the strongest capillary wave propagates at a constant velocity $V_\gamma$, dashed line in figure 3(a). The viscous stress attenuates these waves but does not influence $\theta _c$. Previous studies (Krishnan, Hopfinger & Puthenveettil 2017; Gordillo & Rodríguez-Rodríguez 2019) have found similar results for Newtonian liquids (see Appendix B for more details). The strength of this wave decreases as it propagates down the cavity owing to continuous viscous dissipation. Around $\theta _c \approx {\rm \pi}/2$, the geometry changes leading to flow focusing resulting in an increase in the strength $(\kappa _c)$ of the wave (see figure 3c: $t = 0.2$ to $t = 0.35$). This minimum value of $\kappa _c$ depends nonlinearly on $Oh$ (see Gordillo & Rodríguez-Rodríguez 2019 & Appendix B for details).

As shown in figures 3(a) and 3(b) (and also discussed in § 3.1), the initial changes in $\|\theta _c\|$ and $\|\kappa _c\|$ remain similar to the Newtonian limit, because the highly curved region near the initial rim retraction fully yields the fluid around it. As the flow develops, the plasticity effects become more pronounced as compared to the capillary effects, and the capillary waves no longer follow the path taken by their Newtonian counterpart. The larger the value of $\mathcal {J}$, the sooner the dynamics of the capillary waves deviate from the Newtonian limit, and they become weaker. Eventually, the waves stop at a finite stoppage time, furnishing a finite final $\theta _c$ and $\|k_c\|$ (represented by $\theta _f$ and $\|\kappa _f\|$, respectively). In § 5, we will discuss the variation of these parameters for the final crater shapes.

3.3. Jet formation in the presence of yield stress

Another interesting feature of the bubble bursting process is the Worthington jet formation as the bubble cavity collapses. To characterise this jet, we track the location $(\mathcal {H})$ of the interface at the centre $\mathcal {R} = 0$. figure 4(a) shows the temporal variation of $\mathcal {H}$ for different values of $\mathcal {J}$ at fixed $Oh = 10^{-2}$. As the waves propagate, the cavity begins to collapse, hence the value of $\mathcal {H}$ decreases. A jet forms when the bottom of the cavity crosses the free surface, and $\mathcal {H}$ becomes negative (see figure 4b). For small values of $\mathcal {J}$, the Rayleigh–Plateau instability and a subsequent pinch-off occur, resulting in the kinks shown in figure 4(a). The jets eventually retract, and $\mathcal {H}$ approaches $0$, i.e. a flat final interface. As the value of $\mathcal {J}$ increases, the final value of $\mathcal {H}$ increases, approaching the upper bound of $\mathcal {H} = \mathcal {H}_i = 2$, which is set by the initial condition (twice the bubble radius, i.e. the bottom of the cavity never yields). In fact, for $\mathcal {J} \ge 0.65$, this value remains unchanged, meaning the plug region attached to the bottom of the cavity never yields. Note that, in an intermediate range of $\mathcal {J} \sim 0.35$, the interplay of the capillary waves and the yield stress results in a dimple (underdeveloped jet) that never crosses the free surface (see figure 4c).

Figure 4. Effects of viscoplasticity on the formation of the jet as a result of the collapsing cavity: (a) variation of the depth $\mathcal {H}$ of the cavity at its axis with time. The inset shows the definition of $\mathcal {H}$. Modulus of the deformation tensor $\|\boldsymbol {\mathcal {D}}\|$ for the collapse of the bubble cavity and formation of the jet for (b) $\mathcal {J} = 0.1$ and (c) $\mathcal {J} = 0.3$. Note that each kink in panel (a) is associated with the formation of a drop, as illustrated in the insets of panel (b). For all the cases in this figure, $Oh = 10^{-2}$. Movies 7–8 are available in the supplementary material.

4. What happens to the initial surface energy?

To better understand the bubble bursting dynamics in a viscoplastic medium, we also looked at the energy budgets. The total energy $E$ is the sum of the total kinetic energy of the liquid pool $E_k$, its surface energy $E_s$ and the energy dissipation $E_d$. The latter contains two parts owing to viscous ($E_d^{Oh}$) and yield stress ($E_d^{\mathcal {J}}$) contributions, $E_d = E_d^{Oh} +E_d^{\mathcal {J}}$. Lastly, small energies associated with jet breakup and airflow are summarised in $E_m$. Hence,

(4.1)\begin{equation} E = E_k(t) + E_s(t) + E_d(t) + E_m(t) = E_i, \end{equation}

where, $E_i$ is the initial energy that is purely the surface energy. Readers are referred to Appendix D for details for calculating the energy budget.

Figure 5 shows three representative examples of these energy budgets, normalised by the initial energy $E_i$. In figure 5, the time is normalised by the stoppage time $t_s$. Panel (a) shows the temporal evolution of different modes of the energy transfer for a low $\mathcal {J}$ number. Initially, at $t = 0$, the system's total energy is stored as the bubble-cavity surface energy. As the flow starts, a part of this surface energy converts to the kinetic energy of the flow generated by the travelling capillary waves. The kinetic energy reaches a maximum when the capillary waves focus at the bottom of the cavity, because the focusing process forms a region of high velocity. At the instant of focusing, for $\mathcal {J} = 0.1$ and $Oh = 10^{-1}$, $\sim$60 % of the initial energy is still present in the system as a sum of the kinetic and surface energy of the liquid pool. Subsequently, a Worthington jet forms and high dissipation is observed owing to an increase in the strain rate (see (D3)). The surface energy decreases monotonically throughout the process and reaches a finite near-zero value at the stoppage time $t_s$. This behaviour is different from a Newtonian liquid, where the surface energy would become exactly zero, as $t \to \infty$.

Figure 5. Energy budget for the process of the bubble bursting in a viscoplastic medium: temporal evolution of the different modes of energy transfers for (a) $\mathcal {J} =0.1$ and $Oh =10^{-1}$, (b) $\mathcal {J} = 1.0$ and $Oh =10^{-1}$ and (c) $\mathcal {J} =1.0$ and $Oh =10^{-2}$. (d) Comparison of the energy footprint at the stoppage time, $t = t_s$ for different $\mathcal {J}$ and $Oh$.

For $\mathcal {J} = 1.0$ (figure 5b,c), initially, the surface energy decreases monotonically until it reaches a plateau at $t = t_s$. However, contrary to the example with a small value of $\mathcal {J}$, in these cases, a major part of the cavity never yields. Consequently, more than 70 % for $Oh =10^{-1}$ and over 60 % for $Oh =10^{-2}$ of the initial energy is still stored as the crater's surface energy. In addition, note that in the limit of large $\mathcal {J}$ (and low $Oh$), yield stress is responsible for the majority of energy dissipation, i.e. $E_d^{Oh} \ll E_d^{\mathcal {J}}$.

The energy footprint at $t = t_s$ gives insight into the final static shape of the bubble. Therefore, we compare these energies for different $\mathcal {J}$ and $Oh$ numbers in figure 5(d). For all the conditions, $E_k \to 0$, and only surface energy remains in the system at the stoppage time. The remainder of the energy $(E_i - E_s)$ features as dissipation (except for those cases where drops form and $E_m$ is not negligible).

For low values of $\mathcal {J}$, the final surface energy $E_s$ is close to zero because the final craters are shallow (see § 5 for details). This residual surface energy increases with increasing $\mathcal {J}$ and $Oh$; however, the dependency on $Oh$ is negligible at higher values of $\mathcal {J}$. Finally, the dissipation owing to the yield stress $(E_d^{\mathcal {J}})$ contributes more to the overall dissipation for small $Oh$ numbers.

5. Final crater shapes

The process of bubble bursting in yield-stress fluids results in non-flat final shapes. Figure 6 shows the final crater shapes as observed for different $\mathcal {J}$ and $Oh$ numbers, and figure 7 quantifies the different features of these final shapes by analysing the location $(\theta _f)$ and strength $(\|\kappa _f\|)$ of the strongest capillary wave and the final depth of the crater $(\mathcal {H}_f)$. For the convenience of comparison, we normalise $\mathcal {H}_f$ by its initial value $\mathcal {H}_i \approx 2$ and $\|\kappa _f\|$ by the initial curvature of the bottom of the cavity $\|\kappa _i\| \approx 2$.

Figure 6. Final crater shapes: variation of the final shapes with the $Oh$ at (a) $\mathcal {J} = 0.2$, (b) $\mathcal {J} = 0.4$, (c) $\mathcal {J} = 0.6$, (d) $\mathcal {J} = 1.0$, (e) $\mathcal {J} = 5.0$ and (f) $\mathcal {J} = 10.0$.

Figure 7. Quantifying the characteristics of the final shapes as a function of $\mathcal {J}$ at different $Oh$. (a) Depth $\mathcal {H}_f$ of centreline of the final cavity surface. (b) Location $\theta _f$ and (c) Strength $\|\kappa _f\|$ of the strongest capillary wave in the final crater. The grey dashed lines in panels (b) and (c) are guides to the eye.

At low values of yield stress ($\mathcal {J} \le 0.4$), the final shape of the crater strongly depends on $Oh$ (see figures 6a,b, and 7a). When both $\mathcal {J}$ and $Oh$ are small, a Worthington jet forms (see § 3.3 for detailed discussions) which relaxes back towards the flat surface as $t \to \infty$. This jet relaxation results in shallow final cavities. As $Oh$ increases, viscous dissipation dominates the flow, the capillary waves are damped, and the change in the final cavity height becomes minimal; hence, $\mathcal {H}_f \sim \mathcal {H}_i$. In fact, for $Oh > 10^{-1}$, the capillary wave's amplitude is so close to zero that it becomes impossible to track.

As the value of yield stress ($\mathcal {J}$) increases, the effective viscosity of the domain increases and hence the initial cavity deforms less. For a highly plastic medium, the capillary waves cannot yield the entire cavity. As a result, $\mathcal {H}_f = \mathcal {H}_i$ for $\mathcal {J} \ge 0.65$, independent of the values of $Oh$. For higher $Oh$, this transition is reached (marginally) earlier (e.g. $\mathcal {J} \ge 0.5$ for $Oh \ge 1.0$).

For the cases where the bottom of the cavity never yields (figure 6d–f), we characterise the final crater shapes based on the location and strength of the frozen capillary wave ($\theta _f$ and $\kappa _f$, respectively). The variations of these values are shown in figures 7(b) and 7(c). As $\mathcal {J}$ increases, the final location of the wave is closer to the initial value, $\theta _f \approx \theta _i$. Similarly, the strength of the final wave approaches the value defined by the initial condition, as $\mathcal {J}$ increases. For this regime, the effects of $Oh$ on the final shape seem to be negligible for $\mathcal {J} > 2$.

6. Regime map

In the present study, the two crucial control parameters to describe the process of bursting bubbles in a viscoplastic medium are the plastocapillary number $\mathcal {J}$, and the Ohnesorge number $Oh$. This section uses these dimensionless numbers to summarise the observed features explained in the text, providing a regime map (or phase diagram). Figure 8 shows this map for the bursting bubble process in a viscoplastic medium. Note that we have run more than 750 simulations to arrive at this regime map, but in figure 8 we only show a few representatives at the transition lines.

Figure 8. Regime map in terms of the plastocapillary number $\mathcal {J}$ and the Ohnesorge number $Oh$ showing the transitions between the different categories identified in the current study. The insets show a representative case from each of the four regimes, namely, formation of jet which breaks into droplets (blue), formation of jet without droplets (grey), the entire cavity collapses but the cavity centre never crosses the initial pool free surface (white) and a part of the cavity never yields (red). The symbols represent simulations at the different transition lines.

For Newtonian fluids, the previous studies have found that for $Oh > 0.03$, viscous stresses dominate over the surface tension, such that the Worthington jet does not break up into droplets (San Lee et al. 2011; Ghabache et al. 2014; Walls et al. 2015). In this work, at $\mathcal {J} = 0$, we have reproduced this transition $Oh$ number (see left axis in figure 8). Increasing $\mathcal {J}$ has a similar effect on the jet break up as it manifests itself as increased apparent viscosity of the liquid. Consequently, even when $Oh \to 0$ the capillary waves get severely damped for $\mathcal {J} > 0.3$, and no droplets are formed. The blue area in figure 8 highlights the region in which a Worthington jet forms and disintegrates.

The grey area in figure 8 shows an intermediate regime in which the jet forms and crosses the free-surface line $(\mathcal {Z} = 0)$ but does not break up. This transition, for $\mathcal {J} = 0$, occurs at $Oh \approx 10^{-1}$. For non-zero $\mathcal {J}$, the transition occurs at smaller values of $Oh$ as the jet (if it forms) has less kinetic energy and cannot cross the $\mathcal {Z} = 0$ line. If the jet does not form (beyond the grey area), the collapse of the cavity results in a crater (for $\mathcal {J} \ne 0$).

As discussed in § 3.2, if the surface tension stresses are high enough, the whole cavity yields. Otherwise, for large values of the yield stress $(\mathcal {J} \sim O(1))$, the plug region attached to the bottom of the cavity never yields. As a result, the bottom of the cavity does not move, i.e. $\mathcal {H}_f = \mathcal {H}_i$. This transition from a fully yielded cavity to the cavity with an unyielded bottom is highlighted in figure 8 with the red line.

7. Conclusions

In this work, we have studied the capillary-driven process of bursting bubbles in a viscoplastic medium. As with Newtonian fluids, flow begins when the rim, which connects the bubble cavity to the free surface retracts. Consequently, the fluid is yielded, and a train of capillary waves is generated. The yield stress significantly affects the flow structure inside the pool by making plug regions. The higher the value of the yield stress, the larger is the deviation from the Newtonian counterpart.

Following the analyses of Deike et al. (2018) and Gordillo & Rodríguez-Rodríguez (2019), we provided information on the dynamics of the capillary waves as they travel down the bubble cavity. In liquids with low yield stresses, the cavity collapse leads to a Worthington jet that might break up into drops as a result of a Rayleigh–Plateau instability. However, for liquids with a large yield stress, the jet vanishes. The energy budgets analysis gives insight into the dynamics by showing how the initial surface energy is dissipated. Eventually, in contrast to the Newtonian fluids, where the final state is always a flat film, bubble bursting in a viscoplastic medium results in final crater shapes with high residual surface energy. We have analysed the geometry of these shapes as a function of the governing control parameters, namely, the Ohnesorge and the plastocapillary numbers. Lastly, we use the same numbers to categorise the four different regimes in viscoplastic bubble bursting (see the phase diagram in figure 8).

Our study has direct applications in a range of industrial operations, where bubbles are present at the surface of a yield-stress fluid. The focus of this work is to compare the bursting bubble process in a yield-stress fluid with that of a Newtonian fluid, without the initial shape effects. However, once the exact shape of the bubble at the free surface is known, either from experiment or theory, one can calculate the resulting flow and compare with the present study. Moreover, the current results could be useful in analysing some geophysical flows, such as those in volcanic eruptions. In a broader perspective, the work presents a system in which surface tension and yield stress are the main factors. Such a system is of fundamental interest in design and manufacturing at small scales where capillary action is competing with the yield stress, e.g. in three-dimensional printing and coating with polymeric fluids (Rauzan et al. 2018; Jalaal et al. 2019; Nelson et al. 2019; Jalaal, Stoeber & Balmforth 2021).