# Jets drops ejection in bursting gas bubble process

Sanda-Carmen Georgescu, Jean-Luc Achard, Edouard Canot

► To cite this version:

Sanda-Carmen Georgescu, Jean-Luc Achard, Edouard Canot. Jets drops ejection in bursting gas bubble process. European Journal of Mechanics - B/Fluids, 2002, 21 (2), pp.265-280. 10.1016/S0997-7546(01)01177-3. hal-00182311

HAL Id: hal-00182311

https://hal.science/hal-00182311v1

Submitted on 7 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Jet drops ejection in bursting gas bubble processes

Sanda-Carmen Georgescu<sup>a</sup>, Jean-Luc Achard<sup>b,*</sup>, Édouard Canot<sup>c</sup>

<sup>a</sup> Hydraulic Department, University “Politehnic” of Bucharest, 313 Spl. Independentei, 77206, Bucharest, Romania

<sup>b</sup> Laboratoire des Ecoulements Géophysiques et Industriels, BP 53, 38401, Grenoble cedex 9, France

<sup>c</sup> IRISA-INRIA, Campus de Beaulieu, 35042, Rennes cedex, France

The numerical simulation using a boundary element method is presented for a gas bubble bursting at a free surface in a potential flow with a viscous fluid assumption. Systematic comparisons are given with experimental data on the first “jet drop” size in relation with the parent bubble size, and on the critical bubble radius above which no jet drop forms. The computations were made for different liquids. It is pointed out that an exact description of the jet formation and break up requires the complete Navier–Stokes equations only in the final phase of the evolution.

*Keywords:* Collapsing bubble; Aerosols; Boundary element method

## 1. Introduction

The droplets formation that follows gas bubbles bursting at a free liquid surface plays an important role in mass transfer between phases in stratified two-phase flows. Many applications can be found in chemical and nuclear engineering, as well as in the geophysical field. Processes that involve a wide spectrum of bubble sizes, such as aeration, boiling, degassing, distillation, or evaporation are the most common examples. The considerable amount of salt transfer at sea surface is made through the aerosols ejected in the atmosphere. Moreover, new health and environmental problems appear from the liquid droplets capacity of carrying pathogenic or radioactive substances collected on their surface.

The bursting bubble process leads to the production of two types of droplets:

(i) when rising gas bubbles reach a free liquid surface, the film formed between the bubble's cap and the adjacent gaseous phase is rapidly drained and finally atomised; the ejected tiny film fragments represent the “film droplets”;

(ii) the gaping bubble cavity collapses: a surface wave converges axially and generates an upward high-speed liquid jet; the successive fragmentation of this jet produces some “jet drops”.

Due to the micronic size, Resch et al. [1] performed the first successful visualisations of film droplets only in 1986, by using holographic imaging techniques. The film droplets ejection takes up to 300 µs after the film cap rupture. The number and film droplets size spectrum have been studied by Wu [2].

From the dragged liquid mass point of view, jet drops have been experimentally investigated since 1954 through high-speed photography [3]. It was pointed out that the size of the first released jet drop is about one-tenth of the parent bubble size. There is no significant change in the behaviour of cavity collapsing process over a range of 0.1 to 0.9 mm parent air bubble radius in

<footnote>* Correspondence and reprints.</footnote>

<footnote>E-mail addresses: carmen@hydrop.pub.ro (S.-C. Georgescu), jean-luc.achard@hmg.inpg.fr (J.-L. Achard), edouard.canot@irisa.fr (É. Canot).</footnote>
fresh and seawater, or seawater with oleic acid surface film. The jet drop size distributions were depicted by Garner et al. [4] for air/distilled water, and air/ethanol couples. Typical holograms showing the sequential formation of jet drops were restored by Resch et al. [1] to obtain the true size and three-dimensional location of the drops. The time spanned between film cap rupture and the first jet drop ejection is of the order of milliseconds. The first jet drop velocity rapidly decreases from tens to several m s⁻¹ when bubble size increases from tens of microns to several millimetres. Jet drops number does not exceed 10 and decreases as the bubble diameter increases. For an air/seawater couple, bubbles larger than about 3 mm radius do not eject any jet drop [5,6].

Lately some experimental data are available: Suzuki and Mitachi [7] studied how the number and dimensions of ejected jet
drops vary according to liquid properties and parent bubble size. They also focused on the critical bubble diameter above which
the liquid jet decays without releasing any jet drop.

Due to strong non-linearities, free boundary conditions, and breaking processes, the bursting bubble problem offers
a stimulating challenge to numerical modelling. The first numerical simulations of bubble bursting at a free pure water surface
were made by Boulton-Stone and Blake in 1993 [8], and then improved by Boulton-Stone [9] for a contaminated surface
case. Both studies follow an inviscid fluid model through a boundary element method, completed by boundary-layer effect
consideration. Hence, it follows that surfactants have a significant influence on the studied process, the surface dilatational
viscosity preventing the jet drop formation. The bursting process for large air bubbles in water was modelled by Sussman and
Smereka in 1997 [10], by solving the Navier-Stokes equations in both fluids on a fixed grid, by using a level set method. The
initially submerged spherical bubble was released just below a flat free surface. Due to the computational effort, the bubble
equivalent radius cannot be smaller than 4 mm, and there are no comparisons with the previous numerical studies, which cover
a range of bubble radii from 0.5 to 3 mm.

Our aim is to show that a purely irrotational flow model suffices to reproduce most observed experimental data in the first jet drop geometric and kinetic characteristics, and the critical parent bubble size, for different liquids. To reach this goal, firstly emphasis has been placed on numerical accuracy: a boundary element method has been selected with a second-order time-evolution scheme. Secondly, viscous effects have been incorporated to some extent as it is allowed for potential flows of fluids with constant viscosity [11].

Initial geometric conditions and governing motion equations for the physical model are depicted in Section 2 of the present
paper. Numerical procedure is outlined in Sections 3 and 4. Results and comments are presented in Section 5.

**2. Problem statement**

We focus exclusively on the jet drop formation, which represents the final sequence of a bursting bubble process. The entire
set of phenomena can be divided into three stages:

(i) the approaching stage is defined by the decelerated arrival of the bubble at the free surface simultaneously with the liquid free surface deformation; before bursting occurs, the bubble oscillates vertically with a very small amplitude, and reaches its equilibrium position after oscillations decay [5];

(ii) the thinning stage starts by the liquid film cap drainage that reduces the film thickness down to a certain critical value, and finishes by film disintegration, and film droplets ejection;

(iii) the break up stage includes the bubble cavity collapsing process that creates an unstable liquid jet, which splits up into
several jet drops. Modelling only the behaviour of this last stage implies the difficulty in expressing realistic initial
conditions.

2.1. *Initial geometric conditions*

The coupling between the break up stage of bursting bubble process and the preceding-ones may be simplified. A rising bubble with small terminal velocity once arrived at the liquid free surface, becomes quasi-immobile beneath it and reaches an equilibrium position. The shape of the bubble surface together with the deformed free surface evolves slowly up to the bursting moment: this position can be expressed by a macroscopic balance, in which the film thickness can be neglected with respect to the bubble equivalent radius R₀ (the radius of a sphere with the same volume). After reaching the equilibrium position, the film can be removed artificially in order to dispose of the connected initial interface position requested by the break up stage simulation.

We consider a bounded axisymmetric liquid domain Ω_L, extending largely over the bubble size, in order to obtain
a negligible influence of the phenomena at the periphery. The interfaces shape corresponds to a stationary bubble entrapped
beneath a free surface [12,13]. The resulting bubble cavity surface Σ_b, film equivalent interface Σ_f, and external meniscus of
<figure><img src="image_2.png" /><figcaption>Fig. 1. Air bubble entrapped beneath a free liquid surface: resulting three joined interfaces.</figcaption></figure>

the free liquid surface $\Sigma_m$ are in equilibrium with the surrounding fluid, being joined along the crater line of radius $R_C$ (Fig. 1). At any point of the interface between two phases i and j, there is a positive pressure jump of magnitude

$$p_i - p_j = \sigma_{ij} (1/R_1 + 1/R_2), \qquad (2.1)$$

when passing towards the phase i where the centre of curvature lies; $R_1$ and $R_2$ are the corresponding local principal radii of curvature at the point of interest. When applying the equilibrium condition (2.1) on each of the three previously defined surfaces, $\sigma_{ij}$ equals the surface tension $\sigma$ on $\Sigma_b$ and $\Sigma_m$, and it equals the film tension $\gamma$ on $\Sigma_f$.

At the free liquid surface, the adjacent gas pressure $p_0$ is assumed to be constant since the gas inertia is neglected. Due to the gravity and negligible change of the phase density $\rho_i$ with $z$, the pressure is: $p_i(z) = p_i(0) - \rho_i gz$. When the subscript $i = G$, the relationship refers to the gaseous phase inside the bubble, and when $i = L$, it refers to the liquid phase; $p_i(0)$ represents the corresponding reference pressure at the crater level, where $z = 0$.

The normal force balance at the triple contact point C is written [12] as

$$\gamma \sin \beta_{Cf} = \sigma \sin \beta_{Cb} + \sigma \sin \beta_{Cm}, \qquad (2.2)$$

where the azimuthal angle at the contact point $\beta_{Ck}$ ($k = b, f, m$) is defined on each interface $\Sigma_k$ between the radial unit vector $e_r$ and the corresponding tangent unit vector $t_k$. We add the geometric relation obtained for the radius of curvature of the film assumed to be a part of a sphere [14]:

$$R_f = R_C / \sin \beta_{Cf}. \qquad (2.3)$$

The dimensionless form of the resulting system of Eqs. (2.1)–(2.3) is obtained by taking the bubble equivalent radius $R_0$ as length scale, and is solved numerically according to a shooting method [15], by considering the Froude number $Fr = \sigma/(\rho_L g R_0^2)$ as control parameter. The boundary conditions are matched at the end of the integration through a globally convergent Newton–Raphson method [16]. Bubbles having an equivalent radius less than about 0.5 mm are nearly spherical, while larger bubbles are distorted. For the same size, the bubble equilibrium shape is less distorted in water than in liquids of low surface tension (e.g., ethanol), this being depicted also from photographs [3,7]. We present in Fig. 2 the computed equilibrium position (in dimensionless polar coordinates $r^*$ and $z^*$) of an air bubble of 1.5 mm equivalent radius, entrapped beneath a free surface of ethanol ($Fr = 1.294$), and seawater ($Fr = 3.549$), respectively.

Once we get the equilibrium position of the three interfaces joined at the crater line just before the break up moment, the initial geometric conditions requested by the bursting simulation are well established by a mathematically connected interface $\Sigma_i = \Sigma_b \cup \Sigma_m$, whose shape is defined by the bubble cavity joined only to the external meniscus, the film being removed.

## 2.2. Governing motion equations

The liquid domain $\Omega_L$ has the boundary $\Sigma = \Sigma_i \cup \Sigma_w$, where $\Sigma_w$ is the solid wall and base surface of a cylindrical pool. The depth of the pool equals its radius, and is large enough to neglect any influence of the studied phenomena at this solid boundary. We consider a unit outward normal $n$ on $\Sigma$. It is convenient to adopt $\sigma/R_0$ as pressure scale, $(\sigma/R_0\rho_L)^{1/2}$ as velocity scale, and the ratio between the length scale $R_0$ and velocity as time scale [8]. The dimensionless equations governing the potential fluid flow are classically the Laplace equation for the velocity potential $\phi^*$ and Euler's equation:

$$\nabla^2 \phi^* = 0, \qquad (2.4)$$

$$\frac{\partial(\nabla\phi^*)}{\partial t^*} + \nabla \frac{v^{*2}}{2} = -\nabla p_L^* - \frac{1}{Fr}\nabla z^*. \qquad (2.5)$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Equilibrium position of an air bubble with R₀ = 1.5 mm, entrapped beneath a free surface of (a) ethanol, Fr = 1.294, and (b) seawater, Fr = 3.549.</figcaption></figure>

A potential flow assumption is appropriate because of the impulsive character of the phenomenon. In this type of model, viscous effects can be partially considered through a boundary condition, namely the normal momentum balance at any point of the interface Σᵢ. In this condition appears the normal viscous stress at the interface, 2μ<sub>L</sub>(∂v<sub>n</sub>/∂n), where μ<sub>L</sub> represents the dynamic viscosity of the liquid, and v<sub>n</sub> is the normal component of the velocity. In dimensionless form, this condition becomes:

$$ \frac{1}{We} \left( \frac{1}{R_1^*} + \frac{1}{R_2^*} \right) - p_0^* + p_L^*(z^*) - \frac{2}{Re} \frac{\partial^2 \phi^*}{\partial n^{*2}} = 0, \qquad (2.6) $$

where, according to the choice of scales, the Weber number takes the unit value (We = 1), and the Reynolds number is Re = (ρ<sub>L</sub>σR<sub>0</sub>)<sup>1/2</sup>/μ<sub>L</sub>. A non-penetrability condition at the interface Σᵢ must be added. Combining Eqs. (2.5) and (2.6) to reduce the pressure terms, we obtain the Bernoulli's equation yielding the local time derivative of the velocity potential:

$$ \frac{\partial \phi^*}{\partial t^*} = -\frac{v^{*2}}{2} + \frac{1}{We} \left( \frac{1}{R_1^*} + \frac{1}{R_2^*} \right) - \frac{z^*}{Fr} - \frac{2}{Re} \frac{\partial^2 \phi^*}{\partial n^{*2}}. \qquad (2.7) $$

The expressions of the normal second derivative of the velocity potential ∂²φ*/∂n*², outside of the symmetry axis, and on that axis are defined in Appendix A. The initial kinetic conditions correspond to a velocity field, which obviously is assumed to be zero [8], the interface Σᵢ being frozen till the beginning of the collapsing process.

### 3. Numerical method

The bubble cavity collapse represents a transient free-boundary problem that involves two types of calculations:

(i) The evolution problem is successively divided into tiny time steps Δt*. At a fixed moment t*, we solve the Laplace equation (2.4) to obtain the velocity potential values φ*(t*), hence, the corresponding normal component, and tangential component of the velocity.

(ii) Then, a time-stepping scheme allows the connection of two successive steps to determine the new potential values, and the new interface position at the following moment (t* + Δt*).
The boundary element method (BEM) is well adapted to the first type of calculation: it replaces the worthwhile Laplace's equation extended in the whole liquid domain $\Omega_L$, by a second kind Fredholm integral equation extended only on the boundary $\Sigma$. For this direct BEM computation [17,18], the velocity field is generated by source and normal doublet type singularities spread over $\Sigma$. The integral equation on the boundary $\Sigma$ is written as

$$ -\int_{\Sigma} \phi^*(P) \frac{\partial(1/|MP|)}{\partial n_P^*} d\Sigma_P + \int_{\Sigma} (1/|MP|) \frac{\partial\phi^*(P)}{\partial n_P^*} d\Sigma_P = 2\pi\phi^*(M), \quad M \in \Sigma, \qquad (3.1) $$

where $|MP|$ is the distance between the observation point $M$ and the singularity point $P$. We need to define kinetic conditions of Dirichlet type on the gas-liquid interface $\Sigma_i$, where the velocity potential $\phi^*$ is known, and of Neumann type on the motionless solid surface $\Sigma_w$, where the normal velocity vanishes: $\partial\phi^*/\partial n^* = 0$.

The temporal interface evolution is determined through a Lagrangian description of a variable $n$ number of nodes denoted by $M = M(r^*(t^*), z^*(t^*))$, unevenly redistributed on the boundary $\Sigma(t^*)$ at each time step (see Appendix B), with respect to some criteria like the adaptation at the surface gradients [17]. This leads to a concentration of nodes at places where the interface's curvature is important, or where two portions of the interface approach one another. Thus, the simulation of the interface evolution is performed during the whole collapse process, without applying smoothing techniques that can affect the physics of the phenomenon.

On each boundary element, cubic splines define the geometric variables, and cubic Hermite polynomials approximate the field variables.

Following an explicit numerical scheme [19] improved for the capillary [20], and viscous effects [15,21], the time progression is made through second-order limited Taylor series expansions:

$$ \chi(t^* + \Delta t^*) = \chi(t^*) + \left. \frac{D\chi}{Dt^*} \right|_{t^*} \Delta t^* + \left. \frac{D^2\chi}{Dt^{*2}} \right|_{t^*} \frac{\Delta t^{*2}}{2!} + O(\Delta t^{*3}), \qquad (3.2) $$

where $\chi(t^*)$ denotes $\phi^*(t^*)$, $r^*(t^*)$, and $z^*(t^*)$, respectively. We mention that the matrix of the BEM calculations following (3.1) depends only on geometrical quantities. Thus it has to be computed only once at each time step, and this explains why the temporal scheme (3.2) is very efficient. We add that as emphasised by Machane and Canot [20], “the order of the time-stepping scheme must be consistent with the order of the BEM solver used”. In the case of mixed capillary-gravity waves, a cubic BEM formulation is used ($C^1$ shape basis functions leading $O(\Delta s^3)$ errors). According to the consistency condition, a second order time stepping procedure must be chosen [22].

The material derivative of the velocity potential is: $D\phi^*/Dt^* = \partial\phi^*/\partial t^* + (\nabla\phi^*)^2$, where the local time derivative $\partial\phi^*/\partial t^*$ is defined by the Bernoulli's Eq. (2.7). The corresponding second-order Lagrangian derivative becomes

$$ \frac{D^2\phi^*}{Dt^{*2}} = v_r^* \frac{D^2 r^*}{Dt^{*2}} + v_z^* \frac{D^2 z^*}{Dt^{*2}} + \frac{1}{We} \frac{D}{Dt^*} \left( \frac{1}{R_1^*} + \frac{1}{R_2^*} \right) - \frac{1}{Fr} v_z^* - \frac{2}{Re} \frac{D}{Dt^*} \left( \frac{\partial^2 \phi^*}{\partial n^{*2}} \right). \qquad (3.3) $$

The normal component of the velocity $v_n^* = \partial\phi^*/\partial n^*$, and the tangential-one $v_t^* = \partial\phi^*/\partial s^*$, allow the determination of the radial component of the velocity $v_r^* = (-v_n^* \sin\beta + v_t^* \cos\beta)$, together with the axial component $v_z^* = (v_n^* \cos\beta + v_t^* \sin\beta)$. The angle $\beta$ is defined between the radial unit vector $e_r$ and the tangent unit vector $t$. In a meridian plane, the curvilinear abscissa $s$ starts from the bottom of the bubble cavity, follows the interface $\Sigma_i$, then the solid surface $\Sigma_w$, and ends on the $Oz$ axis. The expressions of Lagrangian derivatives that appear in the right hand term of (3.3) are depicted in Appendix A [15]. All those terms need the calculation of the normal and tangential derivatives of $\partial\phi^*/\partial t^*$. This implies the computation of an extra Laplace equation issued upon application of a time derivative to the Eq. (2.4):

$$ \nabla^2(\partial\phi^*/\partial t^*) = 0. \qquad (3.4) $$

The boundary conditions are: $(\partial\phi^*/\partial t^*)|_{(2.7)}$ on the interface $\Sigma_i$, and $\partial(\partial\phi^*/\partial t^*)/\partial n^* = 0$ on the motionless solid surface $\Sigma_w$.

## 4. Numerical stability and precision

The time step $\Delta t^*$ cannot obviously remain constant during the whole process. It has to be selected at each time through a stability criterion, derived from an eigenvalue problem based on the above temporal scheme [20,22,23]. This latter explicit scheme is stable in the sense of classical numerical analysis, but it diverges at infinity with an exponential-like law. This explains why we have also used, without computational time penalty, an explicit/implicit scheme [23], which provides conditional
stability over a much longer evolution time (the Eulerian part of the scheme – BEM solver – is explicit, but the Lagrangian part can be solved under an implicit form). The time step is defined by a gravity–capillary waves criterion [22]:

$$ \Delta t^* \le 2[\pi^3 We^{-1}(\Delta s_{\min}^*)^{-3} + \pi Fr^{-1}(\Delta s_{\min}^*)^{-1} - 1/2], \quad (4.1) $$

where $\Delta s_{\min}^*$ is the minimum value of the arc length measured between two consecutive points of the boundary $\Sigma(t^*)$.

In the absence of non-linear analytical tests, the evaluation of the numerical code precision is made through a global mechanical energy balance expressed only in surface integrals terms [22]:

$$ \frac{d}{dt^*} \left( \frac{1}{2} \int_{\Sigma} \phi^* \frac{\partial \phi^*}{\partial n^*} dA^* + \frac{1}{We} \int_{\Sigma_i} dA^* + \frac{1}{2Fr} \int_{\Sigma_i} z^{*2} dA^* \right) = -\frac{2}{Re} \int_{\Sigma} \frac{\partial \phi^*}{\partial n^*} \frac{\partial^2 \phi^*}{\partial n^{*2}} dA^*, \quad (4.2) $$

where $dA^*$ is the axisymmetric surface element. Upon an azimuthal integration, the mechanical energy balance deals with line integrals of the terms computed through the BEM, being easily implemented in the computational procedure. Let's observe that the dissipation of the mechanical energy due to shear viscosity depends only on the normal component of the velocity, and on the normal second derivative of the velocity potential $\partial^2\phi^*/\partial n^{*2}$.

The form (4.2) of the mechanical energy balance was successfully used as validation test of the code accuracy in transient free-boundary axisymmetric and bidimensional problems solved through the BEM [24]. More precisely, one has introduced the relative mechanical energy balance deviation $\delta E$, considered to be the ratio (in percent) of the difference between left and right hand terms of (4.2) to the initial mechanic energy, and plotted it versus time in each case. It was shown that the global accuracy is strongly linked to the time step choice: numerical instabilities tend to develop when the time and space steps diminish. In a decelerated process, like the oscillating bidimensional globule of viscous fluid, the $|\delta E|$ values damp from about 2.5% to less than 0.5%.

The collapsing process we are studying in this paper is accelerated. By decreasing the time step, the problem becomes more and more stiff (the time constant related to the capillary oscillation period is very small). It will be seen that when the jet rises, the balance (4.2) is no more accurate, but still preserve an admissible range of $|\delta E|$ that varies typically from zero to less than 5%. Towards the end of the collapse, $|\delta E|$ keeps on increasing especially for the small bubble cases, which involve the greatest velocity values.

## 5. Results and discussion

Numerical simulations of the break up stage of an isolated air bubble bursting at a free surface has been systematically performed for several pure Newtonian liquids (see Table 1) and different parent bubble sizes (from 0.5 mm radius up to the critical bubble radius), in order to compare numerical and experimental data. To begin with, consider typical interface evolutions during the collapse process, started from the equilibrium position, associated to an air bubble in distilled water at 293 K, for an equivalent radius of $R_0 = 0.75$ mm (Fig. 3(a)), $R_0 = 2$ mm (Fig. 4(a)), and $R_0 = 3$ mm (Fig. 4(d)), respectively. The bursting processes in Figs. 3(a) and 4(a) correspond to bubble size smaller than the critical size. In this case, the jet splits up and ejects successive jet drops. The bursting process in Fig. 4(d) corresponds to a bubble size greater than the critical size. For this situation, the jet decays without releasing any jet drop. The corresponding precision is given in the Figs. 3(c), 4(c), and 4(f).

The criterion that defines if a jet drop is going to be ejected is based on a complete jet strangulation (the radial position of another interface node than the central-one becomes zero). Just after the rupture, the surface nodes are repositioned, both on the upper part of the liquid jet, and on the ejected drop, in accordance with the criteria imposed by the regridding procedure [17]. For a bursting bubble computation ending at the first jet drop ejection, we need about 1800–3400 iterations, which are performed in less than 2–4 minutes, on a Pentium III 866 MHz PC (with an average of 14.4–16.2 iterations per second). The velocity of the central node (the one on the Oz axis) increases sharply, and reaches a maximum when the jet starts to form, that instant being linked to a viscous dissipation peak. Further, the velocity magnitude decays to an almost constant level that corresponds to the first jet rupture moment, and approaches the first jet drop ejection speed. Then, the jet suddenly decays and rises again, in accordance with the unstable successive split up process. The temporal variations of the central node dimensionless velocity plotted in Figs. 3(b), 4(b), and 4(e), show that the maximum velocity value decreases with increasing bubble radius. For parent bubble sizes greater than the critical-one, the jet velocity decays slowly after reaching the maximum (Fig. 4(e)).

Now, let us compare various characteristics of our simulations to corresponding data (experimental or numerical) available in the literature. The critical bubble radius will deserve a separate treatment. Finally, the influence of a different type of initial condition (non-equilibrium approaching bubble) will be addressed.
<table><caption>Table 1
Critical bubble radius $R_{cr}$ for different air/liquid couples. The liquid physical properties are considered at 293 K, excepting the chlorofluorocarbon (CFC-113) that is considered at 274 K [7,27]</caption><thead><tr><th rowspan="2">Liquid type</th><th rowspan="2">ρ<sub>L</sub><br>[kg m<sup>-3</sup>]</th><th rowspan="2">μ<sub>L</sub><br>[mPa s]</th><th rowspan="2">σ<br>[N m<sup>-1</sup>]</th><th colspan="2">Experimental data<br>(approximation from plotted points)</th><th rowspan="2">Computed<br>values<br>R<sub>cr</sub><br>[mm]</th></tr><tr><th>Reference</th><th>R<sub>cr</sub><br>[mm]<br>min.-mean-max.</th></tr></thead><tbody><tr><td>Acetic acid</td><td>1046</td><td>1.220</td><td>0.0276</td><td>[7]</td><td>1.90-2.00-2.06</td><td>2.32</td></tr><tr><td>Acetone</td><td>791</td><td>0.324</td><td>0.0233</td><td>[7]</td><td>2.03-2.28-2.37</td><td>2.46</td></tr><tr><td>Aniline</td><td>1022</td><td>4.466</td><td>0.04267</td><td>–</td><td>–</td><td>2.04</td></tr><tr><td>Benzene</td><td>885</td><td>0.647</td><td>0.02888</td><td>–</td><td>–</td><td>2.42</td></tr><tr><td>Carbon tetrachloride</td><td>1594</td><td>0.975</td><td>0.0262</td><td>–</td><td>–</td><td>1.29</td></tr><tr><td>CFC-113</td><td>1617</td><td>0.960</td><td>0.0199</td><td>[7]</td><td>0.77-1.06-1.34</td><td>1.19</td></tr><tr><td>Distilled<br>water</td><td>998.2</td><td>1.002</td><td>0.07274</td><td>[4]</td><td>2.50</td><td>2.87</td></tr><tr><td>Ethanol</td><td>798</td><td>1.194</td><td>0.0228</td><td>[7]</td><td>2.25-2.50-2.70</td><td>2.16</td></tr><tr><td>Ethan.-water<br>2-98 wt%</td><td>993</td><td>1.090</td><td>0.0655</td><td>[7]</td><td>1.84</td><td>2.54</td></tr><tr><td>Glycer.-water<br>85-15 wt%</td><td>1220</td><td>110.8</td><td>0.0657</td><td>–</td><td>–</td><td>1.98</td></tr><tr><td>Glycer.-water<br>79-21 wt%</td><td>1210</td><td>49.7</td><td>0.0664</td><td>–</td><td>–</td><td>2.15</td></tr><tr><td>Glycer.-water<br>70.5-29.5 wt%</td><td>1185</td><td>23.5</td><td>0.0676</td><td>–</td><td>–</td><td>2.34</td></tr><tr><td>Glycer.-water<br>65-35 wt%</td><td>1170</td><td>15.4</td><td>0.0683</td><td>–</td><td>–</td><td>2.60</td></tr><tr><td>Glycer.-water<br>50-50 wt%</td><td>1130</td><td>5.98</td><td>0.0694</td><td>–</td><td>–</td><td>2.72</td></tr><tr><td>i-Propanol</td><td>791</td><td>2.420</td><td>0.0214</td><td>[7]</td><td>1.24-1.31-1.40</td><td>1.42</td></tr><tr><td>Methanol</td><td>804</td><td>0.560</td><td>0.0226</td><td>[7]</td><td>1.97-2.17-2.22</td><td>2.40</td></tr><tr><td>Mineral oil</td><td>840</td><td>17.00</td><td>0.0290</td><td>–</td><td>–</td><td>2.18</td></tr><tr><td>n-Butanol</td><td>816</td><td>2.980</td><td>0.0247</td><td>[7]</td><td>1.00-1.12-1.25</td><td>1.28</td></tr><tr><td>n-Propanol</td><td>805</td><td>2.190</td><td>0.0237</td><td>[7]</td><td>1.47-1.50</td><td>1.87</td></tr><tr><td>Seawater<br>(NaCl 34%)</td><td>1024</td><td>1.126</td><td>0.08023</td><td>[5]</td><td>3.00</td><td>2.94</td></tr></tbody></table>

## 5.1. General characteristics

The computed first jet drop ejection speed seems to approach the real one. The latter decreases sharply from 80 m s⁻¹ for a 35 µm parent bubble radius, to 3.5 m s⁻¹ for an air bubble of 1 mm radius in seawater [25], and approaches 0 m s⁻¹ for bubbles of about 2.7 mm radius, since in experiments this is the upper parent bubble size leading to jet drop ejection for an air/distilled water couple [7]. We also found that the closely related characteristic, i.e. the first drop ejection time, increases with increasing bubble size. Our computed values obtained for an air/distilled water couple at 293 K (e.g., t = {1.702, 2.842, 4.075, 7.458, 10.506} ms for $R_0 = \{0.5, 0.75, 1, 2, 2.5\}$ mm, respectively) are slightly smaller than those of Boulton-Stone and Blake [8] for small parent bubbles, and greater for radii $R_0 \ge 2$ mm. For an air/seawater couple we found t = 7.875 ms for $R_0 = 1.5$ mm. This computed jet drop ejection time is approaching the corresponding experimental data depicted from Resch et al. [1].

The computed dimensionless values of the first jet drop radius $R_d^*$, defined as the radius of a sphere of the same volume, are plotted versus the dimensional equivalent bubble radius $R_0$ for an air/distilled water couple at 293 K (Fig. 5(a)), and air/ethanol couple at 293 K (Fig. 5(b)), together with the available experimental data [3,4,7], and other numerical results [8]. To underline the wide range of the jet drop sizes, we plot in Fig. 5(b) the experimental minima, mean, and maxima radii of droplets obtained by Suzuki and Mitachi [7] for an air/ethanol couple. Our results are slightly overestimated with respect to the corresponding
<figure><img src="image_4.png" /><figcaption>Fig. 3. Interface collapsing evolution starting from an equilibrium position (a), dimensionless velocity of the central node (b), and relative mechanical energy balance deviation $\delta E$ versus the dimensionless time of evolution (c), for $R_0 = 0.75$ mm in distilled water at 293 K: $Fr = 13.21$, $Re = 232.9$ (velocity and time scales: $\bar{v} = 0.312 \text{ m s}^{-1}$, $\bar{t} = 2.406$ ms; range of the total number of grid nodes: $n = 37–79$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Interface collapsing evolution starting from an equilibrium position, dimensionless velocity of the central node, and $\delta E$ versus the dimensionless time of evolution, for an air/distilled water couple at 293 K: (a, b, c) $R_0 = 2$ mm, $Fr = 1.857$, $Re = 380.3$, $\bar{v} = 0.191 \text{ m s}^{-1}$, $\bar{t} = 10.47$ ms, $n = 39–81$, and (d, e, f) $R_0 = 3$ mm, $Fr = 0.826$, $Re = 466.8$, $\bar{v} = 0.156 \text{ m s}^{-1}$, $\bar{t} = 19.25$ ms, $n = 34–54$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. First jet drop dimensionless radius $R_d^*$ against the dimensional equivalent bubble radius $R_0$: (a) air/water couple; and, (b) air/ethanol couple, both at 293 K.</figcaption></figure>

experimental (mean) values. However, the dimension of the first released jet drop depends more on the liquid properties, than on the parent bubble size [15]. When comparing bubbles of the same size, the dimensionless values of the jet drop radius are greater in water than in liquids with lower surface tension (e.g., see the zoomed image of the first drop shape just at the jet rupture moment, for a 0.75 mm radius air bubble bursting in distilled water (Fig. 6(a)), and in ethanol (Fig. 6(b)), both at 293 K, with the surface nodes plotted with *x* marks).

Boulton-Stone and Blake [8] found that there are no distinguishable differences between computations including or not viscous effects. Their boundary-layer calculations show a slightly small dissipative effect of the viscosity only after the jet formation. In order to depict differences between the viscous and the inviscid case, for the same parent bubble size and the same liquid, bursting computations are started from an equilibrium position, with and without the inclusion of viscous effects. We present the interface evolution and the central node velocity profile for a 1.98 mm radius air bubble in 85–15 wt% glycerine-water solution at 293 K, in the inviscid case (Fig. 7(a)), and in the viscous case (Fig. 7(b)). The above parent bubble radius is the critical one in the viscous case. For the same air/liquid couple, we found that:

(i) the dimensionless first jet drop radius and the evolution time till reaching the maximum jet height are slightly greater in the viscous case; and

(ii) the maximum jet speed is almost half the one in the inviscid case.

The temperature influence is examined for the same parent bubble sizes, bursting in the same liquid, over the range of 278–303 K. By increasing the temperature, so by increasing the Reynolds number at an almost constant Froude number, we found that the jet drop size is slightly increased till 293 K (with less than 5% over a change of 5 K in the liquid temperature), and becomes almost constant over the range 293–303 K. The maximum jet speed decreases (with 10% on the whole range of temperatures). The evolution time and the maximum jet height are slightly increased. To conclude, it has been found that geometric and kinetic parameters of the liquid jet and the first ejected drop are influenced more by the surface tension than by the other physical properties (density or viscosity).
<figure><img src="image_7.png" /><figcaption>Fig. 6. First jet drop shape just at the jet rupture moment (zoomed image), for $R_0 = 0.75$ mm: (a) air/distilled water couple, $Fr = 13.21$; $Re = 232.9$, $n = 37–79$, and (b) air/ethanol couple, $Fr = 5.177$, $Re = 97.83$, $n = 39–91$, both at 293 K.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Interface collapsing evolution starting from an equilibrium position, and the central node velocity profile, for $R_0 = 1.98$ mm in 85–15 wt% glycerine-water solution at 293 K ($\bar{v} = 0.165$ m s$^{-1}$, $\bar{t} = 12.00$ ms): (a) inviscid case, $Fr = 1.400$, $n = 39–80$, and (b) viscous case, $Fr = 1.400$, $Re = 3.595$, $n = 38–99$.</figcaption></figure>
## 5.2. Critical bubble radius

As aforementioned, for each air/liquid couple, there is a critical bubble radius $R_{cr}$ for which only a single large jet drop is formed. Above this limit, the liquid jet, thick and short, decays without releasing any jet drop. In order to examine the individual effect of each liquid physical property on the critical bubble radius, our computed values together with the existing experimental values of $R_{cr}$ [4,5,7] are presented in Table 1. Due to the difficulty in measuring accurately the bubble size, the experimental critical bubble radius lies in a rather large range. The minima, maximum likelihood (mean), and maxima values of this parameter are thus presented.

The density effect, rarely investigated, shows a non-linear variation of the critical bubble radius. For liquids lighter than the acetic acid, $R_{cr}$ increases with increasing density, while for liquids heavier than that, $R_{cr}$ decreases with increasing density. This result is in agreement with two natural tendencies [7]:

(i) when increasing density, the inertia of the rising jet becomes larger than the one corresponding to a lighter liquid, thus facilitating the jet rupture;

(ii) larger densities though tend to slowdown the interface motion, the whole collapsing process being decelerated.

The critical bubble size clearly depends on viscous effects: for the same air/liquid couple, $R_{cr}$ is larger in the inviscid case, than in the viscous case. More precisely, it can be shown that it tends to decrease with increasing liquid viscosity. The viscosity slows down the interface motion and delays the collapse. Thus $R_{cr}$ is slightly increased by increasing the liquid temperature. The question arises to know whether there is a limiting viscosity above that $R_{cr}$ does not exists. A very high viscosity may indeed damp any break up phenomenon [26]. Our model is unable to give a firm answer because viscous effects are partially considered in it. More generally, high viscous effects must be depicted on an experimental basis that is lacking at present.

The interface deformation depends strongly on the surface tension. The critical bubble radius increases with increasing surface tension ($R_{cr}$ increases very slowly for $\sigma \in [0.04, 0.06] \text{ N m}^{-1}$, while it increases sharply outside this range).

For the liquids in Table 1, the computed and mean experimental critical bubble radii are plotted against the corresponding Froude number (Fig. 8(a)), and Reynolds number (Fig. 8(b)). The critical radius increases with decreasing Froude number, while

<figure><img src="image_9.png" /><figcaption>Fig. 8. Dimensional critical bubble radius $R_{cr}$ against (a) Froude number, and (b) Reynolds number.</figcaption></figure>
it increases with increasing Reynolds number [15]. Numerical and experimental points exhibit the same trends. A systematic overvaluation appears for the computed $R_{cr}$ values.

### 5.3. Non-equilibrium initial configuration

Finally, we considered the bursting process started from a non-equilibrium position, for a viscous fluid. It was shown that such initial conditions yield to different numerical results [8,10] than those obtained by starting from an equilibrium position. In order to make direct comparison with the above results, computations are performed from a spherical parent bubble, fixed underneath a flat free surface (the liquid film was removed artificially and the resulting connected interface was slightly “smoothed” at the crater by repositioning the surface nodes). It was found that the interface evolution depends strongly on the crater radius: the first jet drop radius, and the evolution time increase with increasing crater radius. The maximum value of the central node velocity decreases sharply with increasing crater radius. These non-equilibrium conditions give critical bubble radius larger than those obtained from an equilibrium position, or provided by experimental data. Choosing the “optimum” crater radius, leading to the “true” first jet drop radius, is difficult and somewhat artificial, though the mechanical energy balance is accurate. It seems that the crater radius approaching the one issued from the equilibrium position leads to a first jet drop size that fits experimental data.

In Fig. 9 we present the interface evolution and the central node velocity profile, for a spherical parent bubble with $R_0 = 4$ mm, underneath a flat free surface of distilled water at 288 K, for 3 different dimensionless crater radii: $R_C^* \in \{0.26, 0.64, 0.91\}$. It must be pointed out that for small crater radius (Fig. 9(a)), the peak of the jet velocity is too important (about $10 \text{ m s}^{-1}$) to be realistic, and the mechanical energy balance deviation $|\delta E|$ reaches 7.64%. Obviously, numerical instabilities develop for punctured parent bubble. For $R_C^* = 0.64$, which approaches the crater radius issued from the equilibrium position, the numerical results become accurate, the deviation $|\delta E|$ does not exceed 2.04%, while the dimensional peak jet velocity is of $1.73 \text{ m s}^{-1}$ (Fig. 9(b)). For large crater radius (Fig. 9(c)), the deviation $|\delta E|$ is low (less than 1.41%), and the peak jet velocity diminishes to $1 \text{ m s}^{-1}$; the jet drop size is not acceptable, since it approaches half the one of the parent bubble size ($R_d^* = 0.415$).

<figure><img src="image_10.png" /><figcaption>Fig. 9. Interface collapsing evolution starting from a non-equilibrium position, and the central node velocity profile, for $R_0 = 4$ mm in distilled water at 288 K, $Fr = 0.464$, $Re = 474$, $\bar{v} = 0.135 \text{ m s}^{-1}$, $\bar{t} = 29.66$ ms, for: (a) $R_C^* = 0.26$ ($R_d^* = 0.142$, $n = 38–88$), (b) $R_C^* = 0.64$ ($R_d^* = 0.289$, $n = 39–78$), and (c) $R_C^* = 0.91$ ($R_d^* = 0.415$, $n = 38–116$).</figcaption></figure>
Sussman and Smereka [10] computed large bubble sizes ($R_0 = 4$, and 5 mm), in distilled water at 288 K, by starting with a submerged spherical bubble released just below a flat free surface. For a dimensionless crater radius $R_C^* = 0.64$, our computed dimensionless first jet drop radius ($R_d^* = 0.289$ for $R_0 = 4$ mm, and $R_d^* = 0.313$ for $R_0 = 5$ mm), and the maximum jet height fit Sussman and Smereka results. Our computed peak jet velocity is almost the same, as well as the associated moment. The time elapsed till the first jet drop ejection is slightly smaller than the one of Sussman and Smereka. It must be noted that the comparison was made by taking into account the Sussman and Smereka frames between the interface position showing a dimensionless crater radius of about 0.64, and the interface position showing a complete jet necking (first jet drop formed, but not released).

So starting computations from a spherical bubble underneath a flat free surface is ambiguous since it is very sensitive to the crater radius, which is arbitrary. It also appears not realistic since the typical computations in Fig. 9 show that jet drops are released. The critical bubble radius issued from this initial condition is thus greater than the experimental-one [7], and equally than the one computed by starting from an equilibrium position. Note that for an air/distilled water couple at 288 K, we found that jet drops are released even for $R_0 = 5$ mm, as Sussman and Smereka [10] did. Even this result is questionable since in this case where bubbles are large, an initially spherical bubble is not a realistic geometric assumption.

## 6. Conclusion

To study the jet drop formation, the final sequence of a bursting bubble process was modelled for several air/liquid couples, using a potential flow of a pure viscous fluid. We focused on collapsing bubble process initiated from an equilibrium position. The numerical simulation of the bubble cavity collapse was made through a boundary element method with an explicit second-order time-evolution scheme, by taking into account capillary and viscous effects. The results mainly point on the first jet drop geometric and kinetic characteristics, and the critical bubble size upon which no jet drop forms.

To prevent numerical instabilities, the progression time step was obtained from a stability criterion linked to the gravity-capillary dispersion equation. The precision of the numerical code was checked through the global mechanical energy balance, which is accurate until the jet formation, and lies in an admissible range during the jet ascension. In the final part of the evolution, interface shapes become highly distorted and a peak of dissipation of mechanical energy due to viscosity occurs. Very small time and space steps have to be introduced, leading to some numerical instabilities.

The comparative analysis made in Section 5 of this paper, for bubbles that burst from an equilibrium position, proves that our computed values are in agreement with available experimental data. It also shows a slightly accelerated behaviour of the phenomena with respect to the real one, especially in the final stage of the bursting process. The observed discrepancies may originate from the partial consideration of the viscous forces, which do not damp enough inertial effects. Thus, only in this final stage of the process, the complete Navier-Stokes equations are worth being used, but their numerical treatment remains far from being easier.

As far as the bursting bubble simulation starts from an equilibrium position, the following main results were found:

(i) the collapsing scenario remains constant irrespective of the considered liquids;

(ii) the dimensionless value of the first released jet drop radius does not vary so much upon the parent bubble radius for the same air/liquid couple, but is closely related to the liquid type. In particular, this radius increases with increasing surface tension;

(iii) the critical bubble radius increases with decreasing Froude number, while it increases with increasing Reynolds number. The computed critical bubble radius is slightly overestimated with respect to the corresponding experimental mean value;

(iv) the dimensionless first jet drop radius is slightly greater in the viscous case, and the maximum jet speed is almost half the one in the inviscid case. For the same air/liquid couple, the critical bubble size is larger in the inviscid case, than in the viscous case.

The first ejected jet drop radius associated to a bursting process started from a non-equilibrium initial position (a spherical parent bubble fixed underneath a flat free surface, with a crater radius value around the one obtained from the equilibrium position), fit the numerical results of Sussman and Smereka [10] that started with a submerged spherical bubble released just below a flat free surface (dynamic initial conditions for the collapsing process). Our computations started from a static non-equilibrium initial position show that the interface evolution depends strongly on the crater radius. It seems that the crater radius approaching the one corresponding to the equilibrium position leads to a jet drop size that fits experimental data. Small or large values of the crater radius should be avoided. When the crater radius increases, the jet drop radius also increases, and the maximum value of the central node velocity decreases sharply. The critical bubble radius is greater than the correspondent experimental value, and greater than the numerical value obtained by starting from the equilibrium position.
Appendix A. Time-stepping procedure formulae

Before depicting the expressions of the Lagrangian derivatives of the velocity potential $\phi^*(t^*)$, of the radial coordinate $r^*(t^*)$ and of the axial coordinate $z^*(t^*)$, that appear in the second-order limited Taylor series expansions (3.2), we point out that there are distinct forms of relations on the axis of symmetry, and outside. The singularities issued at $r^* = 0$ are removed by taking into account that on the $Oz$ axis, the axisymmetric curvature $(1/R_2^* = (\sin\beta)/r^*)$ equals the planar curvature $(1/R_1^*)$, which depends on the variation of the azimuthal angle $\beta$ along the curvilinear abscissa $s^*)$. The total curvature is written $H^* = (\partial\beta/\partial s^* + \sin\beta/r^*)$ if $r^* \neq 0$, and becomes $H^* = 2(\partial\beta/\partial s^*)$ if $r^* = 0$.

In the expression (3.3), we compute the first two terms of the second-order Lagrangian derivatives of the radial and axial coordinates, $D^2r^*/Dt^{*2}$ and $D^2z^*/Dt^{*2}$, through the analytical representations given by Machane and Canot [20].

Outside the $Oz$ axis, the material derivative of the total curvature is written:

$$ \frac{DH^*}{Dt^*} = \left( \left( \frac{\partial \beta}{\partial s^*} \right)^2 + \left( \frac{\sin \beta}{r^*} \right)^2 \right) \frac{\partial \phi^*}{\partial n^*} + \frac{\cos \beta}{r^*} \frac{\partial}{\partial s^*} \left( \frac{\partial \phi^*}{\partial n^*} \right) + \frac{\partial}{\partial s^*} \left( \frac{\partial \beta}{\partial s^*} + \frac{\sin \beta}{r^*} \right) \frac{\partial \phi^*}{\partial s^*} + \frac{\partial^2}{\partial s^{*2}} \left( \frac{\partial \phi^*}{\partial n^*} \right), \quad (A.1) $$

while for $r^* = 0$ it becomes:

$$ \frac{DH^*}{Dt^*} = 2 \left( \frac{\partial \beta}{\partial s^*} \right)^2 \frac{\partial \phi^*}{\partial n^*} + 2 \frac{\partial^2}{\partial s^{*2}} \left( \frac{\partial \phi^*}{\partial n^*} \right). \quad (A.2) $$

When $r^* \neq 0$, the tangential derivative of the axisymmetric curvature is determined by:

$$ \frac{\partial}{\partial s^*} \left( \frac{\sin \beta}{r^*} \right) = \left( \frac{\partial \beta}{\partial s^*} - \frac{\sin \beta}{r^*} \right) \frac{\cos \beta}{r^*}. \quad (A.3) $$

The first and second tangential derivatives of the variables have been obtained directly from a numerical differentiation along the curvilinear abscissa $s^*$ (i.e. $\partial\phi^*/\partial s^*$, and $\partial^2\phi^*/\partial s^{*2}$ for the velocity potential).

Before expressing the last term of Eq. (3.3), we must define the normal second derivative of the velocity potential, $\partial^2\phi^*/\partial n^{*2}$. This quantity is also used in the Bernoulli's equation (2.7), as well as in the mechanical energy balance (4.2). Outside the $Oz$ axis, the corresponding form of $\partial^2\phi^*/\partial n^{*2}$ is:

$$ \frac{\partial^2 \phi^*}{\partial n^{*2}} = -\frac{\partial^2 \phi^*}{\partial s^{*2}} - \frac{\cos \beta}{r^*} \frac{\partial \phi^*}{\partial s^*} + \left( \frac{\partial \beta}{\partial s^*} + \frac{\sin \beta}{r^*} \right) \frac{\partial \phi^*}{\partial n^*}, \quad r^* \neq 0, \quad (A.4) $$

while on the axis of symmetry it becomes:

$$ \frac{\partial^2 \phi^*}{\partial n^{*2}} = -2 \frac{\partial^2 \phi^*}{\partial s^{*2}} + 2 \frac{\partial \beta}{\partial s^*} \frac{\partial \phi^*}{\partial n^*}, \quad r^* = 0. \quad (A.5) $$

From this, by using the Maple software for symbolic computation, the slightly complicated expression of the corresponding material derivative was determined:

$$ \begin{aligned} \frac{D}{Dt^*} \left( \frac{\partial^2 \phi^*}{\partial n^{*2}} \right) = & -2 \left( \frac{\partial^2 \phi^*}{\partial n^* \partial s^*} \right)^2 - \frac{\partial \phi^*}{\partial s^*} \frac{\partial^3 \phi^*}{\partial s^* \partial s^*} + \left( \frac{\partial \phi^*}{\partial s^*} \right)^2 \left( \frac{\cos \beta}{r^*} \right)^2 - \frac{\partial \phi^*}{\partial s^*} \frac{\partial^2 \phi^*}{\partial s^{*2}} \frac{\cos \beta}{r^*} - \frac{\partial^3 \phi^*}{\partial n^* \partial s^{*2}} \frac{\partial \phi^*}{\partial n^*} \\ & - \frac{\partial^2}{\partial s^{*2}} \left( \frac{\partial \phi^*}{\partial t^*} \right) - \frac{\partial \beta}{\partial s^*} \frac{\partial^2 \phi^*}{\partial n^* \partial s^*} \frac{\partial \phi^*}{\partial s^*} - 3 \frac{\partial \beta}{\partial s^*} \frac{\partial^2 \phi^*}{\partial s^{*2}} \frac{\partial \phi^*}{\partial n^*} + \left( \frac{\partial \phi^*}{\partial s^*} \right)^2 \frac{\partial \beta}{\partial s^*} \frac{\sin \beta}{r^*} + \frac{\partial}{\partial n^*} \left( \frac{\partial \phi^*}{\partial t^*} \right) \frac{\sin \beta}{r^*} \\ & + 2 \left( \frac{\partial \beta}{\partial s^*} \right)^2 \left( \frac{\partial \phi^*}{\partial n^*} \right)^2 + 2 \left( \frac{\sin \beta}{r^*} \right)^2 \left( \frac{\partial \phi^*}{\partial n^*} \right)^2 + \frac{\partial}{\partial n^*} \left( \frac{\partial \phi^*}{\partial t^*} \right) \frac{\partial \beta}{\partial s^*} - \frac{\partial \phi^*}{\partial n^*} \frac{\cos \beta}{r^*} \frac{\partial^2 \phi^*}{\partial n^* \partial s^*} \\ & - \frac{\partial^2 \phi^*}{\partial s^{*2}} \frac{\partial \phi^*}{\partial n^*} \frac{\sin \beta}{r^*} + 2 \left( \frac{\partial \phi^*}{\partial n^*} \right)^2 \frac{\partial \beta}{\partial s^*} \frac{\sin \beta}{r^*} + \frac{\partial^2 \phi^*}{\partial n^* \partial s^*} \frac{\partial \phi^*}{\partial s^*} \frac{\sin \beta}{r^*} - \frac{\partial \phi^*}{\partial s^*} \frac{\cos \beta}{r^*} \frac{\partial \phi^*}{\partial n^*} \frac{\partial \beta}{\partial s^*} \\ & - \frac{\cos \beta}{r^*} \frac{\partial}{\partial s^*} \left( \frac{\partial \phi^*}{\partial t^*} \right) - 3 \frac{\partial \phi^*}{\partial s^*} \frac{\cos \beta}{r^*} \frac{\partial \phi^*}{\partial n^*} \frac{\sin \beta}{r^*}, \end{aligned} \quad (A.6) $$

for the boundary nodes in which $r^* \neq 0$, and

$$ \begin{aligned} \frac{D}{Dt^*} \left( \frac{\partial^2 \phi^*}{\partial n^{*2}} \right) = & -2 \frac{\partial^2}{\partial s^{*2}} \left( \frac{\partial \phi^*}{\partial t^*} \right) - 4 \left( \frac{\partial^2 \phi^*}{\partial n^* \partial s^*} \right)^2 - 2 \frac{\partial \phi^*}{\partial n^*} \frac{\partial^3 \phi^*}{\partial n^* \partial s^{*2}} - 2 \frac{\partial \beta}{\partial s^*} \frac{\partial^2 \phi^*}{\partial n^* \partial s^*} \frac{\partial \phi^*}{\partial s^*} - 8 \frac{\partial \beta}{\partial s^*} \frac{\partial \phi^*}{\partial n^*} \frac{\partial^2 \phi^*}{\partial s^{*2}} \\ & + 6 \left( \frac{\partial \beta}{\partial s^*} \right)^2 \left( \frac{\partial \phi^*}{\partial n^*} \right)^2 + 2 \frac{\partial \beta}{\partial s^*} \frac{\partial}{\partial n^*} \left( \frac{\partial \phi^*}{\partial t^*} \right) - 2 \frac{\partial \phi^*}{\partial s^*} \frac{\partial^3 \phi^*}{\partial s^{*3}}, \end{aligned} \quad (A.7) $$

for the $Oz$ axis, where $r^* = 0$.
<figure><img src="image_11.png" /><figcaption>Fig. A1. Total number of grid nodes *n* versus the dimensionless time *t*\*, for the bursting bubble evolution presented in Fig. 3: *R*₀ = 0.75 mm; air/distilled water couple at 293 K; *Fr* = 13.21, *Re* = 232.9.</figcaption></figure>

## Appendix B. Spatial discretization

We examined the numerical code sensitivity at the extension of the solid boundary $Σ_w$, and we established that the external radial limit should be greater than about 4 times the bubble equivalent radius.

During a typical computation, the dimensionless length of a boundary element $\Delta s^*$ is restricted to vary in the range [5 × 10⁻⁵, 0.4]. The ratio between the lengths of two consecutive boundary elements is less than 2. A fixed number of 11 equidistant grid nodes are distributed on each solid surface: lateral wall and base of the cylindrical pool. Initially, about 35 grid nodes are unevenly distributed on the interface $Σ_i$. The regridding (see Section 3) is performed on $Σ_i$ during the whole computation. The total number of grid nodes $n$ distributed on $Σ = Σ_i \cup Σ_w$ varies between 38 and 116. In Fig. A1, the total number of grid nodes is plotted against the dimensionless evolution time, corresponding to the interface collapse started from equilibrium, for a 0.75 mm radius air bubble, in distilled water at 293 K (see for reference Fig. 3).

We established that a grid refinement should be avoided, because it leads to a sharp increase of the peak jet velocity, and an increase of the relative mechanical energy balance deviation. As expected, the computational effort (run time) is increased. The discrepancies are obviously due to numerical instabilities. By decreasing the minimum value of the spatial step $\Delta s^*_{min}$, the time step is decreased as a result of (4.1): $(\Delta t^*)^2 \approx (\Delta s^*_{min})^3$. The error associated with the BEM solver is of the order of $O(\Delta s^3)$. Moreover, in the temporal scheme (3.2), the first-order Lagrangian derivative induces an error of the order of $O(\Delta s^3)$, while the second-order derivative induces the same error multiplied by the matrix coefficient condition number [20], because the latter involves two consecutive Laplace problems: (2.4) and (3.4).

## References

[1] F.J. Resch, J.S. Darrozès, G.M. Afeti, Marine liquid aerosol production from bursting of air bubbles, J. Geophys. Res. 91 (1986) 1019–1029.

[2] J. Wu, Film drops produced by air bubbles bursting at the surface of sea water, J. Geophys. Res. 99 (C8) (1994) 16403–16407.

[3] C.F. Kientzler, A.B. Arons, D.C. Blanchard, A.H. Woodcock, Photographic investigation of the projection of droplets by bubbles bursting at a water surface, Tellus 6 (1) (1954) 1–7.
[4] F.H. Garner, S.R.M. Ellis, J.A. Lacey, The size distribution and entrainment of droplets, Trans. Inst. Chem. Engrg. 32 (1954) 222-235.

[5] F.J. Resch, G.M. Afeti, Gas-liquid conversion at a free interface and its environmental implications, in: J.F. Keffer, R.K. Shah, E.N. Ganic (Eds.), *Experimental Heat Transfer*, Fluid Mechanics and Thermodynamics, Elsevier, 1991, pp. 1066-1076.

[6] G.M. Afeti, F.J. Resch, From breaking waves to breaking bubbles: Production characteristics of the marine droplet aerosol, Trends Geophys. Res. 1 (1992) 113-125.

[7] T. Suzuki, K. Mitachi, Experimental observation of the droplet ejection due to gas bubble bursting at gas liquid interface, in: Proc. 2nd Int. Conf. on Multiphase Flow, Kyoto, Japan, 1995, pp. IP2.3-IP2.9.

[8] J.M. Boulton-Stone, J.R. Blake, Gas bubbles bursting at a free surface, J. Fluid Mech. 254 (1993) 437-466.

[9] J.M. Boulton-Stone, The effect of surfactant on bursting gas bubbles, J. Fluid Mech. 302 (1995) 231-257.

[10] M. Sussman, P. Smereka, Axisymmetric free boundary problems, J. Fluid Mech. 341 (1997) 269-294.

[11] D.D. Joseph, T.Y. Liao, H.H. Hu, Drag and moment in viscous potential flow, Eur. J. Mech. B 12 (1) (1993) 97-106.

[12] I.B. Ivanov, P.A. Kralchevsky, A.D. Nikolov, Film and line tension effects on the attachment of particles to an interface: I. Conditions for mechanical equilibrium of fluid and solid particles at a fluid interface, J. Colloid Interf. Sci. 112 (1) (1986-a) 97-107.

[13] I.B. Ivanov, P.A. Kralchevsky, A.D. Nikolov, Film and line tension effects on the attachment of particles to an interface: II. Shapes of the bubble (drop) and the external meniscus, J. Colloid Interf. Sci. 112 (1) (1986-b) 108-121.

[14] I.B. Ivanov, P.A. Kralchevsky, A.D. Nikolov, Film and line tension effects on the attachment of particles to an interface: IV. Experimental studies with bubbles in solutions of Dodecyl Sodium Sulfate, J. Colloid Interf. Sci. 112 (1) (1986-c) 132-143.

[15] S.-C. Georgescu, Évolution d'une bulle: Formation à partir d'un orifice et éclatement à la traversée d'une surface libre, PhD Thesis, Institut National Polytechnique de Grenoble, Grenoble, 1999.

[16] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd edn., Cambridge University Press, Cambridge, 1992.

[17] E. Canot, J.-L. Achard, An overview of boundary integral formulations for potential flows in fluid-fluid systems, Arch. Mech. 43 (4) (1991) 453-498, Warszawa.

[18] R. Machane, Contributions de la Méthode Intégrale aux Frontières au suivi d'interfaces, PhD Thesis, Université Joseph Fourier, Grenoble 1997.

[19] J.W. Dold, D.H. Peregrine, An efficient boundary-integral method for steep unsteady water waves, in: K.N. Morton, M.J. Baines (Eds.), *Numerical Methods for Fluid Dynamics II*, Oxford University Press, Oxford, 1986, pp. 671-679.

[20] R. Machane, E. Canot, High-order schemes in Boundary Element Methods for transient non-linear free surface problems, Internat. J. Numer. Methods Fluids 24 (1997) 1049-1072.

[21] S.-C. Georgescu, J.-L. Achard, E. Canot, Bubble cavity collapse and liquid jet formation simulation in bursting gas bubble process, in: Book of Abstracts – EUROMECH Colloquium 376, Waves in Two-Phase Flows, Istanbul, 1998, pp. 25-27.

[22] E. Canot, Stability criteria for capillary/gravity free-surface waves in BEM simulations of viscous potential flows, in: M.H. Aliabadi (Ed.), Proc. of the Int. Conf. on Boundary Element Techniques, Queen Mary & Westfield College, London, 1999, pp. 395-404.

[23] M. Wang, A.W. Troesch, Numerical stability analysis for free surface flows, Internat. J. Numer. Methods Fluids 24 (1997) 893-912.

[24] E. Canot, S.-C. Georgescu, Bilan d'énergie mécanique pour le suivi d'interface: Modèle potentiel avec prise en compte des contraintes visqueuses normales, in: Actes du 13ème Congrès Français de Mécanique, Poitiers-Futuroscope, France, 3, 1997, pp. 23-26.

[25] D.C. Blanchard, The electrification of the atmosphere by particles from bubbles in the sea, Prog. Oceanogr. 1 (1963) 71-102.

[26] K. Takayama, S. Hayakawa, Shock wave/bubble interaction in highly viscous liquids, in: Book of Abstracts – EUROMECH Colloquium 376, Waves in Two-Phase Flows, Istanbul, 1998, p. 94.

[27] N.A. Lange, Handbook of Chemistry, 9th edn., McGraw-Hill, New-York, 1956.