# Quasi-steady formation of bubbles and drops viewed as processes that break bifurcation

Jean-Luc Achard, Sanda-Carmen Georgescu

► To cite this version:

Jean-Luc Achard, Sanda-Carmen Georgescu. Quasi-steady formation of bubbles and drops viewed as processes that break bifurcation. Journal of Engineering Mathematics, 2005, 52, pp.147-164. 10.1007/s10665-004-3065-6 . hal-00182329

HAL Id: hal-00182329

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

Submitted on 21 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
# Quasi-steady formation of bubbles and drops viewed as processes that break bifurcation

JEAN-LUC ACHARD¹ and SANDA-CARMEN GEORGESCU²

¹Laboratoire des Écoulements Géophysiques et Industriels, B.P. 53, 38401 Grenoble Cedex 9, France
(E-mail: jean-luc.achard@hmg.inpg.fr); ²Hydraulics Department, University “Politehnica” of Bucharest, 313 Spl.
Independentei, 77206 Bucharest, Romania (E-mail: carmen@hydrop.pub.ro)

**Abstract.** Instabilities occurring during the quasi-steady formation of bubbles and drops at a circular orifice through a thin plate, which separates a cylindrical upper vessel of quiescent liquid from a lower air chamber maintained at constant pressure, are considered in connection with the Rayleigh-Taylor instability of a flat meniscus. The control parameters of the problem are the Eötvös number $E_ö$, and the excess-pressure number $\Delta$, which characterizes the pressure difference between gas and liquid across the orifice. The Rayleigh-Taylor case *i.e.*, $\Delta=0$, can be viewed as a perfect bifurcating problem. A subcritical bifurcation emerges from the critical point, $E_ö_c = 5.783186$, beyond which the flat meniscus is unstable for axisymmetric perturbations. Bubbles ($\Delta > 0$), and drops ($\Delta < 0$) appear as solutions that break bifurcation. When an appropriate measure of their magnitude $\epsilon$ is introduced, it can be shown analytically that the equilibrium surface at the critical point is a cusp; its intermediate sheet is stable, while its two upper and lower sheets are unstable. The analytical bifurcation set onto the control parameters plane is valid only around the critical point.

**Key words:** bifurcation curve, bubbles and drops formation, cusp, meniscus stability, Rayleigh-Taylor

## 1. Introduction

In this paper, we study the well-known problem of the Rayleigh-Taylor instability in connection with the more general problems regarding the formation of emerging bubbles and pendant drops. The Rayleigh-Taylor instability refers to the instability of a plane interface between two immiscible fluids of different density in an external gravitational field. While examining the stability of drops and bubbles, at least in their quasi-steady formation regime, we are led to consider interfaces that are not necessarily flat, but upward (bubbles), or downward (drops) curved. The purpose of this paper is to present a unified theoretical treatment of those two types of instability in a well-defined physical system.

Rayleigh-Taylor instability has generated a large theoretical literature [1], tracing back to Maxwell [2, pp. 56–71], long before Rayleigh [3] and Taylor [4]. As noted by Jacobs and Catton [5], there is a comparable amount of work devoted to the experimental validation of the theory, beginning with the pioneering studies of Duprez [6, Plate 1], [7, Plate 1], and Lewis [8]. The formation of emerging bubbles and pendant drops is not less studied in many cases, both theoretically and experimentally [9]; summaries of the pertinent literature are presented by Tsuge [10, Chapter 9, pp. 191–232], and Sadhal *et al.* [11, Chapter 7]. Facing the huge amount of the above mentioned studies, it is useful to classify them according to some simple criteria, which will characterize either peculiarities of the physical systems, or assumptions of the theoretical models. This classification will help us to specify the framework of our own study.
Physical systems can be distinguished by the following three main conditions: *Control condition* – one of the fluids is maintained either at a constant pressure, or at a constant flow; *Lateral constraints* – the two superposed fluids being contained into a finite apparatus, the interface is attached at some edge (orifice through a thin plate, opening mouth of a nozzle), or it joins the vertical and more or less smooth bounding walls; *Geometrical configuration of the bounding walls* – it may be two-dimensional (parallel and vertical plates separated by a gap), axisymmetric (vertical circular cylinder), or three-dimensional (rectangular vessel).

We are interested in a meniscus attached to the edge of a circular orifice through a thin plate. This plate separates an upper cylindrical vessel of quiescent liquid, from a lower air chamber maintained at constant pressure.

Regarding the Rayleigh-Taylor instabilities, the theoretical models can be distinguished by at least six basic categories of assumptions that form their basis: (I) The interface is of infinite, or finite lateral extent. Many theoretical works, but not all of them [12], have adopted the former condition, leaving aside boundary conditions at the interface. On a smooth, rigid, solid surface, the contact line is free to move, and one is faced with the difficult problem of removing the singularity resulting from the usual hydrodynamic assumptions, i.e., Newtonian fluid behaviour, and no slip on the wall [13]. Fortunately, there are simple situations where the contact line can be considered as prescribed at a fixed location (e.g. sharp edge of an orifice, wire rim); (II) Upper and lower fluids are infinitely thick, or have a finite thickness [14]; (III) Fluids may be inviscid and irrotational, or viscous effects are introduced [15]; (IV) The density of the lighter fluid is negligibly small, or a finite density ratio is allowed [16]; (V) Axisymmetric, two-dimensional, or fully three-dimensional initial disturbances are considered [17]. Oblong apparatuses promote the second case; (VI) The growth rates of instabilities are examined according to linear or nonlinear methods. The former are analytic, being restricted to describe the initial stage of the process, while the latter are generally numerical [18], describing the subsequent stages (spikes, bubbles, drops, ... ). Nevertheless, some analytical studies [12, 19] are devoted to these stages.

In any of these six categories of assumptions, the first option corresponds to the simplest, and to the more manageable assumption, while the second option corresponds to a more realistic, or more general model. Each time, one among the first authors who tried to perform this extension has been quoted.

In our study, we will adopt systematically the simplest option, except for the first category of assumptions where, as stated above, the meniscus is attached to the edge of an orifice. The fixed end-point condition applies on it. As we shall see, there are some loose restrictions on the contact angle at this point.

To some extent our model is modest. Our objective is not to refine the analysis of the Ray-
leigh-Taylor instability, which was already thoroughly studied, but to show that, if we enlarge
the set of control parameters for this phenomenon, we address simultaneously other phenomena
as the slow (quasi-steady) formation of bubbles and drops. This enlargement means basically to
take into account initial and unavoidable pressure steps, which exist at the interface. In order to
provide a new global picture of previously separated instabilities, a simple model is sufficient.

## 2. Formulation

Consider a physical system composed of a cylindrical vessel of liquid over a gas chamber. The vessel is assumed open at the upper side. A thin plate perforated with a single small orifice of radius $R$ separates both fluids. At the orifice, a meniscus forms (Figure 1).

The meniscus may be flat (Rayleigh-Taylor case), upward oriented (emerging bubble), or
downward oriented (pendant drop). The gas density is very small with respect to the liquid
<figure><img src="image_2.png" /><figcaption>Figure 1. Emerging bubble (left), flat meniscus (centre), and pendant drop (right).</figcaption></figure>

density $ρ_L$, so that the gas pressure $p_G$ is uniform, while the vertically oriented gravity $g$ imposes a hydrostatic profile to the liquid pressure $p_L$. The cylindrical vessel containing the liquid is assumed large enough, so that the effect of the sidewalls is negligible. Also, the depth of the liquid is so big that the free surface is ignored. Mathematically, the liquid above the plate will be assumed to extend to infinity.

The steady initial interface shape, as well as its time of evolution are determined by $R$, $ρ_L$, $p_G$, $g$, by the hydrostatic pressure on the plate $p_{L0}$, and also by the surface tension $σ$. It is possible to identify two independent dimensionless parameters, which are supposed to control this evolution. The Eötvös number, which is also known as the Bond number [11], measures the importance of the gravitational force compared to the surface tension force:

$$E_ö = gρ_L R^2 / σ. \qquad (1)$$

The second control parameter is the dimensionless excess pressure:

$$\Delta = (p_G - p_{L0})R/σ, \qquad (2)$$

which characterizes the pressure difference between gas and liquid across the orifice. The Rayleigh-Taylor case, where the meniscus is initially flat, is one situation given by $Δ=0$; by analogy to the bending of an initially straight column, this case will be called the *perfect problem*. Bubble and drop formation is defined by $Δ > 0$, and $Δ < 0$, respectively. They correspond to initially bent column problems, and will be termed *imperfect problems*.

In some literature devoted to equilibrium shapes of drops and bubbles [20], it is quite common to introduce the radius of curvature at the apex of the bubble or drop, instead of the orifice radius in Equation (1). Then the differential equations giving the sought profiles are further modified, by shifting the origin to this apex. These two operations give equations containing a single resulting parameter, the *shape factor*. As mentioned by Rienstra [21], the simplification is only apparent since the distance of the apex to the plate is a part of the solution. Furthermore, the shape factor is not convenient to use when the profile is flat, and the radius of curvature is infinite.

In addition to the set of two parameters {$E_ö, Δ$} controlling the evolution of the interface, we might have introduced some Reynolds number expressing the influence of liquid viscosity. Neglecting such a refinement is justified, provided that the Ohnesorge number is Oh = $μ_L / \sqrt{2ρ_LσR} \ll 1$, where $μ_L$ is the liquid viscosity.

We choose a cylindrical polar coordinate system $(r, θ, z)$, with $r=0$ on the orifice axis, and $z=0$, the horizontal reference plane, at the plate level. Gravity is acting in the negative $z$-direction. Throughout this paper we restrict our attention to axisymmetric configurations with respect to the $z$-axis. For bubble and drop configurations close to the flat meniscus, interfaces have a simple (one-to-one) projection on the reference plane, so that a representation $z = ξ(r, t)$ is convenient.

The appropriate boundary conditions that govern the continuously changing contact line, which at time $t=0$ is attached to the orifice edge, will be simplified due to the following
common hypotheses [22]. First, as the orifice edge is rounded, the contact line can find a location on this edge, respecting the local physics of contact angles. Basically, this physics stipulates that, for one contact location, a region of equilibrium contact angles $\theta_{eq}$ exists [23] with: $\theta_r < \theta_{eq} < \theta_a$, where $\theta_r$ and $\theta_a$ are the limits of the receding and advancing contact angles. Otherwise a complex dynamic process occurs through which the contact line looks for another place of contact, where the above inequalities can be satisfied.

The interface configurations that will be considered in this study are those having a contact line that does not leave the round part of the edge. The corresponding permissible limits can be expressed in terms of the azimuthal angle $\psi_c$ at the contact point. When a bubble emerges over the plate, it can occur that no equilibrium contact angle can be found on the edge as soon as $\psi_c$ decreases to $\theta_r$: the contact line is then on the point of retreating. On the contrary, when a drop hangs under the plate, it can occur that no equilibrium contact angle can be found on the edge as soon as $(\psi_c - \pi)$ increases to $\theta_a$: the contact line is then on the point of advancing.

Secondly, we assume that the radius of curvature of the orifice edge is much smaller than the radius of the orifice. Thus, adopting a macroscopic point of view, the contact line can be considered as attached to the seemingly sharp-edged orifice defined by: $\xi(R, t) = 0$. All the detailed local physics of contact angles can be by-passed, provided that the studied interface configurations are such that $\theta_r < \psi_c < (\pi + \theta_a)$.

A small disturbance is applied to the system, and we wish to study the subsequent motion of the interface, and of the liquid. The imposed axial symmetry of the boundary conditions favours mainly axial disturbances to which our analysis will be restricted.

Now that we have a better definition of our model, we can list the questions that our study will address: (I) What is the stability of the null solution (perfect problem)? This is not more than the classical Rayleigh-Taylor problem, which will be given for sake of completeness. (II) What are the steady solutions that bifurcate from the null solution? More specifically, we shall give the values of Eö (or branch-points) from which at least a bifurcation emanates; a theoretical description of bifurcations (or branches) in the neighbourhood of the branch-points, i.e., number, dependence of the bifurcations on Eö. (III) What is the stability of the bifurcating solution? (IV) What are isolated solutions that break bifurcation? They correspond, of course, to emerging bubbles and pendant drops, and we shall give a theoretical description of these solutions in terms of $\Delta$ in the neighbourhood of the branch-points. (V) What is the stability of isolated solutions?

### 3. Governing equations

The motion in the incompressible inviscid liquid above the interface is assumed to be irrotational initially and remain so, with the velocity potential $\phi$, and velocity $u = \nabla\phi$. The velocity potential $\phi$ satisfies the Laplace equation:

$$ \nabla^2 \phi = 0, \quad \text{for} \quad \begin{cases} 0 \ge r < \infty \\ \xi(r,t) \ge z < \infty, & \text{if } r \le R, \\ 0 \ge z < \infty, & \text{if } r > R \end{cases} \qquad (3) $$

and its evolution is governed by the Bernoulli equation:

$$ \rho_L \left( \frac{\partial \phi}{\partial t} + \frac{1}{2} |\nabla \phi|^2 \right) + p_L + \rho_L g z = C(t), \quad \text{for} \quad \begin{cases} 0 \le r < \infty \\ \xi(r,t) \ge z < \infty, & \text{if } r \le R, \\ 0 \le z < \infty, & \text{if } r > R \end{cases} \qquad (4) $$
The velocity potential must be finite along the z-axis, and it vanishes at infinity for the considered axisymmetric problem:

$$ \frac{\partial \phi}{\partial r} = 0, \quad \text{for } (r=0, \xi(0,t) \le z < \infty), \qquad (5) $$

$$ \phi = 0, \quad \text{for } (r \to \infty, 0 \le z < \infty). \qquad (6) $$

The pressure there becomes the steady hydrostatic pressure. Thus, on the plate it is:

$$ p_L(r, 0, t) = C(t) = p_{L0}. \qquad (7) $$

On the interface, where $[1 + (\partial\xi/\partial r)^2]^{1/2} \mathbf{n} = (\xi_r, -1)$, the so-called kinematic condition applies:

$$ \nabla \phi \cdot \mathbf{n} = -\frac{\partial \xi}{\partial t} \left[ 1 + \left( \frac{\partial \xi}{\partial r} \right)^2 \right]^{-1/2}, \quad \text{or} \quad \frac{\partial \phi}{\partial z} = \frac{\partial \xi}{\partial t} + \frac{\partial \phi}{\partial r} \frac{\partial \xi}{\partial r}, \quad \text{for} \quad \begin{cases} 0 \le r < R \\ z = \xi(r,t), \end{cases} \qquad (8) $$

and the boundary condition on the impermeable plate is given by:

$$ \frac{\partial \phi}{\partial z} = 0, \quad \text{for } (R \le r < \infty, z=0). \qquad (9) $$

Here the momentum balance on the interface reduces to its normal component (the additional pressure jump due to surface tension is the Laplace-Young boundary condition):

$$ p_G - p_L = 2\sigma H, \quad \text{for } (0 \le r \le R, z = \xi(r,t)), \qquad (10) $$

where H is the mean curvature given by:

$$ H = -\frac{1}{2} \nabla \cdot \mathbf{n} = - \frac{r \partial^2 \xi / \partial r^2 + \left[ 1 + \left( \frac{\partial \xi}{\partial r} \right)^2 \right] \partial \xi / \partial r}{2r \left[ 1 + \left( \frac{\partial \xi}{\partial r} \right)^2 \right]^{3/2}}. \qquad (11) $$

Two additional boundary conditions for the interface equation must be specified:

$$ \frac{\partial \xi}{\partial r} = 0 \quad \text{at} \quad (r=0), \quad \text{and} \quad \xi = 0 \quad \text{at} \quad (r=R). \qquad (12) $$

Initial conditions are also prescribed:

$$ \xi(r, 0) = \xi_i(r), \quad \text{and} \quad \frac{\partial \xi}{\partial t}(r, 0) = \dot{\xi}_i(r), \quad \text{for } (0 \le r \le R, t=0). \qquad (13) $$

To get a better idea of the mathematical structure of the above problem, it would be interesting to express it solely in terms of $\xi(r,t)$. One way to do this, is first to solve the potential problem defined by Equations (3), (5), (6), (8), and (9), showing explicitly the dependence of $\phi$ on the various parts of the data. Assuming that the Green function $G_\xi$ exists for this partial problem [24, Volume II, pp. 130–145], and initial conditions for $\phi$ being put equal to 0, we may show that $\phi$ is a highly nonlinear functional of $\xi(r', t)$:

$$ \phi(r, z, t) = -2\pi \int_{0}^{R} r' G_{\xi}(r, z | r', \xi(r', t)) \frac{\partial \xi}{\partial t}(r', t) \left[ 1 + \left( \frac{\partial \xi}{\partial r'} \right)^2 \right]^{-1/2} dr'. \qquad (14) $$

Combining momentum balance and Bernoulli equation evaluated at the interface gives the following equation, written in dimensionless form:

$$ -\mathrm{E} \ddot{\mathrm{o}} \left( \frac{\partial \phi^*}{\partial t^*} + \frac{|\nabla^* \phi^*|^2}{2} \right) = \frac{\partial^2 \xi^*}{\partial r^{*2}} \left[ 1 + \left( \frac{\partial \xi^*}{\partial r^*} \right)^2 \right]^{-\frac{3}{2}} + \frac{1}{r^*} \frac{\partial \xi^*}{\partial r^*} \left[ 1 + \left( \frac{\partial \xi^*}{\partial r^*} \right)^2 \right]^{-\frac{1}{2}} + \Delta + \mathrm{E} \ddot{\mathrm{o}} \xi^*, \qquad (15) $$
where the length scale $R$, and velocity scale $U = \sqrt{gR}$ have been adopted. The dimensionless variables have been taken as:

$$ \xi^* = \xi/R, \quad r^* = r/R, \quad \phi^* = \phi/(UR), \quad \text{and} \quad t^* = Ut/R. \qquad (16) $$

Substituting the dimensionless form of Equation (14) for $\phi^*$ evaluated at $z^* = \xi^*(r^*, t^*)$ in Equation (15), yields a complex nonlinear integro-differential evolution equation governing $\xi^*(r^*, t^*)$. The left-hand side contains the *transient part*. It is well-defined only if the Green function $G_\xi$ is determined explicitly. In fact, such a task cannot be performed analytically, except if we limit ourselves to the study of small disturbances of the interface from a steady flat profile. This will be done in the next section, where the Rayleigh-Taylor stability problem will be addressed. The right-hand side contains the *steady part*; we shall denote this part by $\mathcal{F}(\text{Eö}, \xi^*, \Delta)$. We may observe that the flat profile $\xi^*=0$ (steady solutions are noted with a tilde) is the solution of:

$$ \mathcal{F}(\text{Eö}, \xi^*, 0) = 0. \qquad (17) $$

### 4. Stability of the null solution (Rayleigh-Taylor stability problem)

The particular static interface, for which the stability is tested in this section, is the flat interface $\xi^*=0$; it is given by $\Delta=0$. We linearize the general dynamic problem stated above by considering small perturbations $\hat{\xi}^*, \hat{\phi}^*$ about a given static equilibrium configuration:

$$ \xi^*(r^*, t^*) = \xi^*(r^*) + \hat{\xi}^*(r^*, t^*), \qquad (18) $$

$$ \phi^*(r^*, z^*, t^*) = \hat{\phi}^*(r^*, z^*, t^*). \qquad (19) $$

Following the procedure defined at the end of Section 3, a small disturbance $\hat{\xi}^*$ about $\xi^*$ is shown to satisfy the linear integro-differential equation:

$$ 2\pi \text{Eö} \int_{0}^{1} r'^{*} G_{0}(r^{*}, 0 | r'^{*}, 0) \frac{\partial^{2} \hat{\xi}^{*}}{\partial t'^{*2}}(r'^{*}, t^{*}) dr'^{*} = \mathcal{F}_{\xi}(\text{Eö}, 0, 0 | \hat{\xi}^{*}), \qquad (20) $$

where we have defined the linearized steady operator:

$$ \mathcal{F}_{\xi}(\text{Eö}, 0, 0 | \hat{\xi}^*) = \frac{\partial^2 \hat{\xi}^*}{\partial r^{*2}} + \frac{1}{r^*} \frac{\partial \hat{\xi}^*}{\partial r^*} + \text{Eö} \hat{\xi}^*, \qquad (21) $$

which is the Bessel operator of order 0.

The boundary and initial conditions (12) and (13) become:

$$ \frac{\partial \hat{\xi}^*}{\partial r^*} = 0, \quad \text{and} \quad \hat{\xi}^* \neq \infty, \quad \text{for } (r^* = 0, t^* > 0), \qquad (22) $$

$$ \hat{\xi}^* = 0, \quad \text{for } (r^* = 1, t^* > 0), \qquad (23) $$

$$ \hat{\xi}^* = \hat{\xi}_i^*, \quad \text{and} \quad \frac{\partial \hat{\xi}^*}{\partial t^*} = \hat{\xi}_i^*, \quad \text{for } (0 \le r^* \le 1, t^* = 0), \qquad (24) $$

For convenience asterisks will be suppressed in the remainder of this paper.

The Green function $G_0$ is defined by the following expression (see Appendix A):

$$ G_0(r, z | r', z') = \frac{1}{4\pi} \int_0^\infty J_0(\gamma r) J_0(\gamma r') \left[ e^{-|z-z'|\gamma} + e^{-|z+z'|\gamma} \right] d\gamma, \qquad (25) $$

where $J_0(\gamma r)$ is the Bessel function of first kind and order zero, with parameter $\gamma^2$.
Using an appropriate eigenfunction expansion solves the above stability problem. The corresponding eigenvalue problem [24, Volume I, p. 282] is based on the zero-order Bessel equation. Assuming a uniformly converging series, the solution of Equation (20) can be expanded in a complete set, {$Z_k(r)$}, with coefficients depending parametrically on $t$. Thus,

$$ \hat{\xi}(r,t) = \sum_{k=1}^{\infty} \xi^{(k)}(t) Z_k(r), \quad \text{where } \xi^{(k)}(t) = \int_0^1 r \hat{\xi}(r,t) Z_k(r) dr = \langle \hat{\xi}, Z_k \rangle. \qquad (26) $$

The set of functions:

$$ Z_k(r) = \sqrt{2} J_0(\lambda_k r) / J_1(\lambda_k) \qquad (27) $$

is orthonormal over the range $0<r<1$. Here $\lambda_k$ is the $k$th zero of $J_0$. In our case, the eigenvalue problem is self-adjoint. As a consequence, {$Z_k$} is an orthonormal system of eigenfunctions. Assuming uniform convergence, initial values of $\xi^{(k)}(t)$ result from (24):

$$ \xi^{(k)}(0) = \langle \hat{\xi}_i, Z_k \rangle, \quad \text{and} \quad \frac{\partial \xi^{(k)}(0)}{\partial t} = \langle \hat{\xi}_i, Z_k \rangle. \qquad (28) $$

We multiply both sides of the integro-differential Equation (20) by $(rZ_k(r))$, and integrate over the range $[0, 1]$ to obtain:

$$ \begin{aligned} & 2\pi \text{Eö} \sum_{j=1}^{\infty} \left[ \frac{d^2 \xi^{(j)}(t)}{dt^2} \int_0^1 r' \langle G_0(r, 0|r', 0), Z_k(r) \rangle Z_j(r') dr' \right] \\ & = \int_0^1 \frac{\partial}{\partial r} \left( r \frac{\partial \xi}{\partial r} \right) Z_k dr + \text{Eö} \langle \xi, Z_k \rangle. \end{aligned} \qquad (29) $$

Since

$$ \int_0^1 \frac{\partial}{\partial r} \left( r \frac{\partial \xi}{\partial r} \right) Z_k dr = \int_0^1 \xi \frac{\partial}{\partial r} \left( r \frac{\partial Z_k}{\partial r} \right) dr = -\lambda_k^2 \int_0^1 r \xi Z_k dr = -\lambda_k^2 \xi^{(k)}(t), \qquad (30) $$

and using the definition of Bessel functions, we finally have the differential equation:

$$ \frac{\text{Eö}}{2} [J_1(\lambda_k)]^2 \frac{d^2 \xi^{(k)}(t)}{dt^2} + (\lambda_k^2 - \text{Eö}) \xi^{(k)}(t) = 0, \qquad (31) $$

subject to the initial conditions (28). The coefficient of $(d^2\xi^{(k)}/dt^2)$ is viewed as an *added mass*. For small motions, the upper fluid imposes on the $k$th interface component some inertia, which can be characterized by a single mass coefficient $m_k = (\text{Eö}[J_1(\lambda_k)]^2)/2$. Since the upper fluid has been assumed inviscid, there is no friction term. The coefficient of $\xi^{(k)}$ has the role of a restoring force.

The most unstable component is obtained for $k=1$. Thus, the general condition for loss of stability must be derived from the equation controlling this component:

$$ m_1 \frac{d^2 \xi^{(1)}}{dt^2} - \mu \xi^{(1)} = 0, \qquad (32) $$

where $\mu = \text{Eö} - \lambda_1^2 = \text{Eö} - 5.783186$, which corresponds to the first eigenvalue of the axisymmetric problem; $\lambda_1 = 2.40482556$. Let us define $\xi^{(1)} = e^{\sigma^{(1)}t}x^{(1)}$ as a typical solution of the Equation (32). The corresponding eigenvalues $\sigma^{(1)}$ have to satisfy:

$$ m_1 (\sigma^{(1)})^2 x^{(1)} - \mu x^{(1)} = 0. \qquad (33) $$
<figure><img src="image_3.png" /><figcaption>Figure 2. Behaviour of the eigenvalues around the critical value $\mu=0$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Subcritical bifurcation diagram ($\Delta=0$) from analytical results (Equation 49).</figcaption></figure>

Thus, there are two eigenvalues:

$$ \sigma_{\pm}^{(1)}(\mu) = \pm \sqrt{\mu/m_1}. \qquad (34) $$

They are imaginary when $\mu < 0$, and real when $\mu \ge 0$. In the latter case, there is always a positive eigenvalue, and we conclude that $\tilde{\xi}=0$ is unstable. If $\mu \ge 0$ implies instability, stability may be obtained when $\mu < 0$, i.e., when $\sqrt{\mu}$ is purely imaginary. In fact, in this case, an extra negative real term should be present to dampen oscillations. This term would obviously come from viscous effects, which have been neglected.

At the critical value $\mu=0$, $\sigma_{\pm}^{(1)}=0$ is a double real eigenvalue. We observe that $\hat{\phi}$ and $\xi^{(1)} = \langle \hat{\xi}_i, Z_1 \rangle t + \langle \hat{\xi}_i, Z_1 \rangle$ are given via the Green function as in Equation (14). The behaviour of $\sigma_{\pm}^{(1)}(\mu)$, when $\mu$ increases around 0, is presented in Figure 2. When $\sigma_{\pm}^{(1)}$ coalesce at 0, other eigenvalues $\sigma_{\pm}^{(k)}$ are still on the imaginary axis. They follow the same path as $\sigma_{\pm}^{(1)}$ when $\mu$ keeps increasing.

## 5. Bifurcating solutions

There is a close connection between loss of stability of the linearized problem (i.e., between some critical values of the control parameters), and the condition for existence of a branch-point [25].

So we look for steady bifurcating solutions in the neighbourhood of $E\ddot{o} = \lambda_1^2 = 5.783186$ (or $\mu=0$), $\tilde{\xi}=0$, and $\Delta=0$. Let us expand the nonlinear operator $\mathcal{F}(\lambda_1^2 + \mu, \xi, \Delta)$ around the steady solution $\tilde{\xi}=0$, which satisfies Equation (17):

$$ \begin{aligned} \mathcal{F}(\lambda_1^2 + \mu, \xi, \Delta) = & \mathcal{F}_{\xi}(\lambda_1^2 + \mu, 0, \Delta |\xi) + \frac{1}{2}\mathcal{F}_{\xi\xi}(\lambda_1^2 + \mu, 0, \Delta |\xi |\xi) + \\ & + \frac{1}{6}\mathcal{F}_{\xi\xi\xi}(\lambda_1^2 + \mu, 0, \Delta |\xi |\xi |\xi) + \mathcal{O}(|\xi|^4) = 0. \end{aligned} \qquad (35) $$
To simplify the notation throughout the paper, we shall use $f$ for $\mathcal{F}$, and for all partial derivatives of $\mathcal{F}$ with respect to the indicated subscript, when these quantities are evaluated at the point $(\lambda_1^2, 0, 0)$. We obtain the following relations:

$$
\begin{align}
f_{\Delta} &= \mathcal{F}_{\Delta}(\lambda_1^2, 0, 0) = 1, & f_{\xi}(u) &= \mathcal{F}_{\xi}(\lambda_1^2, 0, 0|u) = \frac{d^2 u}{dr^2} + \frac{1}{r} \frac{du}{dr} + \lambda_1^2 u, \nonumber \\
f_{\xi\mu}(u) &= \mathcal{F}_{\xi\mu}(\lambda_1^2, 0, 0|u) = u, & f_{\xi\Delta}(u) &= \mathcal{F}_{\xi\Delta}(\lambda_1^2, 0, 0|u) = 0, \nonumber \\
f_{\xi\xi}(u|v) &= \mathcal{F}_{\xi\xi}(\lambda_1^2, 0, 0|u|v) = 0, & f_{\xi\xi\mu}(u|v) &= \mathcal{F}_{\xi\xi\mu}(\lambda_1^2, 0, 0|u|v) = 0, \nonumber \\
f_{\xi\xi\Delta}(u|v) &= \mathcal{F}_{\xi\xi\Delta}(\lambda_1^2, 0, 0|u|v) = 0, & f_{\xi\xi\xi}(u|v|w) &= \mathcal{F}_{\xi\xi\xi}(\lambda_1^2, 0, 0|u|v|w) \nonumber \\
&= -3 \left( \frac{d^2 u}{dr^2} \frac{dv}{dr} \frac{dw}{dr} + \frac{d^2 v}{dr^2} \frac{dw}{dr} \frac{du}{dr} + \frac{d^2 w}{dr^2} \frac{du}{dr} \frac{dv}{dr} \right) - \frac{3}{r} \frac{du}{dr} \frac{dv}{dr} \frac{dw}{dr}. \tag{36}
\end{align}
$$

Since $\mathcal{F}$ is linear with respect to $\mu$ and $\Delta$, $f_{\xi\xi\dots\xi\mu} = f_{\xi\xi\dots\xi\Delta} = 0$ for any $n \ge 2$. Moreover, we have $f_{\xi\xi\dots\xi} = 0$. Steady solutions that bifurcate will be constructed as a power series in the solution amplitude $\varepsilon$, that is:

$$
\tilde{\xi}(\varepsilon) = \varepsilon \tilde{\xi}_1 + \frac{\varepsilon^2}{2} \tilde{\xi}_2 + \frac{\varepsilon^3}{6} \tilde{\xi}_3 + O(|\varepsilon|^4), \qquad (37)
$$

$$
\mu(\varepsilon) = \varepsilon\mu_1 + \frac{\varepsilon^2}{2}\mu_2 + \frac{\varepsilon^3}{6}\mu_3 + O(|\varepsilon|^4). \qquad (38)
$$

The amplitude

$$
\varepsilon = \langle \tilde{\xi}(r), Z_1 \rangle = \int_0^1 r \tilde{\xi}(r) Z_1(r) dr = \sqrt{2} \int_0^1 r \tilde{\xi} \frac{J_0(\lambda_1 r)}{J_1(\lambda_1)} dr \quad (39)
$$

is defined by a projection on the eigensubspace associated with the adjoint eigenvector $Z_1(r)$, belonging to the eigenvalue $\mu=0$ of the problem:

$$
f_{\xi}(\tilde{\xi}) + \mu f_{\xi\mu}(\tilde{\xi}) = 0. \qquad (40)
$$

The amplitude definition may be regarded as a normalization, which implies some constraints on $\tilde{\xi}_k$ appearing in Equation (37):

$$
\langle \tilde{\xi}_1, Z_1 \rangle = 1 \quad \text{and} \quad \langle \tilde{\xi}_k, Z_1 \rangle = 0 \quad \text{for } k > 1. \qquad (41)
$$

Physically, $\varepsilon$ can be regarded as a scalar, which characterizes the bubble ($\varepsilon > 0$), or the drop ($\varepsilon < 0$) magnitude. We may try to obtain $\tilde{\xi}_k$, and $\mu_k$ ($k \ge 1$), by expanding the left-hand side of $\mathcal{F}(\lambda_1^2 + \mu(\varepsilon), \tilde{\xi}(\varepsilon), 0) = 0$ in powers of $\varepsilon$, and identifying independent powers of $\varepsilon$, that is:

$$
f_{\xi}(\tilde{\xi}_1) = 0, \qquad (42)
$$

$$
f_{\xi}(\tilde{\xi}_2) + 2\mu_1 f_{\xi\mu}(\tilde{\xi}_1) = 0, \qquad (43)
$$

$$
f_{\xi}(\tilde{\xi}_3) + 3\mu_1 f_{\xi\mu}(\tilde{\xi}_2) + 3\mu_2 f_{\xi\mu}(\tilde{\xi}_1) + f_{\xi\xi\xi}(\tilde{\xi}_1|\tilde{\xi}_1|\tilde{\xi}_1) = 0, \qquad (44)
$$

where relations (36) have been taken into account.

Since $\langle \tilde{\xi}_1, Z_1 \rangle = 1$, the eigenvalue problem (42) has a single solution, namely

$$
\tilde{\xi}_1 = Z_1 = \sqrt{2} J_0(\lambda_1 r) / J_1(\lambda_1). \qquad (45)
$$
The Fredholm alternative applies for the solvability of the other problems (43) and (44). Consider first Equation (43). Since $\mu_1 \langle f_{\xi\mu}(\tilde{\xi}_1), Z_1 \rangle = \mu_1 \langle Z_1, Z_1 \rangle = 0$, and $\langle Z_1, Z_1 \rangle = 1$, we have $\mu_1 = 0$. The resulting equation $f_{\xi}(\tilde{\xi}_2) = 0$ may then be solved. It has also a single solution: $\tilde{\xi}_2 = 0$, since $\langle \tilde{\xi}_2, Z_1 \rangle = 0$. For the resulting problem (44), the compatibility condition gives:

$$ \mu_2 = -\frac{1}{3} \frac{\langle f_{\xi\xi\xi}(Z_1|Z_1|Z_1), Z_1 \rangle}{\langle f_{\xi\mu}(Z_1), Z_1 \rangle} = -73.4685 \quad (46) $$

To obtain this result we have used:

$$ \begin{aligned} \langle f_{\xi\xi\xi}(Z_1|Z_1|Z_1), Z_1 \rangle &= \frac{12\lambda_1^3}{J_1^4(\lambda_1)} \left[ -2 \int_0^1 J_0(\lambda_1 r) J_1^3(\lambda_1 r) dr + 3\lambda_1 \int_0^1 r J_0^2(\lambda_1 r) J_1^2(\lambda_1 r) dr \right] \\ &= 220.4055, \end{aligned} \quad (47) $$

while the denominator is $\langle f_{\xi\mu}(Z_1), Z_1 \rangle = \langle Z_1, Z_1 \rangle = 1$. The associated equation becomes with $\langle \tilde{\xi}_3, Z_1 \rangle = 0$:

$$ f_{\xi}(\tilde{\xi}_3) = -3\mu_2 Z_1 - 9Z_1 \left(\frac{dZ_1}{dr}\right)^2 - \frac{6}{r} \left(\frac{dZ_1}{dr}\right)^3. \quad (48) $$

So, two distinct branches emanate symmetrically from $E_0 = \lambda_1^2$ (or $\mu = 0$):

$$ \mu = \frac{1}{2}\mu_2\epsilon^2 + O(|\epsilon|^4) = -36.73425\epsilon^2 + O(|\epsilon|^4). \quad (49) $$

It is easy to show that only even powers of $\epsilon$ can occur in Equation (49). Conversely, concerning the $\tilde{\xi}$-expansion, only odd powers of $\epsilon$ occur, showing that bubble profiles and drop profiles are symmetrical with respect to the orifice.

As $\mu_2$ is negative, the branching is to the left; as this branching is to the stable side, the system is subcritical. Obviously it is quite difficult to compute further terms in Equation (48) analytically.

A bifurcation diagram restricted to the neighbourhood of $\epsilon = 0$, including the analytical results (Equation 49), is given in Figure 3. The upper branch corresponds to bubbles, the lower to drops. The flat meniscus ($\epsilon = 0$) is stable before the critical point $P_1$, where $\mu < 0$. The solution $\tilde{\xi} = 0$ loses its stability when $\mu = (E_0 - \lambda_1^2)$ is increased beyond zero. New solutions $\mu(\epsilon^2)$, corresponding to curved menisci (bubbles, or drops), undergo double-point bifurcation at the singular point $(\mu, \epsilon) = (0, 0)$ denoted by $P_1$.

## 6. Stability of the bifurcating solution

We study the stability of the steady bifurcating solution we have just found. In Section 4 we showed that, if $\mu$ is increased along the branch $\epsilon = 0$, criticality is reached at $\mu = 0$; a conjugate pair of purely imaginary eigenvalues gives rise to two real eigenvalues, one having a strictly positive part (instability). Instead of increasing $\mu$ from $\mu = 0$, let us increase $|\epsilon|$ along the other two branches, and examine the behaviour of this pair. Such bifurcating solutions $(\mu(\epsilon), \tilde{\xi}(\epsilon))$ satisfy:

$$ \mathcal{F}(\lambda_1^2 + \mu(\epsilon), \tilde{\xi}(\epsilon), 0) = 0. \quad (50) $$

Instead of formulating the fully linearized problem for the small perturbations $\tilde{\xi}, \hat{\phi}$ about the above solutions as in Section 4, let us only consider the associated spectral problem for
$X$ and $F$ defined by:

$$ \xi(\varepsilon; r, t) = \tilde{\xi}(\varepsilon; r) + \hat{\xi}(\varepsilon; r, t) = \tilde{\xi}(\varepsilon; r) + e^{\sigma^{(\varepsilon)t}} X(\varepsilon; r), \quad (51) $$

$$ \phi(\varepsilon; r, z, t) = \hat{\phi}(\varepsilon; r, z, t) = e^{\sigma^{(\varepsilon)t}} F(\varepsilon; r). \quad (52) $$

This spectral (eigenvalue) problem can be decomposed into two parts, one for $F$, the other for $X$. The first part for $F$ is a classical Laplace problem with Neumann boundary conditions; it is complicated by the dependence of $F$ on the shape of the mathematical domain, which is modified by the moving boundary $z = \xi(\varepsilon; r, t)$. The second part, involving $F$ and $X$, results from the interfacial momentum balance incorporating the Bernoulli equation evaluated at the moving boundary. Following a procedure similar to the one defined at the end of Section 3, which collects all the preceding equations, we obtain an extended eigenvalue problem for the single function $X(r)$:

$$ \sigma^2(\varepsilon) K(\varepsilon; X(r')) = \mathcal{F}_{\tilde{\xi}}(\lambda_1^2 + \mu(\varepsilon), \tilde{\xi}(\varepsilon; r), 0 | X(r)), \quad (53) $$

where

$$ K(\varepsilon; X(r')) = \left\langle 2\pi(\lambda_1^2 + \mu(\varepsilon)) G_{\tilde{\xi}}(r, \tilde{\xi}(\varepsilon; r)|r', \tilde{\xi}(\varepsilon; r')) \times \left(1 + \left(\frac{\partial \tilde{\xi}(\varepsilon; r')}{\partial r'}\right)^2\right)^{-\frac{1}{2}} , X(r') \right\rangle. \quad (54) $$

The sign of the real part of $\sigma^{(1)}(\varepsilon)$, which is the eigenvalue with the largest real part, determines the stability of $\tilde{\xi}(\varepsilon; r)$. It can be shown that there is always a positive eigenvalue in the neighbourhood of $\varepsilon=0$, i.e.,

$$ \sigma^{(1)}(\varepsilon) = \pm 9.71\varepsilon + O(|\varepsilon|^2), \quad (55) $$

so we conclude that $\tilde{\xi}(\varepsilon; r)$ is unstable as the solution bifurcates at the critical point.

## 7. Isolated solutions that break bifurcation and their stability

We consider here the imperfect problem, i.e., $\Delta \neq 0$. Using the simplified notations defined in Section 5, we can write the steady equation in the neighbourhood of $(\lambda_1^2, 0, 0)$ as

$$ \begin{aligned} \mathcal{F}(\lambda_1^2 + \mu, \tilde{\xi}, \Delta) = \Delta f_{\Delta} + f_{\tilde{\xi}}(\tilde{\xi}) + \mu f_{\tilde{\xi}\mu}(\tilde{\xi}) &+ \frac{1}{6} f_{\tilde{\xi}\xi\tilde{\xi}}(\tilde{\xi}|\tilde{\xi}|\tilde{\xi}) + \frac{\mu}{6} f_{\tilde{\xi}\xi\tilde{\xi}\mu}(\tilde{\xi}|\tilde{\xi}|\tilde{\xi}) + \\ &+ \mathcal{R}(\mu, \tilde{\xi}) = 0. \end{aligned} \quad (56) $$

In this expansion we have used the Equations (35), and (36). The remainder $\mathcal{R}$ is:

$$ \mathcal{R}(\mu, \tilde{\xi}) = (\lambda_1^2 + \mu) \sum_{n=2}^{\infty} a_n \left( \frac{\partial \tilde{\xi}}{\partial r} \right)^{2n} \tilde{\xi}, \quad (57) $$

where $a_2 = 3/8$, and $a_n = 3(-1)^n(2n-5)/2^nn!$, with $n \ge 3$. We have seen that $f_{\tilde{\xi}}(u) = 0$ at $(\mu, \varepsilon) = (0, 0)$, where $\mathcal{F}(\lambda_1^2 + \mu, \tilde{\xi}, 0) = 0$. Due to the Equation (45), the following key condition
for the remainder of this section is guaranteed:

$$ \langle f_{\Delta}, Z_1 \rangle = \frac{\sqrt{2}}{J_1(\lambda_1)} \int_0^1 r J_0(\lambda_1 r) dr = \frac{\sqrt{2}}{\lambda_1} \neq 0. \qquad (58) $$

The amplitude $\epsilon$ that we have already defined for $\Delta=0$ in Equation (39) is extended here for $\Delta \neq 0$. In that context of imperfection theory, we follow Iooss and Joseph [26, Chapter 6, pp. 87–138], and invoke condition (58), in conjunction with the implicit-function theorem, claim the existence of a smooth solution $\tilde{\xi}(\mu, \epsilon)$, and $\Delta(\mu, \epsilon)$ of:

$$ \mathcal{F}(\lambda_1^2 + \mu, \tilde{\xi}(\mu, \epsilon), \Delta(\mu, \epsilon)) = 0, \qquad (59) $$

$$ \epsilon = \langle \tilde{\xi}(\mu, \epsilon), Z_1 \rangle, \qquad (60) $$

as well as the values of these functions at $\epsilon=0$, i.e., $\tilde{\xi}(\mu, 0)=0$, and $\Delta(\mu, 0)=0$.

We expand $\tilde{\xi}=\tilde{\xi}(\mu, \epsilon)$, and $\Delta=\Delta(\mu, \epsilon)$ as series around $(0, 0)$ up to the third order:

$$ \begin{aligned} \tilde{\xi} = & \tilde{\xi}_{\mu\mu} + \tilde{\xi}_{\epsilon\epsilon} + \frac{1}{2}(\tilde{\xi}_{\mu\mu}\mu^2 + 2\tilde{\xi}_{\mu\epsilon}\mu\epsilon + \tilde{\xi}_{\epsilon\epsilon}\epsilon^2) + \\ & + \frac{1}{6}(\tilde{\xi}_{\mu\mu\mu}\mu^3 + 3\tilde{\xi}_{\mu\mu\epsilon}\mu^2\epsilon + 3\tilde{\xi}_{\mu\epsilon\epsilon}\mu\epsilon^2 + \tilde{\xi}_{\epsilon\epsilon\epsilon}\epsilon^3) + O((|\mu|+|\epsilon|)^4), \end{aligned} \qquad (61) $$

$$ \begin{aligned} \Delta = & \Delta_{\mu\mu} + \Delta_{\epsilon\epsilon} + \frac{1}{2}(\Delta_{\mu\mu}\mu^2 + 2\Delta_{\mu\epsilon}\mu\epsilon + \Delta_{\epsilon\epsilon}\epsilon^2) + \\ & + \frac{1}{6}(\Delta_{\mu\mu\mu}\mu^3 + 3\Delta_{\mu\mu\epsilon}\mu^2\epsilon + 3\Delta_{\mu\epsilon\epsilon}\mu\epsilon^2 + \Delta_{\epsilon\epsilon\epsilon}\epsilon^3) + O((|\mu|+|\epsilon|)^4). \end{aligned} \qquad (62) $$

Owing to the above-mentioned values of these functions at $\epsilon=0$, we can write:

$$ \tilde{\xi}_{\mu}(0, 0) = \tilde{\xi}_{\mu\mu}(0, 0) = \tilde{\xi}_{\mu\mu\mu}(0, 0) = 0, \qquad (63) $$

$$ \Delta_{\mu}(0, 0) = \Delta_{\mu\mu}(0, 0) = \Delta_{\mu\mu\mu}(0, 0) = 0. \qquad (64) $$

Inserting Equations (61), and (62) into (56), and identifying powers of $\epsilon^i$, $\mu^i$, $\epsilon^i\mu^j$, $(i \ge 1, j \ge 1)$, we find to first order:

$$ f_{\xi}(\tilde{\xi}_{\epsilon}) = -\Delta_{\epsilon}, \qquad (65) $$

to second order:

$$ f_{\xi}(\tilde{\xi}_{\mu\epsilon}) = -\Delta_{\mu\epsilon} - f_{\xi\mu}(\tilde{\xi}_{\epsilon}), \qquad (66) $$

$$ f_{\xi}(\tilde{\xi}_{\epsilon\epsilon}) = -\Delta_{\epsilon\epsilon}, \qquad (67) $$

and to third order:

$$ f_{\xi}(\tilde{\xi}_{\mu\mu\epsilon}) = -\Delta_{\mu\mu\epsilon} - 2f_{\xi\mu}(\tilde{\xi}_{\mu\epsilon}), \qquad (68) $$

$$ f_{\xi}(\tilde{\xi}_{\mu\epsilon\epsilon}) = -\Delta_{\mu\epsilon\epsilon} - f_{\xi\mu}(\tilde{\xi}_{\epsilon\epsilon}), \qquad (69) $$

$$ f_{\xi}(\tilde{\xi}_{\epsilon\epsilon\epsilon}) = -\Delta_{\epsilon\epsilon\epsilon} - f_{\xi\xi\xi}(\tilde{\xi}_{\epsilon}|\tilde{\xi}_{\epsilon}|\tilde{\xi}_{\epsilon}). \qquad (70) $$

Again the Fredholm alternative is used to insure the solvability of the Equation (65), and successively of all the other differential Equations (66)–(70). Moreover, the normalization
Equation (39) imposes

$$ \begin{aligned} \langle \tilde{\xi}_{\varepsilon}, Z_1 \rangle &= 1, \quad \langle \tilde{\xi}_{\mu\varepsilon}, Z_1 \rangle &= \langle \tilde{\xi}_{\varepsilon\varepsilon}, Z_1 \rangle = 0, \quad \text{and} \\ \langle \tilde{\xi}_{\mu\mu\varepsilon}, Z_1 \rangle &= \langle \tilde{\xi}_{\mu\varepsilon\varepsilon}, Z_1 \rangle = \langle \tilde{\xi}_{\varepsilon\varepsilon\varepsilon}, Z_1 \rangle = 0. \end{aligned} \quad (71) $$

Thus, considering the first-order Equation (65), we obtain:

$$ \tilde{\xi}_{\varepsilon} = Z_1, \quad \text{and} \quad \Delta_{\varepsilon} = 0. \quad (72) $$

In the second-order Equation (66), we have $\tilde{\xi}_{\mu\varepsilon} = 0$, and the insertion of (72) into the solvability condition gives:

$$ \Delta_{\mu\varepsilon} = - \frac{\langle f_{\tilde{\xi}\mu}(Z_1), Z_1 \rangle}{\langle 1, Z_1 \rangle} = - \frac{\langle Z_1, Z_1 \rangle}{\langle 1, Z_1 \rangle} = - \frac{\lambda_1}{\sqrt{2}} = -1.70046. \quad (73) $$

The solution of this equation may be obtained by the eigenfunction method

$$ \tilde{\xi}_{\mu\varepsilon} = - \sum_{k=2}^{\infty} \frac{\lambda_1 J_1(\lambda_k) Z_k}{\lambda_k(\lambda_k^2 - \lambda_1^2) J_1(\lambda_1)}. \quad (74) $$

In the second-order Equation (67), we observe that $\tilde{\xi}_{\varepsilon\varepsilon} = 0$, and $\Delta_{\varepsilon\varepsilon} = 0$. Now the solvability conditions of the three third-order equations are simply, on account of the preceding results:

$$ \Delta_{\mu\mu\varepsilon} = - \frac{2 \langle f_{\tilde{\xi}\mu}(\tilde{\xi}_{\mu\varepsilon}), Z_1 \rangle}{\langle 1, Z_1 \rangle} = - \frac{2 \langle \tilde{\xi}_{\mu\varepsilon}, Z_1 \rangle}{\langle 1, Z_1 \rangle} = 0, \quad (75) $$

$$ \Delta_{\mu\varepsilon\varepsilon} = - \frac{\langle f_{\tilde{\xi}\mu}(\tilde{\xi}_{\varepsilon\varepsilon}), Z_1 \rangle}{\langle 1, Z_1 \rangle} = - \frac{\langle \tilde{\xi}_{\varepsilon\varepsilon}, Z_1 \rangle}{\langle 1, Z_1 \rangle} = 0, \quad (76) $$

$$ \begin{aligned} \Delta_{\varepsilon\varepsilon\varepsilon} &= - \frac{\langle f_{\tilde{\xi}\xi\xi}(Z_1|Z_1|Z_1), Z_1 \rangle}{\langle 1, Z_1 \rangle} = - \frac{\lambda_1}{\sqrt{2}} \langle f_{\tilde{\xi}\xi\xi}(Z_1|Z_1|Z_1), Z_1 \rangle \\ &= \frac{12\lambda_1^4}{\sqrt{2}J_1^4(\lambda_1)} \left[ 2 \int_0^1 J_0(\lambda_1 r) J_1^3(\lambda_1 r) dr - 3\lambda_1 \int_0^1 r J_0^2(\lambda_1 r) J_1^2(\lambda_1 r) dr \right] \\ &= -374.79173. \end{aligned} \quad (77) $$

Computing $\tilde{\xi}_{\mu\mu\varepsilon}$, $\tilde{\xi}_{\mu\varepsilon\varepsilon}$, and $\tilde{\xi}_{\varepsilon\varepsilon\varepsilon}$ does not present any difficulty, and will not be presented here. What is worth mentioning is the final form of $\Delta$, which reduces to:

$$ \Delta(\mu, \varepsilon) = \Delta_{\mu\varepsilon}\mu\varepsilon + \frac{\Delta_{\varepsilon\varepsilon\varepsilon}}{6}\varepsilon^3 + O\left(|\varepsilon|(|\mu| + |\varepsilon|)^3\right). \quad (78) $$

With (73) and (77), the Equation (78) becomes:

$$ \Delta(\mu, \varepsilon) = -1.70046\mu\varepsilon - 62.46528\varepsilon^3. \quad (79) $$

This analytical steady solution defines the equilibrium 3D surface, which is plotted in Figure 4. In the terminology of Catastrophe Theory we have identified a cusp. The upper cusp sheet corresponds to bubbles, the lower to drops.
<figure><img src="image_5.png" /><figcaption>Figure 4. Cusp near the critical point P₁. Analytical results for μ ∈ [-3.9, 1.2], ε ∈ [-0.32, 0.32].</figcaption></figure>

For some typical μ-values, cross-sections of the equilibrium surface (μ, Δ, ε) are plotted in Figure 5. For μ < 0, there is a stable branch that crosses the μ-axis between the turning point P₁⁻ for drops, and P₁⁺ for bubbles (solid line in Figure 5).

For ε-values larger than the one corresponding to P₁⁺, and for ε-values smaller than the one corresponding to P₁⁻, the branches are unstable (dashed lines in Figure 5). For the particular value μ = 0, the whole bubble branch (ε > 0), and the whole drop branch (ε < 0) are unstable, the cross-section passing beyond the critical point P₁.

The turning points P₁⁺ and P₁⁻ extend into Δ > 0, and Δ < 0 regions as singularity lines (fold curves) of the cusp. So, the intermediate cusp sheet is stable, being delimited by those two singularity lines. The upper and lower cusp sheets are unstable. The bubble and drop configurations, corresponding to the intermediate sheet, are the ones that are usually observed. The intersection of the 3D equilibrium surface with the plane Δ = 0 returns the afore-described Rayleigh-Taylor problem.

The stability study of isolated solutions (Δ≠0), which perturb bifurcating solutions (Δ=0), will be based on the stability study of the latter, presented in Section 6. The relations μ = μ(ε), and ξ̃ = ξ̃(ε) for Δ=0, will be replaced by the relations Δ = Δ(μ, ε), and ξ̃ = ξ̃(μ, ε), in which μ is a fixed parameter such that μ < 0.

The bifurcation curve near the critical point P₁ is given by:

$$ \Delta = \Delta(\mu, \varepsilon), \quad \text{and} \quad 0 = \frac{\partial \Delta}{\partial \varepsilon}(\mu, \varepsilon), \qquad (80) $$

where ε plays the role of a curve parameter.

This curve is composed of two branches through which the number of steady solutions varies when a turning point, P₁⁺ for bubbles, or P₁⁻ for drops, is reached. These two branches merge at the bifurcation point P₁ at μ = 0, ε = 0, and Δ = 0. Each branch is given by the system issued from Equations (79) and (80):

$$ \left\{ \begin{array}{l} \Delta = 124.93056\varepsilon^3 \\ \mu = -110.20302\varepsilon^2 \end{array} \right. . \qquad (81) $$
<figure><img src="image_6.png" /><figcaption>Figure 5. Cross sections of the equilibrium surface for $\mu=0$, $\mu=-1.5$, and $\mu=-3$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Analytical bifurcation set near the critical point $P_1$. The upper branch corresponds to the turning point $P_1^+$ of bubble cusp sheet; the lower branch corresponds to the turning point $P_1^-$ of drop cusp sheet.</figcaption></figure>

The analytical bifurcation curve near the critical point $P_1$ is expanded in Figure 6 on the control parameter plane $(\mu, \Delta)$.

## 8. Conclusions

A systematic theoretical analysis has afforded a unified picture of the Rayleigh-Taylor instability, and of the quasi-steady formation of bubbles and drops, under constant-pressure conditions. Via $\Delta$, the excess pressure number ($\Delta > 0$ for bubbles, and $\Delta < 0$ for drops), these latter processes are viewed in this paper as perturbing the Rayleigh-Taylor physical situation, which appears as a perfect bifurcating problem with $\Delta=0$. In this idealized problem, the flat interface is stretched according to a value measured by the Eötvös number $E_ö$ (or $\mu=E_ö/5.783186$). When $E_ö$ increases, the stretching decreases to some critical value corresponding to $E_ö_c=5.783186$ (or $\mu=0$). Beyond $E_ö_c$, the flat interface is unstable, and in its neighbourhood no steady stable solution is available (a transient process necessarily takes place). For a fixed gas/liquid couple, the critical value $E_ö_c$ has a practical meaning, allowing computation of the maximum orifice radius, which sustains a meniscus in equilibrium.

On the contrary, before reaching the critical value $E_ö_c$, at $\Delta \neq 0$, two sets of steady symmetric deflections $\xi(r)$, one upward (bubble), and one downward (drop), are theoretically available to the interface. These solutions, whose magnitudes can be described by some appropriate scalar $\epsilon = \langle \xi, Z_1 \rangle$, emerge from the critical point $P_1$ to the stable side. They correspond to a subcritical bifurcation.

An explicit picture of the 3D equilibrium surface near the critical point has been given when the excess pressure number $\Delta$ is added to $E_ö$ as an extra control parameter. A cusp has been identified. The turning points $P_1^+$ and $P_1^-$ extend into $\Delta > 0$, and $\Delta < 0$ regions as singularity lines (fold curves) of that cusp. The intermediate sheet of the cusp is stable, while its two upper and lower sheets are unstable. The intersection of the 3D equilibrium surface with the plane $\Delta=0$ returns the afore-described Rayleigh-Taylor problem.

The bubble and drop configurations corresponding to the intermediate cusp sheet are the ones that are usually observed in experiments [27,28]. They are steady, stable curved menisci and represent the first stage of the interface growing evolution [22,29]. Within this first stage, at a given Eötvös number (or $\mu$-value), the gas pressure increases for bubbles from the hydrostatic pressure
$p_{L0}$ (or $\Delta=0$), to its maximum value where $\Delta = \Delta_{max}$. For drops, the gas pressure decreases from $p_{L0}$ to its minimum value where $\Delta = \Delta_{min}$. The gas excess pressure has equal absolute extreme values for bubbles and drops, $\Delta_{max} = |\Delta_{min}|$, at the turning points $P_1^+$ and $P_1^-$ (Figure 5).

The practical value of the proposed description is limited since all the above results are known from the literature. However, it has two important aspects:

- it has an educational value by allowing a unified treatment of separate classical descriptions of interface instabilities. The role of the nondimensional parameters is clearly identified;

- it allows a sound basis for future studies. In fact, it can be shown numerically [30, Chapter 3] that below the critical Eötvös number the upper and lower sheets keep on zigzagging away, and that between the turning points, stable equilibrium states exist, and correspond to bubbles or drops having an undulatory structure. In [31] (in preparation), such computations are discussed, and experimental evidence of the physical reality of some mathematical solutions we have exhibited presented.

## Appendix A. Green's function for the potential of the liquid velocity field

In cylindrical co-ordinates, the Green function $G_0(r, z|r_0, z_0)$ that gives the dependence of $\phi$ on the flat interface velocity evolution, satisfies the equation:

$$ -\frac{\partial}{\partial r} \left( r \frac{\partial G_0}{\partial r} \right) - r \frac{\partial^2 G_0}{\partial z^2} = \frac{1}{2\pi} \delta(r-r_0) \delta(z-z_0), \quad (0<r,r_0<\infty, 0<z,z_0<\infty), \quad (A1) $$

with boundary conditions:

$$ \frac{\partial G_0}{\partial z} = 0, \quad \text{for } (0 < r, r_0 < \infty, z = 0, 0 < z_0 < \infty), \quad (A2) $$

$$ \frac{\partial G_0}{\partial r} = 0, \quad \text{for } (r = 0, 0 < r_0 < \infty, 0 < z, z_0 < \infty), \quad (A3) $$

$$ \lim_{z \to \infty} G_0 = 0, \quad \text{for } (0 < r, r_0 < \infty, 0 < z_0 < \infty), \quad (A4) $$

$$ \lim_{r \to \infty} G_0 = 0, \quad \text{for } (0 < r_0 < \infty, 0 < z, z_0 < \infty). \quad (A5) $$

Multiplying the Equation (A1) by $J_0(\gamma r)$, and integrating from $r=0$ to $\infty$, we obtain:

$$ - \int_{0}^{\infty} J_{0}(\gamma r) \frac{\partial}{\partial r} \left( r \frac{\partial G_{0}}{\partial r} \right) dr - \frac{d^2 G_{0H}}{dz^2} = \frac{1}{2\pi} J_{0}(\gamma r_0) \delta(z - z_0), \quad (A6) $$

where

$$ G_{0H}(\gamma, z|r_0, z_0) = \int_{0}^{\infty} r J_{0}(\gamma r) G_{0}(r, z|r_0, z_0) dr \quad (A7) $$

is the Hankel transform of order zero of $G_0$, with $r, r_0$, and $z_0$ as parameters. Integrating by parts twice, we obtain:

$$ \int_{0}^{\infty} J_{0}(\gamma r) \frac{\partial}{\partial r} \left( r \frac{\partial G_{0}}{\partial r} \right) dr = - \int_{0}^{\infty} r \frac{\partial G_{0}}{\partial r} \frac{dJ_{0}(\gamma r)}{dr} dr = + \int_{0}^{\infty} G_{0} \frac{d}{dr} \left( r \frac{dJ_{0}(\gamma r)}{dr} \right) dr, \quad (A8) $$

where the integrated terms are each time equal to zero.
Since $J_0(\gamma r)$ satisfies the Bessel equation of order zero with parameter $\gamma^2$:

$$ \int_0^\infty J_0(\gamma r) \frac{\partial}{\partial r} \left( r \frac{\partial G_0}{\partial r} \right) dr = -\gamma^2 G_{0H}, \qquad (A9) $$

then (A6) becomes:

$$ \frac{d^2 G_{0H}}{dz^2} - \gamma^2 G_{0H} = -\frac{1}{2\pi} J_0(\gamma r_0) \delta(z - z_0), \quad \text{for } (r_0 < \infty, 0 < z < z_0 < \infty). \qquad (A10) $$

The boundary conditions for this simple differential equation are:

$$ \frac{dG_{0H}}{dz} = 0, \quad \text{for } (0 < r_0 < \infty, z = 0, 0 < z_0 < \infty), \qquad (A11) $$

$$ \lim_{z \to \infty} G_{0H} = 0, \quad \text{for } (0 < r_0 < \infty, 0 < z_0 < \infty), \qquad (A12) $$

according to Equations (A2), (A4), and (A7).

The method of images can be used to construct the transformed Green's function for a semi-infinite positive axis, when the left end of this axis is insulated. Thus, $G_{0H}$ is found by placing a positive image source at $z = -z_0$:

$$ G_{0H}(\gamma, z|r_0, z_0) = \frac{J_0(\gamma r)}{4\pi\gamma} \left[ e^{-|z-z_0|\gamma} + e^{-|z+z_0|\gamma} \right], \qquad (A13) $$

and by the inversion formula we obtain:

$$ G_0(r, z|r_0, z_0) = \frac{1}{4\pi} \int_0^\infty J_0(\gamma r) J_0(\gamma r_0) \left[ e^{-|z-z_0|\gamma} + e^{-|z+z_0|\gamma} \right] d\gamma. \qquad (A14) $$

## References

1. D.H. Sharp, An overview of Rayleigh-Taylor instability. *Physica* D12 (1984) 3–18.

2. J.C. Maxwell, Capillary Action. In: *Encyclopaedia Britannica*. 5, 9th edn., Scientific Papers. II, Cambridge: CUP (1875) 541 pp.

3. J.W.S. Lord Rayleigh, Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density. In: *Scientific Papers*. II, Cambridge: CUP (1900) pp. 200–227.

4. G.I. Taylor, The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. *I. Proc. R. Soc. London A201* (1950), 192–196.

5. J.W. Jacobs and I. Catton, Three-dimensional Rayleigh-Taylor instability. Part 2. Experiment. *J. Fluid Mech.* 187 (1988) 353–371.

6. F. Duprez, Mémoire sur un cas particulier de l'équilibre des liquides. I. *Mém. Acad. R. Sci. Lett. Beaux-Arts Belgique* 26 (1851) 42 pp.

7. F. Duprez, Mémoire sur un cas particulier de l'équilibre des liquides. II. *Mém. Acad. R. Sci. Lett. Beaux-Arts Belgique* 28 (1854) 34 pp.

8. D.J. Lewis, The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. *II. Proc. R. Soc. London A202* (1950) 81–96.

9. D.H. Michael, Meniscus stability. *Ann. Rev. Fluid Mech.* 13 (1981) 189–215.

10. H. Tsuge, Hydrodynamic of bubble formation from submerged orifices. In: N.P. Cheremisinoff (ed.), *Encyclopaedia of Fluid Mechanics: Dynamics of single fluid flows and mixing*. Houston, TX: Gulf Professional Publishing (1986) 1506 pp.

11. S.S. Sadhal, P.S. Ayyaswamy and J.N. Chung, *Transport Phenomena with Drops and Bubbles*. New-York: Springer (1997) 520 pp.

12. G.H. Pimbley, Stationary solutions of the problem of Rayleigh-Taylor instability. *J. Math. Anal. Applics* 55 (1976) 170–206.
13. V.E.B. Dussan and S.H. Davis, Stability in systems with moving contact lines. *J. Fluid Mech.* 173 (1986) 115–130.

14. A. Elgowainy and N. Ashgriz, The Rayleigh-Taylor instability of viscous fluid layers. *Phys. Fluids* 9 (1997) 1635–1649.

15. R. Bellman and R.H. Pennington, Effects of surface tension and viscosity on Taylor instability. *Q. Appl. Maths* 12 (1954) 151–162.

16. S. Chandrasekhar, The character of the equilibrium of an incompressible heavy viscous fluid of variable density. *Proc. Cambridge Philos. Soc.* 51 (1955) 162–178.

17. J.W. Jacobs and I. Catton, Three-dimensional Rayleigh-Taylor instability. Part 1. Weakly nonlinear theory. *J. Fluid Mech.* 187 (1988) 329–352.

18. B.J. Daly, Numerical study of two-fluid Rayleigh-Taylor instability. *Phys. Fluids* 10 (1967) 297–307.

19. M.J. Tan, On the steady solutions of the problem of Rayleigh-Taylor instability. *J. Fluid Mech.* 170 (1986) 339–353.

20. A.K. Chesters, An analytical solution for the profile and volume of a small drop or bubble symmetrical about a vertical axis. *J. Fluid Mech.* 81 (1977) 609–624.

21. S.W. Rienstra, The shape of a sessile drop for small and large surface tension. *J. Eng. Math.* 24 (1990) 193–202.

22. A. Marmur and E. Rubin, Equilibrium shapes and quasi-static formation of bubbles at submerged orifice. *Chem. Engng Sci.* 28 (1973) 1455–1464.

23. V.E.B. Dussan, On the spreading of liquids on solid surfaces: static and dynamic contact lines. *Ann. Rev. Fluid Mech.* II 11 (1979) 371–400.

24. I. Stakgold, *Boundary Value Problems of Mathematical Physics*. I, II. New-York: Macmillan Series in Adv. Math. and Theor. Phys. (1968) 340 pp + 408 pp.

25. I. Stakgold, Branching of solutions of nonlinear equations. *SIAM Rev.* B (1971) 289–332.

26. G. Iooss and D.D. Joseph, *Elementary Stability and Bifurcation Theory*, 2nd edn. New York: Springer (1990) 324 pp.

27. R.R. Hugues, A.E. Handlos, H.D. Evans and R.L. Maycock, The formation of bubbles at simple orifices. *Chem. Engng. Prog.* 51 (1955) 557–563.

28. J.F. Padday, Transient surface properties in liquid bridge and pendant drop menisci in gravity and low gravity. In: A. Steichen (ed.), *Dynamics of Multiphase Flows across Interfaces, Lecture notes in Physics*. Berlin: Springer (1996) pp. 41–55.

29. J.F. Davidson and B.O.G. Schüler, Bubble formation at an orifice in a viscous liquid. *Trans. Inst. Chem. Engrs.* 38 (1960) 144–154.

30. S.-C. Georgescu, *Évolution d'une Bulle: Formation à Partir d'un Orifice et Éclatement à la Traversée d'une Surface Libre*. Ph.D. thesis, Institut National Polytechnique de Grenoble, France (1999) 174 pp.

31. S.-C. Georgescu and J.-L. Achard, Multiple equilibrium states during the quasi-steady formation of bubbles and drops at a circular orifice (to be submitted to J. Eng. Math.).