# Numerical simulation of the buoyancy-driven bouncing of a 2-D bubble at a horizontal wall

Edouard Canot, Laurent Davoust, M. El Hammoumi, D. Lachkar

► To cite this version:

Edouard Canot, Laurent Davoust, M. El Hammoumi, D. Lachkar. Numerical simulation of the buoyancy-driven bouncing of a 2-D bubble at a horizontal wall. Theoretical and Computational Fluid Dynamics, 2003, 17, pp.51-72. 10.1007/s00162-003-0096-y. hal-00203813

HAL Id: hal-00203813

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

Submitted on 18 Aug 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
# Numerical Simulation of the Buoyancy-Driven Bouncing of a 2-D Bubble at a Horizontal Wall

É. Canot

Institut de Recherche en Informatique et Systèmes Aléatoires,
Campus de Beaulieu, 35042 Rennes cedex, France
edouard.canot@irisa.fr

L. Davoust

Laboratoire des Ecoulements Géophysiques et Industriels,
B.P. 53, 38041 Grenoble cedex 09, France
laurent.davoust@hmg.inpg.fr

M. El Hammoumi and D. Lachkar

Faculté des Sciences et Techniques Fès Saïss, Département de Physique,
Laboratoire de Mécanique Appliquée, B.P. 2202, Fès, Maroc

**Abstract.** The rise of a buoyant bubble and its interaction with a target horizontal wall is simulated with a 2-D numerical code based on the Boundary Element Method (BEM). Developed from a viscous potential flow approximation, the BEM takes into account only the part of the energy dissipation related to the normal viscous stresses. Hence, a simple analytical model based on lubrication approximation is coupled to the BEM in order to compute the drainage of the interstitial liquid film filling the gap between the bubble and the near wall. In this way the bubble-wall interaction is fully computed: the approach stage, the bubble deformation stage and, depending on the values of the Reynolds number and the Weber number, the rebound and the bubble oscillations. From computation of both the bubble interface motion and the liquid velocity field, a physical analysis in terms of energy budget is proposed. Though, in the present study, the bubble under consideration is basically supposed to be a 2-D gaseous cylinder, a comparison between our numerical results and the experiments of Tsao and Koch (1997) enlightens interestingly the physics of bouncing.

## 1. Introduction

Bubbly flows play a key role in many chemical engineering processes (gas-liquid extraction, fluidised beds, flotation, sedimentation, electrochemical reactors, etc.) as well as in nature (air entrapment due to breaking waves, chutes and spillways, volcanic activity, etc.).
Investigation of the motion of a deformable bubble at the vicinity of a rigid wall is of particular relevance as it is an essential prerequisite for the development of many Eulerian/Lagrangian codes devoted to the simulation of bubbly flows. Very often, there is actual difficulty in choosing the convenient boundary condition for the dispersed phase at walls; many behaviours are observed: bubble sliding along the wall, bouncing or "sticking" on the wall.

In the literature, most of the theoretical, numerical or experimental studies devoted to the interaction of an inclusion with a wall focus on the case of solid spheres (see, e.g. Brenner, 1961; Dagan et al., 1982; Ambari et al., 1984; Sondergaard et al., 1990; Gondret et al., 1999). The impact and bouncing of droplets on hydrophobic substrates has been recently studied by Richard and Quere (2000). However, here again, like for many other experimental studies involving solid spheres, the interaction with surrounding fluid is disregarded and the rebound is analysed through the reductionist concept of elastic restitution. The other papers, which investigate the behaviour of a deformable bubble or droplet at the vicinity of a wall, are highly motivated by modelling the coalescence of a bubble/droplet with a plane wall or with another bubble/droplet.

To our knowledge, the hydrodynamical interaction of a bubble with a wall has never been computed in its whole, that is, the calculation of both the deformation of a liquid/gas interface and its possible bouncing. Tsao and Koch (1997) are the only authors who analysed experimentally the motion of a bubble within a liquid towards a wall. They observed, before viscosity dissipates completely its initial mechanical energy, that the bubble bounces on the wall in a dynamic process that can be understood in terms of an energy balance including the kinetic energy associated with the liquid motion, the potential energy of the interface and the gravitational potential energy gained or lost by the bubble. They distinguished three convenient times for their physical analysis:

(i) the time when the bubble is rising at its terminal velocity at a sufficient distance from the wall so that the effect of the wall on the bubble shape and fluid velocity can be neglected,

(ii) the time at which the bubble's motion is arrested upon impact with the wall and

(iii) the time when the lowest point of the bounce at which the bubble velocity again goes to zero.

Shopov et al. (1990) present numerical solutions for the unsteady viscous flow induced by a deformable gas bubble approaching or receding from a rigid boundary at moderate values of the bubble Reynolds number and the Bond number. However, the dynamics of the bubble bouncing process was not fully obtained because of the difficulty of simultaneously solving the thin liquid film entrapped between the bubble and the wall and the larger scale outer flow with their finite element solution of the Navier-Stokes equations.

The large number of theoretical or numerical studies dealing with bubble flows usually take the general potential approach: the fluid is supposed to be incompressible and irrotational which implies that the flow is potential. The dynamics of the liquid phase is therefore described by a Laplace equation. Later, using momentum conservation, non-linear equations are derived for the velocity potential at the deformable interface for the motion of the latter. This approximation is valid if the bubble Reynolds number is large and the Weber number is small. In this way Miksis et al. (1982) studied the steady rise of a bubble in an unbounded liquid and accounted for surface tension. A system of integro-differential equations was derived, including the normal component of the free-surface force balance, and treated numerically. The assumption of inviscid external flow used by these authors reduces the practical importance of their results.

In the present study the numerical procedure developed to study the interaction of a bubble with a rigid horizontal wall is based on a potential theory devoted to unsteady viscous flows. A computational model for a 2-D plane flow configuration is implemented; a Boundary Element Method (BEM) is coupled with a fourth-order Runge-Kutta explicit time-stepping technique for the temporal evolution. We thus consider the rise of a deformable gas bubble towards or away from a plane horizontal rigid wall in an unbounded liquid at moderate Reynolds number. The liquid is supposed to be initially motionless and the bubble motion is set up by buoyancy alone. This paper aims at simulating the bouncing dynamics in its whole. This objective is reached by introducing a simple drainage model between the bubble and the rigid wall. This model, coupled with the BEM code, is suitable since it is easily solved numerically (Appendix A). Further physical mechanisms like the Marangoni effect due to interfacial gradients of a surfactant concentration or like the London-Van der Waals forces (Chen and Slattery, 1982) are disregarded here despite their practical importance.
## 2. Mathematical Model

### 2.1. Governing Equations and Boundary Conditions

We consider the irrotational flow of an incompressible viscous liquid generated by the buoyancy-driven rise of a 2-D (cylindrical) bubble to a fixed and rigid horizontal wall of width $2L$ (Figure 1). The initial condition of the configuration is such that the bubble and the horizontal wall are five bubble radii apart. This precision is important since the bouncing process and particularly the height of rebound are expected to be dependent on the height of rise. The liquid flow can be described by the scalar potential, $\varphi$, so that the velocity field $\mathbf{v} = (u, w)$ is given by

$$ \mathbf{v} = \vec{\nabla}\varphi. \qquad (1) $$

Hence, the continuity equation in the fluid domain $\Omega(t)$ with the boundary $\Gamma(t)$ writes as a Laplace equation for $\varphi$:

$$ \nabla^2 \varphi = 0 \quad \text{in } \Omega(t). \qquad (2) $$

At the free surface $\Gamma_f(t)$, two collocated boundary conditions hold. The first one states that the interface is considered as a material surface, i.e. a surface to which the same fluid particles are always attached. In this way the potential $\varphi$ thus satisfies the kinetic boundary condition corresponding to this Lagrangian description of the free surface particles:

$$ \frac{D\mathbf{x}}{Dt} = \left( \frac{\partial}{\partial t} + \mathbf{v} \cdot \vec{\nabla} \right) \mathbf{x} = \vec{\nabla}\varphi \quad \text{at } \Gamma_f(t), \qquad (3) $$

where $\mathbf{x}$ defines the position vector of a free surface fluid particle and $D/Dt$ represents the material derivative following a fluid particle. The second boundary condition on the free surface is a dynamical one:

$$ \frac{\partial \varphi}{\partial t} = -\frac{1}{2}v^2 + \frac{\sigma}{\rho}C_t - \frac{2\mu}{\rho}\frac{\partial v_n}{\partial n} - gz + \frac{p_\infty - p_a}{\rho} \quad \text{at } \Gamma_f(t), \qquad (4) $$

which is derived from the transient Bernoulli’s equation

$$ \frac{\partial \varphi}{\partial t} + \frac{1}{2}v^2 + \frac{p}{\rho} + gz = \text{constant in } \Omega(t), \qquad (5) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. Geometry of the problem.</figcaption></figure>
combined with the normal momentum balance

$$p_a + \tau_{nn} = \sigma C_t, \qquad (6)$$

where $C_t$ is referred to as the sum of the principal curvatures of the liquid/gas interface and $\tau_{nn}$ denotes the normal stress:

$$\tau_{nn} = -p_1 + 2\mu \frac{\partial v_n}{\partial n}. \qquad (7)$$

The non-hydrodynamic variables $p_a$, $p_1$ and $p_\infty$ are referred to as, respectively, the gas pressure, the liquid pressure and the liquid pressure at infinity. The normal viscous effect is taken into account by the introduction of the stress $2\mu \partial v_n / \partial n$ where $\mu$, $\sigma$ and $\rho$ denote, respectively, the dynamic viscosity, the surface tension and the liquid density. This viscous potential model has been extensively described by Joseph et al. (1993).

Due to viscosity, a bulk dissipation occurs throughout the liquid extent and yields the following dissipation rate:

$$\Pi_1 = - \iint_{\Omega(t)} \tau_{ij} \frac{\partial u_i}{\partial x_j} dv.$$

Although a part of the mechanical energy is dissipated throughout the liquid bulk, the potential model cannot take into account the tangential momentum balance at the interface nor the adherence condition at the solid boundary. The result of this is that the viscous drag on the rising bubble is smaller than those obtained in a full Navier-Stokes simulation (Shopov et al., 1990). Equating the rate of dissipation and the energy power due to the buoyancy force which exerts onto a circular cylindrical rising bubble of infinite depth, the terminal velocity of such a 2-D bubble is found to be

$$U = \frac{1}{4} \rho g \frac{R^2}{\mu}. \qquad (8)$$

Subsequently, the associated drag coefficient is:

$$C_D = \frac{8\pi}{Re_p},$$

where $Re_p$ is the particular Reynolds number calculated from the bubble diameter. This drag coefficient is slightly different from the traditional formula $32/Re_p$, which is only valid to 3-D axisymmetric configurations.

To write the mathematical model in a non-dimensional way, the variables of the problem are scaled as follows:

$$t^* = \frac{Ut}{R}, \quad \mathbf{v}^* = \frac{\mathbf{v}}{U}, \quad p^* = \frac{p}{\rho U^2}, \quad \mathbf{x}^* = \frac{\mathbf{x}}{R} \quad \text{and} \quad \varphi^* = \frac{\varphi}{UR}.$$

The characteristic length scale is defined as the equivalent bubble radius, $R = \sqrt{S/\pi}$, where S, the cross-section bubble area, is given constant during all the simulation. For convenience, we drop the * notation. Equation (4) therefore yields

$$\frac{D\varphi}{Dt} = \frac{1}{2}v^2 + \frac{1}{We}C_t - \frac{2}{Re}\frac{\partial^2\varphi}{\partial n^2} - \frac{1}{Fr}z + (p_\infty - p_a) \quad \text{at } \Gamma_f(t), \qquad (9)$$

where $Fr$, $Re$ and $We$ are non-dimensional numbers, namely, the Froude number, the liquid Reynolds number and the Weber number. Only two of them are independent:

$$We = \frac{\rho^3 g^2 R^5}{16 \mu^2 \sigma} \quad \text{and} \quad Re = \frac{\rho^2 g R^3}{4 \mu^2}, \qquad (10)$$
the third one being related to the previous two as follows:

$$Fr = \frac{Re}{4}.$$

Although these numbers are well known, it can be helpful to recall briefly their physical significance:

* We (inertia over interfacial cohesion due to surface tension): it plays a key role in bubble deformation, particularly when a dimple arises or when, after a rebound, the bubble exhibits shape oscillations.

* Re (inertia over viscous forces): it takes into account (i) the viscous stabilising effect of the interfacial stability modes (via the normal stress balance at the gas/liquid interface) when the bubble is submitted to deformation (shape oscillations), and (ii) the thinning kinetics of the interstitial liquid film during the encounter of the bubble with the horizontal wall.

* Fr (inertia over gravity forces): it takes into account buoyancy. It is fairly consistent to get a direct relationship between Re and Fr since the liquid, initially at rest, gets moving via the buoyancy-driven rise of the bubble.

Alternatively, some authors (see, e.g. Shopov et al., 1990) introduce others dimensionless numbers, namely, the Bond (or Eötvös) number (which determines the relative importance of gravity to surface tension), and the Capillary number (relative importance of the viscous shear to surface tension).

The last dimensionless parameter describes the geometry of our problem, this is the aspect ratio,

$$\varepsilon = \frac{R}{L} = \frac{1}{2.5} = 0.4,$$

which is kept constant in our simulations.

Finally, at the rigid and impermeable target wall (here referred to as $\Gamma_r$), the normal velocity is zero:

$$\bar{\nabla}\varphi \cdot \bar{n} = \frac{\partial\varphi}{\partial n} = 0 \quad \text{at } \Gamma_r, \qquad (11)$$

$\bar{n}$ being the unit outward normal vector held by the n-axis. This Neumann condition for the potential $\varphi$ could be easily generalised to moving rigid boundaries.

The term ($p_\infty - p_a$) in (9) is usually constant and taken equal to zero because the potential $\varphi$ is defined up to a constant. However, it is worth keeping this term here in order to introduce the drainage (over)pressure $p_{\text{drainage}}$ due to liquid film thinning which slows down the bubble interface motion to the rigid wall. As a matter of fact, this additional pressure term permits the coupling with the numerical scheme. In Appendix A an analytical model is derived to deal with liquid film thinning between two parallel planes; it allows us to calculate the drainage pressure. In this model a parabolic velocity profile is assumed to hold exactly as for a Poiseuille flow. The transient Bernoulli's equation thus modified can be written as

$$\frac{D\varphi}{Dt} = \frac{1}{2}v^2 + \frac{1}{We}C_t - \frac{2}{Re}\frac{\partial^2\varphi}{\partial n^2} - \frac{1}{Fr}z + p_{\text{drainage}}. \qquad (12)$$

## 2.2. Boundary Element Method

The boundary value problem described by the system of equations (1)–(2) can be transformed applying Green's second identity to the velocity potential $\varphi$ (Stakgold, 1979) into the following boundary integral equation:

$$\alpha(\mathbf{r})\varphi(r) = \int_{r' \in \Gamma(t)} \left( \frac{\partial\varphi}{\partial n'} G(\mathbf{r}, \mathbf{r}') - \varphi \frac{\partial G}{\partial n'}(\mathbf{r}, \mathbf{r}') \right) dr', \qquad (13)$$

where $\mathbf{r}, \mathbf{r}', G$ and $\alpha(\mathbf{r})$ denote, respectively, the field and the source points, the 2-D Green function (free space fundamental solution) for (2),

$$G(\mathbf{r}, \mathbf{r}') = -\frac{1}{2\pi} \ln (|\mathbf{r} - \mathbf{r}'|),$$
and a geometrical boundary coefficient whose value depends on the shape of $\Gamma_f(t) \cup \Gamma_r : \alpha(\mathbf{r}) = \pi$ locally smooth), $\alpha(\mathbf{r}) = \pi/2$ corner), etc. From a physical point of view, (13) can be seen as the potential created by two superposed layers of singularities: a source layer (single layer of density $\partial\varphi/\partial n'$) and a normal doublet layer (double layer of density $-\varphi$).

First, (13) is written on the part of the boundary with a Neumann condition (i.e. the rigid boundary) and then combined with (11) to yield a Fredholm integral equation of the second kind for the unknown $\varphi(r')$. Next, (13) is differentiated with respect to the normal $\mathbf{n}$ at the point $\mathbf{r}$ and transformed as follows:

$$ \alpha(\mathbf{r}) \frac{\partial \varphi}{\partial n} = \int_{r' \in \Gamma(t)} \frac{\partial \varphi}{\partial n'} \frac{\partial G}{\partial n}(\mathbf{r}, \mathbf{r}') dr' + \int_{r' \in \Gamma(t)} \gamma(r') \mathbf{H}(\mathbf{r}, \mathbf{r}') \cdot \mathbf{n}(r) dr', \quad (14) $$

where $\mathbf{H}(\mathbf{r}, \mathbf{r}')$ is the velocity field created by a 2-D vortex point:

$$ \mathbf{H}(\mathbf{r}, \mathbf{r}') = \frac{1}{2\pi} (\mathbf{r} - \mathbf{r}') \wedge \mathbf{k}, $$

$\mathbf{k}$ denoting the third basis vector perpendicular to the 2-D plane, and then is written on the part of the boundary with a Dirichlet condition (i.e. the interface). The second integral in (14) comes from the equivalence between a normal doublet layer and a vortex layer, whose density is related to the former one by

$$ \gamma(r') = \frac{\partial(-\varphi)}{\partial s}, $$

$s$ denoting the boundary arclength.

Thus (14) is another Fredholm integral equation of the second kind for the unknown $\partial\varphi/\partial n'$. The singularity of the kernels $G$ and $\mathbf{H}$ ensures diagonal dominance in the matrix of the corresponding algebraic system and the problem is in general well conditioned (Jaswon and Symm, 1977).

## 2.3. Derivation of the Energy Components

As demonstrated in the following, it is worthwhile studying the physical problem in term of a mechanical energy budget. The classical energy components are derived in a straightforward manner:

$$ \text{Kinetic energy:} \quad E_K = \frac{1}{2} \int_{\Sigma} \varphi \frac{\partial \varphi}{\partial n} dA, $$

$$ \text{Potential energy:} \quad E_P = \frac{1}{We} \int_{\Sigma_i} dA + \frac{1}{2Fr} \int_{\Sigma_i} z^2 dA. $$

The mechanical energy is the sum of the kinetic and the potential energies, the latter being the sum of the surface and gravity energies. The mechanical energy conservation is written as follows:

$$ \frac{d}{dt} (E_K + E_P) = \Pi_1, $$

where $\Pi_1$ is the rate of viscous dissipation:

$$ \Pi_1 = - \frac{2}{Re} \int_{\Sigma} \frac{\partial \varphi}{\partial n} \frac{\partial^2 \varphi}{\partial n^2} dA. $$

However, the term $\Pi_1$ is nothing but a bulk dissipation computed in the framework of a viscous potential model. As a consequence, the viscous dissipation owing to shear stresses (due to the draining stage for instance) is not included in this term. In our computation, the previous balance is rewritten under the equivalent form:

$$ E_K + E_P = \int_0^t \Pi_1 dt. $$
Hence, considering the discretisation errors, the dissipation due to liquid film thinning during the draining stage can be approached from the formula,

$$E_K + E_P - \int_{t=0}^{t} \Pi_1 dt,$$

where $t = 0$ corresponds to the start of the computation. All these energy components involve line integrals of terms which can be directly computed from the BEM and easily implemented within the code.

### 3. Numerical Formulation

The approach to compute the transient motion of the liquid/gas interface of the bubble mixes Eulerian and Lagrangian formulations as described in the following sequence:

1. Eulerian part: the Laplace equation (2) is solved by a BEM solver in relation to the boundary conditions; a Dirichlet condition for the liquid/gas interface, a Neumann condition for the solid boundary.

2. Lagrangian part: (3) and (12) are solved in order to update both the position of the liquid/gas interface and the value of the velocity potential (Figure 2). For this time-stepping procedure, a fourth-order explicit Runge-Kutta algorithm is involved.

#### 3.1. Spatial Discretisation for the BEM

The boundary of the computational domain is described by a number of elements defined by nodal points. On each element, both the boundary geometry and the field functions are mapped on a low-degree polynomial basis. Use of the Gauss-Legendre quadrature is made to perform numerical integration of the discretised form of (13) and (14). To deal with all the singular integrals involved, the classical method of subtraction of singularity is numerically performed (singular integrals of type 1/s are computed in the Cauchy Principal Value). The integral equations thus reduce to a set of linear equations with the unknown $\partial\varphi/\partial n$ on the liquid/gas interface and the unknown $\varphi$ on the rigid boundary. The set of equations is then solved by an LU decomposition method. Useful material for this classical direct approach can be found in Brebbia et al. (1984). In our particular case, however, the boundary geometry is approximated by cubic splines and use is made of a cubic Hermite polynomial approximation for the field functions (Machane and Canot, 1997): this choice has been proved to provide good accuracy in interface-tracking problems.

#### 3.2. Time-Stepping Method

As mentioned in the previous section, the boundary value problem for $\varphi$ is solved by the BEM. This formulation starts from the integral equations (13) and (14) as we have seen in Section 2.2. Provided that initial conditions are known at a given time, i.e. the domain $\Omega(t)$ and the solution of the Laplace problem, the interface position and the velocity potential can be updated to the next time step using (3) and (12), respectively. In time-dependent non-linear interface problems a boundary integral problem is solved at each time step. Since most of the computation time is devoted to that boundary integral problem, a stable and effective solution method is crucial in the time-stepping procedure. In the present case, stability considerations have led us to choose an explicit fourth-order Runge-Kutta scheme which requires information at several intermediate times between $t$ and $t + \Delta t$. As a result, for every time step $\Delta t$ we have to solve several Laplace problems in different geometries. More precisely, four computations of the integral equations are required at each iteration and that can be time-consuming in comparison with an explicit high-order Taylor series expansion<fnref n="1" /> (Machane and Canot, 1997). Here, the Runge-Kutta method suits our problem particularly well since, during the draining stage, the typical mesh size becomes smaller and smaller.

<footnote n="1">¹ In some problems this Taylor series expansion method may be preferred because of its efficiency and despite its poor stability properties as the mesh size becomes small.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Motion of material fluid particles at the bubble surface as sketched.</figcaption></figure>

### 3.3. Stability Analysis

When one deals with non-linear free-surface flows with the help of a mixed Eulerian–Lagrangian formulation, the instability of the free-surface position arises as an important problem. In order to reduce or remove this instability, many authors in the literature develop numerical techniques such as regriding, smoothing, filtering and artificial damping (Yeung, 1982; Wang and Troesch, 1997). Nevertheless, these techniques reduce the physical information considerably. Several methods based on Fourier analysis have also been developed and applied (see, e.g. Baker and Nachbin, 1998). In the present study we prefer to use a linear stability analysis based on the normal modes of a free-surface perturbation to provide a stability criteria useful to our temporal scheme. Some schemes are stable in the sense of classical numerical analysis but they diverge at infinity with an exponential-like law (no absolute stability). Here a linear stability criterion derived in Appendix B, is used, without any computational time penalty. The basic idea is to work in a new frame (s, n) centred on a node, where s and n are referred to as, respectively, the tangential and the normal coordinates of the free surface (see also Canot, 1999).

### 3.4. Nodes Distribution

The temporal interface evolution is determined through a Lagrangian description of a variable number of nodes unevenly redistributed on the boundary $Γ_f(t) ∪ Γ_r$ at each time step with respect to some criteria like the adaptation at the surface gradients of the potential. This leads to a concentration of nodes at places where the interface's curvature is important and when the bubble approaches close to the wall. Consequently, simulation of the interface evolution is performed during the whole bouncing process without applying smoothing techniques that can affect the physics of the phenomenon under consideration. The stability criterion established in Appendix B allows us to take full account of this adaptive meshing strategy. Another strategy advised by Baker and Nachbin (1998) is to keep an even distribution of the nodes along the interface and to develop a stability criterion based on Fourier analysis.

## 4. Results and Discussion

Figure 3 exhibits a series of snapshots obtained from the numerical simulation of the bouncing when $We = 0.735$ and $Re = 210$. These values match the experiments of Tsao and Koch (1997) performed on millimetric bubbles. Figure 3 shows that, during the approach stage, the bubble deforms and its shape progressively becomes elliptic owing to a 2-D stagnation liquid flow at the wall which can no longer be neglected. As a matter of fact, this stagnation flow is driven by the bubble motion itself. The bubble motion seems to stop for a moment before rebound. In agreement with the existing literature (Nakamura and Uchida, 1980; Lin and Slattery, 1982), a dimple progressively grows at the top of the bubble (seventh snapshot). A favourable pressure gradient allows the interstitial liquid film entrapped between the bubble and
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Series of snapshots illustrating the bubble-wall interaction at $Re = 210$ and $We = 0.735$. The sequence is to be considered from left to right, top to bottom. The time is indicated on the first four snapshots; the non-dimensional time delay between the following snapshots is 1.5.</figcaption></figure>

the wall to drain along the wall (typical thickness of about $\frac{1}{100}$ of the bubble radius). The former pressure gradient yields in turn a local overpressure at the bubble apex and is responsible for the formation of the dimple.

After a first rebound, shape oscillations of the bubble are consistently observed. After a second bounce, the bubble position is stabilised around an equilibrium location near the wall and comes to a halt. Figure 4 shows the trajectory of the bubble centroid and its velocity during the full bouncing process. The first bounce occurs at $t = 10.8$ and the second one occurs at $t = 22.5$ (non-dimensional time).

To analyse the bouncing process in its whole, we now distinguish and focus on three typical stages, namely, (i) the approach stage, (ii) the dimpling stage and (iii) the rebound and its attendant bubble shape oscillations.

## 4.1. Approach Stage

In this first stage the simulation of the bubble approaching the wall is compared with the model of O'Neill (1996) which predicts the motion of a rigid sphere to a horizontal wall submitted to a Stokeslet. Under the Stokes approximation, and if the motion of the sphere is bouyancy-driven alone, O'Neill predicts the
<figure><img src="image_5.png" /><figcaption>Figure 4. The trajectory (a) and the velocity (b) of the bubble centroid during bouncing ($Re = 210$ and $We = 0.735$).</figcaption></figure>

modified velocity of the sphere as follows:

$$V(z_0) = 1 - \frac{3}{4}\epsilon \left( \frac{3}{\pi\epsilon z_0} \tan^{-1} \left( \frac{1}{\epsilon z_0} \right) + \frac{3 + 5\epsilon^2 z_0^2}{\pi(1 + \epsilon^2 z_0^2)^2} \right); \quad (15)$$
<figure><img src="image_6.png" /><figcaption>Figure 5. Time-dependence of the distance between the bubble centroid and the wall. Comparison between O’Neill’s analytical prediction and the present study for two computations: $Re = 40$ and $We = 0.01$ and $Re = 40$ and $We = 1$.</figcaption></figure>

$V(z_0)$ denotes the (non-dimensional) velocity of a rigid sphere scaled by its terminal velocity $U$ calculated with the assumption of an unbounded liquid bulk. The scalar $z_0$ denotes the time-dependent distance between the particle centroid and the target wall scaled by the bubble radius. Note that if the aspect ratio $\epsilon$ tends to zero (a situation where the particle is approaching very close to an infinite wall), then (15) simplifies to the well-known model of Brenner (1961):

$$V(z_0) = 1 - \frac{9}{8} \frac{1}{z_0}.$$

Figure 5 exhibits a comparison between (15) above and two numerical computations performed in the present study ($Re = 40$) when the bubble is either rigid ($We = 0.01$) or deformable ($We = 1$). Even if the two-dimensionality of our numerical simulations makes comparison with the 3-D model of O’Neill difficult, a few points are worthy of note. Large discrepancies are made evident which are expected considering the assumptions of the Stokes theory of O’Neill (quasi-steady creeping flow). Despite the fact that our numerical simulations take into account the transient mass-added force, in contrast to O’Neill’s model, they neglect the viscous drag which exerts onto the bubble due to the tangential viscous stresses (viscous potential theory). Consequently, the bubble velocity computed in the present study is consistently larger than the bubble velocity predicted by O’Neill. Moreover, the same goes for bubble acceleration since the mass-added force is not taken into account in O’Neill’s model.

The first simulation ($We = 0.01$) for which the bubble contour remains circular is interestingly compared with the second simulation ($We = 1$) for which the distance $z_0$ can decrease to a value lower than unity since the bubble becomes deformable.

## 4.2. Dimple Formation, Rebound

Bubble motion is dramatically damped at the proximity of the wall due to the induced stagnation liquid flow in a first time and due to the fact that the thin liquid film entrapped between the bubble and the wall has
to drain in a second time. The bubble interface exhibits a dimple on the side closest to the wall as a result of the pressure drop required to drive the fluid out of the gap. Modelling the thinning of the liquid film is the object of a huge literature (see, e.g. Platikanov, 1964; Nakamura and Uchida, 1980; Lin and Slattery, 1982; Chen and Slattery, 1982; Doubliez, 1991; Bazhlekov et al., 2000). Two competitive scenarios hold:

* either the liquid film drains quickly until it becomes thin enough to break out (then the bubble coalesces with the wall),

* or the potential energy due to the interfacial deformation of the bubble is efficiently released and rebound is observed.

In the present study the 2-D bubble under consideration stands as a gaseous cylinder and the dimpling effect is relatively more important than the one we can observe with a spherical (3-D) bubble (see, hereafter, Section 5). The dimple is described by one curvature radius for a 2-D bubble whereas it is described by two equal curvature radii for a 3-D bubble. To store and release the same input kinetic energy by way of the interfacial (potential) energy, the dimpling effect has to be more pronounced in the 2-D case (Laplace's law). Figure 6 represents the effect of the We number on the dimpling effect, whose magnitude consistently increases with We: lowering the surface tension permits a larger deformation of the bubble.

Even in the 2-D configuration, bouncing is not observed for all the points scanned. Figure 7 represents a contour map of the rebound amplitude of the bubble centroid in the (We, Re) coordinates; the Weber number ranges between 0.1 and unity whereas the Reynolds number ranges between 0.05 (Stokes regime) and 200 (inertial regime). Our results show that bouncing occurrence strongly depends on the Reynolds number. The dashed curve in Figure 7 shows the onset of bouncing, arbitrarily defined when the apex bounce is greater than one-tenth of the radius: this curve describes nearly a straight line located at $Re = 15$. For high values of the Reynolds number ($Re \sim 100$), bouncing stands as very vigorous: consecutive rebounds are observed. Oscillations modes are expected to be excited at the bubble surface and could enhance bubble fragmentation. Higher values of the Reynolds number are difficult to simulate because the interstitial film becomes thinner and thinner: the mesh needs to be refined much more. This point is reminiscent of the difficulties encountered by Shopov et al. (1990) with a finite element method when the bubble is approaching very close to the wall. Nevertheless, a dependence of bouncing on the Weber number can be found when the Reynolds number is large enough (larger than 100, typically) or small enough (smaller than 100 typically). Then, when $Re > 100$, it is worth noting that bouncing is strongly damped by lowering the surface tension (large values of the Weber number), which correlates surprisingly with a better ability of the bubble to deform. Anticipating on the following, it is worth emphasizing that the surface potential energy drives the ability of the bubble to rebound or not. Subsequently, when $Re > 100$, the bubble deformation is not large enough to compensate the lowering of the surface tension and of the surface potential energy too.

### 4.3. Rebound and Bubble Oscillations

In the third stage the bubble bounces and shape oscillations are observed. Figure 8 shows different modes of oscillations. These modes are the Fourier-cosine amplitude of the radius $r$ of the markers in polar coordinates defined in a frame attached to the bubble mass centre:

$$r = \sqrt{x^2 + (z - z_G)^2}.$$

As shown in Figure 3, only modes 2–4 seem to be dominant. The higher modes are significantly damped by viscosity.

In Figure 9 the variation of the different energy components during the full bouncing process have been plotted. According to Tsao and Koch (1997), the rebound results from an interchange between (i) the kinetic energy of fluid motion, (ii) the surface (potential) energy of the gas-liquid interface and (iii) the potential energy due to gravity. According to these authors, the energy balance shows that about 41% of the initial energy is converted into gravitational and surface potential energy at the lowest point of the bounce trajectory. Figure 9(a) shows that at time = 16.5, which represents the lowest point of the bounce trajectory (Figure 4(a)), about 60% of the initial energy is converted into potential energy. The deviation with Tsao and Koch's findings is therefore about 20%. To explain this difference, first,
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Thinning of the dimple for three values of We and Re = 210. The profiles, referred to as 1–5, are time-delayed with a step $\delta t = 0.75$ for We = 0.5 (a), $\delta t = 0.735$ for We = 0.735 (b) and $\delta t = 1.12$ for We = 1 (c).</figcaption></figure>

it is worth noting that viscous dissipation remains quasi-permanent except during liquid film thinning (first and second bounces, Figure 9(a)). Moreover, bulk dissipation (through the liquid) is small compared with film dissipation. Consequently, the viscous drag which exerts on the bubble in the experiments of Tsao and Koch, and which is not taken into account in our simulation, contributes in a non-negligible way to energy dissipation. Second, the bubble under consideration here is cylindrical and not spherical as it is obviously the case in the experiments of Tsao and Koch. The transfer of the input energy to surface potential energy is probably much more efficient in our 2-D geometry than in its 3-D counterpart.

Figure 9(b) demonstrates that the surface energy component is a reliable indicator of any variation of bubble shape during bouncing. As a matter of fact, the oscillations made evident on the kinetic energy component after the rebound (Figure 9(a)) shows how the energy exchange between the kinetic component and the surface potential component can be efficient.
<figure><img src="image_8.png" /><figcaption>Figure 7. Contour map of the bouncing magnitude connected to the first rebound $We = 0.1–1$, $Re = 0.05–200$).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Modes of deformation for the bubble surface.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Time-dependence of the energy components (a) and the potential energy components (b) during the bubble-wall interaction ($Re = 210$, $We = 0.735$).</figcaption></figure>
## 5. The Case of a Spherical (3-D) Bubble

The interaction of a spherical bubble at a horizontal wall has been computed (axial-symmetric configuration) using the same conditions as the ones of the former 2-D simulation, $We = 0.735$ and $Re = 210$. To simulate this 3-D case, the reader can find all the numerical changes that have to be made to the model (Green's function, updating of the bubble's location and potential) in a companion paper<fnref n="2" /> (Georgescu et al., 2002).

The simulation shows that the shape of the bubble changes and a dimple consistently arises upon impact. In contrast with the 2-D configuration, the bubble stops at the wall and rebound do not occur. To interpret this, it is worthwhile considering experiments on particle bouncing at a rigid wall. According to Tsao and Koch (1997), a possible physical mechanism responsible for the rebound of a bubble is a phenomenon of boundary layer separation from the rigid wall. During the approach of the bubble, the flow is radially outward, yielding a favourable pressure gradient that is not expected to cause a boundary layer separation. During the rebound, the flow becomes radially inward in the direction of increasing pressure (high pressure at the bubble apex). Then a boundary layer may separate from the rigid wall and lead to a vortical flow behind the bubble that would dissipate much more energy.

Recently, Eames and Dalziel (2000) investigated experimentally the flow around a rigid sphere moving towards or away from a rigid wall. When the sphere approaches the wall and stops on making contact with it, they observed that the wake vortex, initially located at the rear of the sphere, was able to shroud it, to impact on the wall and then to generate a secondary vortex ring. The coherent structure, composed of the wake and the secondary vortices, strikes the wall and pushes liquid, initially adjacent to the wall, to outer.

All these mechanisms are intrinsically related to the viscous diffusion of the vorticity.<fnref n="3" /> They cannot be simulated by our code which assumes the flow around the bubble to be potential. This point, which is clearly the main limitation of the present numerical study, could explain why in the conditions of our 3-D simulation bubble bouncing is not observed.

# Appendix A. The Film Thinning Model

## A. 1. Drainage Pressure

Our objective is to introduce a drainage pressure which exerts on the side of the bubble which impacts the wall. The kinematics of the flow is given by a numerical method. The flow, restricted to the draining film, is considered to be permanent and creeping ($Re \ll 1$); this allows us to neglect inertia. Within the gap between the wall and the bubble, the velocity profile is supposed to be parabolic, $u = U(x, z) \cdot x$, like a Poiseuille flow between parallel plates. On the rigid wall, the boundary condition on the velocity is $u = 0$ which supposes the adherence. A Neumann boundary condition applies at the bubble interface and states that the transversal gradient of velocity vanishes (gas viscosity is neglected). Consequently, we make use of a local relationship between the axial pressure gradient and the mean velocity $\overline{U}$. As in the case of a one-directional flow, the pressure gradient is supposed longitudinal and independent of the (transversal) vertical coordinate $z$ (Figure A. 1). Thus, the pressure is obtained by integration of the gradient along the interface.

## A. 2. Two-Dimensional Plane Model

The velocity profile is given by

$$v_x(x, z) = -\frac{3 \overline{U}(x)}{2 e^2(x)} z(z + 2e(x)). \qquad (A.1)$$

Thus we have

$$\frac{\partial p}{\partial x} = \mu \frac{\partial^2 v_x}{\partial z^2} = -3\mu \frac{\overline{U}(x)}{e^2(x)},$$

<footnote n="2">² The viscous potential BEM method has been also used by Georgescu et al. (2002) to simulate the jet drops ejection in a bursting gas bubble process.</footnote>

<footnote n="3">³ Note that for the particular case of one single particle rising in an unbounded medium, the unsteady diffusion of the vorticity is fairly described by the (history) Basset term of the BBOT equation.</footnote>
<figure><img src="image_12.png" /><figcaption>Figure A. 1. Schematic of the liquid film drainage.</figcaption></figure>

where $\bar{U}$ denotes the velocity calculated from the liquid rate flow, $\mu$ is the dynamical viscosity and $e(x)$ is the $x$-dependent liquid film thickness. The same equation is written in a non-dimensional way as

$$ \frac{\partial p}{\partial x} = - \frac{3 \bar{U}(x)}{Re e^2(x)}. \qquad (A.2) $$

### A. 3. Modification of Bernoulli's Equation

The modified Bernoulli's equation is given as

$$ \frac{D\varphi}{Dt} = \frac{1}{2}v^2 + \frac{1}{We}C_t - \frac{2}{Re}\frac{\partial^2\varphi}{\partial n^2} - \frac{1}{Fr}z + p_{\text{drainage}}. \qquad (A.3) $$

For any point M at the interface, the drainage pressure is derived as

$$ p_{\text{drainage}}(s) = p_R + \int_{s_R}^{s_S} \frac{\partial p}{\partial x} ds. $$

To each point M at the interface, a point $M'$ at the rigid wall can be associated. These points are joined by a parabolic segment issued perpendicularly from the rigid wall. The velocity $\bar{U}(x)$ in the $MM'$ section is computed as the average of the velocities at points M and $M'$. The particular point R (the first point of the integration) is a mesh node of the interface which corresponds to $x_{R'} = 1.5$. This point is considered to be outside the film. Consequently, at this point $x_{R'}$, the following condition holds:

$$ p_R = p_{\text{drainage}}(s = R) = 0. $$

In all cases it is checked that $p_{\text{drainage}}(s) > 0$ if $\bar{U} > 0$ along any cross-section of the liquid film.
# Appendix B. Derivation of the Linear Stability Criterion

A linearized stability method based on the normal modes analysis of free-surface perturbations was presented by Canot (1999). From this method, a stability criterion for the temporal RK4 scheme involved in this paper is made evident via the symbolic computation of two eigenvalues of an amplification matrix. That criterion allows us to choose the optimal time step in a straightforward manner and our BEM calculations are subsequently stable over a large time without any artificial smoothing.

## B. 1. Perturbation's Equations in Linear Theory

Let us suppose that the markers describing the free surface are moved from their original position with a small normal displacement. This leads to the situation illustrated in Figure B. 1. In order to derive an equation for the amplitude $\eta(t)$ of this perturbation, it is worthwhile working in a new frame $(s, n)$ centred at a node, where $s$ and $n$ stand respectively as the tangential and the normal coordinates of the free surface (see Figure B. 2). In this frame the non-dimensional normal acceleration of the reference node and the local Froude number are denoted, respectively, as

$$ \gamma_n = \frac{\Gamma_{en} L}{U^2}, \qquad (B.1) $$

$$ \text{and} \quad \frac{1}{Fr_{\text{loc}}} = \gamma_n + \frac{\cos(\beta)}{Fr}. \qquad (B.2) $$

Consequently, (3) and (9) for the perturbation $\eta(t)$ (see Section 2) can be written

$$ \left\{ \begin{aligned} \frac{\partial \varphi}{\partial t} &= -\frac{1}{Fr_{\text{loc}}} \eta + \frac{1}{We} \kappa + \frac{2}{Re} \frac{\partial^2 \varphi}{\partial s^2}, \\ \frac{\partial \eta}{\partial t} &= \frac{\partial \varphi}{\partial n}, \end{aligned} \right. \qquad (B.3) $$

where the $\frac{1}{2}V^2$ term has been dropped since the perturbation is supposed to be small. By using Fourier decomposition, the problem can be reduced, keeping only the most unstable mode, i.e. the one which has the smallest wavelength. This critical wavelength $\lambda$ is twice the distance $\Delta s$ between two nodes so that we can restrict our study to a solution of the following form:

$$ \eta(s, t) = \eta_0 \sin(ks) H(t), \\ \text{and} \quad \varphi(s, n, t) = \varphi_0 \sin(ks) \exp(kn) \Phi(t), $$

<figure><img src="image_13.png" /><figcaption>Figure B. 1. Displacement of the markers due to interfacial perturbation $\eta(t)$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure B.2. Local frame (s, n) at a node.</figcaption></figure>

where $k = 2\pi/\lambda = \pi/\Delta s$ is the wavelength and $\eta_0, \varphi_0$ are small amplitudes. The symbols $H(t)$ and $\Phi(t)$ are functions we do not need to explicit further. The above solution, injected into (B. 3), leads to the expression

$$ \begin{cases} \frac{\partial \varphi}{\partial t} = -\frac{1}{Fr_{\text{loc}}} \eta - \frac{k^2}{We} \eta - \frac{2k^2}{Re} \varphi, \\ \frac{\partial \eta}{\partial t} = k \varphi, \end{cases} \qquad (B.4) $$

which is equivalent to the following second-order differential equation:

$$ \eta'' + 2\xi\omega_0\eta' + \omega_0^2\eta = 0, \qquad (B.5) $$

$$ \omega_0^2 = \frac{k}{Fr_{\text{loc}}} + \frac{k^3}{We}, \qquad \xi = \frac{k^2}{Re\omega_0}. \qquad (B.6) $$

The symbols $\omega_0$ and $\xi$ stand, respectively, for the undamped frequency and the damping coefficient. The solutions of (B. 5) are well-known damped functions, the explicit form of which (oscillating or not) depends on the value of $\xi$; it follows that possible instabilities occurring during simulation are of a numerical nature, not a physical one. The first term of (B. 6) is the classical dispersion relation for the gravity/capillary waves in the case where there is no viscous term; it shows clearly that capillary phenomena lead to harder problems than gravity phenomena do. In the next section, (B. 5) is involved in deriving a stability criterion since our aim is to solve it by Runge-Kutta methods.

## 2. Derivation of Stability Criteria

In this section we establish a stability criterion which expresses that the time step $\tau$ must be lower than a critical value $\tau_c$. In all cases this critical value $\tau_c$ can also be written as $f(\xi)/\omega_0$ where the function $f(\xi)$ is determined as much as possible in an analytic closed form. Beforehand, the ordinary differential equation (B. 5) is transformed by introducing the new time variable $\omega_0 t$ so that we have to solve a system of two first-order differential equations:

$$ \begin{cases} \eta' = \psi, \\ \psi' = -2\xi\psi - \eta. \end{cases} $$

When Runge-Kutta methods are used to solve the latter system, we obtain a recurrence equation which can be written according to the matrix form:

$$ \begin{pmatrix} \eta_{n+1} \\ \psi_{n+1} \end{pmatrix} = M \begin{pmatrix} \eta_n \\ \psi_n \end{pmatrix}; $$
$M$ is called the amplification matrix. For instance, the explicit first-order Runge-Kutta method leads to

$$ \begin{pmatrix} \eta_{n+1} = \eta_n + \tau \psi_n \\ \psi_{n+1} = \psi_n - \tau (2\xi\psi_n + \eta_n) \end{pmatrix}; $$

the corresponding amplification matrix is written as

$$ M = \begin{bmatrix} 1 & \tau \\ -\tau & 1 - 2\xi\tau \end{bmatrix}. $$

Finite values (non-divergent solution) are obtained if the module of the two eigenvalues of $M$ are both less than unity:

$$ |\lambda_i| < 1, \quad i = 1, 2. \qquad (B.7) $$

Note that the particular structure of the computing sequence presented in Section 2 enables us to employ numerical schemes whose equation in $\psi$ (corresponding to the dynamic equation (9)) is written in an implicit form with respect to the space variable $\eta$ only. For example, the semi-implicit first-order Runge-Kutta method leads to

$$ \begin{cases} \eta_{n+1} = \eta_n + \tau \psi_n, \\ \psi_{n+1} = \psi_n - \tau (2\xi\psi_n + \eta_{n+1}). \end{cases} $$

This kind of semi-implicit scheme greatly improves stability, without iterations, because during the computation we can split the update of the free-surface position $\eta$ and the update of the potential $\varphi$.

We now present analytical results, which have been obtained with the computer algebra system MAPLE, for the fourth-order Runge-Kutta scheme used in this paper. Hereafter, $P(\lambda)$ is the characteristic polynomial

<figure><img src="image_15.png" /><figcaption>Figure B.3. Value of the time step against the parameter $\xi$ (explicit scheme).</figcaption></figure>
of the amplification matrix M.

$$P = \lambda^2 + \lambda \cdot \left( \tau^2 - 2 - \frac{1}{12}\tau^4 - \tau^3\xi - 2\tau^2\xi^2 + \frac{2}{3}\xi^2\tau^4 - \frac{2}{3}\xi^4\tau^4 + 2\tau\xi + \frac{4}{3}\xi^3\tau^3 \right) \\ + 1 + \frac{1}{12}\xi\tau^5 + \frac{1}{576}\tau^8 + \frac{2}{3}\xi^4\tau^4 + 2\tau^2\xi^2 - 2\tau\xi - \frac{4}{3}\xi^3\tau^3 + \frac{1}{12}\xi^2\tau^6 \\ - \frac{1}{3}\xi^3\tau^5 - \frac{1}{72}\xi\tau^7 - \frac{1}{72}\tau^6.$$

Although the two eigenvalues can always be derived analytically, it must be noted that the value of the critical time step must be numerically computed from the condition (B. 7). Time-step values are plotted in Figure B. 3.

### 3. Numerical Implementation

The introduction of the stability criterion in the BEM simulation is done at each time step as follows:

* For each node, compute $k = \pi/\Delta s$ and the local Froude number $Fr_{loc}$ from (B. 2).

* Derive $\omega_0$ and $\xi$ in terms of $k$, $Fr_{loc}$, We and Re from (B. 6).

* Compute the critical time step $\tau_c$ from the above criterion.

* Then retain the smallest value of $\tau_c$ for all the nodes.

### References

* Ambari, A., Gauthier-Manuel, B., and Guyon, E. (1984). Wall effects on a sphere translating at constant velocity. *J. Fluid Mech.*, **149**, 235–253.

* Baker, G., and Nachbin, A. (1998). Stable methods for vortex sheet motion in the presence of surface tension. *SIAM J. Sci. Comput.*, **19**, 1737–1766.

* Bazhlekov, I.B., Chesters, A.K., and Van de Vosse, F.N. (2000). The effect of the dispersed to continuous-phase viscosity ratio on film drainage between interacting drops. *Int. J. Multiphase Flow*, **26**(3), 455–466.

* Brebbia, C.A., Telles, J.C.F., and Wrobel, L.C. (1984). *Boundary Element Techniques. Theory and Applications in Engineering*. Springer-Verlag, New-York.

* Brenner, H. (1961). The slow motion of a sphere through a viscous fluid towards a plane surface. *Chem. Eng. Sci.*, **16**, 242–251.

* Canot, É. (1999). Stability criteria for capillary/gravity waves in BEM simulations of viscous potential flows. *Proc. Int. Conf. on Boundary Element Techniques*, Queen Mary & Westfield College, London, 6–8 July (M.-H. Aliabadi, ed.), pp. 321–330.

* Chen, J.D., and Slattery, J.C. (1982). Effects of London-Van der Waals forces on the thinning if a dimpled liquid as a small drop or bubble approaches a horizontal solid plane. *AIChE. J.*, **28**, 955–962.

* Dagan, Z., Pfeffer, R., and Weinbaum, S. (1982). Axisymmetric stagnation flow of a spherical particle near a finite planar surface at zero Reynolds number. *J. Fluid Mech.*, **122**, 273–294.

* Doubliez, L. (1991). The drainage and rupture of a non-foaming liquid film formed upon bubble impact with a free surface. *Int. J. Multiphase Flow*, **17**(6), 783–803.

* Eames, I., and Dalziel, S.B. (2000). Resuspension of dust by the flow around a sphere impacting a wall. *J. Fluid Mech.*, **403**, 305–328.

* Georgescu, S.-C., Achard, J.-L., and Canot, E. (2002). Jet drops ejection in bursting gas bubble processes. *Eur. J. Mech. B/Fluids*, **21**(2), 265–280.

* Gondret, P., Hallouin E., Lance, M., and Petit, L. (1999). Experiments on the motion of a solid sphere toward a wall: from viscous dissipation to elastohydrodynamic bouncing. *Phys. Fluids.*, **11**(9), 2803–2805.

* Jaswon, M.A., and Symm, G.T. (1977). *Integral Equation Methods in Potential Theory and Elastostatics*, Chapt. 12. Academic Press, New York.

* Joseph, D.D., Liao, T.Y., and Hu, H.H. (1993). Drag and moment in viscous potential flow. *Eur. J. Mech. B/Fluids*, **12**(1), 97–106.

* Lin, C.Y., and Slattery, J.C. (1982). Thinning of a liquid film as a small drop or bubble approaches a solid plane. *AIChE. J.*, **28**, 147–156.

* Machane, R., and Canot, E. (1997). High-order schemes in Boundary Element Methods in fluids. *Int. J. Numer. Methods Fluids*, **24**, 1049–1072.

* Miksis, M.J., Vanden-Broeck, J.M., and Keller, J.B. (1982). Rising bubbles. *J. Fluid Mech.*, **123**, 31–41.

* Nakamura, M., and Uchida, K. (1980). Thinning phenomena of the dimple between a horizontal glass plate and an air bubble in water. *J. Colloid and Interface Sci.*, **78**(2), 479–484.
O'Neill, M.E. (1996). On the modeling of particle-body interaction in stokes flows involving a sphere and circular disc or a torus and circular cylinder using point singularities. *Chem. Eng. Comm.*, **148–150**, 161–182.

Platikanov, D. (1964). Experimental investigation on the dimpling of thin liquid film. *J. Phys. Chem.*, **68**, 3619.

Richard, D., and Quere, D. (2000). Bouncing water drops. *Europhys. Lett.*, **50**(6), 769–775.

Shopov, P.J., Minev, P.D., Bazhlekov, I.B., and Zapryanov, Z.D. (1990). Interaction of a deformable bubble with a rigid wall at moderate Reynolds numbers. *J. Fluid Mech.*, **219**, 241–271.

Sondergaard, R., Channey, K., and Brennen, C.E. (1990). Measurements of solid spheres bouncing off flat plates. *J. Appl. Mech.*, **57**, 694–699.

Stakgold, I. (1979). *Green's Functions and Boundary Value Problems*. Wiley, New York.

Tsao, H.-K., and Koch, D. (1997). Observations of high Reynolds number bubbles interacting with a rigid wall. *Phys. Fluids.*, **9**(1), 44–56.

Wang, M., and Troesch, A.W. (1997). Numerical stability analysis for free surface flows. *Int. J. Numer. Methods Fluids*, **24**, 893–912.

Yeung, R.W. (1982). Numerical methods in free surface flows. *Ann. Rev. Fluid. Mech.*, **74**, 395–442.