Numerical simulation of electrohydrodynamics of a compound drop based on the ternary phase field method

Abstract

Analytical and numerical methods are often used to study the behavior of multiphase fluid under electric field. Compared with analytical methods, numerical methods can simulate the real physical phenomenon of multiphase fluid dynamics in a large deformation range. The finite element method is mainly applied in two-phase fluid currently, although it can be used to analyze the small and large deformation of multiphase fluid under electric field. This article attempts to develop a finite element model of a concentric compound drop immersed in continuous medium under electric field based on the ternary phase field method and simulate the electrohydrodynamics of the compound drop whose core phase, shell phase, and continuous phase are different. The small deformation simulation results of the compound drop under weak electric field are compared with the analytical results of previous researchers from the three aspects, namely, deformation, free charge distribution, and flow pattern. This model is proved to be effective under certain conditions. Based on this premise, the large deformation and breakup of the compound drop under high electric field are further simulated to investigate the mechanism of compound drop breakup preliminarily.

Introduction

The study of drop behavior in electric field plays an important role for the primary theory of electrohydrodynamics (EHD). In recent years, the application of EHD in industrial production such as electrostatic spraying, electrostatic spinning, ink jet printing, oil–water separation, and electrowetting is increasing, which urgently needs the theoretical research to predict and optimize the fluid behavior in electric field. The simplest way is to study the deformation of drop under two parallel electrode plates. Based on the well-known leaky dielectric model of Taylor, ¹ the steady-state EHD of a single drop in weak electric field has been well understood. However, a numerical method has to be used to study the drop behavior in high electric field because the drop fails to reach a balance.

Many numerical methods are used to simulate the EHD of single droplet, such as boundary integral method (BIM), volume of fluid (VOF), level set (LS), and lattice Boltzmann method (LBM). Sherwood ² studied the large deformation and breakup of single droplet subjected to electric field using BIM and found that pointed ends were predicted when the permittivity of droplet was greater than that of the surrounding fluid; however, the droplet experienced division into blobs when the conductivity of the droplet was higher compared with that of the surrounding fluid. Paknemat et al. ³ using the LS method with the ghost fluid method investigated the deformations and breakup modes of three types of droplets in an electric field, and the simulation results were compared with similar numerical and experimental results. The VOF method was utilized by Tomar et al. ⁴ to predict the droplet deformation in an electric field and the electrocoalescence behavior of two droplets. Using the coupled level set and volume of fluid (CLSVOF) method, Huang et al. ⁵ simulated the breakup behavior of the leaky dielectric droplet under direct current electric field. They found that the evolutions of velocity field and pressure field determined the breakup mode transformation. Lin et al. ⁶ used the phase field method to simulate the deformation of a single droplet for leaky and perfect dielectric fluids in an electric field, and the results were found to be in good agreement with the existing theoretical and numerical studies. Based on the three different electric field models, namely, leaky dielectric model, perfect dielectric model, and constant surface charge model, Hua et al. ⁷ used the front tracking method (FTM) to simulate the deformation of single droplet suspended in viscous liquids subjected to an external electric field, and the reasonable prediction was obtained. Singh et al. ⁸ presented a low spurious current lattice Boltzmann model and simulated the interaction of single or multiple droplets suspended in the outer fluid in the presence of an electric field by employing this model. By coupling multi-component LBM with the leaky dielectric model, they also conducted a numerical study of deformation and breakup of a droplet subjected to shear flow under an electric field. ⁹ The numerical studies of drop EHD outlined earlier concentrate on leaky dielectrics, perfect dielectrics, or perfect conductors. The electrical properties are considered to be uniform within each fluid phase in these cases. However, for electrolyte solutions with mobile charge carriers, gradients of ion concentration form in an electric field, causing non-uniform conductance and thereby affecting the operation of microfluidic devices significantly. Therefore, there is a need for a multiphase electrokinetic flow model that is capable of tracking single and multiple species of mobile charge carriers. Based on the CLSVOF technique, Berry et al. ¹⁰ and Pillai et al. ¹¹ developed the first numerical method to incorporate ion transport into a model of electrokinetic two-fluid flow and investigated the deformation and breakup of an electrolyte drop under the action of an electric field. Davidson et al. ¹² extended their model to allow for interfacial charge and simulated the deformation of a hexadecane drop in pressure-driven electrokinetic flow of electrolyte solution through a microchannel by employing this extended model along with a modified pressure formulation.

Compound drop is a fluid particle consisted of an inner drop encapsulated by another immiscible fluid which is itself immersed in a third fluid. It is encountered in a variety of technologically significant processes. Examples of these processes include liquid membranes for selective mass transport, water purification, controlled release of drugs, direct contact heat exchange, and recovery of heavy metals.^13,14 Application of electric field can improve the efficiency of the processes involving compound drops. Hence, the study on the behavior of the compound drop under the influence of electric field is of particular importance to deeply comprehend these processes and provide guidance to practical applications. But compared with single-phase drop, EHD of the compound drop is quite less explored. From the standpoint of experiment, the density of three fluids must be close in order to stabilize the compound drop because of buoyancy caused by the difference in density, which greatly limits the variety of fluid studied. It is only known that Tsukada et al. ¹⁵ performed a deformation experiment of a kind of special compound drop whose core phase and continuous phase are the same in a uniform electric field. Thus, analytical and finite element methods are often used. Ha and Yang ¹⁶ investigated the influence of uniform electric field on the fluid dynamics of a concentric compound drop analytically and successfully characterized the flow patterns by obtaining the approximate solutions of electric and velocity fields and examining the charge distribution at the inner and outer interfaces. They found that the second recirculating flow occurred in the shell phase when the charges at the core–shell and shell–continuous medium interfaces are the same property. Based on the leaky dielectric theory, Behjatian and Esmaeeli^13,14 studied the behavior of a concentric compound drop in a uniform electric field and analyzed the electric field and flow pattern in detail by the analytical method. The results showed that there were four possible flow patterns in the compound drop, and the dielectric constant and conductivity of fluid played a crucial role in determining the steady-state evolution of flow field. However, since the studies of Ha and Yang ¹⁶ and Behjatian and Esmaeeli^13,14 are based on the domain perturbation procedure, they can only accurately predict small deformation of the compound drop in an electric field. Compared with analytical studies, numerical studies can simulate the real physical phenomenon of compound drop dynamics in large deformation range and predict the compound drop breakup accurately. Several numerical studies on a compound drop in electric field have been conducted so far. Soni et al. ¹⁷ used the phase field method to simulate the dynamics of the compound drop in a uniform electric field. The results showed that the simulation results were in good agreement with the experiment results of Tsukada et al. ¹⁵ However, since the core phase of the compound drop studied was the same as the continuous phase, the phase field method employed was still a two-phase flow method. Abbasi et al. ¹⁸ used the LS method to simulate the instability behavior of the compound drop under high electric field. The results suggested that when the core drop was eccentric, the electric field would induce it to be concentric with the shell drop or be close to the shell–continuous phase interface depending on electrical parameters. But the compound drop studied was still a two-phase fluid system in which the core phase was the same as the continuous phase. Moreover, the LS method cannot be extended to three-phase flow. In other words, the phase field and LS methods employed in the aforementioned two studies cannot be employed to simulate the EHD of the compound drop whose core phase, shell phase, and continuous phase are different. To our knowledge, the EHD of the compound drop whose core phase, shell phase, and continuous phase are different has not yet been studied by numerical simulation.

This work mainly includes four parts. First, a theoretical model is proposed for the coupling of electric field and three-phase flow field. Second, a finite element model of a concentric compound drop immersed in continuous medium under electric field is developed based on the ternary phase field method. Then, its feasibility and applicable range are confirmed by comparing the simulation results of small deformation of the compound drop with the analytical results of previous researchers in three aspects, namely, deformation, charge distribution, and flow pattern. Finally, the large deformation and breakup of the compound drop under high electric field are simulated under different electrical parameters to explore the mechanism of compound drop breakup preliminarily.

Theory

Ternary phase field theory

The ternary phase field model proposed by Boyer et al. ¹⁹ is used to track the interfaces among three immiscible fluids. This model is solved based on the following Cahn–Hilliard equations

$$\frac{\partial {\varnothing }_i}{\partial t}+\nabla ·\left(\mathit{u}{\varnothing }_i\right)=\nabla ·(\frac{M_0}{{\xi }_i}\nabla {\eta }_i)$$

(1)

$${\eta }_i=\frac{4{\xi }_{\Gamma }}{\epsilon }\sum _{j\ne i}\left(\frac{1}{{\xi }_j}(\frac{\partial F}{\partial {\varnothing }_i}-\frac{\partial F}{\partial {\varnothing }_j})\right)-\frac{3}{4}\epsilon {\xi }_i{\nabla }^2{\varnothing }_i$$

(2)

$$\sum _{i=A,B,C}{\varnothing }_i=1$$

(3)

$$\rho =\sum _{i=A,B,C}{\rho }_i{\varnothing }_i$$

(4)

$$\mu =\sum _{i=A,B,C}{\mu }_i{\varnothing }_i$$

(5)

where ${\varnothing }_i$ is the phase field variable and represents the volume concentration of phase i, varying from 0 to 1; u is the fluid velocity vector; $M_0$ is the migration; ${\xi }_i$ is the capillary parameter; ${\eta }_i$ is the generalized chemical potential; ${\rho }_i$ is the density; ${\mu }_i$ is the viscosity, $i=A,B,C$ ; F is the free energy of a three-phase system; and $\epsilon$ is the interface thickness. As pointed out by Boyer et al., ¹⁹ one of the phase variables, ${\varnothing }_A$ , ${\varnothing }_B$ , and ${\varnothing }_C$ , will be eliminated a posteriori and the solution must not depend on the choice of the eliminated unknown (equation (3)). Thus, only two Cahn–Hilliard equations should be solved.

$M_0$ can be estimated based on the following empirical formula provided in the help documentation of COMSOL software

$$M_0\approx L_c{U}_c\epsilon \times {10}^{-2}$$

(6)

where $L_c$ is the characteristic length and $U_c$ is the characteristic velocity. Details on the values of $\epsilon$ and $M_0$ are given in section “Parameter set.”

The parameter ${\xi }_{\Gamma }$ in equation (2) is defined as follows

$$\frac{3}{{\xi }_{\Gamma }}=\frac{1}{{\xi }_A}+\frac{1}{{\xi }_B}+\frac{1}{{\xi }_C}$$

(7)

${\xi }_i(i=A,B,C)$ is further defined as

$${\xi }_A={\sigma }_{\mathrm{AB}}+{\sigma }_{\mathrm{AC}}-{\sigma }_{\mathrm{BC}}$$

(8)

$${\xi }_B={\sigma }_{\mathrm{AB}}+{\sigma }_{\mathrm{BC}}-{\sigma }_{\mathrm{AC}}$$

(9)

$${\xi }_C={\sigma }_{\mathrm{BC}}+{\sigma }_{\mathrm{AC}}-{\sigma }_{\mathrm{AB}}$$

(10)

where ${\sigma }_{\mathrm{ij}}$ is the initial surface tension coefficient of phase i and phase j.

F is defined as a function of phase field variable

$$\begin{array}{c}\multicolumn{1}{c}{F={\sigma }_{\mathrm{AB}}{\varnothing }_A^2{\varnothing }_B^2+{\sigma }_{\mathrm{AC}}{\varnothing }_A^2{\varnothing }_C^2+{\sigma }_{\mathrm{BC}}{\varnothing }_B^2{\varnothing }_C^2}\\ +{\varnothing }_A{\varnothing }_B{\varnothing }_C({\xi }_A{\varnothing }_A+{\xi }_B{\varnothing }_B+{\xi }_C{\varnothing }_C)+\mathrm{\Lambda }{\varnothing }_A^2{\varnothing }_B^2{\varnothing }_C^2\hfill \end{array}$$

(11)

where $\mathrm{\Lambda }$ is the free energy outside the system. If there is no special explanation, $\mathrm{\Lambda }$ is 0.

Electric field theory

For leaky dielectric fluids, both free and polarized charges should be considered. So, it can be started with Gauss law

$$Q_f=\nabla ·\mathit{D}$$

(12)

$$\mathit{D}={\epsilon }_0{\epsilon }_r\mathit{E}$$

(13)

where $Q_f$ is the free charge density, D is the electric displacement vector, E is the electric field vector, ${\epsilon }_0$ is the electric permittivity of vacuum, and ${\epsilon }_r$ is the relative dielectric constant.

The electric field can be assumed to be irrotational $(\nabla \times \mathit{E}=0)$ . Therefore, the electric field can be expressed as a gradient function of electric potential $V(\mathit{E}=-\nabla V)$ . The following equation can be achieved by combining the aforementioned formulas

$$Q_f=-\nabla ·({\epsilon }_0{\epsilon }_r\nabla V)$$

(14)

Free charge density $Q_f$ can be obtained by the charge conservation equation

$$\frac{\partial Q_f}{\partial t}=-\mathit{u}·\nabla Q_f-\nabla ·(D_0\nabla Q_f)+\nabla ·(\kappa \nabla V)$$

(15)

where $D_0$ and $\kappa$ represent the charge diffusion coefficient and the fluid conductivity, respectively. The three terms on the right-hand side of equation (15) represent the convection, diffusion, and migration of charge, respectively. The charge diffusion can be safely ignored. But if equation (15) is further simplified, the time scale of flow $\tau \approx L_c/2U_c$ and the charge relaxation time ${\tau }_c={\epsilon }_0{\epsilon }_r/\kappa$ must be considered. There are two cases.

• If each fluid phase satisfies $\tau >>{\tau }_c$ , it means that compared with the flow of fluid, charge transfers to the interface almost instantaneously. Then, charge convection is neglected, and the electric field governing equation is as follows

  $$\nabla ·(\kappa \nabla V)=0$$

  (16)

• If Case a cannot be satisfied, the electric field governing equation is the general form of equation (15) and given as follows

  $$\frac{\partial Q_f}{\partial t}=-\mathit{u}·\nabla Q_f+\nabla ·(\kappa \nabla V)$$

  (17)

Equation (17) becomes the following form by combining equation (14)

$$\nabla ·[{\epsilon }_0{\epsilon }_r\frac{\partial \nabla V}{\partial t}+\kappa \nabla V+\mathit{u}\nabla ·({\epsilon }_0{\epsilon }_r\nabla V)]=0$$

(18)

Finally, for simplicity, the functions in parentheses of equation (16) and square brackets of equation (18) are referred to as J , which represents the current density. Similar to the density and viscosity, $\kappa$ and ${\epsilon }_r$ are also defined as a function of phase field variable

$$\kappa =\sum _{i=A,B,C}{\kappa }_i{\varnothing }_i$$

(19)

$${\epsilon }_r=\sum _{i=A,B,C}{\epsilon }_i{\varnothing }_i$$

(20)

where ${\kappa }_i$ and ${\epsilon }_i(i=A,B,C)$ represent the conductivity and relative dielectric constant of each phase, respectively.

The coupling

Drop deformation generates when electric stress overcomes surface tension and viscous stress. It is assumed that the fluid is incompressible and the viscosity is irrotational. So, the governing equations are as follows

$$\rho \frac{\partial \mathit{u}}{\partial t}+\rho (\mathit{u}·\nabla )\mathit{u}=\nabla ·[-\mathrm{pI}+\mu (\nabla \mathit{u}+{(\nabla \mathit{u})}^{\mathrm{T}})]+{\mathit{F}}_{\mathrm{st}}+{\mathit{F}}_e$$

(21)

$$\rho \nabla ·\left(\mathit{u}\right)=0$$

(22)

where p and I are the pressure and unit matrix, respectively, ${\mathit{F}}_{\mathrm{st}}$ is the surface tension force, and ${\mathit{F}}_e$ is the electric stress per unit volume.

According to the ternary phase field method, the surface tension force ${\mathit{F}}_{\mathrm{st}}$ can be represented in the form of body force by the generalized chemical potential

$${\mathit{F}}_{\mathrm{st}}=\sum _{i=A,B,C}{\eta }_i\nabla {\varnothing }_i$$

(23)

The electric stress per unit volume ${\mathit{F}}_e$ can be obtained by the divergence of Maxwell stress tensor $\overleftrightarrow{\mathit{T}}$

$$\overleftrightarrow{\mathit{T}}=\mathit{E}{\mathit{D}}^T-\frac{1}{2}(\mathit{E}·\mathit{D})\mathit{I}$$

(24)

$${\mathit{F}}_e=\nabla ·\overleftrightarrow{\mathit{T}}=\mathit{E}{Q}_f-\frac{1}{2}{\epsilon }_0{\mathit{E}}^2\nabla {\epsilon }_r$$

(25)

Review of analytical research on a compound drop in electric field

The key factors for EHD of a compound drop are the conductivity ratio $({\mathrm{\kappa }}_{\mathrm{AB}}={\kappa }_A/{\kappa }_B,{\mathrm{\kappa }}_{\mathrm{BC}}={\kappa }_B/{\kappa }_C)$ and the dielectric constant ratio $({\epsilon }_{\mathrm{AB}}={\epsilon }_A/{\epsilon }_B,{\epsilon }_{\mathrm{BC}}={\epsilon }_B/{\epsilon }_C)$ of inner and outer fluids. According to the studies of Ha and Yang ¹⁶ and Behjatian and Esmaeeli, ¹³ there are four relationships among charge distribution, conductivity ratio, and dielectric constant ratio, as shown in Figure 1. When the conductivity ratios of inner and outer fluids are both more or less than the dielectric constant ratios of those, the charge property at the core–shell interface is the same as that at the shell–continuous medium interface (Figure 1(a) and (c)). Otherwise, the charge property at the core–shell interface is contrary to that at the shell–continuous medium interface (Figure 1(b) and (d)).

Figure 1.

Charge distribution of compound drop: (a) ${\kappa }_{\mathrm{AB}}<{\epsilon }_{\mathrm{AB}}{\kappa }_{\mathrm{BC}}<{\epsilon }_{\mathrm{BC}}$ , (b) ${\kappa }_{\mathrm{AB}}>{\epsilon }_{\mathrm{AB}}{\kappa }_{\mathrm{BC}}<{\epsilon }_{\mathrm{BC}}$ , (c) ${\kappa }_{\mathrm{AB}}>{\epsilon }_{\mathrm{AB}}{\kappa }_{\mathrm{BC}}>{\epsilon }_{\mathrm{BC}}$ , and (d) ${\kappa }_{\mathrm{AB}}<{\epsilon }_{\mathrm{AB}}{\kappa }_{\mathrm{BC}}>{\epsilon }_{\mathrm{BC}}$ .

The weak electric field can induce a circulation flow in the compound drop, and the flow pattern also has an important relationship with the electrical properties of fluids. According to the research work of Ha and Yang ¹⁶ and Behjatian and Esmaeeli, ¹³ flow patterns are mainly divided into four cases, as shown in Figure 2. The charge distributions in Figure 1(b) and (d) lead to the flow patterns of Cases a and b in Figure 2, respectively. However, when the charge property at the core–shell interface is the same as that at the shell–continuous medium interface, the flow pattern depends on the magnitude of electric stresses at these two interfaces. If the electric stresses are comparable, the charge distributions in Figure 1(a) and (c) result in the flow patterns of Cases c and d in Figure 2, respectively. If the difference in electric stresses at these two interfaces is much, the flow patterns of Cases a and b in Figure 2 may result from the charge distributions in Figure 1(a) and (c), respectively.

Figure 2.

Flow patterns of steady-state electrohydrodynamics of compound drop: (a) Case a, (b) Case b, (c) Case c, and (d) Case d.

Methodology

Numerical method and geometry

COMSOL Multi-physics 5.3 which is a commercial package based on the finite element method is used to simulate the EHD of compound drop. Since weak form is an integral form and has a low requirement for the continuity of field variable, it is very suitable for nonlinear multi-physics problems. Thus, the Navier–Stokes equations (equations (21) and (22)) and the governing equations of the ternary phase field (equations (1)–(10)) are first transformed into weak form and then solved.

It is assumed that a concentric compound drop is immersed in the infinite continuous medium. Thus, more solution efficiency can be applied to the interest domain as small as possible. As shown in Figure 3, the cylindrical coordinate is used. The compound drop consists of core drop with a radius of $R_0$ and shell drop with a radius of $2R_0$ . Since the charge obeys the anti-symmetric constraint on the other side of OI in a uniform electric field, the solution domain can be further simplified to $6R_0\times 6R_0$ and only the upper half of the MNGH plane, namely IOGH plane, is solved. ²⁰ To simulate an infinite solution domain, the open boundary ( $[-\mathrm{pI}+\mu (\nabla \mathit{u}+(\nabla \mathit{u}{)}^{\mathrm{T}})]n=0$ , where n is the normal gradient operator) is applied on GH and HI. The symmetry condition $(\mathit{u}·n=0\mathrm{\&}K-(K·n)n=0,K=[\mu (\nabla \mathit{u}+(\nabla \mathit{u}{)}^{\mathrm{T}})]n)$ is used on OI. In order to obtain a uniform electric field, it is necessary to add virtual electrodes on GH and OI. Furthermore, the electric potential $V=V_0=6E_0{R}_0$ and $V=0$ are specified on GH and OI, respectively. The boundary HI is electrical insulation $(n·\mathit{J}=0)$ to achieve the definite solution of electric potential.

Figure 3.

Computational domain and boundary condition.

As shown in Figure 4, the uniform triangular element is used in meshing most regions, and the grids are encrypted at the initial interface of the compound drop with a scaling of 0.5. This mesh is applied in all the simulations. The adaptive mesh refinement is performed in order to simulate the EHD of the compound drop more accurately. The quadratic and linearity are used in the discretization of electric field governing equations and flow field governing equations, respectively, so as to ensure the solution accuracy and improve the convergence.

Figure 4.

Mesh used in simulations.

Parameter set

The two key parameters in the ternary phase field method, namely interface thickness $\epsilon$ and migration $M_0$ , need to be clearly specified. Good convergence can be achieved when the interface thickness $\epsilon$ takes the recommended value of COMSOL software, namely $\epsilon =\sqrt{2}{h}_{\max }$ , where $h_{\max }$ is the maximum grid length in computational domain. The flow will be inaccurate if the characteristic velocity $U_c$ in equation (6) is not suitable. So, the migration $M_0$ needs to be determined cautiously. It is found that $U_c={\epsilon }_C{E}_0^2{L}_c/{\mu }_C$ and $L_c=2R_0$ may be appropriate in our simulations. It is proposed that $M_0$ is constrained at the interface for better accuracy. Here, $M_{\mathrm{help}}$ is defined as a function of phase field variable

$$M_{\mathrm{help}}={\varnothing }_A^2{\varnothing }_B^2+{\varnothing }_B^2{\varnothing }_C^2+{\varnothing }_A^2{\varnothing }_C^2$$

(26)

Due to the fact that $\varnothing =0.5$ represents the interface, the initial value of $M_{\mathrm{help}}$ at the interface will become $(1/2^2)\times (1/2^2)=(1/16)$ . Therefore, the migration $M_0$ may be calculated as follows

$$M_0=M_{\mathrm{help}}{L}_c{U}_c\epsilon \times {10}^{-2}\times 16$$

(27)

Simulation conditions

Six leaky dielectric systems from the study by Behjatian and Esmaeeli ¹³ are used in the small deformation simulations of the compound drop under weak electric field. The key physical properties of six leaky dielectric systems are extracted from Behjatian and Esmaeeli ¹³ and given in Table 1.

Table 1.

Properties of leaky dielectric systems extracted from the study of Behjatian and Esmaeeli. ¹³

System	Phase A–B–C	${\kappa }_{\mathrm{AB}}$	${\kappa }_{\mathrm{BC}}$	${\epsilon }_{\mathrm{AB}}$	${\epsilon }_{\mathrm{BC}}$
1	M2–T–B1	0.073	0.0891	0.46	0.72
2	N–M3–N	30.03	0.033	2.17	0.46
3	B1–A1–A2	481.13	3.97	2.6	1.22
4	M3–N–M3	0.033	30.3	0.46	2.17
5	B2–A1–M4	0.5	0.32	0.86	1.17
6	M2–A1–B2	3.14	2.0	0.85	1.16

For comparison, the definition of time scale of flow $\tau$ is the same as that of Behjatian and Esmaeeli, ¹⁴ as listed as follows

$$\tau =t_s\frac{(19\tilde{\mu }+16)(2\tilde{\mu }+3)}{40(\tilde{\mu }+1)}$$

(28)

where $t_s={\mu }_{C}2R_0/{\sigma }_{\mathrm{BC}}$ , $\tilde{\mu }={\mu }_B/{\mu }_C$ or $t_s={\mu }_B{R}_0/{\sigma }_{\mathrm{AB}}$ , $\tilde{\mu }={\mu }_A/{\mu }_B$ . According to Behjatian and Esmaeeli ¹³ and equation (28), the ratio of conductivity ratio to dielectric constant ratio, charge relaxation time, and time scale of flow for six leaky dielectric systems are calculated and provided in Table 2. It can be seen from Table 2 that the charge relaxation time ${\tau }_c$ of at least one phase is not much less than the time scale of flow $\tau$ in the three phases of the compound drop. Hence, the electric field governing equation is equation (18).

Table 2.

Other data of systems 1–6 in Table 1 obtained by calculation.

System	${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}$	${\kappa }_{\mathrm{BC}}/{\epsilon }_{\mathrm{BC}}$	${\tau }_{\mathrm{c}}$ (s)	$\tau$ (s)
1	0.16 < 1	0.12 < 1	0.015–0.74	5.52
2	13.96 > 1	0.07 < 1	0.053–0.73	4.68
3	222.75 > 1	3.25 > 1	0.014–8.82	0.37
4	0.07 < 1	13.96 > 1	0.053–0.73	4.68
5	0.58 < 1	0.27 < 1	0.74–4.67	20
6	3.69 > 1	1.72 > 1	0.73–4.67	5.80

The small deformation simulation results are compared with the analytical results of Behjatian and Esmaeeli^13,14 in terms of deformation, free charge distribution, and flow pattern to evaluate the feasibility and applicable range of the proposed model. After that, the large deformation and breakup of the compound drop are simulated under high electric field to show the topological advantage of the proposed model. In these simulations, the magnitude of conductivity and relative dielectric constant of three fluids is about 10⁻³ S/m and 10, respectively. Thus, the charge relaxation time ${\tau }_c$ and the time scale of flow $\tau$ are about 10⁻⁷ s and 10⁻³ s, respectively. So the electric field governing equation used is equation (16) due to the fact that ${\tau }_c$ is much less than $\tau$ .

Results and discussion

Small deformation of the compound drop under weak electric field

Under small deformation condition, a dimensionless parameter D is usually used to assess the drop deformation. It is defined as

$$D=\frac{b-a}{b+a}$$

(29)

where b and a are the length of drop axis parallel and perpendicular to the electric field, respectively. Due to the presence of core and shell for the compound drop, core drop deformation and shell drop deformation are defined as $D_{12}$ and $D_{23}$ , respectively.

Figure 5 shows the variation in core drop deformation $D_{12}$ and shell drop deformation $D_{23}$ with time scale $t/\tau$ for systems 1–6 and their comparison with the analytical results of Behjatian and Esmaeeli.^13,14 It is pointed out that in the analytical results of Behjatian and Esmaeeli,^13,14 half of the six systems (Figure 5(b), (d), and (f)) contain the time-dependent deformation history, while the other half (Figure 5(a), (c), and (e)) show only the final steady-state deformation. It can be seen from Figure 5 that for system 3 (Figure 5(c)), the core drop deformation $D_{12}$ and shell drop deformation $D_{23}$ deviate obviously from the analytical values of Behjatian and Esmaeeli, ¹³ especially the core drop deformation $D_{12}$ . The deformation of compound drop for the other five systems is in good agreement with the analytical results of Behjatian and Esmaeeli.^13,14 The free charge distribution is further examined in order to find out the reason why a big error generates in the deformation simulation of system 3.

Figure 5.

Comparison of core drop deformation $D_{12}$ and shell drop deformation $D_{23}$ with the analytical results of Behjatian and Esmaeeli for systems 1–6: (a) system 1, (b) system 2, (c) system 3, (d) system 4, (e) system 5, and (f) system 6.

Figure 6 shows the distribution of free charge and electric stress for systems 1–6. It can be found that the distribution of positive and negative charges for systems 1–6 is in accord with the analytical results of Behjatian and Esmaeeli. ¹³ As can be seen in Figure 6, except system 3, the charges for the other five systems distribute very closely to the interfaces. The free charges in system 3 (Figure 6(c)) enter into the shell phase, which results in the relatively large deviation in the deformation of the compound drop from the analytical results of Behjatian and Esmaeeli (Figure 5(c)). It can be seen from Table 1 that the conductivity ratio of the core phase to the shell phase for system 3 $({\kappa }_{\mathrm{AB}}=481.13)$ is far more than that for the other five systems. This may be the reason why the free charge distribution for system 3 is different from that for the other five systems. Here, the influence of dielectric constant ratio is excluded because the dielectric constant only affects the magnitude of free charge, not making the free charge deviate from the interface.

Figure 6.

The distribution of free charge and electric stress for systems 1–6: (a) system 1, (b) system 2, (c) system 3, (d) system 4, (e) system 5, and (f) system 6.

By observing Figures 5 and 6, it can be found that although the charge distribution for system 1 is similar to that for system 5, the tendency of deformation of the compound drop is different. The same phenomenon also exists in systems 3 and 6. Thus, it is necessary to further examine the flow pattern. Figure 7 shows the flow patterns of systems 1–6. It can be found that the flow patterns of systems 2, 4, 5, and 6 are exactly consistent with the analytical results of Behjatian and Esmaeeli, ¹³ but there is a deviation for those of systems 1 and 3. The flow direction obtained by Behjatian and Esmaeeli ¹³ is marked using red rotating arrow in systems 1 and 3 (Figure 7(a) and (c)). The abnormal distribution of free charge (Figure 6(c)) causes a deviation in the flow pattern of system 3. However, the deviation in the flow pattern of system 1 surprised us. In the research work of Behjatian and Esmaeeli, ¹³ the flow pattern of system 1 is Case a in Figure 2, which generates when the difference in electric stresses at the core–shell and shell–continuous medium interfaces is very big. But in our simulation, there is not much difference in the electric stresses at these two interfaces (Figure 6(a)), consequently bringing about a deviation of flow pattern of system 1. As can be seen in Figure 7(a), there is a large velocity in the core drop close to the symmetry axis, which is roughly equal to the velocity at the shell–continuous medium interface. This velocity drops by two orders of magnitude at the core–shell interface. It is believed that the flow field outside the core–shell interface is not sufficient to overcome the internal flow field of core drop and change the flow direction in core drop. The secondary flow circulation is infinitely small at the core–shell interface, making the flow at the interface relatively stationary (Figure 7(a)). Therefore, this phenomenon can be understood as the transition from Case c to Case a in Figure 2.

Figure 7.

Velocity stream lines and direction arrows for systems 1–6: (a) system 1, (b) system 2, (c) system 3, (d) system 4, (e) system 5, and (f) system 6.

In conclusion, the conductivity ratio of two adjacent fluids for the compound drop is found to be the only parameter that significantly restricts the accuracy of the proposed model. Figure 8 shows the relationship between conductivity ratio of two adjacent fluids ${\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$ (in for inner fluid and out for outer fluid) and ratio of D to D_BE for systems 1–6. D and D_BE in Figure 8 represent the deformations of the compound drop obtained by the authors and Behjatian and Esmaeeli,^13,14 respectively. Thus, the ratio of D to D_BE reflects the accuracy of the proposed model. It can be seen from Figure 8 that when ${\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$ ranges from 0.1 to 10, good accuracy can be achieved. However, the accuracy sharply declines as ${\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$ is more than 100. Table 3 gives the conductivity ratio range of leaky dielectric fluids used in prior simulations with different interface tracking methods. As can be seen in Table 3, regardless of the interface tracking method used, the conductivity ratio of two adjacent fluids ${\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$ is generally not more than 100 or less than 0.01. Furthermore, its range is usually 0.1–10 when the simulation results are in good agreement with the analytical results. Therefore, to ensure the accuracy, the conductivity ratio of two adjacent fluids ${\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$ for the compound drop is set between 0.1 and 10 in the following simulations.

Figure 8.

The relationship between ${\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$ and D/D_BE for systems 1–6.

Table 3.

Conductivity ratio of leaky dielectric fluids used in different ITMs.

	Authors
	Hua et al. 7	Paknemat et al. 3	Lin et al. 6	Lima et al. 21	Abbasi et al. 18
${\kappa }_r={\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$	0.01–70	0.1–10	0.1–10	2–30	0.1–100
ITM	FTM	GFM and LS	PF	VOF	LS
Agreement a	0.1–10	0.1–10	0.1–10	<10	0.1–10

ITM: interface tracking method; FTM: front tracking method; GFM: ghost fluid method; LS: level set method; PF: phase field method; VOF: volume of fluid method.

a

The best agreement range compared to the analytical results.

Large deformation and breakup of the compound drop under high electric field

Figure 9 shows the variation in core–shell and shell–continuous medium interfaces of the compound drop with different electrical parameters under high electric field. When the ratio of ${\kappa }_{\mathrm{AB}}$ to ${\epsilon }_{\mathrm{AB}}$ is equal to that of ${\kappa }_{\mathrm{BC}}$ to ${\epsilon }_{\mathrm{BC}}$ (Figure 9(a)), the charge property at the core–shell interface is the same as that at the shell–continuous medium interface. As can be seen in Figure 9(a), when the ratio is more than 1, the core–shell and shell–continuous medium interfaces are both elongated. As the ratio is less than 1, these two interfaces are both flattened. Moreover, these tendencies become more obvious as the ratio is far away from 1. When the ratio of ${\kappa }_{\mathrm{AB}}$ to ${\epsilon }_{\mathrm{AB}}$ is equal to that of ${\epsilon }_{\mathrm{BC}}$ to ${\kappa }_{\mathrm{BC}}$ (Figure 9(b)), the charge property at the core–shell interface is contrary to that at the shell–continuous medium interface. As shown in Figure 9(b), as the ratio is less than 1, the core–shell and shell–continuous medium interfaces are flattened and elongated, respectively. However, it is contrary as the ratio is more than 1. Similarly, when the ratio is far away from 1, these changes in core–shell and shell–continuous medium interfaces also become more evident.

Figure 9.

Variation in core–shell and shell–continuous medium interfaces of compound drop with electrical parameters under high electric field. (a) The color black, magenta, blue, red, cyan correspond to ${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}={\kappa }_{\mathrm{BC}}/{\epsilon }_{\mathrm{BC}}$  = 4, 2, 1, 0.5, and 0.25, respectively. (b) The color black, red, blue, magenta, cyan correspond to ${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}={\epsilon }_{\mathrm{BC}}/{\kappa }_{\mathrm{BC}}$  = 0.25, 0.5, 1, 2, and 16, respectively (E₀ = 1.5 MV/m, ${\rho }_A={\rho }_B={\rho }_C=980\mathrm{kg}/{\mathrm{m}}^3$ , ${\mu }_A={\mu }_B={\mu }_C=0.01\mathrm{Pa}\mathrm{s}$ , ${\sigma }_A={\sigma }_B={\sigma }_C=0.3\mathrm{N}/\mathrm{m}$ , and $R_0=1\mathrm{mm}$ ): (a) ${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}={\kappa }_{\mathrm{BC}}/{\epsilon }_{\mathrm{BC}}$ and (b) ${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}={\epsilon }_{\mathrm{BC}}/{\kappa }_{\mathrm{BC}}$ .

Figure 10 shows the breakup processes of compound drops for different electrical parameters. As shown in Figure 10, the breakup first occurs in the core drops and near the symmetry axis. Figure 11 shows the distribution of electric field, charge, velocity, and pressure when compound drops nearly break up. As can be seen in Figure 11(a), under the action of electric field, the core–shell and shell–continuous medium interfaces are driven toward the opposite directions due to the positive charges at the shell–continuous medium interface and the negative charges at the core–shell interface. When the two interfaces contact, part of positive and negative charges cancel each other, leading to a sharp decline in electrical force. The charges concentrate near the largest curvature of the core–shell interface (Figure 11(a)), thus forcing the core drop to break up into two drops. One attaches to the shell–continuous medium interface (Figure 10(a)) due to the fact that the charges are neutralized and then electric force cannot overcome the viscous force, the other reaches a balance under the action of electric force and surface tension (Figure 10(a)). As shown in Figure 11(b), since the charges at the core–shell interface are much more than those at the shell–continuous medium interface, the movement of compound drop is almost determined by the core–shell interface. This means that the breakup of shell–continuous medium interface may be driven by the core–shell interface, and consequently a new compound drop forms.

Figure 10.

The breakup of compound drop for different electrical parameters under high electric field (E₀ = 3 MV/m, ${\rho }_A={\rho }_B={\rho }_C=980\mathrm{kg}/{\mathrm{m}}^3$ , ${\mu }_B=0.02\mathrm{Pa}\mathrm{s}$ , ${\mu }_A={\mu }_C=0.01\mathrm{Pa}\mathrm{s}$ , ${\sigma }_A={\sigma }_B={\sigma }_C=0.3\mathrm{N}/\mathrm{m}$ , and $R_0=1\mathrm{mm}$ ): (a) ${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}={\epsilon }_{\mathrm{BC}}/{\kappa }_{\mathrm{BC}}=5$ and (b) ${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}=60{\kappa }_{\mathrm{BC}}/{\epsilon }_{\mathrm{BC}}=1.2$ .

Figure 11.

The distribution of electric field, charge, velocity and pressure for the forthcoming breakup of compound drop in Figure 10. (E₀ = 3 MV/m, ${\rho }_A={\rho }_B={\rho }_C=980\mathrm{kg}/{\mathrm{m}}^3$ , ${\mu }_B=0.02\mathrm{Pa}\mathrm{s}$ , ${\mu }_A={\mu }_C=0.01\mathrm{Pa}\mathrm{s}$ , ${\sigma }_A={\sigma }_B={\sigma }_C=0.3\mathrm{N}/\mathrm{m}$ , and $R_0=1\mathrm{mm}$ ): (a) ${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}={\epsilon }_{\mathrm{BC}}/{\kappa }_{\mathrm{BC}}=5$ and (b) ${\kappa }_{\mathrm{AB}}/{\epsilon }_{\mathrm{AB}}=60{\kappa }_{\mathrm{BC}}/{\epsilon }_{\mathrm{BC}}=1.2$ .

Conclusion

Based on the ternary phase field method, a finite element model is developed to simulate the EHD of the compound drop whose core phase, shell phase, and continuous phase are different in this article. The simulation results of small deformation of the compound drop under weak electric field show that the proposed model may have a good accuracy when the conductivity ratio of two adjacent fluids ${\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$ for the compound drop is in the range of 0.1–10. It is found from the simulations of large deformation and breakup of the compound drop under high electric field that the ratio of electrical parameter of two adjacent fluids ${\kappa }_{\mathrm{io}}/{\epsilon }_{\mathrm{io}}$ ( ${\kappa }_{\mathrm{io}}={\kappa }_{\mathrm{in}}/{\kappa }_{\mathrm{out}}$ , ${\epsilon }_{\mathrm{io}}={\epsilon }_{\mathrm{in}}/{\epsilon }_{\mathrm{out}}$ ) for the compound drop significantly affects the amount of charge at the corresponding interface. The difference in charge distribution is the key to the prediction of breakup form of the compound drop. This study provides the basic theory for predicting the behavior of three immiscible fluids under electric field and contributes to exploring new processes. In further work, a series of simulations will be conducted with varied electrical capillary number over a wide range of physical parameters and at different volume fractions of inner droplet to predict the qualitative and quantitative behavior of large deformation and breakup of the compound drop in the presence of an electric field.