EVOLUTION, STABILITY AND EQUILIBRIUM SHAPES
OF ROTATING DROPS WHICH ARE CHARGED OR
SUBJECT TO ELECTRIC FIELDS
by V. J. GARCÍA-GARRIDO, M. A. FONTELOS
†
(Instituto de Ciencias Matemáticas (CSIC - UAM - UC3M - UCM). C/Nicolás Cabrera,
13-15, Campus de Cantoblanco, 28049 Madrid, Spain)
and
U. KINDELÁN
(Departamento de Matemática Aplicada y Met. Inf., Universidad Politécnica de Madrid,
Alenza 4, 28003 Madrid, Spain)
[Received 2 August 2012. Revise 27 May 2013. Accepted 20 July 2013]
Summary
In this article, we study by means of the boundary element method the effect that rotation at
constant angular momentum L has on the evolution of a conducting and viscous drop when
it holds an amount of charge Q on its surface or is immersed in an external electric field of
magnitude E
∞
acting in the direction of the rotation axis. This droplet is considered to be
contained in another viscous and insulating fluid. Our numerical simulations and stability analysis
show that the Rayleigh fissibility ratio χ at which charged drops become unstable decreases with
angular momentum. For neutral drops subject to an electric field, the critical value of the field
which destabilizes the drop increases with rotation. Concerning equilibrium shapes, approximate
spheroids and ellipsoids are obtained and the transition values between these two families of
solutions is described. When the drop becomes unstable, a two-lobed structure forms where a
pinch-off occurs in finite time or dynamic Taylor cones (in the sense of [Betelú et al., Phys.
Fluids. 18 (2006)]) develop, whose semiangle, for small L, remains the same as if there was no
rotation in the system.
1. Introduction
Electrohydrodynamics is an area of fluid mechanics which has become fundamental to many
industrial and technological applications in recent years. The development of techniques such as
electrospinning, electrospraying and the design of Field Emission Electric Propulsion (FEEP) colloid
thrusters for space vehicles and satellites (1) is clear evidence of its importance. We should also
mention its contributions to electrowetting and electronic paper (2) and the promising microfluidic
chips, where electric fields are used to control chemicals inside a small device with very thin channels
(3, 4). Concerning physical applications, the understanding of coalescence and fission processes for
charged droplets is crucial to study how thunderstorm clouds are formed (5).
Since the pioneering works in 1843 by J. Plateau (6), where an oil drop immersed in a neutral
buoyancy tank was rotated by turning a shaft, many scientists have become interested in studying
†
<marco.fontelos@icmat.es>
Q. Jl Mech. Appl. Math, Vol. 66. No. 4 © The Author, 2013. Published by Oxford University Press;
all rights reserved. For Permissions, please email: journals.permissions@oup.com
Advance Access publication 5 November 2013. doi:10.1093/qjmam/hbt015
Downloaded from https://academic.oup.com/qjmam/article/66/4/489/1943361 by guest on 21 August 2022
490 V. J. GARCÍA-GARRIDO et al.
the behaviour of droplets subject to different forces such as rotation, gravity, surface tension and
electromagnetic fields. The rotating drop problem for example was used by Chandrasekhar in the
field of astrophysics to determine the shape of self-gravitating masses (7). Adetailed numerical study
of equilibrium configurations and their stability for rotating drops was first undertaken by Brown
and Scriven in (8) and improved by Heine (9). From the point of view of the evolution problem,
a detailed analysis with boundary element method (BEM) that considers also rotation at constant
agular momentum has recently been presented in (10).
The effects of charge on droplets was first investigated by Lord Rayleigh in 1882 (11). With an
energy stability analysis he showed that, for a spherical, isolated and conducting drop of radius R
and surface tension γ surrounded by an insulating medium of permittivity ε
0
, a loss of stability
occurs at critical values of the charge Q
n
c
= 4π
ε
0
γ R
3
(
n +2
)
to shape perturbations given by the
n-th order Legendre polynomial (number of lobes on the perturbed shape). At each critical value
Q
n
c
, the sphere is neutrally stable and a family of n-lobed shapes branch. These families of n-lobed
shapes were described numerically in 1989 by Basaran and Scriven using finite elements (12).
Regarding conducting drops under the influence of electric fields, Taylor obtained in 1964 a family
of approximate prolate solutions and identified a critical value of the electric field for which these
configurations become unstable and develop Taylor cones (13, 14). Numerical evolution with BEM
of charged or neutral viscous and conducting drops subject to uniform electric fields is addressed
in (15).
The aim of this article is to describe the effects that rotation at constant angular momentum has
on the evolution of charged or neutral viscous and conducting drops subject to an electric field. This
is done by means of BEM, which is a suitable choice for tracking interfaces in Stokes regimes. An
interesting aspect is that our algorithm is implemented to be adaptive, since a good resolution is
needed in regions of the mesh where singularities develop. The rotating and charged drop model
became relevant when Bohr and Wheeler proposed it in 1939as a simplified version of the mechanism
of nuclear fission (16). Our main motivation for this work comes from the existent discrepancy in
the values between theory (14) and laboratory experiments (5) for the measured semiangle of Taylor
cones and Rayleigh’s critical charge. We will establish how the critical charge or electric field a drop
can sustain before becoming unstable varies with angular mometum. We also present several methods
to obtain approximate equilibrium solutions and compare them with those introduced by Rosenkilde
and Randall (17, 18), who used an appropriate extension of Chandrasekhar’s virial method when
rotation takes place at constant angular velocity.
We start by introducing in Section 2 the mathematical equations that govern the physical system
and dedicate Section 3 to briefly describe the numerical method implemented to solve the evolution
problem. We follow with Section 4, where a theoretical and numerical approach to compute
axisymmetric solutions of charged rotating droplets is presented. Section 5 discusses axisymmetric
rotating drops in electric fields and Section 6 studies the effects of rotation on dynamic Taylor cones.
To finish, in Sections 7 and 8, a 3D stability analysis is conducted for charged rotating drops or
rotating drops immersed in an electric field.
2. Mathematical formulation
We are interested in the evolution, stability and equilibrium configurations of a viscous and
conducting drop surrounded by another viscous and insulating fluid. Both fluids rotate about
a common axis (e.g. the z axis) with constant angular momentum L and are assumed to be
incompressible. We consider the following two situations: in the first, the droplet has an amount
Downloaded from https://academic.oup.com/qjmam/article/66/4/489/1943361 by guest on 21 August 2022
EVOLUTION, STABILITY AND EQUILIBRIUM SHAPES OF ROTATING DROPS 491
of charge Q distributed over its surface, and in the second, the drop is neutral and subject to a
uniform external electric field of magnitude E
∞
along its axis of rotation. Working in a rotating
frame of reference, Navier-Stokes equations for a fluid with viscosity μ
i
, density
i
, pressure p
(i)
and velocity field u
(i)
have the form (19):
⎧
⎪
⎨
⎪
⎩
i
∂u
∂t
(i)
+ u
(i)
·∇u
(i)
=−∇p
(i)
+ μ
i
u
(i)
− 2
i
ω × u
(i)
−
i
ω ×
(
ω × r
)
, in D
i
(
t
)
∇·u
(i)
= 0, in D
i
(
t
)
,
(2.1)
where ω ≡ (0, 0,ω) is the angular velocity and D
1
(
t
)
is the region enclosed by the droplet and
D
2
(
t
)
that of the surrounding fluid. Observe that, since we will work in regimes where the Reynolds
number is small, there will be time for diffusion of vorticity to occur provided that deformation is
sufficiently slow, thus leading to an almost solid-body rotation. These allows to neglect the ficticious
force arising from the variable rate of rotation of the frame of reference. The terms ω × u
(i)
and
ω ×
(
ω × r
)
represent the Coriolis and centrifugal forces respectively. We can write the centrifugal
force as:
−ω ×
(
ω × r
)
= ω
2
r
axis
e
r
=∇
1
2
ω
2
r
2
axis
, (2.2)
where r
axis
is the orthogonal distance from a point at the surface of the drop to the axis of rotation.
Suppose now that the density of the surrounding fluid,
2
, is very small with respect to the density
of the droplet,
1
, and set
2
= 0. If we define reduced pressures as:
(1)
= p
(1)
−
1
L
2
2I
2
r
2
axis
,
(2)
= p
(2)
, (2.3)
and the drop’s moment of inertia:
I =
1
D
1
(
t
)
r
2
axis
dV, (2.4)
then (2.1) becomes:
⎧
⎪
⎨
⎪
⎩
i
∂u
∂t
(i)
+ u
(i)
·∇u
(i)
=∇−
(
i
)
+ μ
i
u
(
i
)
− 2
i
ω × u
(i)
, in D
i
(
t
)
.
∇·u
(
i
)
= 0, in D
i
(
t
)
(2.5)
In order to non-dimensionalize (2.5) we introduce a characteristic length l, a typical velocity U and
a characteristic time scale τ = l/U so that:
u
(i)
=
u
(i)
U
,
r =
r
l
,
t =
t
τ
,
(i)
=
l
μ
1
U
(i)
, ω = ω
z,λ=
μ
2
μ
1
,ζ=
ρ
2
ρ
1
.
Omitting overbars to simplify notation, the non-dimensional problem is:
⎧
⎪
⎨
⎪
⎩
ζ
i−1
Re
∂u
∂t
(i)
+ u
(i)
·∇u
(i)
=−∇
(i)
+ λ
i−1
u
(i)
− 2
ζ
i−1
Ek
z × u
(i)
, in D
i
(
t
)
.
∇·u
(i)
= 0, in D
i
(
t
)
(2.6)
Downloaded from https://academic.oup.com/qjmam/article/66/4/489/1943361 by guest on 21 August 2022
492 V. J. GARCÍA-GARRIDO et al.
Two dimensionless parameters arise: Re is the Reynolds number (measures the relative importance
between inertial and viscous forces) and Ek is the Ekman number (characterizing the relation between
Coriolis and viscous forces). They are defined as:
Re =
1
Ul
μ
1
, Ek =
μ
1
1
ωl
2
, (2.7)
In this article we will study the case where viscous forces dominate inertial and Coriolis forces. This
regime is known as Stokes flow and is characterized by Re ≪ 1 and Ek ≫ 1. Thus we have Stokes
system:
−∇
(i)
+ λ
i−1
u
(i)
= 0, in D
i
(
t
)
∇·u
(i)
= 0, in D
i
(
t
)
, (2.8)
to which we impose that the normal component of the velocity is continuous across the boundary:
u
(1)
· n = u
(2)
· n ≡ u · n, (2.9)
and the kinematic condition:
v
n
= u · n, on ∂D
(
t
)
, (2.10)
with v
n
being the normal velocity of the free boundary.
Now, since the drop is considered an ideal conductor, the potential V must be constant inside and
at the drop’s surface, and all the charge is located at the boundary. The electric potential satisfies the
Laplace equation:
⎧
⎪
⎨
⎪
⎩
V = 0, in D
2
(
t
)
V = V
0
, in ∂D
(
t
)
V →−E
∞
z + O(
|
r
|
−1
), as
|
r
|
→∞
, (2.11)
and V
0
has to be chosen so that the total charge is Q. At the boundary of the drop the surface charge
density σ is given by the normal derivative of the potential, σ =−ε
0
∂V
∂n
, with ε
0
the permittivity of
the surrounding fluid and n is the outward unit normal to the surface of the droplet. At the surface
of a conductor, the repulsive electrostatic force per unit area is:
F
e
=
ε
0
2
∂V
∂n
2
n =
σ
2
2ε
0
n.
To solve (2.8) we set a balance between viscous stresses and capillary, electrostatic and centrifugal
forces at the interface of both fluids:
T
(2)
− T
(1)
n =
2γ H −
1
L
2
2I
2
r
2
axis
−
σ
2
2ε
0
n, on ∂D
1
(
t
)
, (2.12)
where γ is the surface tension and H the mean curvature and, for Newtonian fluids:
T
(k)
ij
=−
(k)
δ
ij
+ μ
k
⎛
⎝
∂u
(k)
i
∂x
j
+
∂u
(k)
j
∂x
i
⎞
⎠
, k = 1, 2. (2.13)
Downloaded from https://academic.oup.com/qjmam/article/66/4/489/1943361 by guest on 21 August 2022
EVOLUTION, STABILITY AND EQUILIBRIUM SHAPES OF ROTATING DROPS 493
To make the boundary condition (2.12) dimensionless we introduce the characteristic scales:
l =
3
√
V, U =
γ
μ
1
,τ=
l
U
, (2.14)
and the quantities:
T
(
i
)
=
l
γ
T
(
i
)
,
r
axis
=
r
axis
l
,
H = lH, I =
I
l
5
1
,
σ =
l
γε
0
σ, (2.15)
where V is the volume occupied by the drop. By defining:
L =
L
1
l
7
γ
,χ=
Q
2
48π
=
Q
2
48πγǫ
0
l
3
, E
∞
=
ε
0
l
γ
E
∞
, (2.16)
we get for the boundary condition, ommiting overbars:
T
(2)
− T
(1)
n =
2H −
L
2
2I
2
r
2
axis
−
σ
2
2
n, on ∂D
(
t
)
. (2.17)
In what follows we will describe the results in terms of the parameters in (2.16) where χ is known
as Rayleigh’s fissibility ratio. Notice that one can easily recover the physical values by solving
for L, Q and E
∞
Concerning the stationary problem (where u
(i)
≡ 0), equilibrium solutions can be calculated by
solving the modified Young–Laplace equation:
δ = 2γ H −
1
L
2
2I
2
r
2
axis
−
σ
2
2ε
0
, on ∂D, (2.18)
where D is the region enclosed by the drop and δ =
(1)
−
(2)
is the reduced pressure difference
across the drop’s surface. It is important to note that this equilibrium equation can be seen as
the Euler–Lagrange equation of a particular energy functional. In this variational framework, our
problem reduces to the minimization of the total energy for the closed system, which we can write as:
E
total
= E
area
+ E
kinetic
+ E
electrostatic
− δ
(
V − V
0
)
, (2.19)
where δ
(
V − V
0
)
is the constraint (δ plays the role of a Lagrange multiplier) so the volume of
the drop is V = V
0
, E
area
is the energy due to surface area, E
kinetic
is the rotational kinetic energy
for constant angular momentum rotation and E
electrostatic
is the electrical energy. These energies can
be obtained from:
E
area
= γ Area
(
∂D
)
, E
kinetic
=
L
2
2I
, E
electrostatic
=
ε
0
2
R
3
\D
|
E
|
2
dV. (2.20)
Since the electric field verifies E =−V, after integrating by parts one can write the electrostatic
energy as:
E
electrostatic
=
ε
0
2
R
3
\D
|
E
|
2
dV =−
ε
0
2
∂D
V
∂V
∂n
dS =
1
2
∂D
V
0
σ dV =
1
2
QV
0
=
Q
2
2C
, (2.21)
Downloaded from https://academic.oup.com/qjmam/article/66/4/489/1943361 by guest on 21 August 2022