Dynamics of double emulsion interfaces under the combined effects of electric field and shear flow

https://doi.org/10.1007/s00466-021-02045-x O R I G I N A L P A P E R

Dynamics of double emulsion interfaces under the combined effects of

electric field and shear flow

Roozbeh Saghatchi1,2,3· Murat Ozbulut4· Mehmet Yildiz1,2,3 Received: 2 April 2021 / Accepted: 5 June 2021

© The Author(s), under exclusive licence to Springer-Verlag GmbH Germany, part of Springer Nature 2021

Abstract

In this paper, the dynamics of a 2D double emulsion under the combined effects of the electric field and shear flow are studied by using an incompressible smoothed particle hydrodynamics (ISPH) method. Six different systems are used, each corresponding to different electrical properties. The effects of capillary and electrical capillary numbers, and core to shell radius ratio on the deformation and orientation angle of core and shell droplets are discussed thoroughly. It is shown that the deformation is highly dependent on these values, as well as the electrical properties. Electric force components and hydrodynamic stresses on the core-shell and shell-medium interfaces are calculated and discussed in detail. It is demonstrated that in some systems, a breakup occurs, which can be circumvented by changing the capillary and electrical capillary numbers as well as the core to shell droplet radius ratio. Finally, different breakup forms are investigated comprehensively.

Keywords Double emulsion· Electric field · Shear flow · Multiphase flow · Incompressible smoothed particle hydrodynamics method· Breakup

1 Introduction

Double emulsions (compound droplet) are structures con-sisting of two immiscible fluids, namely, a core fluid fully encapsulated by a shell fluid, which are dispersed in a third fluid medium. The distinctive features of double emulsions are due to the presence of a shell fluid that forms an encap-sulating layer around the core fluid. This makes double emulsions highly desirable in a vast number of applications,

B

1 _{Faculty of Engineering and Natural Sciences, Sabanci}

University, 34956 Tuzla, Istanbul, Turkey

2 _{Integrated Manufacturing Technology Research and}

Application Center, Sabanci University, 34956 Tuzla, Istanbul, Turkey

3 _{Composite Technologies Center of Excellence, Sabanci}

University-Kordsa, 34906 Pendik, Istanbul, Turkey

4 _{Naval Architecture and Marine Engineering Department,}

Faculty of Engineering, Piri Reis University, Istanbul, Turkey

including drug delivery [1], food production [2], cosmetics [3], wastewater treatment [4], amongst others.

External forces such as electric and shear forces play a vital role in the preparation and manipulation of emulsions, and both of these phenomena are of an important relevance to microfluidic applications and lab-on-a-chip technologies. Therefore, understanding the concurrent effect of electric and shear forces on the dynamics of emulsions is critical in terms of having a better control on the emulsification process.

Electrohydrodynamics (EHD) of a single emulsion, which consists of a suspended single droplet in a fluid medium, have been widely studied by many researchers. Taylor conducted a theoretical study based on an extensive simplification where both droplet and medium fluids were treated either perfect dielectrics (insulators) or perfect conductors, which causes the droplet to deform only in the direction of the applied elec-tric field (prolate deformation) [5]. Later, Allan and Mason experimentally observed that the droplet could also elongate in the direction perpendicular to the electric field (oblate deformation) [6]. Accordingly, Taylor corrected his model and introduced the leaky dielectric model considering the fluids to be slightly conductive that allows for the accumu-lation of free electric charges on the interface [7]. Based on the leaky dielectric model, various studies have been con-ducted considering different aspects of single emulsion EHD

[8–11]. Dynamics of viscous single emulsions subjected to a linear shear in a Couette device have also been analyzed quite extensively using various methods [12–17]. In this problem, surface tension and inertia forces play significant roles in the dynamics of the emulsions.

Unlike a single emulsion droplet, relatively few studies have been conducted on the dynamics of double emulsions under the effects of electric and shear forces. For example, Tsukada et al. [18] conducted both numerical and experimen-tal studies on the EHD of a double emulsion and observed that the single and double emulsions exhibit different behav-ior at the same condition. They also investigated the effects of the electric field strength and the volume ratio on the cir-culations induced in each phase and the drop deformations. Spasic et al. [19] studied the EHD effect on the rupture of double emulsions. They focused on the electrocoalescence process that occurs during the breaking of double emulsions. Santra et al. [20] performed a numerical study on the EHD of a double emulsion consisting of leaky dielectric and perfect dielectric fluids in the confined microfluidic domain. They evaluated the various parameters, such as electrical capillary, confinement ratio, viscosity ratio on the steady-state and tran-sient behavior of the droplet. Abbasi et al. [21] conducted an experimental study on the EHD deformation and breakup of double emulsion. They concluded that the breakup modes are dependent on the viscosity, electrical properties, and volume fraction of core/shell fluids.

Hua et al. [22] numerically investigated the effect of shear flow on the dynamics of a double emulsion and showed the impacts of droplet radius ratios, surface tension ratio, and core droplet eccentricity on the deformation and inclination angle of the double emulsion. Chen et al. [23] performed a numerical study on the deformation and breakup of double emulsions located in a shear flow and discussed the various breakup modes that occurred during their simulations. Luo et al. [24] discussed the mechanisms of the core and shell deformation of a double emulsion subjected to a shear flow. They also investigated the effect of capillary number and cores droplet size on the deformation of double emulsions.

Although many researchers have studied the behaviors of single emulsions under the combined effects of electric and shear force [25–27], the dynamics of double emulsions have only been investigated under the simultaneous influence of either pressure-driven (Poiseuille) flow [28] or extensional flow and [29,30] electric field. Effort on scrutinizing the concurrent effect of electric field and shear flow on the dou-ble emulsions are limited to very few recent papers [31,32]. Study of Santra et al. [31] investigated the deformation dynamics of a double emulsion in the presence of confined shear flow while the detailed electric forces analysis on the emulsion interfaces is conspicuously absent. Borthakur et al. [32] have discussed this important issue superficially while considering only the limited combination of

shell-core fluids. Moreover, both studies utilized the mesh-based methods, in other words, Finite Element based commercial COMSOL Multiphysiscs software [31] and Finite Volume based BASILISK open-source code [32].

Due to the dependency of mesh-based methods on the domain grid, it can be troublesome and arduous to deal with the complex geometries and deforming boundaries. To circumvent such computational difficulties, one may use mesh-free methods. Smoothed particle hydrodynamics (SPH) is a very powerful method for solving complex phys-ical problems in the field of fluid mechanics [33]. Due to its mesh independency and Lagrangian nature, SPH can be applied to simulations involving complicated geometries [34], and flow structures with extensive deformations [35].

Motivated by these considerations, in the present study, the dynamics of a 2D Leaky dielectric double emulsion is thoroughly investigated under the combined effect of applied electric field and shear flow using mesh-free incompress-ible smoothed particle hydrodynamics methods (ISPH). To address the impacts of electrical property ratios of the core, shell, and medium, six different fluid configurations are numerically handled. All the components of the electrical and hydrodynamics forces are calculated on the droplet inter-faces and within the flow domain correspondingly and their influence on the droplet deformation are discussed in detail considering both their magnitude and sign. Effects of cap-illary and electrical capcap-illary numbers, droplet radius ratio, electrical conductivity and permittivity ratios on the emul-sion deformation and angular orientation are explored.

The organization of the paper is as follows. In Sect.2, the governing equation of the problem and its non-dimensional parameters are mentioned. In Sect.3, the numerical method, ISPH, is introduced together with the discretization of the governing equations with the the ISPH framework. In Sect.4, numerical problems are defined and the accuracy of our numerical method is investigated by verifying the obtained results with those available in literature other studies. In Sect. 5, numerical results are presented for droplet behav-ior in the confined domain, and finally concluding remarks are provided in Sect.6.

2 Governing equations

Assuming that the shell, core and background medium are incompressible and Newtonian fluids and subjected to an external electric field, the flows inside and outside the double emulsion are governed by the continuity and conservation of linear momentum equations as follows [36]:

∇ · u = 0, (1)

ρDu

where u, p,ρ, t, andτ are the velocity vector, pressure, density, time, and viscous stress tensor defined as τ =

μ(∇u + (∇u)†_{) with μ being a viscosity, respectively. The} dagger symbol † as a superscript indicates transpose opera-tion. fsand feare surface tension and electric force densities

in the given order. D/Dt is the material time derivative and is defined as D/Dt = ∂/∂t +u ·∇, and ∇ is the nabla operator. It should be noted that the gravitational force is neglected.

Surface tension force density fs on the interface is

for-mulated through utilizing continuum surface force (CSF) method [37],

fs = γ κnδ, (3)

in whichγ represents the surface tension coefficient, κ is the local interface curvature obtained from the divergence of the unit surface normal vector n, andδ is the Dirac delta function. To be able to write the relevant set of governing equations together with the associated jump or interface conditions for the electrostatic phenomena to be investigated in this study, several important assumptions should be made in place. The multiphase fluid system is polarizable and non-magnetizable. Furthermore, the system under investigation can be accurately described by a quasi-static electric field model. The latter assumption presumes that the time varia-tion of electric field (also referred to as dynamic current) in Ampere’s law is negligible whereby magnetic induction can be ignored. Therefore, the Ampere’s law is dropped off the list of governing equations. Moreover, neglecting magnetic induction, Faraday’s law,∇ × E = ∂B/∂t, requires that the curl of electric field vector E be zero thereby leading to the fact that gradient of electric field vector is sym-metric, ∇E = (∇E)†. Upon using the vector identity of ∇ × (∇φ) = 0 for any arbitrary scalar field, the electric field vector can be expressed as the gradient of the electric potential as [38],

E= −∇φ. (4)

The total volume current is introduced as J= qvu+j [39]. Herein, q_vis the volume-charge density of free charges, the term q_vu is the convection current due to the fluid induced motion free charges, and j is the volume conduction cur-rent density, ohmic curcur-rent, that is coupled with the electric field vector through the relation j = σ E or equivalently j= −σ ∇φ with σ being the electrical conductivity. Remem-bering the assumption that the mediums handled in this study are non-polarizable, there are no bound charges in the domain of interest. The only motion of charges here is due to the free charges.

The Gauss’s law for electricity can be written in terms of the electric displacement vector, D= εE as [39]:

∇ · D = qv, (5)

whereε is the electric permitivity.

The conservation of charge can be readily formulated through taking the divergence of the differential form of Ampere’s law, and subsequently employing the vector rela-tion ∇ · (∇ × B) = 0 (the divergence of the curl of any arbitrary vector field, in this case magnetic induction vector, B, is equal to zero) together with the Gauss’s law as [38],

Dq_v

Dt = −∇ · j, (6)

The interface or jump conditions for Faraday’s law, Gauss’s law for electricity and the conservation of charge can be respectively formulated as n× [E] = 0, n · [D] = qs, and

¯δqs/δt+∇s·K+n·[J − q_vv]= 0 [39,40]. Here, the symbol

[·] indicates the jump of the enclosed quantities across the discontinuity surfaceξ. As such, [β] = β+− β−withβ+ andβ−being respectively the values of any arbitrary fieldβ on the positive and negative sides of the discontinuity surface

ξ. ¯δqs/δt= ∂/∂t+ (v · n)(n · ∇) is the total time derivative

in following the motion of the discontinuity surface along its unit normal vector n, v is the velocity of the discontinuity surface, K= k+qsv is the total surface current composed of

the surface conduction and convection currents, respectively,

qsis the surface-charge density of free charge, and∇s is the

surface gradient operator. The former two jump conditions in the above given order state that tangential component of the electric field vector is continuous across the discontinuity surface while the normal component of electric displacement is discontinuous therein and balanced by the surface-charge density of free charge.

Electrohydrodynamics of the fluid system is investigated through coupling the electrostatics and hydrodynamics phe-nomena by using the divergence of the Maxwell stress tensor in the linear momentum balance equation in Eq. (2). As such, the volumetric electrical force vector can be written as fe = ∇ · te. In literature, one may find various forms

of and derivation methodology for the Maxwell stress ten-sor [10,36,40–42]. The one used in this study is of the form excluding the contribution from the induced magnetic field such that t = D ⊗ E − 0.5(D · E) [40]. In what follows, one can write,

fe= qvE−

1

2E· E∇ε, (7)

Here, the first term on the right hand side of this equation represents the Coulomb force due to the interaction of free

charges with the electric field while the second term corre-sponds to the polarization force because of pairs of charges. In this study, the multi-phase fluid system is treated as conducting-conducting (leaky dielectric medium) due to the assumption that the viscous time scale of the fluid motion is much larger than the bulk relaxation time, namely, t_μ  t_e Here, t_μ = ρL2/μ where L is the characteristic length scale and t_e = ε/σ [40]. Due to the leaky dielectric assump-tion, both volume and surface charges can reach steady state whereby Dq_v/Dt = 0 and ¯δqs/δt = 0. Moreover, if the

charge conservation equation for the discontinuity surface is written in a non-dimensional form, recalling that t_μ  t_e, one can prudently neglect the contribution of the surface current to the physics of problems. Furthermore, if the discontinuity surface is regarded as the material interface, then, u = v. Under all these assumptions, one write the conversation of charge in the volume and on the discontinuity surface, respec-tively, as [10]:

∇ · (σ∇φ) = 0, n· [σ∇φ] = 0 (8)

Moreover, the final version of the Gauss’s law for elec-tricity for the volume and the discontinuity surface follows as [10];

∇ · (ε∇φ) = qv, n· [ε∇φ] = 0, (9)

The ratio of physical properties in the domain is defined as:

Dmn= ρm/ρn, Vmn= μm/μn, Pmn= εm/εn,

Cmn= σm/σn, (10)

whereas subscripts m and n refer to mth and nth fluids, respectively (m, n = 1, 2, 3). Subscripts 1, 2 and 3 denote the core (inner droplet), shell (intermediate droplet) and back-ground or medium fluid phases, respectively, whereas the double subscript mn represents the interface between mth and nth fluids. Important non-dimensional parameters of the current study are Reynolds number (Re = ρ2Ud2

μ2 ), Weber

number (We = ρ2U2d2

γ23 ), Electro-Weber number (Ew =

ρ2U2

ε2E_∞2 ), capillary number (Ca =

We Re = μ2

U

γ23), and

electri-cal capillary number (Ec=We_Ew = ε2E_∞2d2

γ23 ), where E∞is the

magnitude of the electric field vector, U is the velocity of top and bottom plates and d is the droplet diameter. We will also use dimensionless time t in our simulation which is defined as t= t_2HU where 2H is the domain width.

3 Numerical method

We have chosen smoothed particle hydrodynamics (SPH) method to solve the governing equations of this study. SPH had been widely used by many researchers to simulate the broad range of physical problems including free surface [43, 44], turbulent [45,46], multi-phase [47], heat transfer [48,49], biological problems [50]. This section looks introduces the SPH method used in the current study.

In the present SPH method, a color function ˆc is utilized to track the interface between different phases. For a two-phase problem, for example, a value of zero is assigned to one phase and unity for the other and these values remain constant during the entire simulation. To achieve a smooth and finite transition region which can improve the accuracy of com-puted interface features such as interface unit normal vector and curvature, thereby enhancing the convergence behavior and the robustness and fidelity of numerical simulations, at each time step, the initial color function is smoothed out for particles in the vicinity of interface by using the following equation: ci= Jn  j=1 ˆcjWij ψi , (11) in which ψi = Jn

j₌₁Wij is the particle number density and Wij is the interpolation kernel which is a function of

h, the smoothing length, and the magnitude of distance

vec-tor, rij= ri− rj, between the positions of particle of interest i and its neighboring particles j [51,52], and the summation is performed over all neighbors Jn of particle i. It is noted

that the two-dimensional quintic spline kernel is used in this study [10].

Unit surface normal vector n = ∇c/ |∇c| and Dirac delta function δ  |∇c| in Eq. (3) are calculated using the smoothed color function through applying a constraint |∇ci| ≥ α/h to avoid the possible erroneous normal vector in the calculation of surface tension force [53].α is a numer-ical constant with the value ofα = 0.08 to obtain reliable and accurate results [54].

Thermodynamic properties (i.e., ρ) and transport coef-ficients (i.e., σ, ε and μ) can experience a sharp jump across the interface between core and shell or shell and medium, thereby adversely affecting the accuracy and robust-ness of the numerical simulations and consequently leading to erroneous results. Thus, an averaging scheme is used to interpolate the phase properties and eliminate the sharp jump across the interface. In the current study, a weighted arith-metic mean (WAM) approach is used and defined as:

whereχ stands for any hydrodynamic or electrical fluid prop-erties (i.e.,ρ, μ, σ and ε).

Temporal discretization of the governing equations is accomplished by predictor–corrector scheme. The first-order Euler approach is implemented for the time discretization along with the Courant criterium to determine the variable time step size,Δt = ζ h/umax, where umaxis the largest parti-cle velocity magnitude andζ is taken to be equal to 0.25 [54]. The time marching proceeds to a temporary or intermediate stage using: r∗_i = r(n)_i + u(n)_i Δt + δr_i(n), (13) u∗i = u(n)i + 1 ρi(n)  ∇ · τi+ fsi+ fei _(n) Δt, (14) while the densities are calculated using the following rela-tions: ψi∗= Jn  j=1 W_ij∗, ρ_i∗= miψi∗. (15) Here, starred entities represent the intermediate values of the given variables and superscript(n) denotes values of vari-ables at the nth time step. Here, it should be noticed that the linear momentum balance equation at the intermediate stage given in Eq. (14) does not include the pressure gradient term, hence leading to an intermediate velocity field which is not divergence-free. To prevent the particle clustering, artificial particle displacement (APD) method is used [44,55–57]. In this method, the APD vector is implemented throughδr(n)_i as, δri(n)= α ⎡ ⎣umax Jn  j=1  rij r_ij3r 2 avg,i ⎤ ⎦ (n) Δt. (16)

ravg,i =_jJ₌₁n rij/Jnis used to find the average cut-off

dis-tance for a given particle.αis a coefficient that controls the magnitude of the APD vector. It is found that theαequal to 0.06 is sufficient to render a uniform particle distribu-tion through eliminating the possible non-physical particle clustering without forfeiting the physics of the problem in question. Here, it is important to note that the values of APD for each particle at each time step are much smaller than the physical movement of these particles. Furthermore, APD is an odd function and therefore is only operational for flow regions with particle clustering whereas in the regions with uniform particle distribution, it takes zero value [47].

In the ISPH method, the pressure distribution is found by solving the pressure Poisson equation with a source term being the divergence of the intermediate velocity as given in Eq. (17). As a result, the pressure values can be obtained

in a way that incompressibility conditions is enforced. To complete the time marching scheme, the velocities and then the positions of the particles are corrected via Eqs. (18) and (19), respectively. ∇ ·  1 ρ∗ i ∇ pi(n+1) =∇ · u∗i Δt , (17) u(n+1)_i = u∗i − 1 ρ∗ i ∇ p(n+1)i Δt, (18) r(n+1)_i = r(n)_i +1 2  u(n)_i + u(n+1)_i  Δt+ δr(n)i . (19) First derivative and Laplace operator of vector functions are approximated as [58]: ∂ fm i ∂xk i a_ikl = j 1 ψj  f_jm− f_im ∂Wij ∂xl i , (20) ∂ ∂xk i  ϕi ∂ fm i ∂xk i a_iml=8 Jn  j=1 2ϕiϕj ϕi+ϕj 1 ψj  f_im− f_jm r m ij r2 ij ∂Wij ∂xl i . (21) Here, akl_i = j rk ji ψj ∂Wij ∂xl i

is a corrective second rank ten-sor that eliminates particle inconsistencies [59] whileϕ may denotev, μ, ε, or σ, where v denotes the specific volume and it is equivalent toρ−1. Equation (8) for the volume and the left hand sides of Eq. (17) as well as the left hand side of Eq. (5) are discretized as:

∂ ∂xk i  ϕi∂ fi ∂xk i  2+ akk_i  = 8 Jn  j=1 2ϕiϕj ϕi+ ϕj 1 ψj  fi− fj  rk ij r_ij2 ∂Wij ∂xk i . (22) Multiple boundary tangent (MBT) method is used to apply the boundary conditions on the walls [60].

4 Numerical consistency and accuracy

studies

4.1 Problem setup

The schematic of the physical system is shown in Fig.1. As can be seen from the figure, a 2D double emulsion droplet consists of concentric core and shell droplets with initial radii of r1and r2, respectively, and is suspended in a third medium fluid subjected to a simple shear flow in a square domain with dimensions of W = 2H. A uniform electric field E_∞ is applied in the downward direction. It is assumed that there is no relative motion between droplets and ambient fluid

Fig. 1 Schematic representation of double emulsion droplet in a

Cou-ette device under the effect of applied constant electric field

initially, and no gravitational force is present. No-slip bound-ary condition is applied on top and bottom boundaries with constant velocities U and−U, respectively, while periodic boundary condition is enforced on the right and left bound-aries of the domain to reproduce the shear-driven (Couette) flow. As for the electric potential, Dirichlet and Neumann boundary conditions are employed for horizontal and verti-cal walls, respectively. All particles inside and outside the double emulsion are arranged in the form of a uniformly spaced Cartesian grid.

4.2 Particle resolution independency and numerical

code accuracy

Particle resolution independency test is carried out by simu-lating a single emulsion droplet (r1/r2= 0) under the effect of pure shear flow (E_∞= 0) whereby the deformation and angular orientation of the droplet is analyzed. Considering different particle resolution X/d = 80, 110, 160, and 200 where X/d is the number of particles per unit of droplet’s initial diameter, numerical simulations are performed for Re 1 and Ca = 0.094, which imitate the conditions used in the experimental study of Bruijn [12]. It should be noted that the density and viscosity ratios are considered as

Di j = Vi j = 1. The extent of the deformation for the single

emulsion droplet is calculated as,

D= l1− l2 l1+ l2,

(23) where l1and l2are the longest and shortest diameters of the deformed droplet, respectively. In Fig.2, the results of these simulations are given both qualitatively and qualitatively for

t= 0.25 . As seen from the figure, upon increasing the

parti-cle resolution from x/d = 80 to 110 and 160, the numerical

accuracy also gets improved. However, further increase in the particle resolution does not bring about significant change in the obtained results. Thus, x/d = 160 is selected as a refer-ence particle resolution.

The fidelity and the reliability of our current in-house code has been proven in various cases including the EHD simu-lation of a droplet [10], and other multiphase flow problems [47,61]. However, for the additional validation and verifi-cation tests, the deformation of a single emulsion droplet is simulated under two distinct conditions: pure shear flow, and pure EHD. For the former test case, the same simulation condition is considered as in the case of the particle resolu-tion dependency analysis. SPH results are compared with the available results from previous numerical and experi-mental studies. To further demonstrate the accuracy of the SPH method, another simulation is also performed using the open-source Finite Volume code, OpenFOAM [62] based on interFOAM multiphase solver. This solver is based on a modified volume of fluid (VOF) approach [63]. Having performed a mesh and time step independency analysis, a uniform Cartesian mesh composed of 120 by 120 cells and 10−4time step size are selected for all the simulations with OpenFOAM, while the CFL condition is utilized to adjust the time step size throughout the simulation. Gauss upwind scheme is used for the discretization of convection term and the Euler scheme is used for the time marching with PISO algorithm for the correction. Transient droplet deformation is calculated for different values of Ca number and the results are shown in Fig.3a, b. Reynolds number is set to Re= 1 and Re 1 for Fig.3a, b, respectively. As shown in these fig-ures, a close agreement is observed between SPH and other numerical/experimental results as well as the values obtained by OpenFOAM. A qualitative comparison is also realized between SPH and OpenFOAM result and shown in Fig. 4 for Ca = 0.194 at t = 1. As can be readily inferred from this figure, both numerical results exhibit similar patterns of streamlines and contour of normalized velocity magnitude.

For the pure EHD case, SPH results are compared with the analytical results of Feng and Scott [64]. Considering the leaky dielectric model, the droplet may deform into two distinct configurations forming prolate or oblate shapes. The prolate shape is achieved when the droplet elongates in the direction of the applied electric field whereas the transverse elongation of the droplet is known as oblate deformation. Here, the deformation is calculated using the same equation as (23). However, here, l1and l2are considered as the two main axes of the deformed droplet in the directions parallel and perpendicular to the applied electric field, respectively. So, positive and negative values calculated by Eq. (23) repre-sent prolate and oblate deformations, respectively. Feng and Scott [64] used the following equation for calculating the deformation of the droplet as:

Fig. 2 Effect of particle resolution on the droplet deformation in pure shear flow. a Quantitative, and b qualitative results

D= fdE

2

∞εd(d/2)

3(1 + C )2_Pγ, (24) where fdis the discriminating function defined as:

fd= C2+ C + 1 − 3P. (25)

The comparison of the results for the SPH method and analytical solution is rendered in Fig.3c forDi j = Vi j = 1, Pi j = 0.5, Ci j = 2 and Pi j = 5, Ci j = 0.2 for prolate

and oblate cases, respectively. As shown in this figure, SPH results perfectly match with the analytical data at lower val-ues of Ec whereas they start deviating from the analytical results for a comparatively higher values of Ec. This is due to the assumptions made in the theory which considers the droplet to remain almost circular even during the deforma-tion. Thus, the analytical predictions are only valid for small deformations [10,36,65].

5 Numerical results

In this section, the dynamics of a double emulsion droplet under the combined effects of electric field and shear flow is presented by considering six different systems of fluids, as tabulated in Table1. The first three systems are selected from experiments in micro- and bio-fluidics applications. The properties of the remaining systems are chosen hypotheti-cally but carefully to scrutinize all possible forms of double emulsion deformation. In all simulations, core and medium fluids are identical, which is the condition that has also been utilized in some previous numerical and experimental

stud-ies [18,21,66,67]. Furthermore, in all cases, viscosity and density ratios are considered to be unity,Dmn = Vmn = 1.

Double emulsion dimensions are chosen as r2/H = 0.25 and core to shell droplet radius ratio is r1/r2 = 0.5, unless stated otherwise.

Figure5represents the streamlines inside the domain for all systems at t = 0.4. All sub-figures demonstrate the deformed double emulsions at Ca = 0.4 and Ec = 0.4. As shown in this figure, the deformed shape of the double emulsions is closely related to the conductivity ratios of the fluids. For systems I, II, III, and VI, which have relatively large values of conductivity ratio (C21 > O(10−1)), double emulsions are deformed and oriented without breakup. How-ever, in systems IV and V with smaller conductivity ratios (C21 = O(10−2)), breakup occurs. For the double emul-sions with no breakup, there is a high strength vortex inside the core, and two weak vortices at both sides of it. For sys-tems IV and V, however, the pattern of streamlines is more complex due to the breakup of the droplet.

To see the effect of electric field and shear flow separately, a streamline of the domain is shown in Fig.6. In the sub-figures (a) and (b), the dynamics of the double emulsions are presented in the presence of pure electric field and pure shear flow, respectively. As it is shown in Fig.6b, all systems exhibit a similar behavior under the effect of pure shear flow regardless of their different electrical properties as it should be. Since in the pure EHD, the flow field is symmetric in both x- and y- directions, only half of the domain is shown in Fig.6a. It should be noted that all figures are plotted for the same instant of the time (i.e., t = 0.4) as in Fig.5. Three different vortex patterns can be observed in the flow field

Fig. 3 The comparison of the single emulsion droplet deformation

under the a, b pure shear flow and c pure EHD between SPH results and previous studies. For the pure shear flow, numerical results of Chinyoka et al. [14], Mahlmann and Papageorgiou [15], and Sheth and

Pozrikidis [13], and experimental results of Bruijn [12] are used here. Results of Finite Volume simulation with OpenFOAM are compared as well. For the EHD case, analytical data of Feng and Scott [64] is used

Table 1 Dielectric properties that are used in this study

System Ambient and core fluids σ1, σ3_mS ε1, ε3F_m Shell fluid σ2_mS ε2_mF C21 P21 References

I Silicon oil 1 2.67 × 10−12 2.66 Corn oil 1.06 × 10−11 3.24 3.9700 1.2180 [68]

II Medium 1.8 × 10−2 80 Cytoplasm 0.5 75 27.780 0.9375 [69]

III Corn oil 1.06 × 10−11 3.24 Silicon oil 1 2.67 × 10−12 2.66 0.2519 0.8210 [68]

IV Hypothetical 27.78 0.821 Hypothetical 1 1 0.0360 1.2180 –

V Hypothetical 50 10 Hypothetical 1 1 0.020 0.10 –

Fig. 4 The comparison of normalized velocity magnitude contour and streamlines for the single emulsion droplet under the pure shear flow between

SPH (left) and OpenFOAM (right) results at Ca= 0.194, and t = 1

under the pure effect of the electric field. In systems II, V, and VI, there is only one vortex in each quarter. However, for the other systems, there is a pair of vortices in each quarter, one in the ambient fluid, and the other in either shell (systems III and IV) or core (system I) droplets. The direction of the ambient vortex depends on the conductivity and permittivity ratios such that the velocities run from the top and the bottom towards the left and the right of the droplet for the systems I, II, V and VI while it rotates in the opposite direction for the others. We have shown in our previous study that for the single emulsion, velocity in medium fluid runs from poles to the equators whenC > P and vice versa [70]. For the double emulsion, however, this is only true when the deviation ofP from unity is not too high. Thus, for the systems I, II, III, and IV, the circulation direction is similar to the single emulsion studied in our previous publication [70]. For the systems V and VI, circulation behaves reversely due to the significant difference betweenP and unity.

To shed light on why double emulsion breaks up and also vortex location changes, we studied the individual effect of electrical properties of droplet, electric field strength and the shear flow characteristic on the deformation values of core/shell droplet as well as the droplet orientation. The deformation D and orientationθ of the core and shell droplets are obtained for each system at t = 0.4 and tabulated in Table2for different cases including pure EHD, pure shear flow, and combined EHD and shear flow. The orientation

angleθ = 90 and θ = 0 correspond to prolate and oblate elongations, respectively. It should be noted that in the pure EHD flow, the necessary condition for the double emulsion breakup is the large difference in deformation values between core and shell droplets (Dcore  Dshell). The breakup will

be further facilitated if the core and shell droplets deform into different shapes and have the opposite deformation sign (i.e., one assumesθ < 45 and the other one has θ > 45). These conditions are realized for systems I, III, IV, and V. However, breakup occurs only in the system V since the strength of the electric force is large enough to elongate the core droplet adequately. As a result, the core droplet splits the shell droplet and thereby lead to the breakup of the double emulsion. In the system V, both core and shell droplets are deformed into the prolate shape with much higher deforma-tion for the core. For example, the deformadeforma-tion values before the breakup at t = 0.2 are D = 0.281 and D = 0.055 for core and shell droplets, respectively. Applying a shear force in addition to the EHD would intensify the breakup tendency of these systems. As can be recalled from Fig.5, the applied shear force would cause the double emulsion for the sys-tem IV to undergo a breakup. Another interesting finding is related to the orientation angle of core and shell droplets. As was remarked by Hua et al. [22], the orientation angle of the core is always greater than or equal to the shell’s for the double emulsion under the pure shear flow. Our numerical results also confirm the claim of Hua et al. for pure shear

Fig. 5 Streamlines inside and outside of the deformed double emulsions under the combined effect of electric field and shear flow for Ca= 0.4

and Ec= 0.4 at the non-dimensional time t = 0.4

flow as well as the combined EHD and shear flow. Recall-ing from Fig.6, there is a pair of vortices in each quarter of the systems I, III, and IV whereas there is only one in the systems II, V, and VI. This can be reasoned out through con-sidering the orientation angles given in Table2for the pure EHD case. The difference between the orientation angle of core and shell droplets in the systems with a pair of vortices, isΔθ = 90◦whereas for the systems with only one vortex is

Δθ = 0◦_{. Moreover, in the systems with a pair of vortices (I,} III, and IV ), only for the system I, there are vortices inside the core droplet. As shown in Table2for the system I, the core and shell droplets deform into the oblate and prolate shapes, respectively unlike the system III and IV, in which the core and shell droplets deform into the prolate and oblate shapes, respectively.

Figure 7 present the variation of the absolute value of double emulsion deformation as a function of the non-dimensional time for both core and shell droplets under the effect of combined EHD and shear flow. In this figure, core and shell deformations are shown with solid and dashed lines, respectively. As mentioned previously, the deforma-tion of the core droplet is much higher for systems IV and V which undergo a breakup. Moreover, the core deformation

for the system V is much greater than that for the system IV. Consequently, the system V breaks up before the system IV. Nonetheless, the difference between core and shell deforma-tions are not significantly high for the rest of the systems.

To investigate the core and shell droplet deformations, interface forces should be evaluated in detail. As described in Sect.2, one of the important interfacial forces is the elec-trical force composed of two important components, namely, Coulomb and polarization forces. Recalling from Eq. (7), it can be easily concluded that the Coulomb force acts along the direction of the electric field while the direction of the polar-ization force is normal to the droplet interface in the opposite direction of the electrical permittivity gradient. These two forces are demonstrated in Fig.8for all systems at t = 0.1. The contours and directions of the Coulomb and polarization forces are shown in Fig.8a, b, respectively. In the systems I and II, the Coulomb force is exerted on the top and bot-tom poles of the shell droplet in the outward direction while the polarization forces act in the outward and inward direc-tions for the systems I and II, respectively. By comparing the magnitude of forces for the system II, the resultant elec-tric force is in the outward direction similar to the system I, which elongates both shell droplets into prolate shape. For

Fig. 6 Streamlines inside and outside of the deformed double

emul-sions under a the pure electric field with Ec= 0.4 and b the pure shear flow with Ca= 0.4 at the non-dimensional time t = 0.4. Due to the

symmetry in the y- direction of the flow in the presence of the pure electric field, only the half of the domain is shown

Table 2 Deformation and orientation of core and shell droplets in all six systems at t= 0.4 for a double emulsion under the effect of pure electric

field, shear flow and combined EHD and shear

Case System I System II System III System IV System V System VI

Core Shell Core Shell Core Shell Core Shell Core Shell Core Shell

Pure EHD D 0.001 0.015 0.043 0.046 0.121 0.004 0.345 0.016 – – 0.006 0.004 θ 0 90 90 90 90 0 90 0 – – 90 90 Pure shear D 0.044 0.179 0.044 0.179 0.044 0.179 0.044 0.179 0.044 0.179 0.044 0.179 θ 47.90 36.38 47.90 36.38 47.90 36.38 47.90 36.38 47.90 36.38 47.90 36.38 Combined effect D 0.049 0.190 0.108 0.250 0.134 0.192 – – – – 0.047 0.181 θ 51.63 42.09 63.67 51.72 75.60 31.20 – – – – 52.00 37.80

the other systems, Coulomb forces act in the inward direction whereas the polarization forces act toward the center of the systems III and II unlike the systems IV and VI. By compar-ing the magnitude of the Coulomb and polarization forces, the total electric forces lead to the oblate deformation of the shell droplet for the systems III, IV, and VI. For the system V, however, the magnitude of the electric force is higher at the sides of the shell droplet, thus deforming it into the pro-late shape. A similar analysis may be conducted for the core

droplet. Such investigation reveals that the direction of both Coulomb and polarization forces are toward the center for the system I only, and the resultant total electric force acts on the left and right sides of the core droplet and stretches it to deform into the oblate shape. For the other systems, however, the force direction and its point of application are such that the total electric force deforms the core droplet into a prolate shape.

Fig. 7 Time evolution of core (solid lines) and shell (dashed lines)

droplets deformation in Ca= 0.4 and Ec = 0.4

In addition to the electric forces, hydrodynamic interac-tion plays an important role in the double emulsion dynamics [24]. These interaction can be studied considering the hydro-dynamic stress tensor. Lou et al. [24] claimed that the hydrodynamic interaction mainly depends on the distribution of the stresses inside the shell droplet, which is significantly affected by the core size and the deformations of the shell and core droplets. Hence, pressure and viscous stresses are shown in Fig.9. Since the stress distribution of double emulsions in the systems I, II, III, and VI are almost similar, for the sake of brevity, only system III is selected as their representative and shown alongside the system IV and V. As shown in Fig.9, the pressure force has a significant contribution to the hydrody-namic interactions while the viscous force is trivial in these systems since the viscosity ratio is unity. Therefore, we only consider the pressure term in the following. The pressure value is high in the core droplet of the system III while the maximum pressure occurs at the shell droplet for the system V. The pressure distribution determines the curvature of the droplet interface and subsequently, its deformation value. As presented in Fig.9a for the systems IV and V, pressure in the core droplet has its lower value near the tips (high curva-ture) while it reaches its higher values at the equator of the core droplet (low curvature). Within the shell droplet, regions with maximum and minimum pressure values coincide with the different locations for the systems IV and V. As repre-sented in Fig.9a, there is high- and low-pressure regions at the equator and tips for the systems IV, respectively. For the system V, however, tips and equator correspond to the high- and low-pressure regions, respectively. The pressure is distributed uniformly inside the core and shell droplets for the system III, which is commensurate with the uniform

curvature of these droplets, thereby resulting in the low defor-mation value as represented in Fig.7. The pressure gradient inside the droplets for the system V is higher than that of the system IV, which leads to the higher value of deformation in the system V.

To particularize the effect of shell droplet thickness on the double emulsion dynamics and deformation, two additional core to shell droplet radius ratios, i.e. r1/r2= 0.25 and 0.75 are also simulated. The results of these simulations are pre-sented in Fig.10a in terms of the deformation and the angular orientation values of the double emulsion. A detailed evalua-tion of this figure reveals that the variaevalua-tion of the deformaevalua-tion with respect to r1/r2for the core droplet has an opposite trend in comparison to that for the shell droplet. Namely, as the deformation increases for the core droplet, it decreases for the shell one. This behavior can be explained through referred to Eq. (24) given that the deformation of droplet is directly linked to its diameter. Note that the increase in r1/r2can be controlled by either increasing or decreasing the diameter of core or shell droplet, respectively. This can explain the oppo-site behavior of core and shell deformations. Droplets with the higher extent of prolate deformation would experience a higher inertial force under the shear flow. This force will compel the droplet to orient itself in a horizontal position. Thus, the orientation angle of the core droplet decreases by the increase in r1/r2. Remembering Fig.8, the direction of the electric force on the core droplet for the system III is such that it elongates the core, hence causing prolate deformation. Because of this elongation, the pressure force inside the shell would be higher at the poles and causes the shell droplet to elongate and to acquire a more prolate shape accordingly. The larger value of r1/r2corresponds to a thinner shell droplet, in which there is not enough space for the core droplet to deform freely, resulting in an increase in the pressure force at the shell poles. Thus, at higher r1/r2, the shell droplet is forced to elongate in the vertical direction and in turn has a higher orientation angle value.

To investigate the effect of electrical capillary number on the double emulsion dynamics, the capillary number is maintained constant (Ca= 0.4) while the electrical capillary number is changed for each system. The deformation value and the angular orientation are calculated at t = 0.4 and presented in Fig.10b. Regarding the shell droplet, the defor-mation value is observed to be directly related to the electrical capillary number for all the systems except the system VI in which there is no significant change in the deformation value. Similarly, the orientation angle is also affected by the elec-trical capillary number for all systems except the system III. These correlations can be explained considering the electric forces on the interface of the shell droplet as has been shown in Fig.8. It can be seen that the resultant electric force exerted on the top and bottom regions of the shell droplet is high for the system III and acts inwardly to the interface. Shear flow

Fig. 8 Comparison of a Coulomb force and b polarization force on the interface of core and shell droplets at t= 0.1 for all systems. Note that

Fig. 9 Comparison of a p, bμdu/dx, c (μ/2)(du/dy + dv/dx), and d μdv/dy in the double emulsion under the effect of combined electric field

induced angular orientation of the droplet leads to the occur-rence of a couple force that tends to rotate the droplet in the clockwise direction, hence decreasing its orientation angle. However, for the other systems, this couple force acts such that it rotates the shell droplet in the counterclockwise direc-tion thereby increasing the orientadirec-tion angle. Moreover, in the system VI, Coulomb and polarization forces on the shell droplet interface have a relatively similar order of magnitude and act in the opposite directions, thus resulting in a negligi-bly small net electric force. As a consequence, the variation in electrical capillary number does not have a considerable effect on the shell deformation for the system VI. It can be concluded that for those systems in which the shell droplet undergoes an oblate deformation in the absence of shear flow, their orientation angle decreases in mixed EHD and shear flow with an increase in the electrical capillary number. In the system IV, one can see a minimum value in the deforma-tion and orientadeforma-tion angle at Ec= 0.2. It should be noted that applying electric field on this system first reduces the defor-mation and orientation angle of the shell due to the pure shear flow. However, on further increasing the electrical capillary number, electrical forces on the interface become dominant and affect the droplet behavior. Likewise, core droplet behav-ior can be explained by investigating the electrical forces on its interface. As a result, the core droplet experiences a high deformation and orientation angle at high electrical capillary numbers. For the double emulsion in which breakup occurs, the variation of deformation with electric capillary number is significant. For instance, in the system IV, the deformation of the core droplet at Ec= 0.3 becomes 9.1 times of its value for the pure shear flow (Ec= 0).

The effect of capillary number (Ca) on the double emul-sion dynamics is shown in Fig.10c for constant electrical capillary number, (Ec = 0.4). The interplay between elec-trical and hydrodynamics forces determines the droplet deformation and orientation angle. As elaborated previously, core and shell droplets have a significant impact on each other. For instance, the deformation value for the shell droplet of the system III first increases up to Ca= 0.1. The further increase in the capillary number decreases the deformation value as the inertial force becomes dominant. At high capil-lary numbers, however, the shell droplet thickness decreases at the pole regions, which leads to an increase in the pressure near the poles (Fig.9a). As a consequence, the shell droplet elongates and hence acquires a larger deformation value. For the system II, there is an extremum at Ca= 0.1 which acts as a maximum point for the shell deformation. For the sys-tems I and VI, the variation of deformation as a function of the capillary number has an increasing, which points out that the hydrodynamics forces are the dominant in these systems. The orientation angle of shell droplet increases for the sys-tem III whereas it decreases for the syssys-tems I, II, and VI. Core droplet behavior can be explained in a similar manner.

It can be shown that the capillary number may enhance the likelihood of the breakup of double emulsion in the system IV. For the system V, however, breakup occurs even without the shear flow (Ca = 0) due to the high value of resultant electric force.

Transient behavior of double emulsions for the systems IV and V are demonstrated in Fig.11. As mentioned previously, due to the large deformation of the core droplets in the sys-tems IV and V, a breakup occurs. However, the behavior of double emulsions is quite different for these two systems. As such, the core droplet in the system IV deforms into an “S” shape whereas the “S” shape like deformation of the core in the system V is insignificant. Accordingly, after the pinch-off, the shapes of two daughter droplets are different for these two systems. As shown in Fig.11, the daughter droplets in the sys-tem IV are “onion-shaped” whereas they are “bean-shaped” in the system V. It should be noted that in the presence of shear flow, two daughter droplets exhibit radial symmetry, while, in the absence of shear flow, there is a reflective sym-metry between them. After a complete breakup, two daughter droplets move away from each other due to the shear forces and each act as a single emulsion. As the time progress, the surface tension force removes the sharp edges of the daughter droplets, thus leading to a smoother interface shape.

6 Conclusion

In this paper, the behavior of a double emulsion is studied under the combined effect of the electric field and shear flow using the multi-phase ISPH method. Six different systems are chosen based on electrical properties. Deformation and orientation angles of the core and the shell droplets are calcu-lated under pure shear, pure EHD, and combined EHD-shear flow. Flow streamline patterns are utilized to discuss the results. It is shown that the flow vortex direction depends on the conductivity and permittivity ratios. Electric force components are calculated on the interfaces and their effects on the double emulsion dynamics are discussed comprehen-sively. In addition to the EHD forces, hydrodynamics forces are also computed within the the entire flow domain. It is pointed out that in the double emulsion with a viscosity ratio of unity, pressure has a significant contribution to hydro-dynamic forces. The effects of capillary, electrical capillary numbers, and core to shell radius ratio on the double emul-sion dynamics are investigated as well. It is indicated that the deformation of core/shell droplet increases/decreases with an increase in radius ratio. Double emulsion exhibits relatively complex behavior in response to the variation of capillary and electrical capillary numbers which has been explained con-sidering the interaction of hydrodynamics and EHD forces. Finally, the breakup phenomenon is elaborated in double emulsions and the necessary condition for the breakup is

explained. Different breakup patterns are identified in the double emulsion for the specific electrical properties.

References

1. Pontrelli G, Carr EJ, Tiribocchi A, Succi S (2020) Modeling drug delivery from multiple emulsions. Phys Rev E 102:023114 2. Paximada P, Howarth M, Dubey BN (2021) Double emulsions

for-tified with plant and milk proteins as fat replacers in cheese. J Food Eng 288:110229

3. Tarnowska M, Briançon S, Resende de Azevedo J, Chevalier Y, Arquier D, Barratier C, Bolzinger M-A (2020) The effect of vehicle on skin absorption of Mg2+and Ca2+from thermal spring water. Int J Cosmet Sci 42(3):248–258

4. Raghuraman B, Tirmizi N, Wiencek J (1994) Emulsion liquid mem-branes for wastewater treatment: Equilibrium models for some typical metal-extractant systems. Environ Sci Technol 28(6):1090– 1098 PMID: 22176234

5. Taylor GI (1964) Disintegration of water drops in an electric field. Proc R Soc Lond Ser A Math Phys Sci 280(1382):383–397 6. Allan RS, Mason SG, Marion LE (1962) Particle behaviour in shear

and electric fields I. Deformation and burst of fluid drops. Proc R Soc Lond Ser A Math Phys Sci 267(1328):45–61

7. Taylor GI, McEwan AD, de Jong LNJ (1966) Studies in electrohy-drodynamics. I. The circulation produced in a drop by an electric field. Proc R Soc Lond Ser A Math Phys Sci 291(1425):159–166 8. Lac E, Homsy GM (2007) Axisymmetric deformation and stability

of a viscous drop in a steady electric field. J Fluid Mech 590:239– 264

9. Marrivada Nanchara Reddy and Asghar Esmaeeli (2009) The EHD-driven fluid flow and deformation of a liquid jet by a transverse electric field. Int J Multiph Flow 35(11):1051–1065

10. Shadloo MS, Rahmat A, Yildiz M (2013) A smoothed particle hydrodynamics study on the electrohydrodynamic deformation of a droplet suspended in a neutrally buoyant Newtonian fluid. Comput Mech 52(3):693–707

11. Xu J-J, Shi W, Hu W-F, Huang J-J (2020) A level-set immersed interface method for simulating the electrohydrodynamics. J Com-put Phys 400:108956

12. de Bruijn RA (1989) Deformation and breakup of drops in simple shear flows. PhD thesis, TUE, Department of Applied Physics 13. Sheth KS, Pozrikidis C (1995) Effects of inertia on the deformation

of liquid drops in simple shear flow. Comput Fluids 24(2):101–119 14. Chinyoka T, Renardy YY, Renardy M, Khismatullin DB (2005) Two-dimensional study of drop deformation under simple shear for Oldroyd-B liquids. J Non-Newtonian Fluid Mech 130(1):45– 56

15. Mählmann S, Papageorgiou DT (2009) Numerical study of electric field effects on the deformation of two-dimensional liquid drops in simple shear flow at arbitrary Reynolds number. J Fluid Mech 626:367–393

16. Zainali A, Tofighi N, Shadloo MS, Yildiz M (2013) Numerical investigation of Newtonian and non-Newtonian multiphase flows using ISPH method. Comput Methods Appl Mech Eng 254:99–113 17. Amani A, Balcázar N, Castro J, Oliva A (2019) Numerical study of droplet deformation in shear flow using a conservative level-set method. Chem Eng Sci 207:153–171

18. Tsukada T, Mayama J, Sato M, Hozawa M (1997) Theoretical and experimental studies on the behavior of a compound drop under a uniform DC electric field. J Chem Eng Jpn 30(2):215–222 19. Spasic AM, Jovanovic JM, Manojlovic V, Jovanovic M (2016)

Breaking of double emulsions based on electrohydrodynamics principles. J Colloid Interface Sci 479:165–172

20. Santra S, Das S, Chakraborty S (2020) Electrically modulated dynamics of a compound droplet in a confined microfluidic envi-ronment. J Fluid Mech 882:A23

21. Abbasi MS, Song R, Kim H, Lee J (2019) Multimodal breakup of a double emulsion droplet under an electric field. Soft Matter 15:2292–2300

22. Hua H, Shin J, Kim J (2014) Dynamics of a compound droplet in shear flow. Int J Heat Fluid Flow 50:63–71

23. Chen Y, Liu X, Shi M (2013) Hydrodynamics of double emulsion droplet in shear flow. Appl Phys Lett 102(5):051609

24. Luo ZY, He L, Bai BF (2015) Deformation of spherical compound capsules in simple shear flow. J Fluid Mech 775:77–104 25. Almasi F, Shadloo MS, Hadjadj A, Ozbulut M, Tofighi N,

Yildiz M (2019) Numerical simulations of multi-phase electro-hydrodynamics flows using a simple incompressible smoothed particle hydrodynamics method. Comput Math Appl

26. Santra S, Mandal S, Chakraborty S (2019) Confinement effect on electrically induced dynamics of a droplet in shear flow. Phys Rev E 100:033101

27. Nath B, Borthakur MP, Biswas G (2020) Electric field induced dynamics of viscoplastic droplets in shear flow. Phys Fluids 32(9):092110

28. Santra S, Das S, Chakraborty S (2019) Electric field-induced pinch-off of a compound droplet in Poiseuille flow. Phys Fluids 31(6):062004

29. Santra S, Panigrahi DP, Das S, Chakraborty S (2020) Shape evo-lution of compound droplet in combined presence of electric field and extensional flow. Phys Rev Fluids 5:063602

30. Lee HM, Choi SB, Kim JH, Lee JS (2020) Interfacial behavior of surfactant-covered double emulsion in extensional flow. Phys Rev E 102:053104

31. Santra S, Jana A, Chakraborty S (2020) Electric field modulated deformation dynamics of a compound drop in the presence of con-fined shear flow. Phys Fluids 32(12):122006

32. Borthakur MP, Nath B, Biswas G (2021) Dynamics of a compound droplet under the combined influence of electric field and shear flow. Phys Rev Fluids 6:023603

33. Shadloo MS, Oger G, Le Touzé D (2016) Smoothed particle hydro-dynamics method for fluid flows, towards industrial applications: motivations, current state, and challenges. Comput Fluids 136:11– 34

34. Fürstenau JP, Weißenfels C, Wriggers P (2020) Free surface ten-sion in incompressible smoothed particle hydrodynamcis (ISPH). Comput Mech 65:487–502

35. Huang C, Liu MB (2020) Modeling hydrate-bearing sediment with a mixed smoothed particle hydrodynamics. Comput Mech 66:877– 891

36. Roghair I, Musterd M, van den Ende D, Kleijn C, Kreutzer M, Mugele F (2015) A numerical technique to simulate display pixels based on electrowetting. Microfluid Nanofluidics 19(2):465–482 37. Brackbill JU, Kothe DB, Zemach C (1992) A continuum method

for modeling surface tension. J Comput Phys 100(2):335–354 38. Castellanos A (1998) Electrohydrodynamics, vol 380. Springer,

Berlin

39. Eringen AC, Maugin GA (2012) Electrodynamics of continua II: fluids and complex media. Springer, Berlin

40. Saville DA (1997) Electrohydrodynamics: the Taylor–Melcher leaky dielectric model. Ann Rev Fluid Mech 29(1):27–64 41. Reich FA, Rickert W, Müller WH (2018) An investigation into

elec-tromagnetic force models: differences in global and local effects demonstrated by selected problems. Contin Mech Thermodyn 30:233–266

42. Das D, Saintillan D (2021) A three-dimensional small-deformation theory for electrohydrodynamics of dielectric drops. J Fluid Mech 914:A22

43. Saghatchi R, Ghazanfarian J, Gorji-Bandpy M (2014) Numerical simulation of water-entry and sedimentation of an elliptic cylinder using smoothed-particle hydrodynamics method. J Offshore Mech Arctic Eng 136(3):031801

44. Ozbulut M, Tofighi N, Goren O, Yildiz M (2017) Investigation of wave characteristics in oscillatory motion of partially filled rectan-gular tanks. J Fluids Eng 140(4):041204

45. Ghazanfarian J, Saghatchi R, Gorji-Bandpy M (2015) Turbulent fluid-structure interaction of water-entry/exit of a rotating circular cylinder using SPH method. Int J Mod Phys C 26(08):1550088 46. Ghazanfarian J, Saghatchi R, Gorji-Bandpy M (2016) SPH

sim-ulation of turbulent flow past a high-frequency in-line oscillating cylinder near free-surface. Int J Mod Phys C 27(12):1650152 47. Zainali A, Tofighi N, Shadloo MS, Yildiz M (2013) Numerical

investigation of Newtonian and non-Newtonian multiphase flows using ISPH method. Comput Methods Appl Mech Eng 254:99–113 48. Rook R, Yildiz M, Dost S (2007) Modeling transient heat transfer using SPH and implicit time integration. Numer Heat Transf Part B Fundam 51(1):1–23

49. Saghatchi R, Ghazanfarian J (2015) A novel SPH method for the solution of dual-phase-lag model with temperature-jump boundary condition in nanoscale. Appl Math Model 39(3):1063–1073 50. Ghazanfarian J, Saghatchi R, Patil DV (2015) Implementation of

smoothed-particle hydrodynamics for non-linear Pennes’ bioheat transfer equation. Appl Math Comput 259:21–31

51. Monaghan JJ, Lattanzio JC (1985) A refined particle method for astrophysical problems. Astron Astrophys 149(1):135–143 52. Monaghan JJ, Kocharyan A (1995) SPH simulation of multiphase

flow. Comput Phys Commun 87(1–2):225–235

53. Morris JP (2000) Simulating surface tension with smoothed particle hydrodynamics. Int J Numer Methods Fluids 33(3):333–353 54. Tofighi N, Yildiz M (2013) Numerical simulation of single droplet

dynamics in three-phase flows using ISPH. Comput Math Appl 66(4):525–536

55. Ozbulut M, Yildiz M, Goren O (2014) A numerical investigation into the correction algorithms for SPH method in modeling violent free surface flows. Int J Mech Sci 79:56–65

56. Ozbulut M, Ramezanzadeh S, Yildiz M, Goren O (2020) Modelling of wave generation in a numerical tank by SPH method. J Ocean Eng Mar Energy 6:2198–6452

57. Kolukisa DC, Ozbulut M, Pesman E, Yildiz M (2020) Develop-ment of computationally efficient augDevelop-mented Lagrangian SPH for incompressible flows and its quantitative comparison with WCSPH simulating flow past a circular cylinder. Int J Numer Methods Eng 121(18):4187–4207

58. Shadloo MS, Zainali A, Sadek SH, Yildiz M (2011) Improved incompressible smoothed particle hydrodynamics method for sim-ulating flow around bluff bodies. Comput Methods Appl Mech Eng 200(9):1008–1020

59. Tofighi N, Ozbulut M, Rahmat A, Feng JJ, Yildiz M (2015) An incompressible smoothed particle hydrodynamics method for the motion of rigid bodies in fluids. J Comput Phys 297:207–220 60. Yildiz M, Rook RA, Suleman A (2009) SPH with the multiple

boundary tangent method. Int J Numer Methods Eng 77(10):1416– 1438

61. Tofighi N, Yildiz M (2013) Numerical simulation of single droplet dynamics in three-phase flows using ISPH. Comput Math Appl 66(4):525–536

62. Weller HG, Tabor G, Jasak H, Fureby C (1998) A tensorial approach to computational continuum mechanics using object-oriented tech-niques. Comput Phys 12(6):620–631

63. Deshpande SS, Anumolu L, Trujillo MF (2012) Evaluating the performance of the two-phase flow solver interFoam. Comput Sci Discov 5(1):014016

64. Feng JQ, Scott TC (1996) A computational analysis of electrohy-drodynamics of a leaky dielectric drop in an electric field. J Fluid Mech 311:289–326

65. Zhang J, Kwok DY (2005) A 2d lattice Boltzmann study on electro-hydrodynamic drop deformation with the leaky dielectric theory. J Comput Phys 206(1):150–161

66. Song Y, Shum HC (2012) Monodisperse w/w/w double emul-sion induced by phase separation. Langmuir 28(33):12054–12059 (PMID: 22849828)

67. Opalski AS, Makuch K, Derzsi L, Garstecki P (2020) Split or slip– passive generation of monodisperse double emulsions with cores of varying viscosity in microfluidic tandem step emulsification sys-tem. RSC Adv 10:23058–23065

68. Tsukada T, Yamamoto Y, Katayama T, Hozawa M (1994) Effect of an electric field on the behavior of a drop moving in a quiescent liquid. J Chem Eng Jpn 27(5):662–666

69. Hunt TP, Issadore D, Brown KA, Lee H, Westervelt RM (2009) Integrated circuit/microfluidic chips for dielectric manipulation. Caister Academic Press, Wymondham

70. Saghatchi R, Rahmat A, Yildiz M (2020) Electrohydrodynamics of a droplet in a highly confined domain: a numerical study. Phys Fluids 32(12):123305

Publisher’s Note Springer Nature remains neutral with regard to