# Numerical simulation of uncharged droplets in a uniform electrical field

Adel Mustapha Benselama, Jean-Luc Achard, Pascale Pham

► To cite this version:

Adel Mustapha Benselama, Jean-Luc Achard, Pascale Pham. Numerical simulation of uncharged droplets in a uniform electrical field. Journal of Electrostatics, 2006, 64 (7-9), pp.562-658. 10.1016/j.elstat.2005.10.011. hal-00182345

HAL Id: hal-00182345

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

Submitted on 24 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
Numerical simulation of an uncharged droplet in a uniform electric field

Adel M. Benselama<sup>a,b,*</sup>, Jean-Luc Achard<sup>c</sup>, Pascale Pham<sup>b</sup>

<sup>a</sup>LEMD, CNRS–UJF, 25 Rue des Martyrs, BP166, 38042 Grenoble Cedex 9, France

<sup>b</sup>LETI (CEA–DRT)–DTBS, 17 Avenue des Martyrs, 38054 Grenoble Cedex 9, France

<sup>c</sup>LEGI, CNRS, 1025 rue de la Piscine, Domaine Universitaire, BP 53, 38041 Grenoble Cedex, France

The deformations of a real dielectric droplet subjected to a uniform alternating electric field and immersed in an insulating fluid are numerically studied by the Boundary Element Method. The alternating electric field time scale is taken much smaller than the hydrodynamic time scale of the droplet shape deformation. The influence of the frequency and the conductivity of the droplet upon the critical electric field, beyond which instabilities develop, are compared with the experimental measurements. Numerical results well account for experimental data while an unexpected good correspondence with Taylor’s theory is found.

Keywords: Dielectric droplet deformation; Electric field effects; Boundary Element Method

# 1. Introduction

## 1.1. Objective

Applying an electric field onto water droplets immersed in oil is an efficient process to make them coalesce and thus extract water from crude oil for instance. In that case, water droplets can be supposed to behave like electric dipoles and attract each other to form larger drops. However, when the electric field reaches a critical value, it just may cause the opposite effect; instead of coalescing, droplets may break into smaller ones and the emulsion may become more difficult to separate. To prevent such a behavior, many studies on droplet’s deformations were performed, specially by experimental means [1]. However, only a few numerical simulations studies involving general models, have been made to follow droplet’s deformations. In fact, almost all proposed analysis are *local*: linear stability is used or droplets shapes are restricted to spheroids that are close to spherical forms.

Hence, we have developed a numerical model that simulates the behavior of a droplet subjected to a uniform electric field. This model is based on the Boundary Element Method (BEM) formulation of the hydrodynamic’s equations of irrotational fluids in a steady harmonic electrokinetic regime. In order to validate our model, the simulation results are compared to Berg et al. experimental data [1], where more or less conductive droplets are subjected to an alternating current regime. The 2D-axisymmetric test configuration consists of a conductive droplet (water) immersed in a more insulating liquid (mineral oil) and subjected to a uniform sinusoidal electric field generated between two plane electrodes (see Fig. 1).

In response to the applied electric field, Berg et al. found two specific behaviors: (i) at low electric field magnitude, spherical droplets distort either to prolate or oblate spheroids, (ii) for electric field magnitude larger than a threshold (that we call the critical electric field $E_c$), instabilities develop (sub-droplet ejection, thin jet, cloud of fine droplets, etc.) at the poles following the appearance of Taylor’s cones [2]. The adjustable parameters in Berg’s experiments are the magnitude, the frequency and the waveform of the applied electric field, in one hand, and the droplet electric conductivity, in the other.

*Corresponding author. LEMD, CNRS–UJF, 25 Rue des Martyrs, BP166, 38042 Grenoble Cedex 9, France. Tel.: +33 4 38 78 24 12; fax: +33 4 38 78 57 87.

E-mail address: adel.benselama@cea.fr (A.M. Benselama).
<figure><img src="image_2.png" /><figcaption>Fig. 1. Configuration for numerical computations.</figcaption></figure>

In another viewpoint, a comparison is made to distinguish the effect, on droplet breakup, of a DC electric field (via Taylor's theory [2]) from an AC electric field.

The numerical model we developed is expected to provide a less patchy sight of uncharged droplet instability phenomena, as spotlighted all this paper along.

## 1.2. *Studied configuration*

The studied configuration (see Fig. 1) is likely to computationally reproduce Berg's experiments [1]. A water droplet $\Omega^{(1)}$, of radius $r$, is placed between two plane electrodes ($\Gamma_1$ and $\Gamma_2$ boundaries). The droplet is immersed in oil $\Omega^{(2)}$. The boundary $\Gamma_3$ is an insulating material. The electric field applied by the electrodes, is assumed to be uniform far from the droplet. Hence, the electrode-to-electrode distance $L$ is chosen at least one order larger than the droplet radius $r$. The 2D-axisymmetric configuration details are shown in Fig. 1 below.

An alternating electric potential is imposed to the left electrode

$$ \tilde{V}(t) = (1 - e^{-t/\tau})\mathrm{Re}(V_{\max}e^{i\omega t}), \quad (1) $$

where $\omega$ is the angular frequency ($\omega = 2\pi f$ and $f$ the frequency), $\tau$ is the specific electric field rise time. In any case, $\tau$ must be larger than the electric relaxation time $\tau_e$ of any interface. For a tap water droplet immersed in oil, $\tau_e = 0.1$ µs. The right electrode is grounded.

The time scale of hydrodynamic phenomena, and in particular of the interface motion, will be used to describe the whole process. The weak coupling assumption, which permits to treat the electric field as noninfluenced by the droplet time evolution, is valid when this time scale $t$ is much larger than the applied electric field period. As interface movement is concerned, two main frequencies are involved; the capillary and the viscous frequencies. For a water droplet in oil, both are presented in Fig. 2.

Because of high oil viscosity, we do consider that the droplet dynamics are fundamentally governed by longitudinal viscous stresses, so we take the viscous frequency as a reference for the electric frequency, which has to satisfy the condition

$$ f \gg f_v \quad (2) $$

in order that the weak coupling hypothesis is insured.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Capillary, $f_v$, and viscous, $f_v$, typical frequency scales with respect to the droplet radius ($v$ is the oil kinematic viscosity).</figcaption></figure>

# 2. Mathematical model

The equations, which govern the motion of a multiphase flow subjected to electric forces and the numerical model we used to solve this problem, are presented below. Each phase is assumed to consist in a Newtonian isotropic fluid with uniform isotropic mechanical and electrical properties. The fluids are immiscible.

## 2.1. *The electrokinetic model*

### 2.1.1. *General equations*

Under the assumption of weak coupling between electric and hydrodynamic effects, the linearity of the electric system implies that its response is also sinusoidal, i.e.

$$ \mathbf{E}^{(k)}(\mathbf{x}, t) = \mathrm{Re} \left[ \bar{\mathbf{E}}^{(k)}(\mathbf{x}) e^{i\omega t} \right]. \quad (3) $$

Also, the quasi-static approximation, where the coupling between the electric field and the magnetic one is neglected, allows us to assume that the electric field derives from an electric potential $\bar{\mathbf{V}}^{(k)}$:

$$ \bar{\mathbf{E}}^{(k)} = -\nabla \bar{\mathbf{V}}^{(k)} \quad (4) $$

which implies [3]

$$ \nabla \cdot (\bar{\varepsilon}^{(k)} \nabla \bar{\mathbf{V}}^{(k)}) = 0, \quad k = 1, 2, \quad (5) $$

where $\bar{\varepsilon}^{(k)}$ is the complex permittivity of the fluid $k$, defined by

$$ \bar{\varepsilon}^{(k)} = \varepsilon^{(k)} + \frac{\sigma^{(k)}}{i\omega} \quad (6) $$

$\varepsilon^{(k)}$ and $\sigma^{(k)}$ are, respectively, its permittivity and electric conductivity. Imposed potential at the left electrode $\Gamma_1$ implies that

$$ \bar{\mathbf{V}}_{1}^{(2)} = (1 - e^{-t/\tau})V_{\max} \quad (7) $$

while at the grounded one, $\Gamma_2$, we have

$$ \bar{\mathbf{V}}_{2}^{(2)} = 0. \quad (8) $$
The insulating condition at the boundary $\Gamma_3$ writes down as

$$ \partial \tilde{V}^{(2)}/\partial n \Big|_{\Gamma_3} = 0. \qquad (9) $$

$\mathbf{n}$ is the outward normal to the boundary $\Gamma_3$.

The cross condition through the interface writes as

$$ (\tilde{\epsilon}^{(1)}\tilde{\mathbf{E}}^{(1)} - \tilde{\epsilon}^{(2)}\tilde{\mathbf{E}}^{(2)}) \cdot \mathbf{n} \Big|_{S} = 0. \qquad (10) $$

### 2.1.2. Boundary element applied to electrokinetics

Since the electric properties are assumed uniform, for each domain $\Omega^{(k)}$, Eq. (5) writes

$$ \nabla^2 \tilde{V}^{(k)} = 0. \qquad (11) $$

By using the Green second identity, Eq. (11) implies the following integral equation:

$$ \alpha^{(k)}(\mathbf{x}) \tilde{V}^{(k)}(\mathbf{x}) = \int_{\partial \Omega^{(k)}} \left( G \partial \tilde{V}^{(k)} / \partial n'^{(k)} - \tilde{V}^{(k)} \partial G / n'^{(k)} \right) dS', \qquad (12) $$

where $\partial\Omega^{(k)}$ is the boundary of $\Omega^{(k)}$ (when considering Fig. 1, $\partial\Omega^{(2)} = \Gamma_1 + \Gamma_2 + \Gamma_3 + S$ and $\partial\Omega^{(1)} = S$). $\mathbf{n}^{(k)}$ is the outward normal to $\partial\Omega^{(k)}$ and $\alpha^{(k)}$ is given by

$$ \alpha^{(k)}(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \in \Omega^{(k)}, \\ \frac{1}{2} & \text{if } \mathbf{x} \in \partial\Omega^{(k)}, \\ 0 & \text{if } \mathbf{x} \notin \Omega^{(k)}. \end{cases} \qquad (13) $$

$G$ is the Green's kernel defined by

$$ G(\mathbf{x}, \mathbf{x}') = 1/4\pi \|\mathbf{x} - \mathbf{x}'\|. \qquad (14) $$

Differentiating Eq. (12), with respect to the normal to the boundary $\partial\Omega^{(k)}$, gives

$$ \alpha^{(k)} \partial \tilde{V}^{(k)} / \partial n \Big|_{\partial \Omega^{(k)}} = - \mathbf{n} \cdot \int_{\partial \Omega^{(k)}} \nabla' G \partial \tilde{V}^{(k)} / \partial n'^{(k)} dS' \\ + \mathbf{n} \cdot \int_{\partial \Omega^{(k)}} \tilde{V}^{(k)} \nabla' \partial G / \partial n'^{(k)} dS'. \qquad (15) $$

Let us denote the jump of the normal component, to the interface, of the electric field by $\tilde{\epsilon} = \mathbf{n} \cdot (\tilde{\mathbf{E}}^{(1)} - \tilde{\mathbf{E}}^{(2)})|_S$, and $\tilde{\delta}$ the normal component $E_n^{(2)}$, to the boundaries $\Gamma_1$ and $\Gamma_2$, of the electric field. By using condition (10), the sum, side by side, of the product of Eq. (15) by $\tilde{\epsilon}^{(k)}$, for $k=1,2$, writes, when considering boundary conditions, as

$$ \frac{1}{2} \tilde{\epsilon} = (\tilde{\epsilon}^{(1)} - \tilde{\epsilon}^{(2)}) / (\tilde{\epsilon}^{(1)} + \tilde{\epsilon}^{(2)}) \times \left[ \int_S (\mathbf{n} \cdot \nabla' G) \tilde{\epsilon} \, dS' - \int_{\Gamma_1} (\mathbf{n} \cdot \nabla' G) \tilde{\delta} + V_1^{(2)} (\mathbf{n} \cdot \nabla' \partial G / \partial n') \, dS' - \int_{\Gamma_2} \tilde{\delta} (\mathbf{n} \cdot \nabla' G) \, dS' - \int_{\Gamma_3} \tilde{V}^{(2)} (\mathbf{n} \cdot \nabla' \partial G / \partial n') \, dS' \right]. \qquad (16) $$

Once the electric field solved, we can get the mean electric stress contribution, by integrating the Maxwell

stress tensor over one period of the alternating electric field, which gives

$$ \langle T^{(\text{Maxwell})} \rangle = \frac{1}{2} \epsilon \text{Re} [\tilde{\mathbf{E}} \cdot \tilde{\mathbf{E}}^* - \frac{1}{2} (\tilde{\mathbf{E}} \cdot \tilde{\mathbf{E}}^*)^2] I. \qquad (17) $$

## 2.2. Mechanical model

The Newtonian fluid mechanics equations are described hereafter. Fluid flows are supposed to be irrotational in the whole volume they occupy. Yet, it remains possible to take into account a partial effect of the viscosity; at the interface $S$, the balance of the normal momentum may consider the dissipative phenomenon due to the longitudinal deformation of the underlying fluids. To this effect, is added the surface tension effect as well as the electric one.

### 2.2.1. Basic equations

The equation of conservation of mass and the Navier-Stokes equations, for a Newtonian incompressible fluid $k$ ($k=1,2$), write as

$$ \tilde{\nabla} \cdot \tilde{\mathbf{u}}^{(k)} = 0, \qquad (18) $$

$$ \partial \tilde{\mathbf{u}}^{(k)} / \partial \tilde{t} + \frac{1}{2} \tilde{\nabla} (\tilde{\mathbf{u}}^{(k)})^2 = \tilde{\nabla} \tilde{\psi}^{(k)} - 1 / \rho^{(k)} \times \tilde{\nabla} \tilde{\rho}^{(k)} + \mu^{(k)} / \rho^{(k)} \times \tilde{\nabla}^2 \tilde{\mathbf{u}}, \qquad (19) $$

where $\tilde{\mathbf{u}}^{(k)}, \tilde{\rho}^{(k)}, \tilde{\psi}^{(k)}, \rho^{(k)}$ and $\mu^{(k)}$ stand for the velocity, the pressure, the volume force potential, the density and the dynamic viscosity of the fluid $k$, respectively.

At the electrodes, $\Gamma_1$ and $\Gamma_2$, and the insulating boundary, $\Gamma_3$, a no slip condition is assigned. That is

$$ \mathbf{u}|_{\Gamma_1 + \Gamma_2 + \Gamma_3} = 0. \qquad (20) $$

The impermeability condition at the interface $S$ implies the continuity of the normal velocity across it

$$ u_n^{(1)}|_S = u_n^{(2)}|_S. \qquad (21) $$

If we suppose the flow irrotational, then the velocity field derives from a hydrodynamic potential $\tilde{\phi}^{(k)}$

$$ \tilde{\mathbf{u}}^{(k)} = \tilde{\nabla} \tilde{\phi}^{(k)}. \qquad (22) $$

Thus, the traditional Bernoulli equation is obtained

$$ \partial \tilde{\phi}^{(k)} / \partial \tilde{t} = - \frac{1}{2} (\tilde{\mathbf{u}}^{(k)})^2 + \tilde{\psi}^{(k)} - \tilde{\rho}^{(k)} / \rho^{(k)}. \qquad (23) $$

In the same manner, introducing (22) into (18) gives

$$ \tilde{\nabla}^2 \tilde{\phi}^{(k)} = 0. \qquad (24) $$

The irrotationality hypothesis does not allow to consider the tangential momentum balance especially on the boundaries $\Gamma_1, \Gamma_2$ and $\Gamma_3$. Rather than (20), the condition to be satisfied at these boundaries will be

$$ u_n|_{\Gamma_1 + \Gamma_2 + \Gamma_3} = (\mathbf{u} \cdot \mathbf{n})|_{\Gamma_1 + \Gamma_2 + \Gamma_3} = 0. \qquad (25) $$

Let us denote $U$ the velocity scale, that will be determined later. The resulting nondimensional quantities
are $\mathbf{u}^{(k)} = \tilde{\mathbf{u}}^{(k)}/U$, $p^{(k)} = \tilde{p}^{(k)}/\rho^{(k)}U^2$, $\psi^{(k)} = \tilde{\psi}^{(k)}/U^2$, $t = U\tilde{t}/L$ and $\phi^{(k)} = \tilde{\phi}^{(k)}/UL$.

Hence, the nondimensional form of Eqs. (24) and (23) will be

$$\nabla^2 \phi^{(k)} = 0, \qquad (26)$$

$$\partial_t \phi^{(k)} / \partial t = -\frac{1}{2} (\mathbf{u}^{(k)})^2 + \psi^{(k)} - p^{(k)} / \rho^{(k)}. \qquad (27)$$

Let us also suppose that both fluids have the same density $\rho$. By subtracting, side by side, Eq. (23), we get

$$\partial \Phi / \partial t = -\frac{1}{2} (\mathbf{u}^{(1)})^2 + \frac{1}{2} (\mathbf{u}^{(2)})^2 - (p^{(1)} - p^{(2)}) / \rho, \qquad (28)$$

where $\Phi = (\phi^{(1)} - \phi^{(2)})|_S$ is called the generalized hydrodynamic potential.

The normal momentum quantity balance at each point in the interface $S$ gives [4]

$$[\mathbf{n} \cdot \mathbf{T}^{\text{dyn}}] \cdot \mathbf{n}^{(1)} + \gamma \nabla_{(s)} \cdot \mathbf{n}^{(1)} = -[\mathbf{n} \cdot \langle \mathbf{T}^{\text{Maxwell}} \rangle] \cdot \mathbf{n}^{(1)}, \qquad (29)$$

where the symbol $[\mathbf{n} \cdot \Phi] = \sum_{k=1,2} \mathbf{n}^{(k)} \cdot \Phi^{(k)}$ is used to represent the effect of the underlying media to the interface $S$, $\mathbf{n}$ is the normal vector to the interface and $\gamma$ its surface tension. $2\tilde{C}_M = \nabla_{(s)} \cdot \mathbf{n}^{(1)}$, where $\tilde{C}_M$ is the mean curvature of the interface $S$. Tensors $T^{\text{dyn}}$ represents the hydrodynamic tensor. For a Newtonian fluid in quasi-static regime, this tensor is given by

$$\mathbf{T}^{(\text{dyn}, k)} = -\tilde{p}^{(k)} \mathbf{I} + \mu^{(k)} (\tilde{\nabla} \tilde{\mathbf{u}}^{(k)} + (\tilde{\nabla} \tilde{\mathbf{u}}^{(k)})^{\dagger}). \qquad (30)$$

The superscript $\dagger$ stands for the transpose matrix symbol and $\mathbf{I}$ is the identity tensor. The dissipation contribution from fluid volume elongation is retained at the interface (although the fluids are assumed irrotational) in order to bring in a dumping effect to the fluid underlying layers to the interface. Moore [5], who considers the motion of a spherical air bubble in a fluid, shows that the irrotational solution provides, as the Reynolds number $Re = \rho UL/\mu$ increases, a uniform approximation for normal component of $\mathbf{n} \cdot \mathbf{T}^{\text{dyn}}$.

If we use Eq. (29) in order to exclude the pressure jump in Eq. (28), we will obtain the equation

$$\partial \Phi / \partial t = \frac{1}{2} (\mathbf{u}^{(2)})^2 - \frac{1}{2} (\mathbf{u}^{(1)})^2 \\ - 2C_M / We - \left[ Ne^{(2)} \left( (E_n^{(2)})^2 - (E_t^{(2)})^2 \right) \right. \\ - Ne^{(1)} \left( (E_n^{(1)})^2 - (E_t^{(1)})^2 \right) \\ \left. + 2 / \text{Re}^{(2)} \times \partial u_n^{(2)} / \partial n - 2 / \text{Re}^{(1)} \times \partial u_n^{(1)} / \partial n \right], \qquad (31)$$

where $We = \rho U^2 L/\gamma$, $Re^{(k)} = \rho UL/\mu^{(k)}$, $Ne^{(k)} = \varepsilon^{(k)} V_{\text{max}}^2 / \rho U^2 L^2$, $k = 1,2$, and $C_M = \tilde{C}_M \times L$. We remind that, consequently from (17), we have

$$(\mathbf{E}^{(k)})^2 = \frac{1}{2} \text{Re} [\tilde{\mathbf{E}}^{(k)} (\tilde{\mathbf{E}}^{(k)})^*], \quad k = 1, 2. \qquad (32)$$

The velocity scale $U$ can be determined either by setting $We = 1$ or $Re^{(2)} = 1$. As aforementioned, the latter has been selected.

### 2.2.2. Boundary element applied to mechanics

The Green's second identity applied to Eq. (26) gives the following integral equation:

$$\alpha^{(k)}(\mathbf{x})\phi^{(k)}(\mathbf{x}) = \int_{\partial\Omega^{(k)}} G \frac{\partial\phi^{(k)}}{\partial n^{(k)}} dS' - \int_{\partial\Omega^{(k)}} \frac{\partial G}{\partial n^{(k)}} \phi^{(k)} dS'. \qquad (33)$$

If we take the sum, side by side, of Eq. (33), we obtain, for $\mathbf{x} \in S$,

$$\frac{1}{2} (\phi^{(2)} + \phi^{(1)})(\mathbf{x})|_S = \int_S \frac{\partial G}{\partial n'} \Phi \, dS' \\ - \int_{\Gamma_1 + \Gamma_2 + \Gamma_3} G \frac{\partial \phi^{(2)}}{\partial n'} \, dS' \\ + \int_{\Gamma_1 + \Gamma_2 + \Gamma_3} \frac{\partial G}{\partial n'} \phi^{(2)} \, dS'. \qquad (34)$$

$\mathbf{n}'$ is directed from $\Omega^{(1)}$ to $\Omega^{(2)}$.

The same sum (of Eq. (33)) gives, for $\mathbf{x} \in \Gamma_1 + \Gamma_2 + \Gamma_3$,

$$\alpha(\mathbf{x})\phi^{(2)}(\mathbf{x})|_{\Gamma_1+\Gamma_2+\Gamma_3} - \int_{\Gamma_1+\Gamma_2+\Gamma_3} \frac{\partial G}{\partial n'} \phi^{(2)} \, dS' \\ - \int_S \frac{\partial G}{\partial n'} \Phi \, dS' = - \int_{\Gamma_1+\Gamma_2+\Gamma_3} G \frac{\partial \phi^{(2)}}{\partial n'} \, dS'. \qquad (35)$$

If we assume that the generalized hydrodynamic potential is known, from Eq. (31), then the system, which consists of Eqs. (34) and (35), is closed. After resolution, both hydrodynamics potentials $\phi^{(1)}$ and $\phi^{(2)}$ will be available.

Differentiating Eq. (34), with respect to the normal $\mathbf{n}$ to the interface $S$, permits to get the normal velocity on this interface

$$u_n(\mathbf{x})|_S = \int_S \frac{\partial^2 G}{\partial n \partial n'} \Phi \, dS' - \int_{\Gamma_1 + \Gamma_2 + \Gamma_3} \frac{\partial G}{\partial n} \frac{\partial \phi^{(2)}}{\partial n'} \, dS' \\ + \int_{\Gamma_1 + \Gamma_2 + \Gamma_3} \frac{\partial^2 G}{\partial n \partial n'} \phi^{(2)} \, dS'. \qquad (36)$$

The tangential velocities are deduced by differentiating the hydrodynamic potentials $\phi^{(k)}$ ($k = 1,2$) with respect to curvilinear coordinates

$$\mathbf{u}_t^{(k)} = \nabla_{(s)} \phi^{(k)}, \quad k = 1, 2. \qquad (37)$$

Once the electric field and the velocity are known, the generalized hydrodynamic potential $\Phi$ is updated, by time integrating of Eq. (31), as well as the shape of the interface $S$. Consequently, a new resolution is possible for the next time step.

## 3. Results: comparison with experiments and Taylor's theory

The simulation results are compared with the experimental data of Berg et al. [1]. The influence of the
frequency *f* and the droplet conductivity on the magnitude *E*<sub>c</sub> of the critical electric field, beyond which the droplet becomes unstable, is studied. Transient deformation of the droplet shape is also presented. Among frequencies studied by Berg et al. we have retained two of them: 50 and 2000 Hz. Although Berg et al. have obtained different results according to the electric applied waveform, we study only the case where the waveform is sinusoidal. The water droplet is immersed in a very slightly conductive oil. Berg added 5 and 50 g/l of NaCl to the droplet in order to modify its conductivity: these concentrations, respectively correspond (according to Kohlrausch law) to a conductivity of 10<sup>5</sup> and 10<sup>4</sup> µS/cm. Berg showed that at this order of magnitude, the droplet conductivity has no effect on its behavior. This point strengthens the droplet infinite conductivity hypothesis. To see whether or not this hypothesis still holds true, we also simulate pure water droplet to make comparison with tap water conductivity

effect. The different parameters values used in the simulations are summarized in Table 1. Comparison is extended finally to the traditional Taylor's theory of uncharged droplet stability.

### 3.1. Frequency effect

Fig. 3(a) illustrates the variations of *E*<sub>c</sub> according to the droplet radius at *f* = 2000 Hz. Agreement with experiment is reasonably good.

When *f* = 50 Hz, since hydrodynamic low-frequency perturbations are excluded from our simulations, the latter are expected to only give a sufficient condition of droplet instability. In fact, as shown in Fig. 3(b), the numerical stability limit is always higher than the stability obtained from experiment.

Numerical results show, in both cases, the same tendency as the experiments: as the droplet gets smaller, the required electric field to deform and make the droplet unstable gets larger. This behavior is due to the fact that the capillarity effect increases as the droplet size decreases. As the frequency increases, small droplets become more stable, whereas no significant changes are noticeable for bigger ones. The mechanism that renders droplets more stable as the frequency gets larger seems to be included in the range of 50<f<2000 Hz for the smallest droplets whilst it stands at lower frequencies for the biggest ones.

### 3.2. Effect of droplet conductivity

Fig. 4 displays the variation of *E*<sub>c</sub> for two droplet conductivities (σ = 1 and 100 µS/cm). As mentioned by Berg et al. [1], we can see, on this figure, that the droplet conductivity has no significant influence upon *E*<sub>c</sub>. Therefore, the droplet may be considered as infinitely conductive.

### 3.3. Equilibrium shapes

In practice, various deformation regimes are observed, according to the droplet size, frequency and intensity of the

<table><caption>Table 1
Values of the main parameters used in simulations</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>Droplet relative permittivity ε<sup>(1)</sup></td><td>78.8</td><td>—</td></tr><tr><td>External fluid relative permittivity ε<sup>(2)</sup></td><td>2.55</td><td>—</td></tr><tr><td>Droplet conductivity, σ<sup>(1)</sup> (tap water)</td><td>100</td><td>μS/cm</td></tr><tr><td>Droplet conductivity, σ<sup>(1)</sup> (very pure water)</td><td>1</td><td>μS/cm</td></tr><tr><td>External fluid conductivity σ<sup>(2)</sup> (mineral oil)</td><td>3 × 10<sup>-10</sup></td><td>μS/cm</td></tr><tr><td>Surface tension γ of the interface (water-oil)</td><td>40</td><td>mN/m</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Critical electric field versus droplet radius, (a) *f* = 2000 Hz and (b) *f* = 50 Hz.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Effect of conductivity upon the critical electric field for various droplet radii (*f* = 2000 Hz).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. From left to right, sequence of droplet deformation ($r = 500~\mu\text{m}$ and $f = 2000~\text{Hz}$) calculated after applying the electric field. The electric field is gradually applied.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. From left to right, sequence of droplet deformation ($r = 500~\mu\text{m}$, $f = 2000~\text{Hz}$ and $E_0 = 12~\text{kV/cm}$) calculated after applying the electric field. (a) Very slowly applied electric field, $\tau = 0.1~\text{s}$, and (b) very rapidly applied electric field $\tau = 1~\text{µs}$.</figcaption></figure>

applied electric field. Fig. 5 points out the droplet forms obtained from simulations at 2000 Hz and for 500 µm droplet radius. The obtained forms and those recorded by the experiments [1] match up (even if no time evolution is examined in these experiments). If the electric field is applied progressively from 0 to $E_0$, where $E_0 < E_c$, ($E_{\text{max}} = E_0(1 - e^{-t/\tau})$, see Eq. (1), with $\tau \gg 1/f_v$), then the droplet deforms in a sequence of near-equilibrium prolate spheroidal forms.

### 3.4. Instability regimes

At the same applied frequency and for $E_0 > E_c$ and $\tau \gg 1/f_v$, the droplet starts to stretch and suddenly the geometry changes from a prolate to an oblate spheroid (Fig. 6(a)). This behavior has been observed in Berg et al. experiments [1]. If $\tau$ is reduced, but larger than $\tau_e$, and $E_0$ reaches the critical value $E_c$ then the droplet is deformed until a cone appears at each droplet's pole (Fig. 6(b)).

### 3.5. Taylor's instabilities

The nonlinear stability analysis (involving high deformation rates from equilibrium shapes) being presented in Section 3.1, we focus our attention now on linear stability,

<figure><img src="image_8.png" /><figcaption>Fig. 7. Equilibrium shape laws given by Taylor and the computations in terms of applied electric field root-mean-square. Equivalent droplet radius $r = 500~\mu\text{m}$ and $f = 2000~\text{Hz}$.</figcaption></figure>

which corresponds to small deformations from equilibrium. In Taylor's theory of droplets instability [2], the applied electric field $E_s$ is assumed to be continuous. Whilst, in our simulations, the applied electric field is time harmonic with an amplitude $E_{\text{max}} = E_s$.

We choose to work at $r = 500~\mu\text{m}$ and $f = 2000~\text{Hz}$. At this frequency, simulations have been compared with Taylor theory on the evolution of a spheroidal droplet characterized by its ellipticity with respect to the applied electric field (Fig. 7). We try to estimate how far the two stability curves are apart from each other. According to Taylor's [2], this law is two-valued. Surprisingly enough, the AC electric field law that we found is very close to the DC Taylor results, on the condition that the electric field root-mean-square is used. As in Taylor's theory, the upper curve is shown to be unstable and the lower curve, stable. These curves merge into a turning point, which corresponds to very close critical electric field values, in both cases.

## 4. Conclusions

Numerical simulations on dielectric droplets, more or less conductive, immersed in a highly insulating fluid, has been performed. An irrotational two-phase flow model, subjected to a harmonic electric field, has been solved by the BEM. The model presented is applied to high electric frequencies so the electric behavior is weakly coupled with the hydrodynamic one. The first part of our work consists in a comparison between the computational results and the experimental data. The main stability feature is correctly depicted in terms of variation of the critical electric field according to the droplet initial radius. In other respects, the critical electric field increases with respect to frequency (in the tested frequency range). As water droplet immersed in oil is considered, droplet conductivity seems to have no significant effect on the critical electric field as emphasized by experiment. Some instability regimes may develop according to the way the electric field is initially applied;
for slowly applied electric field, oblate droplet forms
appear, while for sharply applied electric field, cones
develop at the poles. The second part performs comparison
between numerical results and Taylor’s theory. Although it
was unexpected, close correspondence between both results
was found.

Acknowledgments

The authors are indebted to G. Berg for providing them
with his experimental data.

References

[1] G. Berg, L.E. Lundgaard, M. Becidan, R.S. Sigmond, Instability of electrically stressed water droplets in oil, Proceedings of the 14th International Conference on Dielectric Liquids, Graz, Austria, 2002, pp. 220–224.

[2] G.I. Taylor, Disintegration of water drops in an electric field, Proc. Roy. Soc. 280 (1964) 383-397.

[3] A. Castellanos, Electrohydrodynamics, Springer Series, Springer, Berlin, 1998.

[4] Y. Fouillet et, J.-L. Achard, Microfluidique discrète et biotechnologie, C. R. Phys. Acad. Française Sci 5 (2004) 577-588.

[5] D.W. Moore, The boundary layer on a spherical gas bubble, J. Fluid Mech. 16 (1963) 161-176.