# Bursting air bubble at a free surface: reggridding influence on the interface evolution

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

► To cite this version:

Edouard Canot, Sanda-Carmen Georgescu, Jean-Luc Achard. Bursting air bubble at a free surface: reggridding influence on the interface evolution. Workshop on Numerical Simulation for Fluid Mechanics and Magnetic liquids, 2001, Timisoara, Romania. pp.121-132. hal-00197477

HAL Id: hal-00197477

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

Submitted on 3 Apr 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
# BURSTING AIR BUBBLE AT A FREE SURFACE: REGRIDDING INFLUENCE ON THE INTERFACE EVOLUTION

Édouard CANOT (1), Sanda-Carmen GEORGESCU (2), Jean-Luc ACHARD (3)

**Abstract:** The numerical simulation of an air bubble bursting at a free surface was made in a potential flow of a viscous fluid assumption, through the Boundary Element Method. The evaluation of the numerical code precision was allowed trough a global mechanical energy balance, expressed only in surface integrals terms. The bubble cavity collapse behaviour, finalised by the jet formation and its successive break up, depends on the regridding performed during the whole computation.

## 1. INTRODUCTION

When a rising air bubble reaches a free surface, the liquid film formed between the bubble cap and the adjacent gaseous phase is rapidly drained, and finally atomised into tiny “film droplets”. The gaping bubble cavity collapses then: a surface wave converges axially and generates an upward high-speed liquid jet that splits up into several “jet drops”. The droplets formation following air bubbles bursting at a free liquid surface plays an important role in mass transfer between phases in stratified two-phase flows.

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 simulation of the cavity collapsing process for an air bubble bursting at a free pure water surface was made by Boulton-Stone and Blake [1]; it was improved by Boulton-Stone [2] for a contaminated surface. Both studies follow an inviscid fluid model through a Boundary Element Method completed by boundary-layer effect consideration. The bursting process for large air bubbles in water was modelled by Sussman and Smereka [3], by solving the Navier-Stokes equations in both fluids on a fixed grid, through the level set method; because of the important computational effort, the bubble equivalent radius was greater than 4 mm. There are no comparisons with the previous numerical studies, which cover a range of bubble radius from 0.5 to 3 mm.

We simulated the interface evolution during the cavity collapsing process of a bursting air bubble, by using a purely irrotational flow model, where viscous effects have been incorporated to some extent as it was allowed for potential flows of fluids with constant viscosity [4]. For different liquids, we reproduced most observed experimental data, namely the first jet drop geometric and kinetic characteristics, and the critical bubble diameter upon which the liquid jet decays without releasing any jet drop [5], [6]. Emphasis has been placed on numerical accuracy: a Boundary Element Method has been selected with a second-order time-evolution scheme. In the absence of analytical tests, the evaluation of numerical code precision was allowed trough a global mechanical energy balance expressed only in surface integrals terms. In this paper we point on the regridding requirements allowing the simulation of the interface evolution during the whole collapsing process.

## 2. PROBLEM STATEMENT

We focus exclusively on the “break up stage” of a bursting bubble process: after the liquid film disintegration, the bubble cavity collapses, being finally solved in an unstable liquid jet that splits up into several jet drops. Modelling only this stage stands with the difficulty to express realistic initial conditions.

We consider a bounded axisymmetric liquid domain $\Omega_L$, extending largely over the bubble size, in order to obtain a negligible influence of phenomena at the periphery. The interfaces shape corresponds to a stationary bubble entrapped beneath a free surface [7]. The resulting bubble cavity surface $\Sigma_b$, film equivalent interface $\Sigma_f$, and external meniscus $\Sigma_m$ of the free surface are in equilibrium with the surrounding fluid, being joined along the crater line (figure 1). At any point of the interface between two phases i and j, there is a positive pressure jump

$$ (1) \qquad p_i - p_j = \sigma_{ij} \left( \frac{1}{R_1} + \frac{1}{R_2} \right), $$

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; $\sigma_{ij}$ equals the surface tension $\sigma$ on $\Sigma_b$ and $\Sigma_m$, and it equals the film tension $\gamma$ on $\Sigma_f$. The adjacent gas pressure $p_0$ is assumed to be constant. The gas pressure $p_G$ inside the bubble and the liquid pressure $p_L$ depend on $z$. The system derived from (1), completed with the normal force balance at the triple contact point C, and with the geometric relation issued for the radius of curvature of the film assumed to be a part of a sphere, is nondimensionalized by taking the bubble equivalent radius $R_0$ as length scale (the radius of a sphere with the same volume), and is solved through a globally convergent Newton Method [5], by considering the Froude number
$Fr = \sigma/(\rho_L g R_0^2)$ as control parameter. Once we get the equilibrium position of the three joined interfaces, the initial geometric conditions requested by the bursting simulation are established by the connected interface $\Sigma_i = \Sigma_b \cup \Sigma_m$, the film being artificially removed.

<figure><img src="image_2.png" /><figcaption>Figure 1. Resulting three joined interfaces (meridian plane $rOz$)</figcaption></figure>

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, being fourth times greater than the bubble radius. We take a unit outward normal $\mathbf{n}$ on $\Sigma$. A potential flow assumption is appropriate because of the impulsive character of the phenomenon. We will take into account the viscous effects only by the normal viscous stress at the interface: $2\mu_L (\partial v_n / \partial n)$, where $\mu_L$ represents the dynamic viscosity of the liquid, and $v_n$ is the normal component of the velocity.

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 length scale $R_0$ and velocity as time scale [1]. The dimensionless equations governing the potential fluid flow are classically the Laplace equation for the velocity potential $\phi^*$ and Euler's equation. The viscous effects are considered through boundary conditions, namely the normal momentum balance at any point of the interface $\Sigma_i$. Combining those last two equations to reduce the pressure terms, we obtain the Bernoulli's equation:

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

where, according to the choice of scales, the Weber number is $We = 1$, and the Reynolds number is $Re = (\rho_L \sigma R_0)^{1/2} / \mu_L$. We add also a non-penetrability condition at the interface $\Sigma_i$. The expressions of the normal gradient of normal velocity $\partial^2\phi^* / \partial n^{*2}$ are defined in [5]. The initial kinetic conditions correspond to a velocity field assumed to be zero [1], $\Sigma_i$ being frozen till the beginning of the cavity collapse.

### 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 $\Delta t^*$. At a fixed instant $t^*$, we solve the Laplace equation, $\nabla^2\phi^* = 0$, to obtain the velocity potential values $\phi^*(t^*)$, hence, the corresponding normal component, and tangential component of the velocity. (ii) A time-stepping scheme allows the connection of two successive steps to determine the new potential values and interface position at the following instant $(t^* + \Delta t^*)$.

The Boundary Element Method (BEM) is well adapted to the first type of calculation: it replaces 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 calculations [8], [9], the velocity field is generated by source and normal doublet type singularities spread over $\Sigma$. The integral equation on boundary $\Sigma$ is written

$$ (3) \qquad -\int_{\Sigma} \phi^*(\mathbf{P}) \frac{\partial(1/|\mathbf{MP}|)}{\partial n_{\mathbf{P}}^*} d\Sigma_{\mathbf{P}} + \int_{\Sigma} (1/|\mathbf{MP}|) \frac{\partial\phi^*(\mathbf{P})}{\partial n_{\mathbf{P}}^*} d\Sigma_{\mathbf{P}} = 2\pi\phi^*(\mathbf{M}), \quad \mathbf{M} \in \Sigma, $$
where $|\mathbf{MP}|$ is the distance between the observation point $\mathbf{M}$ and the singularity point $\mathbf{P}$. We define the 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 immobile 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 $M = M(r^*(t^*), z^*(t^*)$, unevenly redistributed on the boundary $\Sigma(t^*)$ at each time step, with respect to some criteria like the adaptation at surface gradients [8]: that leads to a concentration of nodes at places where the interface curvature is important, or where two portions of the interface approach one another (see Appendix). On each boundary element, cubic splines define the geometric variables, and cubic Hermite polynomials approximate the field variables. Following an explicit numerical scheme [10] improved for the
capillary [11], and viscous effects [12], the time progression is made through second-order limited Taylor series expansions of a function $\chi(t^*)$

$$ (4) \qquad \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}), $$

where $\chi(t^*)$ denotes $\phi^*(t^*)$, $r^*(t^*)$ and $z^*(t^*)$ respectively. We mention that the matrix of the BEM calculations following (3) depends only on geometrical quantities. Thus it has to be computed only once for each time step, and this explains why the temporal scheme (4) is very efficient. The order of the time stepping scheme must be consistent with the order of the BEM solver used [11]. In the case of mixed capillary-gravity waves, a cubic BEM formulation is used; according to the consistency condition, a second order time stepping procedure must be chosen [13]. The material derivative of the velocity potential is:

$$ \frac{D\phi^*/Dt^*}{D^2\phi^*/Dt^{*2}} = \frac{\partial\phi^*/\partial t^*}{\nabla(\phi^*)^2}, \text{ where the local time derivative } \frac{\partial\phi^*/\partial t^*}{} \text{ is defined by (2). The corresponding second-order Lagrangian derivative becomes:} $$

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

The normal component of the velocity $v_n^* = \partial\phi^*/\partial n^*$, and the tangential-one

$$ v_t^* = \frac{\partial\phi^*}{\partial s^*}, \text{ allow the determination of the radial component of the velocity} $$

$$ v_r^* = (-v_n^* \sin\beta + v_t^* \cos\beta), \text{ together with the axial-one } 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$, than the solid surface $\Sigma_w$, and ends on the Oz axis. The expressions of Lagrangian derivatives that appear in the right hand of (5) are depicted in [5]. All those terms need the calculations of the normal and tangential derivatives of $\partial\phi^*/\partial t^*$. This involves the computation of an extra Laplace equation: $\nabla^2(\partial\phi^*/\partial t^*) = 0$. The boundary conditions are: $\partial\phi^*/\partial t^*$ issued from (2) on $\Sigma_i$, and $\partial(\partial\phi^*/\partial t^*)/\partial n^* = 0$ on $\Sigma_w$.

The variable time step $\Delta t^*$ is selected at each time through a stability criterion derived from an eigenvalue problem based on the above temporal scheme [11], [13]. 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 [14], 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 [13]

$$ (6) \qquad \Delta t^* \le 2 \left( \frac{\pi^3}{\text{We}} \left( \Delta s_{\min}^* \right)^3 + \frac{\pi}{\text{Fr}} \Delta s_{\min}^* \right)^{-1/2}, $$

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 numerical code precision is allowed through the global mechanical energy balance [6], [13]:

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

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. The dissipation of mechanical energy due to shear viscosity depends only on the normal component of the velocity, and on its normal derivative, $\partial^2\phi^*/\partial n^{*2}$.

## 4. RESULTS

The numerical simulation of the break up stage of an isolated air bubble bursting at a free surface has been systematically made for 21 pure Newtonian liquids and different parent bubble sizes [5], starting from 0.5-mm-radius up to the critical bubble radius, in order to compare numerical and available experimental data. To exemplify, we present the interface evolution during the collapse process of an air bubble in distilled water at 20°C, for an equivalent radius of $R_0 = 0.5$ mm (figure 2). After the first jet drop ejection, the liquid jet still rises and ejects drops. Finally, because of the energy decrease, the jet (thick and short) falls down.

The velocity of the central node (on the Oz axis) increases rapidly, and reaches a maximum when the jet starts to form (that instant being linked to a viscous dissipation peak); then, the velocity magnitude decays sharply to successive levels that correspond to each jet rupture moment. After ejecting the last drop, the velocity decays slowly. The maximum value of the velocity decreases with increasing bubble radius. The first drop ejection time increases with increasing bubble size. According to the velocity variation, the mechanical energy balance is affected through the rates of the kinetic energy and dissipation of energy due to viscosity: when the jet rises, the balance (7) is no more well verified, but still preserves an admissible range of variation relative to the initial mechanical energy. This discrepancy that we have imputed above to a viscous energy dissipation peak,
is also due to numerical instabilities, which appear when steps diminish in the final part of the collapsing process.

<figure><img src="image_3.png" /><figcaption>Figure 2. Bursting bubble at a free surface: $R_0 = 0.5$ mm, $Fr = 29.7$, $Re = 190$</figcaption></figure>

By decreasing the time step, the problem becomes more and more stiff. Towards the end of the collapse, the mechanical energy balance is less verified for the small bubble case (the relative mechanical energy balance deviation $\delta E$, considered to be the ratio of the difference between left and right hand terms of (7) to the initial mechanic energy is less than 5%), and is verified for large bubbles ($\delta E < 1%$), because of the greatest velocity values in the former case. The form (7) of the mechanical energy balance was successfully used as validation test of code accuracy in transient free-boundary axisymmetric and 2D problems solved through the BEM [15].

In accordance with the regridding procedure described in the Appendix, in figures 3 and 4 we present two examples of the interface discretization, namely zoomed images of the upper part of the liquid jet, both at $Fr = 29.7$ and $Re = 190$, but for different values of the constant $e$ in (A.3). The sequences are taken prior to the moment of the first jet pinch-off and ejection of the formed “jet drop”.

<figure><img src="image_4.png" /><figcaption>Figure 3. Meshing of the upper part of the jet: $e = 0.75$, $t^* = 0.938$, $N = 77$</figcaption></figure>

One can see clearly the concentration of nodes (plotted by “+” marks) at places where the curvature is important (e.g. the inferior part of the formed drop, at the connection with the jet neck - figure 3), or where two portions of the interface get closer (e.g. the long tiny jet – figure 4). Only the computation in figure 3 can be trustworthy, the balance (7) being well verified. The computation in figure 4 cannot be validated because it leads to $\delta E$ values exceeding the admissible range; besides, the phenomenon seems to be retarded, the interface evolution behaviour being not realistic since the jet height exceeds 4-5 times the parent bubble radius, and the first jet rupture comes two times slower than in the previous case. We found that values of the constant $e$ around the unity lead to a higher concentration of nodes when approaching the axis of symmetry, which involve a faster pinch-off after the “jet drop” is formed.

## 5. 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 started 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, considering capillary and viscous effects. To prevent numerical instabilities, the progression time step was issued 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, being very well verified until the jet formation, and lying in an admissible range during the jet ascension.

<figure><img src="image_5.png" /><figcaption>Figure 4. Meshing of the upper part of the jet: $e = 2$, $t^* = 1.841$, $N = 91$</figcaption></figure>

In order to avoid numerical instabilities, a variable number of nodes are unevenly redistributed on the boundary at each time step, in accordance with the criteria imposed by the regridding procedure. 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.

# 6. APPENDIX: REGRIDDING PROCEDURE

The regridding procedure [8] generates the distribution of boundary elements on $\Sigma(t)$, in the meridian plane $(rOz)$ – we will suppress here the asterisk related to nondimensional variables. A moving surface $\Sigma(t)$ can be viewed as a family of surfaces, one for each time $t$, which can be characterized either by the equation $f(M,t) = 0$, where $f$ is smooth with respect to $t$, or by the smooth mapping valid at least on $\Sigma(t): M' = M'(r, z, t)$. We need a coordinate system $x^j = x^j(r, z, t)$, ($j = 1,2$) to give the usual parametric representation of $\Sigma(t): x'^j = x'^j(r, z, t)$, ($j = 1,2$). The basic ideas of the meshing method will not be obscured by the complex formalism of the surface representation. At a fixed time $t$, the vectorial form of the curve equation is: $\mathbf{x} = \mathbf{x}(u)$, the curve coordinate $u$ being assumed to be positively oriented. To make calculations easier, a preliminary transformation must be done from $u$ to $s$, the arc-length parameter of the curve. The local increasing transformation $\bar{u} = \bar{u}(s)$ meets two requirements: local smoothness and adaptation of the boundary elements to the variation of a function $F$. The variable $\bar{u}$ allows the numbering of the surface nodes, by attaching them an integer value: $\bar{u} = 1,2,3,...$ (the node $\bar{u} = 1$ is on the axis Oz, and the last node $\bar{u} = N$ is on the solid wall $\Sigma_w$). The vector $\Delta\mathbf{x}(\bar{u}) = \mathbf{x}(\bar{u} + \Delta\bar{u}) - \mathbf{x}(\bar{u})$ representing a boundary element is given by

$$ (A.1) \qquad \Delta\mathbf{x}(\bar{u}) = \frac{ds}{du} \Delta\bar{u} \mathbf{t} + O(\Delta\bar{u}), $$

where $\Delta\bar{u} = 1$, and $\mathbf{t} = \frac{d\mathbf{x}}{ds}$ is the unit tangent vector. The difference between two successive boundary elements is given by

$$ (A.2) \qquad \Delta\mathbf{x}(\bar{u} + 1) - \Delta\mathbf{x}(\bar{u}) = \left[ \frac{d^2 s}{d\bar{u}^2} \mathbf{t} + \left( \frac{ds}{d\bar{u}} \right)^2 \mathbf{\kappa} \mathbf{n} \right] + O(1), $$

where $\mathbf{\kappa}$ is the normal curvature and $\mathbf{n}$ is the unit normal. A mesh will be considered as locally smooth if two subsequent boundary elements have [8]: (i) two orientations differing by less than a given amount $\theta_{max}$, (ii) two component magnitude along $\mathbf{t}$ differing by less than a given maximum amount, and (iii) if the product of the strictly positive weight function $F(s)$ by the boundary element length is smaller than a constant $c$. The function $F$ reflects the fact that the nodes must be concentrated at some place: firstly, to avoid great discrepancies, the segment length must be adapted to the gradient of the velocity potential $\phi$; secondly, some self-crossing of the surface can occur if not enough points are used when different parts of the surface approach one another. We add a fourth constraint (iv) related to the proximity to the axis of symmetry: $(ds/d\bar{u})\mathbf{\kappa}_{axi}(s) \le e$, where $\mathbf{\kappa}_{axi}(s)$ is the axisymmetric curvature; it diminishes the distance between two successive nodes when approaching the Oz axis. The function
$$
(A.3) \quad \zeta(s) = \frac{ds}{d\bar{u}} = \min\left[\frac{\theta_{\max}}{|\kappa|}, \frac{c}{F(s)}, \frac{\epsilon}{\kappa_{\text{axi}}(s)}\right]
$$

verifies the constraints (i), (iii) and (iv), and it is possible to correct this function in
order to respect the second constraint. Thus the unknown function $\bar{u}(s)$ is:

$$
(A.4) \qquad \bar{u}(s) = \int_0^s \frac{dl}{\zeta(l)},
$$

which gives also the inverse function $s(\bar{u})$. However, the value $\bar{u}_{\max} = \bar{u}(s_{\max})$ is
not an integer, but we must consider a node at the final extremity of the boundary.
There are two options: (a) the first-one consists in imposing $N = \text{int}(\bar{u}_{\max}) + 1$, the
curvilinear abscissa of each node being

$$
(A.5) \quad s_j = (\bar{u})^{-1} \left[ \frac{(j-1)\bar{u}_{\max}}{N-1} \right], \quad j=1,\dots,N.
$$

Since the maximum boundary length $s_{\max}$ varies at each time $t$, the node distribution
is not continuous at each regridding. The second option, (b), consists in inverting
the function $\bar{u}(s)$ only for the first ($N-1$) nodes, which gives

$$
(A.6) \quad s_j = (\bar{u})^{-1} (j-1), \quad j=1, \dots, N-1.
$$

The last node is positioned then exactly at the end of the curve, at $s_{\max}$. That last
option leads to a continuous regridding for the first ($N-1$) nodes, and it is preferred
to the first one, because it gives 20% improvement on the global mechanical energy
balance, the interpolation errors being diminished (at each iteration, the nodes are
relatively very close to their precedent position).

7. REFERENCES

1. Boulton-Stone, J. M., Blake, J. R., 1993 "Gas bubbles bursting at a free surface", J. Fluid Mech., Vol.254, 437-466.

2. Boulton-Stone, J. M., 1995 "The effect of surfactant on bursting gas bubbles", J. Fluid Mech., Vol.302, 231-257.

3. Sussman, M., Smereka, P., 1997 "Axisymmetric free boundary problems", J. Fluid Mech., Vol.341, 269-294.

4. Joseph, D. D., Liao, T. Y., Hu, H. H., 1993 "Drag and moment in viscous potential flow", Eur. J. Mech., B / Fluids, Vol.12, No.1, 97-106.

5. Georgescu, S.-C., 1999 "Evolution d'une bulle: Formation à partir d'un orifice et éclatement à la traversée d'une surface libre", PhD Thesis, Institut National Polytechnique de Grenoble, France.

6. Georgescu, S.-C., Achard, J.-L., Canot, E., "Jet drops ejection in bursting gas bubble processes", Eur. J. Mech., B / Fluids (submitted; last revision Jan. 2001)

7. Ivanov, I. B., Kralchevsky, P. A., Nikolov, A. D., 1986 "Film and line tension effects on the attachment of particles to an interface", J. Colloid Interf. Sci., Vol.112, No.1, 97-143.

8. Canot, E., Achard, J.-L., 1991 "An overview of boundary integral formulations for potential flows in fluid-fluid systems", Arch. Mech., Vol.43, No.4, Warszawa, 453-498.

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

10. Dold, J. W., Peregrine, D. H., 1986 "An efficient boundary-integral method for steep unsteady water waves", In: Numerical Methods for Fluid Dynamics II, eds. K. N. Morton & M. J. Baines, pp. 671-679, Oxford University Press, Oxford.

11. Machane, R., Canot, E., 1997 "High-order schemes in Boundary Element Methods for transient non-linear free surface problems", Int. J. Numer. Methods Fluids, Vol.24, 1049-1072.

12. Georgescu, S.-C., Achard, J.-L., Canot, E., 1998 "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», pp. 25-27, Istanbul, Turkey.

13. Canot, E., 1999 "Stability criteria for capillary/ gravity free-surface waves in BEM simulations of viscous potential flows", In: Proc. of the Int. Conf. on Boundary Element Techniques, ed. M. H. Aliabadi, pp. 395-404, London, U.K.

14. Wang, M., Troesch, A. W., 1997 "Numerical stability analysis for free surface flows", Int. J. Numer. Methods Fluids, Vol.24, 893-912.

15. Canot, E., Georgescu, S.-C., 1997 "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, Vol.3, pp. 23-26, Poitiers-Futuroscope, France.

ADDRESSES

(1) Chargé de Recherche au C.N.R.S. dr. ing.; IRISA-INRIA, Campus de Beaulieu, 35042 Rennes cedex; France; E-mail: Edouard.Canot@irisa.fr

(2) Asistent dr. ing.; Universitatea "Politehnica" Bucuresti, Spl. Independentei, No. 313, sector 6, 77206, Bucharest; Romania; E-mail: carmen@hydrop.pub.ro

(3) Directeur des Recherches au C.N.R.S. dr. ing.; Laboratoire des Écoulements Géophysiques et Industriels de Grenoble; B.P. 53, 38401, Grenoble Cedex 9; France; E-mail: Jean-Luc.Achard@hmg.inpg.fr