Electrohydrodynamic Rayleigh-Taylor instability in leaky dielectric fluids

International Journal of Heat and Mass Transfer 109 (2017) 690–704

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Electrohydrodynamic Rayleigh-Taylor instability in leaky dielectric fluids Qingzhen Yang a,b, Ben Q. Li c,⇑, Feng Xu a,b,⇑ a The Key Laboratory of Biomedical Information Engineering of Ministry of Education, School of Life Science and Technology, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, PR China b Bioinspired Engineering and Biomechanics Center (BEBC), Xi’an Jiaotong University, Xi’an, Shaanxi 710049, PR China c Department of Mechanical Engineering, University of Michigan–Dearborn, Dearborn, MI 48128, USA

a r t i c l e

i n f o

Article history: Received 24 August 2016 Received in revised form 16 February 2017 Accepted 17 February 2017

a b s t r a c t In this paper, an analysis is presented of the influence of electric field on the Rayleigh-Taylor instability (RTI) of two leaky dielectric fluids. A numerical model is developed based on the coupled solution of the governing equations of the electric field (i.e., Poisson equation) and the fluid flow (i.e., Navier-Stokes equation), where the leaky dielectric model is employed to take into account the free charges in the fluids. These equations are formulated within the framework of phase field and solved simultaneously. The numerical model, after validated with existing data, is applied to study the effect of the electric field on interfacial morphology associated with the RTI. Different from the perfect dielectric model, the free charges originating from the conductivity in a leaky dielectric fluid plays an important role in the evolution of the RTI under the influence of an applied field. Various interfacial morphologies have been observed in the electrohydrodynamic RTI, suggesting that the electric field has a significant influence on the RTI. An electric field (horizontal or vertical) may suppress or aggravate the instability depending on the system parameters the permittivity ratio and the conductivity ratio of the two leaky dielectric fluids. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Rayleigh-Taylor instability (RTI) is a phenomenon that takes place at the interface of two stratified immiscible fluids when the denser fluid placed over the lighter one under the action of gravity, where gravity force breaks the equilibrium at the interface [1,2]. The nature remains the same when the problem is posed as a lighter fluid accreting into a heavier one as pointed by Taylor [3]. As a fundamental type of instability, the RTI occurs in nature and also has found applications in a wide range of areas, such as geophysics [4], astrophysics [5] and inertial confinement fusion [6]. In many of these applications, the RTI is influenced by the presence of other external fields (or forces), such as an external electric field [7–11] or a magnetic field [12–15]. In this study, we focus our interest on the RTI in an electric field, which sometimes is referred to as the electrohydrodynamic RTI. The subject of electrohydrodynamics, which is concerned the effect of electric field on fluid motion, was first considered by ⇑ Corresponding authors at: The Key Laboratory of Biomedical Information Engineering of Ministry of Education, School of Life Science and Technology, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, PR China (F. Xu). E-mail addresses: [email protected] (B.Q. Li), [email protected] (F. Xu). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.02.049 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

Rayleigh [16] in the 19th century. The interest in the subject appears to wade quickly until more than half a century later when it was rekindled by a series of papers published by Taylor and Melcher in 1960s [3,17–19], owing to the technological advances demanding a better understanding of physics governing the electrohydrodynamic phenomena. In the theory of electrodynamics, materials are in general categorized into two types, conductors and insulator [20]. However, this viewpoint was contested by the experimental results of the droplet deformation immersed in another host fluid with the presence of an electric field [21]. The data showed that the electrical conductivity of these poorly conducting fluids (e.g., polymers), albeit being infinitesimal, plays a crucial role in the electrohydrodynamic process. In order to account for the behavior of conductivity, Taylor introduced the concept of leaky dielectrics, by which he was able to explain the discrepancy between the experimental observations and theoretical analysis [18]. In leaky dielectric model, there are free charges (including both positive and negative ones) in the material. Without the electric field, the charges are distributed randomly and the material is electrically neutral. However, once the electric field is applied, the positive charges will move along the electric field and the negative charges in the opposite direction. Thus the free

691

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

charges manifest themselves [22,23]. Thus, with the framework of the electrohydrodynamics of leaky dielectric fluids, an electric field affects fluid motion in two ways. One comes from the electrical force caused by the polarized charges (or the KortewegHelmholtz force), that is, the dielectric effect; and the other from the Coulomb force associated with the free charges, that is, the electrically conducting effect. Considerable amount of literature has appeared on the subject of leaky dielectrics since the concept was introduced in 1960s. Taylor and McEwan studied the stability of conducting and/or nonconducting fluids subjected to a normal electric field [19]. The effect of a tangential electric field on the stability of the fluid-fluid interface was investigated by Melcher and Schwartz employing a leaky dielectric model [17]. Mohamed and Shehawey analyzed the nonlinear effect of electrohydrodynamic RTI when there are preset free charges at the interface [10,11]. Recently, the electrohydrodynamic instability on a thin liquid film confined between two electrodes has been exploited to produce the micro/nano-scaled patterning in the liquid [24] or cause defects of thin film coating [25]. Kumar studied the electrohydrodynamic instabilities in thin liquid films under an AC voltage [22], which can be extended to the RTI with ease. The above works on the subject of electrohydrodynamic RTI have employed the classical analytical perturbation method. Although capable of predicting the onset of the stability, the method is inadequate to provide information on the interfacial morphology and flow pattern, especially when the interface deformation evolves dynamically with fully nonlinear effects. A reliable and accurate, fully nonlinear numerical model is required to analyze the electrohydrodynamic RTI in the nonlinear evolution regime. To model the two phase fluid flow encountered in the RTI, the phase field method is adopted to account the flow and free interface deformation in this study. Instead of assuming a sharp interface with zero thickness, phase field treats the interface as a thin diffusive layer separating the two fluids. In this approach, the two bulk fluids and their interface can be treated in a unified manner through an introduction of the phase field parameter [26]. A major merit of phase field is remeshing-free during calculation and hence is employed broadly in multiphase flow problems [27,28]. In this sense, the phase field approach falls into the category of fixed grid methods. Moreover, the phase field model of a two-phase fluid interface in general is considered more physically realistic for small scale interfacial problems, yet a sharp interface employed in the free-moving grids methods represents a mathematical idealization of the interface. As such, it has been employed to simulate the conventional RTI by Ding et al. [29], Celani et al. [30], and Jacqmin [26]. The phase field methodology is also extended to the study of the electrohydrodynamic multiphase flow by incorporating the calculation of electric field, see a 2D modeling by Lin et al. [31] and a 3D case by Yang et al. [32]. Studies have shown that for the perfect dielectric fluids, a vertical electric field aggravates the RTI while a horizontal field suppresses it, which have been studied by linear perturbation analysis [7,8] and by numerical modeling [9]. In this paper, the interest is focused on the phase field modeling of electrohydrodynamic RTI of leaky dielectric fluids, with an emphasis on the fully nonlinear effects during the dynamic development of fluid-fluid interface in the presence of an electric field. To the best of our knowledge, such a study does not appear to have been carried out elsewhere. The paper is organized as follows. First the mathematical formulation of electric field, fluid flow and phase field for leaky dielectric fluids is given, which is followed by the computational details. The numerical phase field model is tested and validated against available theoretical results and linear perturbation analysis. After validation, the model is employed to simulate the electrohydrodynamic RTI. Extensive computations are carried out

for different parameters, such as the direction of electric field, electrical permittivity ratio of two fluids e2/e1 and also conductivity ratio r2/r1 (1 refers to the bottom fluid, 2 the top one). It is found that the electric field plays a dual role on the RTI, it may suppress or aggravate the instability depending on the detailed parameters. When the electrical permittivity ratio is equivalent to the conductivity ratio, the leaky dielectric model reduces to a perfect dielectric model, free charge is eliminated and only polarized charge effect exists. 2. Problem statement Consider a two layer liquid confined between a pair of horizontal infinite plates separated by distance l, subject to an external electric field, as shown in Fig. 1. Two fluids each are assumed immiscible and incompressible, with a difference in their properties, such as density, viscosity, permittivity and conductivity (see Fig. 1), which are denoted correspondingly by qi, li, ei and ri (i = 1, 2 for fluids 1 and 2, resepectively). The fluid-fluid interface is also subject to the action of interface surface tension c. As usual, for instability studies, a sinusoidal perturbation is imposed upon an initially flat surface. Our interest here centers on the electrohydrodynamic RTI for the case with q1 < q2, for which the interface is unstable due to the downward gravity. Since the electric field in different directions produces different effects, the effect of horizontal electric field (Ex) and vertical one (Ey) shall be discussed separately. 3. Mathematical formulation and numerical solution The above electrohydrodynamic RTI represents a coupled system of electric field and fluid flow with a moving interface demarking the two fluids. The dynamic development of the free surface is modeled by phase field. The mathematical formulation of the problem is described by a set of governing equations as follows. 3.1. The phase field equation Different from a sharp interface method, the phase field approach treats the fluid-fluid interface as a thin diffusive layer. We define a phase parameter u to distinct these two phases (u = 1 represents one phase and u = 0 for another one). Other physical parameters (e.g., electrical permittivity, conductivity, viscosity and density) are represented by a function of u across the interface (a linear function is employed in this study). Thus, the bulk fluids together with the diffusive interface can be treated as a single domain with continuous physical properties. For the system under consideration, the Ginzburg-Landau form of free energy density ffree is applicable, which takes the form of [33]

f free ðuÞ ¼

1 1k 2 kjruj2 þ u ð1  uÞ2 2 4n

2XT

ð1Þ

where n is the capillary width representative of liquid-liquid interface thickness and the sharp interface recovers when n ? 0. pffiffiffi At equilibrium, k satisfies the relationship k ¼ 6 2nc, where c is the surface tension coefficient. In the above equation, X is the computational domain, and T the time duration. The total free energy is an integration of the free-energy density over the entire domain X, which can be expressed as follows,

Z 

Z F free ¼

X

f free ðuÞdX ¼

X

 1 1k 2 kjruj2 þ u ð1  uÞ2 dX 2 4n

ð2Þ

A variational procedure [26,34] leads to the bulk chemical potential G which is the differential of the free energy density with respect to the phase parameter u, viz,

692

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

Fig. 1. Schematic sketch of the two-fluid system for the study of electrohydrodynamic RTI. The interface separating two fluids is represented by f(x). Both horizontal electric field (Ex) and vertical one (Ey) are considered.

G¼

  df free ðu  1Þu ¼ k r2 u þ du n2

ð3Þ

We trace the variation of phase parameter u which could reveal the interfacial evolution. The phase parameter is governed by the convective Cahn-Hilliard equation [26], !

ut ¼  u ru þ Mr2 G

ð4Þ

!

where u denotes the fluid velocity and M represents the mobility determining the diffusivity. Essentially, Eq. (4) describes the conservation of phase parameter. On the right hand side, the first term stands for the convection effect and the second one diffusion. 3.2. The electric field equation In a leaky dielectric model, both the polarized charges and free charges are considered. For the electrohydrodynamic RTI, there are two bulk fluids with distinctive properties. According to the Gauss’s law, the governing equation is [31,32], !

e0 r  ðer ðuÞ E Þ ¼ qe

ð5Þ

where e0 is the electrical permittivity of vacuum, er(u) the dielectric !

constant of the fluid, qe the free charge density and E electric field intensity. The analysis here is restricted to the DC fields and the free charge motion is usually slow. Hence the magnetic effect can be neglected and the electric field is assumed to be irrotational

is thought to be satisfied here and more discussions on the charge relaxation can be found elsewhere [38,39]. Thus, the transient and convection terms are expected to be of minor importance and may be neglected in the analysis. The governing equations of the electric field therefore are simplified,

r  ðrðuÞrVÞ ¼ 0

ð8aÞ

qe ¼ r  ðe0 er ðuÞrVÞ

ð8bÞ

Here, Eq. (8a) describes the electric potential and Eq. (8b) the calculation of free charges.A special case is worth noting. If the conductivity ratio of two fluids equals the permittivity ratio (r2/r1 = e2/ e1), then Eq. (8a) is equivalent to r  ðer ðuÞrVÞ ¼ 0, which leads to the condition that qe = r  (e0er(u)rV) = 0. For this case the leaky dielectric model reduces to the perfect dielectric model and free charge becomes ignorable. 3.3. The fluid flow field equation The fluid flow is governed by the laws of mass and momentum conservation. Due to the presence of a diffuse interface, the Navier-Stokes equation needs to be changed accordingly. As to the RTI system, both fluids are considered incompressible [29], and the external source consists of gravity, electrical force and also surface tension. The governing equations are of the form, !

ru ¼0

ð9aÞ

!

(r  E ¼ 0), allowing us to approximate the electric field as a gradi!

ent of the electrical potential E ¼ rV. Thus, Eq. (5) can be simplified to the Poisson equation,

e0 r  ðer ðuÞrVÞ ¼ q

e

ð6Þ

where V represents the electric potential. As for the free charge density qe, it could be determined by the law of charge conservation, which is expressed as, !

qet ¼  u rqe þ r  ðrðuÞrVÞ !

!

!

!

!

!

!

!

ð9bÞ where



!

fe ¼ r 

!

!

1 2



!

e0 er E  E  e0 er j E j2 di;j

!

¼ qe E 

ð7Þ

with u being the fluid velocity and r(u) the conductivity. The first term on the right side stands for the convection effect and, the second one accounts for the electromigration. It is worth pointing out that the diffusion effect is ignored here (Saville, 1997). To obtain the free charge distribution from Eq. (7), the initial value for qe is needed; in this study it is assumed zero everywhere. Further simplification also may be possible. For many applications [31,35–37], the time scale of the flow t c ¼ Lc =U c (Lc is the length scale, and Uc the characteristic velocity) is much larger than the time scale of charge relaxation t r ¼ e0 er =r. This condition also

!

qðuÞðu t þ u r u Þ ¼ rp þ r½lðuÞðr u þr uT Þ þ qðuÞ g þ f e þ f c

1 ! 2 e0 j E j rer ðuÞ 2

ð10Þ

!

f c ¼ Gru

ð11Þ

!

with f e denoting the total electrical force in the system. The term !

qe E is the Coulomb force due to the free charges and the one !

e0 j E j2 rer ðuÞ=2 represents the dielectric force due to the polar!

ized charges. The body force f c is the diffuse-interface equivalent !

of the interfacial tension [40]. Also, g represents the gravity, di,j identity matrix, p the pressure and l stands for the viscosity. These

693

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

equations are applicable to the bulk region of each fluid and also the diffusion interface. 3.4. The dimensionless parameters To ease the analysis, we define some dimensionless parameters. In the RTI, one of the key factors is the density difference, represented by the Atwood ratio A ¼ ðq2  q1 Þ=ðq2 þ q1 Þ. Some characpffiffiffiffiffiffiffiffiffiffi teristic parameters are chosen, length Lc = w, velocity U c ¼ Agw , pffiffiffiffiffiffiffiffiffiffiffiffi time t c ¼ Lc =U c ¼ w=Ag and free charge qec ¼ e0 e2 V 0 =w2 . The aforementioned governing equations are normalized by these parameters. The following dimensionless parameters arise during the normalization of the governing equations,

pffiffiffiffiffiffiffiffiffi

Re ¼ Pe ¼

q2 w3 g c e0 e2 V 20 ; We ¼ ; Eb ¼ ; 2 q2 gw l2 q2 gw3

pffiffiffiffiffiffiffiffiffi n w3 g ; Mc

Cn ¼

n : w

The properties of the denser fluid (i.e., fluid 2 in Fig. 1) are chosen as the denominators. As gravity plays an important role in RTI, the importance of other forces is defined as parameters by comparing against gravitational force. For instance, the electrical Weber number (Eb) reflects the relative magnitude of gravitational force and electrical force. As to the case of vertical field, 2

Eb ¼ e0 e2 V 20 =ðq2 gwl Þ. In phase field model, two more dimensionless parameters are needed; the Peclet number (Pe) denotes the ratio of the convective over diffusive mass transport and the Cahn number (Cn) defines the relative magnitude of the diffusive interface thickness compared to the characteristic length. The viscosity, electrical permittivity and conductivity ratios of the two fluids are non-dimensionalized as kl = l2/l1, ke = e2/e1 and kr = r2/r1, respectively. Employing these dimensionless parameters, the above governing equations can be recast as !

C t  þ u   r C 

1 2  r G ¼0 Pe

e







q ¼ r  ðer ðuÞr V Þ 

!

r u ¼0 !

Agw



ð16Þ

c

1 au2 ð1  uÞ , b is a shape factor, b = 1 for f ðuÞ ¼ 12 ajr uj Cn þ 4Cn horizontal electric field and, b = l/w for vertical electric field. The dimensionless density, viscosity, electrical permittivity and conductivity are 2

!

evolve starting from a quiescent state, u  ¼ 0. The interface is assumed y⁄(x⁄) = 2 + 0.025cos(2px⁄) which represents a planar interface superimposed by a perturbation of wave number k = 1 with an amplitude 0.025. The phase parameter is set u = 0 below the interface and u = 1 above it. 3.6. The numerical methodology

4. Linear stability analysis

To differentiate with the dimensional parameters, the dimensionless parameters are marked by the superscription *. The aforementioned dimensionless parameters are defined as qffiffiffiffi ! ! !  e e 2   u t  ¼ ttc ¼ t Ag , u  ¼ Uuc ¼ pffiffiffiffiffiffi , V  ¼ VV 0 , qe ¼ qqe ¼ eq0 e2wV 0 , G ¼ dfdu and w 

!

the side boundaries, @u  =@n ¼ 0. For pressure, the derivatives in the normal direction all are set to zero, @p =@n ¼ 0. If the electric field is oriented horizontally, the Dirichlet conditions are applied at the two side boundaries V⁄ = 1 at x⁄ = 1 and, V⁄ = 0 at x⁄ = 0 and the Neumann boundary conditions are imposed at the top and bottom boundary@V  =@n ¼ 0. For a vertical electric field, the Dirichlet boundary conditions are imposed at the top and bottom boundaries V⁄ = 1 at y⁄ = l/w and, V⁄ = 0 at y = 0 and the Neumann conditions at the two side boundaries @V=@n ¼ 0. About the initial conditions, the entire system is assumed to

ð14Þ

 !

Aq ðuÞ½u t þ ðu  r Þu  pffiffiffi ! ! A    ¼ r p þ r  ½l ðuÞðr u  þ r u T Þ Re ! 1  þ q ðuÞ j Eb b2 qe r V   Eb b2 jr V  j2 r er ðuÞ 2 þ WeG r u

!

applied, u  ¼ 0. The Neumann boundary conditions are applied at

ð13Þ

ð15Þ

!

The boundary and initial conditions are necessary for the numerical solution of the governing equations. For phase field, the Neumann boundary conditions are imposed at all boundaries for phase parameter and chemical potential, @ u=@n ¼ 0, and @G =@n ¼ 0, where n denotes the normal derivative; this guarantees the conservation of total mass. For fluid flow, the top and bottom boundaries are solid walls, and hence no slip conditions are

All the above governing equations are solved numerically by modifying the in-house FORTRAN code program. Details of the numerical scheme are described elsewhere [32]. In brief, the equations are discretized and solved iteratively for each time step by the finite difference method. The explicit scheme is employed to discretize the temporal terms, and the first order terms (such as the convection terms) are discretized by the upwind scheme to ensure the numerical stability and central difference scheme is employed for the second order terms.

ð12Þ

r  ðr ðuÞr V  Þ ¼ 0

3.5. The boundary conditions and initial conditions

2

q ðuÞ ¼ ð2Au þ 1  AÞ=ð1 þ AÞ

ð17aÞ

l ðuÞ ¼ u þ ð1  uÞ=kl

ð17bÞ

er ðuÞ ¼ u þ ð1  uÞ=ke

ð17cÞ

r ðuÞ ¼ u þ ð1  uÞ=kr

ð17dÞ

Prior to presenting the detailed numerical analysis, a classical linear stability analysis is performed to gain a fundamental understanding of the RTI of leaky dielectric fluids in electric field. The linear analysis also serves as a platform against which numerical computations are compared. For the configuration in Fig. 1, the governing equations are the Laplace equation for electric field, the continuity equation for incompressible fluids and the N-S equation for fluid flows, which take the following forms [17],

r2 V i ¼ 0

ð18Þ

!

r  ui ¼ 0 !

!

ð19Þ !

!

!

qðu i;t þ u i  ru i Þ ¼ rpi þ li r2 ui þ qi g

ð20Þ

In the linear stability analysis, the sharp interface model is employed. The electric field and fluid flow are sufficient to describe the system and hence the phase field is trivial. The linear stability analysis is carried out for both the horizontal and vertical electric fields. 4.1. Linear stability analysis of RTI in a horizontal electric field The analysis starts with an initially flat interface, the zero flow field and a uniformly imposed electric field everywhere. These base states are written as

694

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

V 1 ¼ V 2 ¼ E0 x

ð21Þ

! u1

ð22Þ

!

¼ u2 ¼ 0

p1 ¼ q1 gy þ ðq1  q2 ÞgH0 p2 ¼ q2 gy 

ð23Þ

1 e0 ðe2  e1 ÞE20 2

ð24Þ ð25Þ

f ðxÞ ¼ H0

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    2 @f @f n ¼  ;1 þ 1; @x @x

!

With the above base states and perturbed states substituted into the governing equations and the boundary conditions applied, a system of equations is obtained. The condition for non-trivial solutions to this system leads to the following dispersion equation,

" # 3 2 gkða2  a1 Þ k c C 3 k e0 E20 ðe2  e1 Þðr2  r1 Þ  2  C3  s2 s2 ðr2 þ r1 Þðq1 þ q2 Þ s ðq1 þ q2 Þ  ðC 1 q2 a1 þ C 2 q1 a2  kÞ  4C 3 ka1 a2

To the above base state, a small perturbation is added,

^ 1 ðyÞestjkx V 1 ¼ E0 x þ V

ð26Þ

^ 2 ðyÞestjkx V 2 ¼ E0 x þ V

ð27Þ

! u1

! ^ ¼ u 1 ðyÞestjkx

ð28Þ

! u2

! ^ ¼ u 2 ðyÞestjkx

ð29Þ

^1 ðyÞestjkx  q1 gy þ ðq1  q2 ÞgH0 p1 ¼ p ^2 ðyÞestjkx  q2 gy  p2 ¼ p

ð30Þ

1 e0 ðe2  e1 ÞE20 2

f ¼ H0 þ ^f estjkx

¼0

@V 2 ¼ 0; @y

! u2

¼ 0 at y ¼ l

at y ¼ 0

ð33Þ ð34Þ

@f ! @f þ u 1;x @t @x

ð35Þ

u 2;y ¼

@f ! @f þ u 2;x @t @x

ð36Þ

For the electric field, the tangential component is continuous across the interface. In the normal direction, the continuity of rrV is imposed, which can be derived from Eq. (8a). These conditions are stated as !

!

n ðr1 rV 1 Þ ¼ n ðr2 rV 2 Þ

!

ð37Þ

!

t rV 1 ¼ t rV 2

ð38Þ

The normal stress at the interface needs to be balanced, namely, !

!

!

!

ð2Þ rð1Þ ij  n rij  n ¼ cðrs  n Þ  n

where

C 1 ¼ tanhðkH0 Þ

The first one 

¼

!

where

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 q2 ¼ k þ sq2 =l2 ,

ð32Þ

The fluids are immiscible, that is, neither of the fluids is able to cross the interface; this geometrical constraint leads to the two kinematic conditions at the interface y ¼ f ðxÞ  H0 ! u 1;y

ð40Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 q1 ¼ k þ sq1 =l1 ,

ð31Þ

where the symbol indicates the amplitude of the perturbation, s is the growth rate and k = 2p/k (k the wave length) is the wave number. The top and bottom boundaries of the domain are confined by solid walls and, no slip boundary conditions are imposed for the fluid flow. For the electric field, the Neumann conditions are applied [9], ! u1

2

4C 3 k ða2 m2  a1 m1 Þ½C 1 q2 a1  C 2 q1 a2 þ C 3 kða2  a1 Þ s 3 4C 3 k þ 2 ða2 m2  a1 m1 Þ2 ðq2 C 1  kÞðq1 C 2  kÞ s 2 k e0 E20 ðe1 r2  e2 r1 Þ þ 4C 3 2 ½C 1 q2 a1  C 2 q1 a2 þ C 3 kða2  a1 Þ ¼ 0 s ðr2 þ r1 Þðq1 þ q2 Þ

þ

cothðq2 H0 Þ, C 2 ¼ tanhðkH0 Þ cothðq1 H0 Þ,C 3 ¼ cothðkH0 Þ, a1 ¼ q1 = ðq1 þ q2 Þ and a2 ¼ q2 =ðq1 þ q2 Þ. It is worth pointing out that l = 2H0 has been set in this case. Inspection of the above equation indicates that the electric field affects the instability via two terms.

^

@V 1 ¼ 0; @y

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    2 @f @f t ¼ 1; þ 1: @x @x

!

ð39Þ ! !

rijð1Þ and rijð2Þ are the stress tensors, rs ¼ r  n ðn rÞ is the

surface divergence operator and, the normal and tangential vectors of the interface are given by

C 3 k2 e0 E20 ðe2 e1 Þðr2 r1 Þ ðC 1 q2 s2 ðr2 þr1 Þðq1 þq2 Þ

a1 þ C 2 q1 a2  kÞ represents the electrical force normal to the interface. When e1 = e2 or r2 = r1, the electrical force becomes uniform along the interface and thus has no influence on the instability. The second term 4C 3

k2 e0 E20 ðe1 r2 e2 r1 Þ ½C 1 q2 s2 ðr2 þr1 Þðq1 þq2 Þ

a1  C 2 q1 a2 þ C 3 kða2  a1 Þ originates from

the tangential fluid flow due to the tangential electrical force. This term in essence describes the effect resulting from the free charges accumulated at the interface. It is important to note that when r2/r1 = e2/e1, the free charges vanish and no tangential flow occurs. For the case of e2/e1 = r2/r1, the leaky model reduces to the perfect dielectric model and Eq. (40) simplifies to the following form,

" # 3 2 gkða2  a1 Þ k c C 3 k E20 e0 ðe2  e1 Þ2  2  C3  2 s2 s ðq1 þ q2 Þ s ðe2 þ e1 Þðq1 þ q2 Þ  ðC 1 q2 a1 þ C 2 q1 a2  kÞ  4C 3 ka1 a2 2

4C 3 k ða2 m2  a1 m1 Þ½C 1 q2 a1  C 2 q1 a2 þ C 3 kða2  a1 Þ s 3 4C 3 k þ 2 ða2 m2  a1 m1 Þ2 ðq2 C 1  kÞðq1 C 2  kÞ ¼ 0 s

þ

ð41Þ

For the RTI in a horizontal electric field, Melcher & Schwarz developed similar linear stability analysis for the leaky dielectric fluids. In their model, the governing equation is the charge conservation (similar with equations 6 and 7) and the vertical depth is extended to infinity (see equation 34 in [17]). In Melcher & Schwarz (MS) model, the governing equations of electric field are,

(

r  ðrEÞ ¼ 0 !

!

!

!

!

!

qt þ u rs q  q n ðn rÞ u þkr E k  n ¼ 0 !

!

with the associated boundary conditions, ke E k  n ¼ q. In this paper, we assume the time scale of the flow t c ¼ Lc =U c is much larger than the time scale of charge relaxation t r ¼ e0 er =r, e.g., r  e0 er Uc =Lc . The transient and convection terms are negligible !

!

in the transportation equation which is simplified as kr E k  n ¼ 0. Thus the governing equation of electric field is simplified to !

!

r  (rE) = 0 with the boundary condition kr E k  n ¼ 0. The original

695

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704 !

!

boundary condition ke E k  n ¼ q becomes redundant.Here we shall compare the results of the present leaky dielectric model with those of the MS model and also of the perfect dielectric model. For this purpose, the parameters are chosen as q1 = 1000 kg/m3, q2 = 3000 kg/m3, l1 = l2 = 0.01 Pas, e1 = 1, e2 = 3, r2/r1 = 1.5, E0 = 2  106 V/m, and the depth is assumed as l ! 1, which yields C1 = C2 = C3 = 1. In the MS model, both r1 and r2 influence the final results. However, in the present case, the conductivity ratio (r2/r1) affects the stability. Thus, to facilitate comparison, r2/r1 is fixed with r1 changing. The results are illustrated in Fig. 2, which shows that the present leaky model is a good approximation of the MS model when the conductivity is high. Specifically, the MS model becomes essentially the same as the present leaky model for r1 > 107 S/m; and the MS model also reduces to the perfect dielectric model when r1 < 1011 S/m. 4.2. Linear stability analysis on RTI in a vertical electric field For the case of imposed vertical electric field, the system also evolves starting from quiescent state (i.e., the flow is zero) with a flat interface. The electric field is uniform everywhere. These base states can be written as,

V1 ¼ V2 ¼ !

D0

r1 D0

r2

ð42Þ

y

yþ

ðr2  r1 ÞV 0 r2 þ r1

ð43Þ

!

u1 ¼ u2 ¼ 0

ð44Þ

p1 ¼ q1 gy þ ðq1  q2 ÞgH0

ð45Þ

p2 ¼ q2 gy 

1 e0 2





e1 e2 2 D  r21 r22 0

f ðxÞ ¼ H0

ð46Þ ð47Þ

the base state, the following expressions are obtained,

D0

r1

^ 1 ðyÞestjkx yþV

D0

r2

yþ

ðr2  r1 ÞV 0 ^ þ V 2 ðyÞestjkx r2 þ r1

ð49Þ

! u1

! ^ ¼ u 1 ðyÞestjkx

ð50Þ

! u2

! ^ ¼ u 2 ðyÞestjkx

ð51Þ

^1 ðyÞestjkx  q1 gy þ ðq1  q2 ÞgH0 p1 ¼ p   1 e e ^2 ðyÞestjkx  q2 gy  e0 12  22 D20 p2 ¼ p 2 r1 r2

ð52Þ

f ¼ H0 þ ^f estjkx

ð54Þ

ð53Þ

where k = 2p/k is the wave number. The upper and bottom boundaries are considered as walls and, no slip boundary conditions are imposed for the fluid flow. For the vertical electric field, the Dirichlet conditions are imposed at the walls,

V 1 ¼ 0;

! u1

¼ 0 at y ¼ 0

ð55Þ

u 2 ¼ 0 at y ¼ l

ð56Þ

!

V 2 ¼ V 0;

The kinematic conditions, stress balance equations and the boundary conditions for the electric field along the interface,y ¼ f ðxÞ  H0 , are the same as stated above for the horizontal electric field. For electric field, the same boundary conditions are imposed at the interface. Also the force is balanced at the interface. With the above perturbations and boundary conditions substituted, one has the dispersion equation for instability,

"  # 3 2 gkða2  a1 Þ k c C 4 k D20 e0 ðr2  r1 Þ e2 e1  2   C3 þ 2 s2 s ðq1 þ q2 Þðr2 þ r1 Þ r22 r21 s ðq1 þ q2 Þ  ðC 1 q2 a1 þ C 2 q1 a2  kÞ  4ka1 a2 2

4C 3 k ða2 m2  a1 m1 Þ½C 1 q2 a1  C 2 q1 a2 þ kða2  a1 Þ s 3 4C 3 k þ 2 ða2 m2  a1 m1 Þ2 ðq2 C 1  kÞðq1 C 2  kÞ s 2 4C 3 k e0 ðe2 r1  e1 r2 Þ D20 þ 2 ðq1 þ q2 Þðr1 þ r2 Þ r1 r2 s

þ

0 where D0 ¼ ðrr11þrr22VÞH . With an infinitesimal perturbation added to 0

V1 ¼

V2 ¼

ð48Þ

 ½C 1 q2 a1  C 2 q1 a2 þ kða2  a1 Þ ¼ 0 where

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 q1 ¼ k þ sq1 =l1 ,

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 q2 ¼ k þ sq2 =l2 ,

ð57Þ C 1 ¼ tanhðkH0 Þ

cothðq2 H0 Þ, C 2 ¼ tanhðkH0 Þ cothðq1 H0 Þ,C 3 ¼ cothðkH0 Þ, a1 ¼ q1 = ðq1 þ q2 Þ and a2 ¼ q2 =ðq1 þ q2 Þ. Also, the assumption l = 2H0 has been used in this case. Once again, when e2/e1 = r2/r1, the leaky model becomes the perfect dielectric model, and Eq. (57) simplifies to

" # 3 2 gkða2 a1 Þ k c C 4 k D20 e0 ðe2  e1 Þ2  2  C3 þ s2 s ðq1 þ q2 Þ e1 e2 s2 ðq1 þ q2 Þðe2 þ e1 Þ  ðC 1 q2 a1 þ C 2 q1 a2  kÞ  4ka1 a2 2

4C 3 k ða2 m2  a1 m1 Þ½C 1 q2 a1  C 2 q1 a2 þ kða2  a1 Þ s 3 4C 3 k þ 2 ða2 m2  a1 m1 Þ2 ðq2 C 1  kÞðq1 C 2  kÞ ¼ 0 s

þ

ð58Þ

5. Numerical results 5.1. Model validation Fig. 2. Comparison of Melcher & Schwarz’s model (solid lines) with the present leaky model and the perfect dielectric model.

In order to test the numerical model, we considered the classical problem of a droplet deformation with the presence of

696

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

an electric field, the solution to which has been studied analytically by Taylor [18] and, is widely employed as the benchmark test. The deformation of a droplet in an electric field is characterized by [31,32],

D¼

LB 9 Eb f d ðrd =rm ; ed =em ; ld =lm Þ ¼   L þ B 16 We ð2 þ rd =rm Þ2

ð59Þ

where L and B are the lengths of the parallel and perpendicular axes of the droplet with respect to the imposed electric field and, subscripts m and d denote the surrounding host fluid and the droplet, respectively. The discriminating function, f d ðrd =rm ; ed =em ; ld =lm Þ, is expressed as follows,

 fd





rd ed ld rd ¼ ; ; rm em lm rm

2



ed



þ12 em   3 rd ed ð2 þ 3ld =lm Þ þ  5 rm em ð1 þ ld =lm Þ

ð60Þ

where the droplet takes a spherical shape for fd = 0, an oblate shape for fd > 0 and a prolate shape for fd < 0. The present numerical and analytical results (equation 60) are compared in Fig. 3. In the numerical simulation, a mesh of 201  201  201 is chosen and the computational domain is fixed as (6R  6R  6R) where R is the radius of the droplet. Apparently, they agree well with each other except when the droplet deformation becomes severe. The reason lies in that Taylor’s analytical solution is based on the small deformation assumption, which no longer holds true for a large deformation [31,32].Linear stability analysis presented in the last section is used here to further validate the present numerical model. It is worth noting that the numerical model we developed in Section 3 is a 3D general model. To compare with the linear stability analysis we would simplified the model to 2D. The geometry of the interface is given by f(x, t) = estcos(kx), where the real part of s represents the growth rate of perturbation. A value of s < 0 indicates that the perturbation vanishes with time and the system is stabilized. A positive s stands for a destabilized system. The electrohydrodynamic RTI in the horizontal electric field is chosen as the testing case with gravity and surface tension ignored. The numerical phase field model presented above is also capable of determining the growth rate. The phase field computation starts with a cosine perturbation added to the flat interface. The two phase interface is then captured as the contour u = 0.5 and monitored during the numerical solution. It is worth mentioning that, the initial amplitude of perturbation Dh should be sufficiently small to fall

Fig. 3. Comparison of analytical and numerical solutions for deformation of a leaky dielectric droplet in an electric field as a function of ed/em: Eb/We = 0.2; rd/rm = 3.0; ld/lm = 1.0.

in the linear phase (i.e., Dh k, k being the wavelength) and, Dh is sufficiently large compared with mixing interface width (i.e., Dh  n, n being the interfacial thickness) [30]. A growth rate was obtained by the best fit of function f(x, t) and the results obtained by linear stability analysis and numerical method are shown in Fig. 4. Clearly, agreement is gratifying, thereby further validating the numerical model. For the numerical simulation in Fig. 4, the mesh is chosen as 241  961. 5.2. Stability analysis Having been validated as discussed above, the phase field numerical model is now employed to study the RTI of the leaky fluids in electric fields. The results for the horizontal electric field are discussed first, followed by those for the vertical electric field. A set of default parameters used for calculations presented in this paper are given in Table 1. In this study, the Peclet number (Pe) depends pffiffiffiffiffiffiffiffiffi on the surface tension, and the term n w3 g =M is taken 2.5, leading to Pe = 2.5/c. To emphasize the effect of the electric field on the RTI, we would conduct the 2D numerical simulation as the first step and the 3D effect shall be studied in the future work. 5.2.1. Electrohydrodynamic RTI in a horizontal electric field It is known that for perfect dielectric fluids the horizontal electric field suppresses the RTI [7,8], for which only polarized charges induced at the interface are responsible. For leaky dielectric fluids, the effect of electric field becomes more involved as both polarized charges and free charges now are present, the latter being a result of considering electrical conductivity. The free and polarized charges generate the electrical force at the interface and consequently influence the instability and the development of interfacial morphology. Emphasis here is then on the effect of the electrical conductivity on the electrohydrodynamic RTI. To initiate the numerical calculations, a cosine wave perturbation is added to the flat interface and the initial perturbation is given by y⁄(x⁄) = 2 + 0.025cos(2px⁄). Then the bottom of falling fluid (i.e., falling tip, see Fig. 1) and the top of rising fluid (i.e., rising tip) are monitored during the time development of the interface deformation. The structure height Dh is used to characterize the instability,

Dh ¼ yr  yf

ð61Þ

where yr⁄ and yf⁄ represent the positions of rising and falling tips, respectively.

Fig. 4. Growth rate of perturbation in horizontal electric field calculated by linear stability analysis and numerical model. The gravity and surface tension are not considered. Other parameters are listed in Table 1.

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704 Table 1 Parameters or constants used in simulation. Parameters

Value

Permittivity of vacuum (e0) Atwood ratio (A) Reynolds number (Re) Cahn number (Cn) Weber number (We) Ratio of viscosity (kl) Ratio of electrical Permittivity (ke) Ratio of electrical conductivity (kr) Ratio of height/width (l/w)

8.85  10–12 F/m 0.50 1000 0.03 103 1 3 5 4

Table 2 The mesh sensitivity of the RTI.a

a

Mesh

Dh

Error (%)

31  121 61  241 121  481 241  961 481  1921

0.1187 0.1228 0.1243 0.1252 0.1269

6.45 3.20 2.07 1.34

The case 481  1921 is chosen as the base for comparison.

In order to test the dependence of numerical results on the mesh, the RTI without electric field is simulated with different meshes. The parameter are chosen as in Table 1 and the structure height Dh at dimensionless time t⁄ = 0.7 are listed in Table 2. As can be seen from this table, the numerical result is of high accuracy when the mesh is 241  961 and thus such a mesh would be employed hereafter in this study. We characterized the stability of the system in terms of structure height Dh at dimensionless time t⁄ = 0.7, Fig. 5(a). The permittivity ratio is fixed at ke = e2/e1 = 3, with kr varied. In Fig. 5(a), a large Dh indicates that the system is unstable and a small Dh means the system is stable. If Dh(t⁄ = 0.7) is less than 0.05 which is the initial value at t⁄=0, the instability is suppressed by the electric field and the RTI would decay eventually. Compared with the case of electric field free (Eb = 0), the instability is aggravated by the electric field when kr is small (kr < 1.7), but is suppressed when kr is large (kr > 1.7). When kr is large enough (kr > 2.5), the electric field produces a strong stabilizing effect and suppresses the instability completely, i.e., the system is stabilized. In this system, the suppression or aggravation effect of electric field originates from the electrical force along the liquid-liquid interface. For this case, the Coulomb force originating from the free charges is mainly in the horizontal direction. The force caused by the polarized charges however orientates vertically and is respon-

697

sible for the influence of electric field on the RTI. To investigate the mechanism governing the deformation behavior of the interface, the electrical force in the y direction along the interface are calculated for different kr and, the forces, normalized by

e0 je2  e1 jV 20 =ð2w2 Þ, are depicted in Fig. 5(b). In phase field model, the interface is defined as the zone of 0.01 u 0.99 (u represents the phase parameter), the electrical force in Fig. 5(b) is an integration of electrical force density over the entire diffusive layer. For the interface instability under consideration here, it is the distribution of the force along the interface, rather than the absolute value of the force, that determines the effect on the surface deformation. Since both fluids are incompressible, the mass is conserved for each phase. The total electrical force is balanced by the hydrostatic pressure and, only the non-uniformity of the force affects the interfacial deformation. Although the interfacial morphology is evolving in time, the electrical force distribution at the initial interface would reveal the essential mechanism by which the interface responds to the applied field. Analyzing the force distribution for different kr in Fig. 5(b), one can conclude that the electric field enhances the RTI if kr is small (see kr = 0.5). On the other hand, an increase of kr leads to the effect of suppression. This is consistent with the observation of structure height shown in Fig. 5(a). To gain a perspective of the effect of electrical permittivity, computations were carried out with the ratio ke inverted, that is, ke = e2/e1 = 1/3, and the results are shown in Fig. 6(a). For this case, the electrical conductivity produces a different effect, where a smaller kr results in a stronger suppression effect and the instability is completely suppressed when kr is small enough. However, with increasing kr, the electric field would aggravate the RTI. Summarizing the results from these two cases (ke = 3 and ke = 1/3), one can conclude that a horizontal electric field plays a dual role in the RTI of leaky dielectric fluids. More specifically, the effect can be either suppressing or aggravating depends on the permittivity ratio ke and the conductivity ratio kr. In contrast, for pure dielectric fluids, there will be present no free charges and as a result a horizontal electric field produces only a suppressing effect on the RTI in these fluids. The effect of electric field strength, which is represented by the electrical Weber number Eb, is studied for two cases and the computed results are shown in Fig. 7. For the first case (ke = 3 and kr = 5), the horizontal electric field tends to suppress the instability and, this suppression effect becomes stronger with an increase of the electric field strength; in particular, the instability is completely suppressed when Eb > 0.22. For the second case (ke = 3 and kr = 0.5), the electric field generates a strong aggravation effect. This aggravation effect is enhanced by a high electrical Weber number. Thus, the effect of the electric field on the RTI (i.e., suppression or aggravation) is determined by ke and kr. A

Fig. 5. Electrohydrodynamic RTI in a horizontal electric field with electrical strength Eb = 0.35, electrical permittivity ratio ke = 3, and the conductivity ratio kr is varied to investigate its effect on the RTI. Other parameters are listed in Table 1. (a) Structure height at dimensionless time t⁄ = 0.7. (b) The electrical force at the fluid-fluid interface at t⁄ = 0.

698

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

Fig. 6. Electrohydrodynamic RTI in a horizontal electric field with electric field strength Eb = 0.35. The electrical permittivities of the two fluids are interchanged or ke = 1/3. Other parameters are listed in Table 1. (a) The structure height at dimensionless time t⁄ = 0.7. (b) The electrical force along the interface at t⁄ = 0.

at one point (kc = r2/r1 = 1.7). Compared with the case without an electric field (Eb = 0), the electric field tends to aggravate the instability when the conductivity ratio is smaller than the critical value of kc = 1.7. If the ratio is larger than kc, the instability would be suppressed. This critical value corresponds to the point where the electric field has no influence on the instability in despite of the electric field strength. This phenomenon may also be revealed through the analysis of the dispersion Eq. (40). The electric field affects the instability through two terms. The term A¼

C 3 k2 e0 E20 ðe2 e1 Þðr2 r1 Þ ðC 1 q2 s2 ðr2 þr1 Þðq1 þq2 Þ

a1 þ C 2 q1 a2  kÞ represents the effect

of electrical force normal to the interface that tends to induce the interfacial deformation, leading to the instability, whereas the

Fig. 7. The structure height of electrohydrodynamic RTI in horizontal electric field at dimensionless time t⁄ = 0.7 for two case kr = 5 and kr = 0.5. The electrical permittivity ratio is fixed at ke = 3 and Eb is varied to investigate the effect of electric field strength. Other parameters are listed in Table 1.

strong electric field would enhance this specific effect, be it suppression or aggravation. It is important to note that this differs from the perfect dielectric case where a horizontal electric field always suppresses the instability. To further validate the above numerical results, the linear stability analysis presented in Section 4.1 above is employed. For the case of e2/e1 = 3, the growth rate decreases with the conductivity ratio kr = r2/r1. As shown in Fig. 8(a), all the curves converge

term

B ¼ 4C 3

k2 e0 E20 ðe1 r2 e2 r1 Þ ½C 1 q2 s2 ðr2 þr1 Þðq1 þq2 Þ

a1  C 2 q1 a2 þ C 3 kða2  a1 Þ

stands for the tangential electrical force that causes the tangential flow and thus affect the instability. If we add these two terms together and let it be zero (i.e., A + B = 0), then the electric field becomes noneffective. For a set of given parameters, the equation A + B = 0 can be numerically solved and a critical kr = r2/r1 can be obtained. These two terms counterbalance each other when kc = kr. For the permittivity ratio of e2/e1 = 1/3, the growth rate increases with the conductivity ratio r2/r1. There also exists a critical conductivity ratio kc = 2.3, below which the electric field suppresses, and above which the field aggravates, the instability. This is consistent with the results from the numerical computations.

5.2.2. Electrohydrodynamic RTI in a vertical electric field Early studies on the RTI indicate that the direction of electric field can strongly affect the instability. For perfect dielectric fluids, it has been found that the horizontal electric field suppresses the

Fig. 8. The effect of electrical strength and conductivity ratio r2/r1 on the RTI by linear stability analysis: (a) Growth rate vs conductivity ratio r2/r1 in a horizontal electric field with e2/e1 = 3 and (b) e2/e1 = 1/3.

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

699

Fig. 9. The effect of conductivity ratio kr on the electrohydrodynamic RTI in a vertical electric field (ke = 3 and Eb = 0.05): (a) The structure height at dimensionless time t⁄ = 1.7 for different electrical permittivity ratios, where the dashed red line indicates the electric field free case; and (b) Distribution of electrical force along the interface at t⁄ = 0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. The effect of conductivity ratio kr on the electrohydrodynamic RTI in vertical electric field (ke = 1/3 and Eb = 0.05): (a) The structure height at dimensionless time t⁄ = 1.7 for different conductivity ratios, the dashed red line the E-free case; and (b) The distribution of electrical force along the interface at t⁄ = 0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

RTI, whereas the vertical electric field intensifies it. The effect of a horizontal field on the RTI of leaky dielectric fluids has been analyzed above; we may now turn our attention to the effect of a vertical electric field. Two cases are selected for discussion. For the first case, the permittivity ratio is set at ke = 3 and electric field strength at Eb = 0.05. To better exhibit the instability, the structure height at dimensionless time t⁄ = 1.7 are plotted with respect to the conductivity ratio, and results are shown in Fig. 9(a). Apparently, the conductivity ratio kr affects the instability significantly. The case kr = ke = 3 corresponds to a perfect dielectric case, wherein the free charges are zero, as stated before. As can be seen, when kr increases, the structure height decreases, attains a minimum at kr = 1.2 and then increases. A small structure height represents a weak aggravation effect, whereas a large one indicates a strong effect. Thus, the aggravation effect becomes weaker as kr increases up to 1.2; the effect of electric field may actually be considered suppressing in comparison with the field-free case (see the dashed red line in Fig. 9a). On the other hand, the aggravation effect increases with a further increase of kr above 1.2. Different from the case of the horizontal electric field discussed above, the aggravation effect of a vertical electric field is not a monotonic function of the conductivity ratio kr. To better understand the effect of the vertical field, the electrical force along the interface is calculated and plotted in Fig. 9(b), 2

where the force is normalized by the factor e0 je2  e1 jV 20 =ð2l Þ . As stated before, the hydrostatic pressure would counteract the total electrical force, and the non-uniform part of the electrical force acts on the interface to affect the instability. From Fig. 9(b)

it is clear that at kr = 0.5, the electrical force is positive (i.e., pointing upwards) but the magnitude is small at the falling point and large at the rising point. Consequently, the rising point is pushed upwards by the force while the falling tip goes downwards to satisfy the requirement of mass conservation, thereby aggravating the instability. With an increase of kr, the non-uniformity becomes smaller and thus the aggravation effect is reduced (see the curves of kr = 1, 2). A further increase of kr, however, leads to a severe non-uniformity of the electrical force distribution (see the curves of kr = 3, 5 and 10), and thus the aggravation effect is enhanced. For the second case, the electrical permittivity ratio was set at ke = 1/3 and the electrical strength Eb = 0.05. To show the dependence of structure height on the conductivity, the structure height at dimensionless time t⁄=1.7 is plotted in Fig. 10(a). Once again, the structure height decreases first and then increases, as kr increases. The minimum occurs at kr = 0.6. This is similar to the case given in Fig. 9. One different phenomenon in Fig. 10(a) is that the electric field always aggravates the RTI regardless of the conductivity ratio. The linear stability analysis may also be made of the effect of the vertical electric field case and the results are shown in Fig. 11. Clearly, the growth rate decreases first and then increases with increasing the conductivity ratio kr. Comparing with the electric field free case, at e2/e1 = 3, two critical conductivity ratios kc1 and kc2 are found at which the electric field becomes noneffective. If kr < kc1 or kr > kc2, then the electric field intensifies the instability. However, the field tends to suppress the instability when kc1 < kr < kc2. Similar to the horizontal electric field, the effect of a vertical field originates from two terms. The first term  C 4 k2 D20 e0 ðr2 r1 Þ e2  re12 ðC 1 q2 a1 þ C 2 q1 a2  kÞ represents the effect s2 ðq þq Þðr þr Þ r2 1

2

2

1

2

1

700

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

Fig. 11. Dependence of growth rate on the conductivity ratio in the presence of a vertical electric field by linear stability analysis: (a) e2/e1 = 3. (b) e2/e1 = 1/3.

of electrical force normal to the interface, and the second term 4C 3 k2 e0 ðe2 r1 e1 r2 Þ ðq1 þq2 Þðr1 þr2 Þ s2

5.3. Interfacial morphology of electrohydrodynamic RTI

D20

r1 r2 ½C 1 q2 a1  C 2 q1 a2 þ kða2  a1 Þ for the tangen-

tial electrical force which would induce the tangential fluid flow and affect the instability. Balancing these two terms yields the value of critical conductivity ratio where the electric field becomes noneffective. On the other hand, for the permittivity ratio of e2/ e1 = 1/3, no critical conductivity ratio is found (i.e., the electric field always aggravates the instability). The above results are consistent with the numerical calculations.

The stability of the interface under the action of electric field has been discussed in the previous section, in which the stability is characterized by the structure height Dh. Another important issue of the stability is the interfacial morphology. The interface could be of different morphology even with the same structure height. For the case without electric field, the gravity drags the falling tip downwards and forming a ‘‘spike” like structure. With the

Fig. 12. The interfacial morphology of the electrohydrodynamic RTI of leaky fluids and electrical force distribution along the interface in a horizontal electric field (Eb = 0.2, ke = 3) as a function of the conductivity ratio kr: (a) Snapshots of the interfacial morphology in the horizontal electric field at dimensionless time t⁄ = 2.5. The conductivity ratio kr, from left to right, are 1, 2, 3, 4 and 6. The corresponding structure heights Dh are 1.89, 1.10, 0.62, 0.35, and 0.12. (b) The vertical component of electrical force along the interface for different conductivity ratios. (c) The horizontal component of electrical force along the interface.

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

evolution of interfacial deformation, roll-ups appear at the two sides of the falling tip due to the nonlinear effect [1]. Daly [41] studied the interfacial morphology of the RTI for different density ratio. It was found that when the density ratio is approximately 1 (i.e., Atwood ratio A  0), the falling tip and rising tip evolve in the same pattern. However when the density ratio goes to infinity (i.e., Atwood ratio A  1), a narrow and long falling ‘‘spike” was observed and no roll-ups appear. This occurs because the drag applied on the fluid is proportional to the density of fluid that it goes into [41]. When the two fluids have approximately the same density, the drags on each fluid are the same. Thus the falling tip and rising tip evolve in the same pattern. However, if the density ratio goes to infinity, the drag on the falling fluid is ignorable. Hence the falling tip goes much faster than the rising tip, resulting in a long and thin falling tip due to the requirement of mass conservation. In this section, the interfacial morphology in an electric field is studied, with the focus on the effect of conductivity ratio kr (or leaky effect) on the interfacial deformation. For this purpose, the density ratio and electrical strength are fixed. 5.3.1. Interfacial morphology in horizontal electric field The interfacial morphology in the horizontal electric field is studied first. Fig. 12(a) plots the shapes of the interfacial morphology with different conductivity ratios at the same time as the perturbation sets in (see also Supplemental Movie 1 for the dynamic evolution). The electrical permittivity ratio is fixed at ke = 3 and

701

electrical strength at Eb = 0.2. Clearly, the structure height is large for kr = 1 and then decrease with an increase of kr, indicating that the instability is suppressed by the electric field and this effect increases with increasing kr. These results are consistent with the stability analysis in the previous section. Another phenomenon evident in Fig. 12(a) is that the falling ‘‘spike” is narrow for kr = 1, while it is broadening with increasing kr. For a better understanding of the role of the electric field in the deformation of the interfacial morphology, the electrical force distribution along the interface is calculated and the results for the horizontal and vertical components are discussed in Section 5.2.1 (see Fig. 5). In Fig. 5, the interface shape is taken to be the perturbation, y⁄(x⁄) = 2 + 0.025cos(2px⁄), at dimensionless time t⁄ = 0. The force distribution in Fig. 5(b) follows cosine functions, this Fig. helps understanding the stability at the initial time. With the development of the deformation, nonlinear effects come into play and the interfacial morphology do not follow cosine functions anymore. To explain the interfacial morphology shown in Fig. 12(a), we plotted the electrical force along the interface with a large amplitude, where the interface is assumed as y⁄(x⁄) = 2 + 0.25cos (2px⁄). Although this is not the real interface profile, the electrical force along this interface could provide us some clues. The results are depicted in Fig. 12(b) and (c), corresponding to the vertical component and horizontal component respectively. Fig. 12(b) illustrates that the force generated by the vertical electric field is negative (i.e., pointing downwards) and assumes a larger value at

Fig. 13. The interfacial morphology of the leaky electrohydrodynamic RTI and electrical force with ke = 1/3 and Eb = 0.2 in a horizontal field: (a) Snapshots of the interface deformation at dimensionless time t⁄ = 2.5, and he conductivity ratio kr, from left to right, are 1/6, 1/4, 1/3, 1/2 and 1. The corresponding structure heights Dh are 0.48, 0.90, 0.92, 1.10, and 1.20. (b) The vertical component of electrical force along the interface y⁄(x⁄) = 2 + 0.25cos(2px⁄); and (c) The horizontal component of electrical force along the interface.

702

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

Fig. 14. The interfacial morphology in vertical electric field. The electrical permittivity ratio is ke = 3 and electrical strength Eb = 0.05. (a) Snapshots of interface deformation at dimensionless time t⁄ = 2.5. The conductivity ratio kr, from left to right, are 0.5, 1, 3, 5 and 10. The corresponding structure heights Dh are 1.31, 1.22, 1.49, 1.67, and 1.98. (b) The vertical electrical force along the interface, where the interface is assumed as y⁄(x⁄) = 2 + 0.025cos(2px⁄).

Fig. 15. The interfacial morphology after interchanging the electrical permittivity, the electrical permittivity ratio is ke = 1/3 and electrical strength Eb = 0.05. (a) The snapshots of interfacial morphology for different conductivity ratios. The conductivity ratio kr, from left to right, are 1/10, 1/5, 1/3, 1 and 2. The corresponding structure heights Dh are 1.66, 1.53, 1.43, 1.52, and 1.97. (b) The electrical force in vertical direction along the interface y⁄(x⁄) = 2 + 0.025cos(2px⁄).

the rising tip. This force enforces the drag at the rising tip and therefore slows down its upward motion. At the same time, the falling ‘‘spike” is stretched long and thin to satisfy the constraint of mass conservation. The electrical force in the horizontal direction also affects the interfacial morphology, when kr = 1 the horizontal electrical force tends to squeeze the ‘‘spike”. However if kr = 6, the ‘‘spike” would be stretched in width. Thus with the increase of conductivity ratio kr, the ‘‘spike” becomes wider. The electrical permittivity affects the interfacial morphology as well. To better understand this effect, numerical simulations were carried with the electrical permittivities of the two fluids interchanged, that is, ke = 1/3. The snapshots of interfacial morphology at dimensionless time t⁄=2.5 are captured and depicted in Fig. 13 (a) for different values of kr (see Supplemental Movie 2 for the dynamic evolution). At kr = 1/6, the falling ‘‘spike” is short and thin. With an increase of kr, the ‘‘spike” becomes long (i.e., large Dh) and fat (see kr = 1). The vertical component which determines the motion of interface in the vertical direction is plotted in Fig. 13 (b) and the horizontal component is shown in Fig. 13(c). The electrical force is positive (i.e., pointing upwards) and attains a maximum at the falling tip. This suggests that during the evolution of interface deformation, the electrical force adds on to the drag at the falling tip, thereby slowing down the falling speed. The rising tip develops relatively fast than the falling tip, and hence the requirement of mass conservation causes the formation of a fat ‘‘spike”. The role of the horizontal electrical force seems to be com-

plicated. As shown in Fig. 13(c), the horizontal electrical force tends to squeeze the ‘‘spike” for small kr (=1/6, 1/4, 1/3 and 1/2) and helps to form a thin ‘‘falling” spike (see also the snapshots of kr = 1/6 and 1/4). For a large kr = 1, the horizontal force tends to stretch the ‘‘spike”, causing it to flatten. 5.3.2. Interfacial morphology in vertical electric field In a vertical electric field, the conductivity ratio kr can have a significant effect on the development of interfacial morphology. This is shown in Fig. 14(a), which plots the snapshots of interfacial morphology with different conductivity ratio at dimensionless time t⁄ = 2.5 (see also Supplemental Movie 3 for the dynamic process). At kr = 0.5, the falling ‘‘spike” is fat and short. With increasing kr, the ‘‘spike” becomes thin and long. The effect of conductivity ratio on the structure height (i.e., the length of the ‘‘spike”) has been explained in a previous section (Section 5.2.2). A detailed explanation of the interfacial morphology is provided here. Towards this end, the electrical force along the interface y⁄(x⁄) = 2 + 0.025cos(2px⁄) is plotted in Fig. 14(b). Clearly, for kr = 0.5, the electrical force points upwards (i.e., positive) with the maximum at the rising tip. Thus the force counters the action of the drag at the rising tip. This interaction speeds up the rising tip and renders it long and thin, while meantime broadening the falling ‘‘spike”. For the case of kr = 10, the electrical force points downwards (i.e., negative) and the largest magnitude occurs at the falling tip. This electrical force counterbalances the drag at

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

the falling tip. Thus the tip falls faster, leading to a long and thin ‘‘spike”. For the other cases (kr = 1, 3 and 5), the interfacial morphology are in between the case of kr = 0.5 and that of kr = 10. Further numerical results on vertical field effect are plotted in Fig. 15 for ke = 1/3, that is, the electrical permittivities of the two fluids are interchanged (see Supplemental Movie 4 for the dynamic process). The snapshots of interface deformation shows qualitatively the similar behavior as with the case of ke = 3; specifically, the ‘‘spike” is wide for a small kr and thin for a large kr. This is true for the pattern of the electrical force distribution also, though there are some differences in detail. As a result, the development of the interface morphology shown in Figs. 14 and 15 is dictated by the same mechanism.

703

Science Foundation (2015M570826), the Fundamental Research Funds for the Central Universities (xjj2016074), the National Engineering Laboratory for Highway Maintenance Equipment (Chang’an University) (310825161103) and the Foundation of Shaanxi Postdoctoral Science.

Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.02.049.

References 6. Conclusions A phase field numerical model has been developed to analyze the RTI of leaky dielectric fluids in an electric field. The model has been validated by the available data and linear analysis. Extensive numerical simulations were carried out. Results indicate that the effect of electric field on the RTI of leak dielectric fluids is more complicated than of pure dielectric fluids; on the electrical permittivity ratio ke, conductivity ratio kr and the direction of electric field. Also, the conductivity of the leak dielectric fluids, howbeit small, has a significant effect on the RTI. Main findings may be summarized below. (1) For a horizontal electric field with the permittivity ratio greater than 1 (ke > 1), the applied field suppresses the instability when the conductivity ratio larger than a critical value kr > kc and this suppressing effect becomes stronger with increasing conductivity ratio. However, it may aggravate the RTI when kr < kc. The falling ‘‘spike” is long and thin when conductivity ratio is small, while short and broad for a large conductivity ratio. (2) For a horizontal electric field with ke < 1, the electric field would suppress the RTI if the conductivity ratio smaller than a critical value kr < kc, and the suppressing effect is weakened as the conductivity ratio increases. The applied field starts to enhance the RTI when kr > kc. The falling ‘‘spike” is short and thin for a small conductivity ratio, and the ‘‘spike” becomes fat and long with increasing the conductivity ratio. (3) The effect of a vertical electric field is a non-monotonous function of the conductivity ratio. The length of the ‘‘spike” (i.e. structure height) decreases first and then increases with the conductivity ratio. A stationary point exists which corresponds to the minimum of the structure height. For different electrical permittivity ratios, the stationary points occur at different conductivity ratios. If ke > 1, two critical conductivity ratio kc1 and kc2 are found at which the electric field has no effect on the RTI. The electric field suppresses the instability when kc1 < kr < kc2 but aggravates it when kr < kc1 or kc2 < kr. On the other hand, if ke < 1, a critical conductivity ratio cannot be found. As for the interfacial morphology in a vertical electric field, the falling ‘‘spike” is broad when the conductivity ratio is small and becomes long and thin for a large conductivity ratio.

Acknowledgments Financial support of this work by the following sources is gratefully acknowledged: the National Natural Science Foundation of China (51605377, 11522219, 11372243), the China Postdoctoral

[1] D.H. Sharp, An overview of Rayleigh-Taylor instability, Physica D 12 (1–3) (1984) 3–18. [2] H.J. Kull, Theory of the Rayleigh-Taylor instability, Phys. Rep. 206 (5) (1991) 197–325. [3] G. Taylor, The Instability of liquid surfaces when accelerated in a direction perpendicular to their planes 1, Proc. Roy. Soc. Lon. Ser.-A 201 (1065) (1950) 192–196. [4] H. Michioka, I. Sumita, Rayleigh-Taylor instability of a particle packed viscous fluid: Implications for a solidifying magma, Geophys. Res. Lett. 32 (3) (2005) L03309. [5] X. Ribeyre, V.T. Tikhonchuk, S. Bouquet, Compressible Rayleigh-Taylor instabilities in supernova remnants, Phys. Fluids 16 (12) (2004) 4661–4670. [6] J.D. Kilkenny, S.G. Glendinning, S.W. Haan, B.A. Hammel, J.D. Lindl, D. Munro, B. A. Remington, S.V. Weber, J.P. Knauer, C.P. Verdon, A Review of the ablative stabilization of the Rayleigh-Taylor instability in regimes relevant to inertial confinement fusion, Phys. Plasmas 1 (5) (1994) 1379–1389. [7] L.L. Barannyk, D.T. Papageorgiou, P.G. Petropoulos, Suppression of RayleighTaylor instability using electric fields, Math. Comput. Simulat. 82 (6) (2012) 1008–1016. [8] V.M. Korovin, Effect of tangential electric field on the evolution of the RayleighTaylor instability of a dielectric liquid film, Tech. Phys.+ 56 (10) (2011) 1390– 1397. [9] R. Cimpeanu, D.T. Papageorgiou, P.G. Petropoulos, On the control and suppression of the Rayleigh-Taylor instability using electric fields, Phys. Fluids 26 (2) (2014) 022109. [10] A. Elmagd, A. Mohamed, E.S.F. Elshehawey, Non-linear electrohydrodynamic Rayleigh-Taylor instability 1. A perpendicular field in the absence of surfacecharges, J. Fluid Mech. 129 (Apr) (1983) 473–494. [11] A. Elmagd, A. Mohamed, E.F. Elshehawey, Non-linear electrohydrodynamic Rayleigh-Taylor instability. 2. A Perpendicular field producing surface-charge, Phys. Fluids 26 (7) (1983) 1724–1730. [12] J. White, J. Oakley, M. Anderson, R. Bonazza, Experimental measurements of the nonlinear Rayleigh-Taylor instability using a magnetorheological fluid, Phys. Rev. E 81 (2) (2010) 026303. [13] M.C. Renoult, R.G. Petschek, C. Rosenblatt, P. Carles, Deforming static fluid interfaces with magnetic fields: application to the Rayleigh-Taylor instability, Exp. Fluids 51 (4) (2011) 1073–1083. [14] J.M. Stone, T. Gardiner, The magnetic Rayleigh-Taylor instability in three dimensions, Astrophys. J. 671 (2) (2007) 1726–1735. [15] O. Porth, S.S. Komissarov, R. Keppens, Rayleigh-Taylor instability in magnetohydrodynamic simulations of the Crab nebula, Mon. Not. R. Astron. Soc. 443 (1) (2014) 547–558. [16] L. Rayleigh, Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density, Proceed. Lond. Math. Soc. 14 (1883) 170–177. [17] J.R. Melcher, W.J. Schwarz, Interfacial relaxation overstability in a tangential electric field, Phys. Fluids 11 (12) (1968) 2604–2616. [18] G. Taylor, Studies in electrohydrodynamics. I. Circulation produced in a drop by an electric field, Proc. Roy. Soc. Lon. Ser.-A 291 (1425) (1966) 159–166. [19] G.I. Taylor, A.D. Mcewan, Stability of a horizontal fluid interface in a vertical electric field, J. Fluid Mech. 22 (1) (1965) 1–15. [20] J.D. Jackson, Classical Electrodynamics, Wiley, 1999. [21] R.S. Allan, S.G. Mason, Particle behaviour in shear and electric fields. 1. Deformation and burst of fluid drops, Proc. Roy. Soc. Lon. Ser.-A 267 (1328) (1962) 45–61. [22] S.A. Roberts, S. Kumar, AC electrohydrodynamic instabilities in thin liquid films, J. Fluid Mech. 631 (2009) 255–279. [23] A. Corbett, S. Kumar, Spreading of thin droplets of perfect and leaky dielectric liquids on inclined surfaces, Langmuir 32 (26) (2016) 6606–6617. [24] N. Wu, W.B. Russel, Microand nano-patterns created via electrohydrodynamic instabilities, Nano Today 4 (2) (2009) 180–192. [25] A. Ramkrishnan, S. Kumar, Electrohydrodynamic effects in the leveling of coatings, Chem. Eng. Sci. 101 (2013) 785–799. [26] D. Jacqmin, Calculation of two-phase Navier-Stokes flows using phase-field modeling, J. Comput. Phys. 155 (1) (1999) 96–127.

704

Q. Yang et al. / International Journal of Heat and Mass Transfer 109 (2017) 690–704

[27] X.Y. He, S.Y. Chen, R.Y. Zhang, A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability, J. Comput. Phys. 152 (2) (1999) 642–663. [28] A. Celani, A. Mazzino, P. Muratore-Ginnaneschi, L. Vozella, Rayleigh-Taylor instability in two dimensions and phase-field method, Springer Proc. Phys. 132 (2009) 169–172. [29] H. Ding, P.D.M. Spelt, C. Shu, Diffuse interface model for incompressible twophase flows with large density ratios, J. Comput. Phys. 226 (2) (2007) 2078– 2095. [30] A. Celani, A. Mazzino, P. Muratore-Ginanneschi, L. Vozella, Phase-field model for the Rayleigh-Taylor instability of immiscible fluids, J. Fluid Mech. 622 (2009) 115–134. [31] Y. Lin, P. Skjetne, A. Carlson, A phase field model for multiphase electrohydrodynamic flow, Int. J. Multiphas Flow 45 (2012) 1–11. [32] Q. Yang, B.Q. Li, Y. Ding, 3D phase field modeling of electrohydrodynamic multiphase flows, Int. J. Multiphas Flow 57 (2013) 1–9. [33] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy, J. Chem. Phys. 28 (2) (1958) 258–267. [34] T. Qian, X.-P. Wang, P. Sheng, A variational approach to moving contact line hydrodynamics, J. Fluid Mech. 564 (2006) 333–360.

[35] J.S. Hua, L.K. Lim, C.H. Wang, Numerical simulation of deformation/motion of a drop suspended in viscous liquids under influence of steady electric fields, Phys. Fluids 20 (11) (2008) 113302. [36] D.A. Saville, Electrohydrodynamics: the Taylor-Melcher leaky dielectric model, Annu. Rev. Fluid Mech. 29 (1997) 27–64. [37] G. Tomar, D. Gerlach, G. Biswas, N. Alleborn, A. Sharma, F. Durst, S.W.J. Welch, A. Delgado, Two-phase electrohydrodynamic simulations using a volume-offluid approach, J. Comput. Phys. 227 (2) (2007) 1267–1285. [38] J.A. Lanauze, L.M. Walker, A.S. Khair, The influence of inertia and charge relaxation on electrohydrodynamic drop deformation, Phys. Fluids (1994present) 25 (11) (2013) 112101. [39] J.A. Lanauze, L.M. Walker, A.S. Khair, Nonlinear electrohydrodynamics of slightly deformed oblate drops, J. Fluid Mech. 774 (2015) 245–266. [40] P. Yue, C. Zhou, J.J. Feng, C.F. Ollivier-Gooch, H.H. Hu, Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing, J. Comput. Phys. 219 (1) (2006) 47–67. [41] B.J. Daly, Numerical Study of two Fluid Rayleigh-Taylor Instability, Phys. Fluids 10 (2) (1967) 297–307.