An immersed interface method for axisymmetric electrohydrodynamic simulations in Stokes flow

http://journals.cambridge.org/CPH

Additional services for Communications in Computational Physics:

Email alerts: Click here Subscriptions: Click here Commercial reprints: Click here Terms of use : Click here

An Immersed Interface Method for Axisymmetric Electrohydrodynamic Simulations in Stokes ow

H. Nganguia, Y.-N. Young, A. T. Layton, W.-F. Hu and M.-C. Lai

Communications in Computational Physics / Volume 18 / Issue 02 / August 2015, pp 429 - 449 DOI: 10.4208/cicp.171014.270315a, Published online: 30 July 2015

Link to this article: http://journals.cambridge.org/abstract_S1815240615000651

How to cite this article:

H. Nganguia, Y.-N. Young, A. T. Layton, W.-F. Hu and M.-C. Lai (2015). An Immersed Interface Method for Axisymmetric Electrohydrodynamic Simulations in Stokes ow. Communications in Computational Physics, 18, pp 429-449 doi:10.4208/

cicp.171014.270315a

Request Permissions : Click here

Downloaded from http://journals.cambridge.org/CPH, IP address: 178.251.81.110 on 31 Jul 2015

An Immersed Interface Method for Axisymmetric Electrohydrodynamic Simulations in Stokes flow

H. Nganguia^1,2, Y.-N. Young^2,∗, A. T. Layton³, W.-F. Hu⁴ and M.-C. Lai⁴

1Department of Biological Sciences, New Jersey Institute of Technology, Newark, NJ 07102, USA.

2Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102, USA.

3 Department of Mathematics, Duke University, Box 90320, Durham, NC 27708, USA.

4Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsingchu 300, Taiwan, R.O.C.

Received 17 October 2014; Accepted (in revised version) 27 March 2015

Abstract. A numerical scheme based on the immersed interface method (IIM) is devel- oped to simulate the dynamics of an axisymmetric viscous drop under an electric field.

In this work, the IIM is used to solve both the fluid velocity field and the electric po- tential field. Detailed numerical studies on the numerical scheme show a second-order convergence. Moreover, our numerical scheme is validated by the good agreement with previous analytical models [1, 31, 39], and numerical results from the boundary integral simulations [17]. Our method can be extended to Navier-Stokes fluid flow with nonlinear inertia effects.

AMS subject classifications: 76W05, 65M06, 65M12

Key words: Immersed interface method, least squares interpolation, electrohydrodynamics.

1 Introduction

Under a direct current (DC) electric field an electrohydrodynamic flow is induced both inside and around a viscous drop [28, 37, 39]. Depending on the electrical properties of the fluids the viscous drop can deform into a prolate (oblate) spheroid with its major axis parallel (perpendicular) to the electric field. Such electrohydrodynamics of a viscous drop has a wide range of applications in micro-fluidic systems [36, 45], from drop ma- nipulation by electro-wetting [6] to enhanced mixing inside and on a drop in a Stokes

∗Corresponding author. Email addresses: hn9@njit.edu (H. Nganguia), yyoung@njit.edu (Y.-N. Young), alayton@math.duke.edu(A. T. Layton), weifanhu.am95g@g2.nctu.edu.tw (W.-F. Hu),

mclai@math.nctu.edu.tw(M.-C. Lai)

http://www.global-sci.com/ 429 2015 Global-Science Pressc

flow [2, 11, 32, 44]. Several theoretical models and numerical schemes have been devel- oped to study the drop dynamics under a DC electric field. Theoretical studies are often restricted to nearly spherical [1, 39, 42] or spheroidal shapes [4, 31, 46]. Consequently these models can only capture the drop electrohydrodynamics up to a moderate elec- tric field, and cannot predict the extreme drop elongation, lobed steady drop shape and several break-up modes under a strong DC field observed in experiments [13, 14] and di- rect numerical simulations. Various numerical schemes have been developed to simulate electrohydrodynamic (EHD) flow with certain limitations. For example, the boundary integral method [7, 17, 38] cannot be easily extended to different configurations such as solving the full Navier-Stokes equations, or accounting for charge effects in the bulk fluid.

Sharp methods such as the ghost-fluid method [33, 35] are generally first-order accurate.

Tomar et al. [41] used the volume of fluid (VOF) method, where jumps in electrical and fluid properties across the drop interface are smoothed out in a transition region around the moving interface. A direct consequence of the smoothing is the reduction in the order of accuracy. The finite element method [10] is too computationally expensive for deform- ing/moving interfaces.

To investigate the electrodeformation of a viscous drop with inertia effects (finite Reynolds number), recently Hu et al. [15] developed a hybrid numerical scheme, using the immersed boundary (IB) method for the incompressible Navier-Stokes flow (finite Reynolds number) and the augmented immersed interface method (AIIM) for the elec- tric potential. The authors treated the electric Maxwell stress as an interfacial force to fa- cilitate a unified immersed boundary framework, and illustrated first-order convergence for the electric force calculations. In addition this code can efficiently simulate large drop deformation under both an electric field and a shear flow.

Based on results in [15], we have developed an immersed interface method (IIM) for axisymmetric electrohydrodynamics of a viscous drop in the present study. Both the electric and the fluid equations are solved within the immersed interface framework with axisymmetry. The fluid solver is based on an axisymmetric immersed interface algorithm [24]. The electric solver is an AIIM modified from [15] to account for the axisymmetry of the drop interface.

The IIM was initially developed to solve elliptic equations with discontinuous co- efficients and/or singular sources [20]. The basic idea behind the IIM is to modify the finite difference discretization of the discretized governing equations near irregular grid points by adding correction terms that depend on the boundary conditions at an inter- face. The IIM differs from the immersed boundary method (IB method) [34], another popular method for solving similar problems, in that the IIM does not regularize the sin- gular forces at the interface. Instead, it captures the interfacial discontinuity in a sharp manner, and yields a second-order convergence rate. While second-order IB methods exist [5, 12, 19, 30], the second-order accuracy is reduced to first order for problems that involve an infinitely sharp interface [29], such as the drop interface that we consider here.

Therefore in this work we use IIM with consistent order of accuracy (second-order) for both the electric field and the fluid flow. In addition, our IIM algorithm is developed for

an axisymmetric viscous drop in three dimensions.

The paper is organized as follows: The problem formulation and governing equations are presented in Section 2. The numerical methods for both the fluid and electric compo- nents are described in Section 3. The IIM to solve the electric potential in axisymmetric coordinates is developed in Section 3.1, and the implementation of the IIM for the flow field is presented in Section 3.2. Finally, the convergence analysis and comparison against analytical models (such as the spheroidal model [31]) and boundary integral numerical simulation results [17] are discussed in Section 4.

2 Problem formulation

A schematic diagram of the problem formulation is shown in Fig. 1: an initially spherical viscous drop under a DC field E0, whose direction is parallel to the axis of symmetry. The interface Γ separates the exterior fluid (superscript ‘+’ for the exterior domain Ω⁺) from the interior fluid (superscript ‘-’ for the interior domain Ω⁻).

¡^ï!"m^ï!"!^ï 1^ï

#$

K

¡^%!"m^%!"!^%

1^%

Figure 1: Sketch of the problem: A leaky dielectric viscous drop (in domain Ω⁻) immersed in another dielectric fluid (in domain Ω⁺), with an external electric field E₀ in the z direction. Subscript ‘+’ and ‘-’ denote the exterior and interior fluids, respectively.

Within each subdomain, the fluid flow is described by the incompressible Stokes equation

−∇p^j+µ^j∇²u^j=0, (2.1)

∇·u^j=0, (2.2)

where j=‘+^′for the exterior fluid and j=‘−^′ for the interior fluid. p, u= (u,w)and the parameter µ denote the pressure, flow field, and fluid viscosity, respectively. In the bulk the electric permittivity ε^jand conductivity σ^j may be different between the exterior (‘+’) and interior (‘-’) of the drop. The ratios of permittivities and conductivities are defined

as

εr=^ε

−

ε⁺≡Q, σr=^σ

−

σ⁺≡¹

R, (2.3)

where Q and R are notations used in [17]. For convenient comparison to boundary inte- gral results in [17], we use (Q, R) instead of (εr, σr) when we present our IIM simulation results in Section 4.

The drop interface Γ is given in parametric form X(s,t) = (ρ(s,t),Z(s,t)), where the parameter 0≤s≤2π. (ρ(s,t),Z(s,t)) are the polar coordinates of the interface Γ in the r and z directions, respectively. Under an external electric field, bulk charges neutral- ize instantaneously in the leaky dielectric framework [28], and consequently the electric potential φ satisfies the Laplace equation

∇·(ε∇φ) =0, (2.4)

with the far-field boundary condition∇φ⁺= −E₀ˆz, and the boundary conditions on the drop interface

JφK=0, Jσ∇φ·nK=^{d ˜q}

dt, (2.5)

where n is the unit outward normal, ˜q=_Jε∇φ·nK represents the surface charge density, and J·K denotes the jump between outside and inside quantities. For ”weakly” con- ducting dielectric fluids, the charging time tc≡ε^j/σ^j is much faster than the time scale tEHD=µ^j/ε^jE²₀ of the electrohydrodynamic flow in most applications: tc≪tEHD. There- fore the charge convection term on the right hand side of Eq. (2.5) is often ignored in the leaky dielectric framework [9], and the jump condition on the normal electric field reduces toJσ∇φ·nK=0.

The electric force FE is related to the Maxwell stress M as FE= ∇·M, and M is com- puted from the electric field E as

M=ε

 EE−¹

2E²I



, (2.6)

where I is the identity tensor. In the leaky dielectric formulation, the electric force is important only when there is a gradient or jump in the electrical conductivity σ and/or permittivity ε. As these electrical properties are assumed to be piecewise constant with a possible jump across the drop interface in our formulation (see Fig. 2), the electric force FEcan be treated as an interfacial force given by the jump in Maxwell stress in the normal direction, as done in the boundary integral method [27].

At the drop interface the stress balance gives

−p⁺+p⁻
n+_JT^hd_K·n+f=0, (2.7) where(T^hd)_ij≡µ(∂_iu_j+∂_ju_i)is the ij-th component of the viscous stress tensor, andJ·_K denotes the difference between exterior and interior fluids. The traction f consists of the

Ω⁻ ε⁻,σ⁻

BC 2 r

Ω⁺ ε⁺,σ⁺ BC 4 z

BC 1

BC 3 Equator Pole

Figure 2: Computational domain on the(r,z)-plane. See text for the boundary conditions. On the walls BC 1, 2, 3, and 4 denote the boundary conditions defined in Eqs. (3.1)-(3.4).

surface tension [8, 18] and electric force, and is given by

f=fγ+f_E=γ ρsZss−ρssZs

|X_s|³ +^Zs ρ



n+_JM·nK, (2.8)

where γ is the surface tension along the drop interface X= (ρ,Z), and |X_s| =pρ²_s+Z²_s, with the subscript ‘s’ denoting partial derivatives with respect to the Lagrangian param- eter s.

3 Numerical methods

The governing equations in Section 2 are cast in cylindrical coordinates(r,z)for our ax- isymmetric system (Fig. 2), and the velocity field u= (u,w). These equations are solved over the r≥0 half-plane and then extended to the left half-plane r<0 by symmetry. We denote ∆t as the time step size and tn=n∆t as the nth time-level with n=0,1,2,···. At a time level tn, the position of the drop interface is given by a set of(NB+1)Lagrangian markers, X= (ρⁿ_k,Zⁿ_k)for k=0,1,2,···, NB, with k=0 at the north pole and k=NB/2 at the south pole. The drop is immersed in the domain Ω=[−L, L]×[−L, L], discretized as ri=ih and zj=jh along the r- and z-axis, respectively, with i, j=−N,−N+1,···,0,···, N−1, N and the mesh width h=L/N.

Fig. 2 shows the computational domain. For the electric potential, Dirichlet boundary conditions are imposed at the top (BC4) and bottom (BC3),

φ⁺= −E0L

2 at z=L, φ⁺=E0L

2 at z= −L. (3.1)

Figure 3: Grid geometry at an irregular point(i, j)from [20].

Neumann boundary condition for the electric field is imposed along the axis of symmetry (r=0, BC1) and the side wall (r=L, BC2),

∂φ

∂r =0 at r=0 and r=L. (3.2)

For the Stokes equations, boundary conditions for pressure and velocity are specified at r=0 (BC1)

∂p

∂r =0, ∂w

∂r =0, u=0. (3.3)

Dirichlet boundary conditions for p, u, and w on the other three sides (BC2-BC4) of the domain are determined from boundary integrals [3]:

p=

Z Z

− ¹

4π f·∇^{ 1} r



dS, (3.4)

u=

Z Z 1 8πµ

δij

r +^x

i′x^′_j r³

!

fjdS(x^′). (3.5)

In the IIM, grid points are classified as regular or irregular points based on their prox- imity to the interface. A point(i, j)is irregular if its stencil consists of neighbors on either side of the interface, while a regular grid point has its four neighbors on the same side of the interface (see Fig. 3).

3.1 Augmented IIM for the electric potential with axial symmetry

The IIM for a Poisson problem has frequently been utilized to solve Eq. (2.4). For the electric potential we adopt the following general IIM formulation for cases where the net bulk electric charge density (denoted as c in the following equation) is not zero

∇·(ε∇φ) =c, (3.6)

with boundary conditions at Γ

JφK=α, JσφnK=β. (3.7)

In our formulation we consider ε to be piecewise constant, therefore Eq. (3.6) can be recast as ∆φ=c/ε in each domain. Moreover, we can use the AIIM (augmented immersed interface method) to solve for the potential, with JφnK=g as the augmented variable.

First developed as a fast algorithm for elliptic interface problems with discontinuous coefficients [25], the AIIM introduces an unknown function (the augmented variable) to the original set of interface boundary conditions. This augmented variable is solved simultaneously in addition to solution of the elliptic equation. Recently AIIM has been applied to various problems [16, 26], and various modifications have been made to solve the Navier-Stokes equations with piecewise constant fluid viscosity.

The augmented method we use for the electric potential is based on [15], and here we extend it from two-dimensions to three-dimensions with axisymmetry. We use a fixed staggered grid, where the grid points xij and potential φij are defined at the cell center.

For the electric potential, the system of Eqs. (3.6)-(3.7) becomes

∆φ=

 c⁻/ε⁻ in Ω⁻,

c⁺/ε⁺ in Ω⁺, (3.8)

with boundary conditions

JφK=α, JφnK=g, JσφnK=β. (3.9)

We use the standard second-order difference stencil to discretize the Poisson equation on either side of the interface. The finite difference discretization near the interface is modified by the correction term Cij

1

h² φi+1,j+φi−1,j+φi,j+1+φi,j−1−4φi,j

+ ¹

2hri φi+1,j−φi−1,j

+^Cij

h² =¯cij, (3.10) or in matrix form,

AΦ+EΨ=T. (3.11)

Here, Φ and Ψ correspond to φij andJφnKx^∗_k, respectively, and x^∗_k is the orthogonal pro- jection of x at the interface; A is the difference operator, E is the operator acting on the augmented variable, and T is the vector of known solution in the right-hand side. The correction term Cij is non-zero only at irregular grid points (see Fig. 3), and is obtained from performing a Taylor’s series expansion about an irregular point on either side of the interface:

Cij= [φ]_x^∗_k+d[φn]_x^∗_k+^d

2

2



[¯cij]−κ[φn]− ¹

|Xs|

∂

∂s

 1

|Xs|

∂[φ]

∂s



x^∗_k

, (3.12)

where d is the signed distance between the irregular grid point and its orthogonal projec- tion x^∗_k, and κ is the local curvature of the interface. We refer the reader to [15] for detailed derivation of the correction terms, Cij.

At the interface, the relation betweenJφnK and JσφnK is expressed by φ⁺_n+^σ

−

JσKJφnK=^JσφnK

JσK , for σ⁻>_σ⁺_{, (3.13)} or

φ⁻_n+^σ

+

JσKJφnK=^JσφnK

JσK , for σ⁺>_σ⁻_{. (3.14)} The choice of Eq. (3.13) or Eq. (3.14) depends on the value of R (=1/σr). This helps ensure that coefficients of σ^∓/JσKI are large enough, so that the matrix in Eq. (3.16) is diagonally dominant. We cast Eqs. (3.13)-(3.14) in matrix form, B^±Φ+σ^∓/JσKΨ=G, where B⁺and B⁻denote the matrices resulting from the one-sided difference approximation along the normal directions, respectively. The overall linear system becomes

"

A E

B^{± σ}_[σ^∓_]I

# Φ Ψ



=

 T G



. (3.15)

We solve Eq. (3.15) by first eliminating Φ to obtain



B^±A⁻¹E−^σ

∓

JσKI



Ψ=B^±(A⁻¹T)−G. (3.16) Eq. (3.16) is a(NB+1)×(NB+1)system. Also, it is one dimension less compared with Eq. (3.15). The electric potential Φ is obtained in three steps: (1) Solve the system AΦ^∗=T for Φ^∗; (2) Use GMRES to calculate Ψ from Eq. (3.16); (3) Compute Φ from AΦ=T−EΨ.

We use a fast Poisson solver in steps (1) and (3).

Once we know the potential, we can compute the jump in the Maxwell stress. We use the one-sided least square interpolation [25, 43] to calculate the interfacial potential and electrostatic field. Here, we briefly summarize the one-sided least square interpolation scheme, modified for our axisymmetric case, and the interested readers are referred to [43] for more details. Denoting q the number of grid points used for the interpolation, we derive a system of q(≥6)equations from the second-order Taylor expansion near the interface

φ(r,z)_q=S1+S2(r−ρ_k)_q+S3(z−Z_k)_q+S4(r−ρ_k)²_q+S5(z−Z_k)²_q+S6(r−ρ_k)_q(z−Z_k)_q. (3.17) The problem is cast in matrix form,

ζs=Y, (3.18)

where s= [S1,S2,···,S6]^T, Y=φ(r,z)_l, and ζ is the coefficient matrix of relative distances.

The unknown coefficients s is solved using the least square minimization method, where the matrix ζ is expressed as ζ=UΣV^∗ by Singular Value Decomposition (SVD).

Uq×qand V₆^∗_×6are unitary matrices, and

Σ=





















Σ_1,1 0 ... 0 0 Σ_2,2 ... 0 ... ... . .. ... 0 0 ... Σ_6,6

... ... ...

0 0 ... 0





















. (3.19)

Eq. (3.18) is then recast as UΣV^∗s=Y, and we obtain the solution s=VΣ⁻¹U^∗Y. From this calculation, the electric potential and field at(ρk,Zk)are immediately obtained:

φ(ρk,Zk) ≈S1, (3.20)

φr(ρ_k,Z_k) ≈S2, (3.21)

φz(ρ_k,Z_k) ≈S3. (3.22)

3.2 IIM for the fluid flow field with axial symmetry

For the fluid flow field we use the IIM and adopt the approach developed by Li et al. [24]

for the incompressible Stokes equations

−∇p+µ∆u+F=0, (3.23)

∇·u=0. (3.24)

Eqs. (3.23)-(3.24) are recast as a set of three Poisson equations for the Stokes pressure p, and for components of the Stokes velocity u and w. The divergence of Eq. (3.23) gives

∇²p= ∇·F, JpK=fn,

s∂p

∂n {

=¹ r

∂

∂s(fτr), (3.25) where fnand fτ are the components of the boundary force normal and tangential to the interface, respectively. In practice, the singular source term is eliminated from the right- hand side. The pressure satisfies the Laplace equation,

∇²p=0, (3.26)

with specified jump conditions. The standard second-order finite difference discretiza- tion gives

1

h² pi+1,j+pi−1,j+pi,j+1+pi,j−1−4pi,j

+ ¹

2hri pi+1,j−pi−1,j

=C_i,j^p, (3.27)

where the correction term C_i,j^p is determined using the same methodology as in [20]. We solve the resulting system of equations for p, using fast Fourier transforms (FFT). Note that C_i,j^p is non-zero only at irregular grid points.

We then compute the velocity field u= (u,w), which is governed by 1

µ∇p= ∇²u, (3.28)

with jump conditions

JuK=0, µ s∂u

∂n {

=fτsinθ, µ s∂w

∂n {

= −fτcosθ. (3.29) We discretize the r-component of velocity to get

1

h² ui+1,j+ui−1,j+ui,j+1+ui,j−1−4ui,j

+ ¹

2hri ui+1,j−ui−1,j

=πi,j+C^u_i,j, (3.30)

where πi,j= pi+1,j−pi−1,j
/(2hµ). A similar discretization is obtained for the z-component of u. As in the case with p, we obtain solutions for the velocity by FFT [22, 24].

The position of the drop interface is updated according to dX

dt (s,t) = (UΓ,WΓ), (3.31) where

UΓ=γ₁ui1,j1+γ2ui2,j2+γ3ui3,j3−C^U, (3.32) W^Γ=δ1wi1,j1+δ2wi2,j2+δ3wi3,j3−C^W, (3.33) result from interpolating the Eulerian velocity, u using a second-order interpolation scheme that incorporates the jumps across the boundary [21]. The coefficients γ’s and δ’s must be determined. In the case of U^Γ, we first choose the three closest points,(ri1,zj1),(ri2,zj2), and(ri3,zj3)to the point(ρk,Zk)on the interface. Then the coefficients γk are,

γ2=(zj1−Zk)(ri3−ri1)−(ri1−ρk)(zj3−zj1)

(ri2−ri1)(zj3−zj1)−(ri3−ri2)(zj2−zj1)^{, (3.34)} γ3=(zj3−zj1)(ri1−ρ_k)−(ri2−ri1)(zj1−Z_k)

(ri2−ri1)(zj3−zj1)−(ri3−ri2)(zj2−zj1)^{, (3.35)}

γ1=−(γ2+γ3), (3.36)

and the correction term

C^U= −(a2[u]+a4[ur]+a6[uz]), (3.37)

where aiare linear combinations of the γ_k. u is the r-component of velocity. We use a sim- ilar approach to determine the interpolated velocity W^Γfrom the z-component of veloc- ity, w. Finally, we advance the markers according to the second-order Adams-Bashforth method

Xⁿ⁺¹=Xⁿ+^∆t 3 2UΓⁿ−¹

2U_Γⁿ⁻¹



. (3.38)

3.3 Summary of the numerical algorithm

Prior to solving the governing equations, the system is non-dimensionalized using the following scale: x=r0x^∗, p=_r^γ

0p^∗, and u=Udu^∗. Here r0 is the initial radius of the drop, and Ud is the characteristic velocity. The dimensionless governing equations become (after dropping the^∗)

−∇p+Ca∆u+

Z _2π

0 (f_γ+CaEf_E)δ²(x−X(s,t))ds=0, (3.39)

∇·u=0, (3.40)

∇·(ε∇φ) =0, JφK=0, JσφnK=0, (3.41)

E= −∇φ, M=ε

 EE−¹

2E²I



, f_E=_JMK·n and fγ=^{ ρsZss}−ρssZs

|X_s|³ +^Zs ρ



n. (3.42) The electrohydrodynamics thus depends on the capillary number Ca=µU/γ, the elec- tric capillary number CaE=ε⁺E₀²r0/γ, and the two ratios (Q, R). The capillary number represents the ratio of viscous force to surface tension, and the electric capillary number reflects the strength of the electric field.

Given the position of the drop interface Xⁿat time tn, the numerical implementation of the equations proceeds as follows:

1. Compute the electric potential φⁿusing the augmented IIM described in Section 3.1.

With φⁿ known, use the one-sided least squares interpolation scheme to calculate the interfacial electric field from (3.21) and (3.22). Then use the interfacial electric field to calculate the interfacial electric force, F_Eⁿfrom the Maxwell stresses M⁺_E and M⁻_E.

2. Compute the interfacial tension force, fγ from (3.42).

3. Using the IIM described in Section 3.2, solve the sequence of three Poisson prob- lems (3.26)-(3.28), for the Stokes pressure and velocity components at time tⁿ⁺¹. 4. Interpolate the fluid velocity to the boundary markers, and update the boundary to

its new position Xⁿ⁺¹.

5. Shift the drop center to the origin(r,z) = (0,0)to off-set the drift of a viscous drop under a DC electric field [10].

As noted in [10] a drop between two parallel plate capacitor (with a fixed jump in electric potential) exhibits a constant drift towards one of the plates. The drift of a viscous drop under a steady electric field is cancelled by an added hydrostatic pressure [10].

Others take advantage of the symmetry with respect to the equatorial plane to reduce the computational domain by half [33,40] (obtaining the solution in only one quadrant of the space). In this work, the drop center is shifted back to the origin at every time step to avoid the drift.

4 Results

4.1 Convergence study

In this subsection, we present the convergence results of our numerical schemes. The computational domain is(r,z) ∈ [0,3]×[−3,3], and the mesh size h=3/N, where N is the grid number. The permittivity and conductivity ratios are(Q, R) = (50,0.04). We used a time step of ∆t=0.001, and ran the simulations from t=0 to t=0.2. The number NB of Lagrangian markers representing the interface scales as NB=N/2.

Table 1 shows the convergence results as the ratio of the solution errors to h². We use the solution computed with N=512 as the true solution, Yexact, and define E^∞= kYN− Yexactk∞/h²and E1= kYN−Yexactk₁/h². The results show that the ratio is approximately constant and non-increasing, which indicates second-order spatial accuracy (similar to results in [23, 24]). Another way to determine the accuracy of the numerical method is to define the ratio

ratio=log₂

ku2N−uNk^∞ ku4N−u2Nk∞



. (4.1)

A ratio around 2 implies second-order accuracy. Table 2 shows the convergence analysis.

Table 1: Mesh refinement results for the velocity components,(u,w, p), and electric potential, φ at grid points of the computational domain.

N E^∞(u) E^∞(w) E^∞(p) E^∞(φ) E₁(u) E₁(w) E₁(p) E₁(φ) 32 4.77×10² 1.59×10³ 2.3×10³ 4.19×10³ 7.27×10² 8.80×10² 1.33×10³ 3.79×10³ 64 3.22×10² 1.05×10³ 1.86×10³ 4.03×10³ 5.49×10² 5.7×10² 8.5×10² 3.69×10³ 128 1.86×10² 6.62×10² 1.41×10³ 2.3×10³ 2.38×10² 3.80×10² 5.58×10² 2.13×10³ 256 1.12×10² 1.7×10² 8.8×10² 1.16×10³ 2.34×10² 1.55×10² 4.43×10² 1.17×10³

Table 2: Mesh refinement results for the velocity components,(u,w, p), and interface position, X.

N ku_2N−u_Nk^∞ ratio kw_2N−w_Nk^∞ ratio kp_2N−p_Nk^∞ ratio kX_2N−X_Nk^∞ ratio

32 3.87×10⁻¹ - 1.29 - 2.00 - 1.99 -

64 6.75×10⁻² 2.52 2.16×10⁻¹ 2.58 3.68×10⁻¹ 2.45 7.66×10⁻⁴ 11.34 128 1.01×10⁻² 2.74 3.82×10⁻² 2.5 7.26×10⁻² 2.34 7.44×10⁻⁵ 3.36 256 1.71×10⁻³ 2.57 2.59×10⁻³ 3.88 1.34×10⁻² 2.44 1.06×10⁻⁵ 2.81

Table 3: Mesh refinement results for the electric potential φ at all grid points of the computational domain, and for the interfacial electric potential, φ⁻, φ⁺. Results for φ⁻, φ⁺ are obtained using the one-sided least squares interpolation.

N kφ_2N−φ_Nk^∞ ratio kφ⁻_2N−φ⁻_Nk^∞ ratio kφ_2N⁺ −φ⁺_Nk^∞ ratio

32 3.12 - 1.27×10⁻¹ - 1.25×10⁻¹ -

64 8.48×10⁻¹ 1.88 3.38×10⁻² 1.91 3.51×10⁻² 1.84 128 1.29×10⁻¹ 2.72 5.35×10⁻³ 2.66 5.55×10⁻³ 2.66 256 1.77×10⁻² 2.86 7.12×10⁻⁴ 2.91 7.39×10⁻⁴ 2.91

For the interface position, X= (ρ,Z), the errors and rates are averages between compo- nents ρ and Z. One can immediately observe that the rate of convergence is second-order for u, w, p and the interface position, X.

Table 3 shows the convergence results for the the electric potential values at grid points, and on the interface. We also determine the accuracy of the interfacial electric force by plotting successive errors in the electric force as a function of the grid number, N. Fig. 4 shows the errors between successive grid refinement,k(F_E)_2N−(F_E)_Nk, for the interfacial electric force. We observe that the error convergence rate is second-order.

32 64 128 256

10⁻⁴ 10⁻³ 10⁻² 10⁻¹ 10⁰

Grid number, N

Successive error

r−direction z−direction

Figure 4: Successive errors for the electric interfacial force (F_E) in r (solid curve) and z-direction (dashed curve), as a function of the grid number N. The rate of convergence is nearly second order.

4.2 Equilibrium deformation as a function of conductivity ratioR

In this section we validate our model by comparing to analytical and numerical results in the literature. The equilibrium drop deformation under a DC electric field is quantified

by the deformation number

D=^L−B

L+B, (4.2)

where L and B are the drop size along the major and minor axes, respectively. Depending on the combination of electrical properties of fluids (Q and R, see appendix in [31] and references therein), the equilibrium drop shape under a DC electric field could be either prolate (the major axis is parallel to the electric field and D>0), spherical (D=0) or oblate (the major axis is orthogonal to the electric field and D<0). Moreover there exist two different modes for the prolate shape: one is prolate A with a circulation from the equator to the pole, and the other is prolate B with a circulation from the pole to the equator.

In the limit of small-deformation (|D| ≪1), a first-order approximation is given by Taylor [39]:

D≈ ^9CaE 16(2+1/R)²

 1

R²+1−2Q+3 1

R−Q  2+3µr

5+5µr



, (4.3)

where the ratio of viscosities µr=µ⁻/µ⁺. To validate our IIM solver, we perform simula- tions with varying conductivity ratio (R) and compare our equilibrium drop deformation D against models and and the VOF simulation results as follows.

We use parameters from Tomar et al. [41] with Q=10. The electric capillary num- ber (equivalent to the Bond number in [41]) is set at CaE=0.18. Moreover, we fixed N=128. The step size is h=5/N, and the time step ∆t=h/12.5 (for prolate) and ∆t=h/25 (for oblate). Fig. 5(a) shows a comparison of the numerical simulations with the predic- tions from Taylor’s first-order theory from Eq. (4.3), the spheroidal model of Nganguia et al. [31], and results using the volume of fluid (VOF) method in [41]. We note a transition from prolate (D>0) to oblate (D<0) shapes with increasing R. Results from our IIM sim- ulations are in good agreement with Eq. (4.3) for a small range of R, whereas we observe

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

−0.25

−0.2

−0.15

−0.1

−0.05 0 0.05 0.1 0.15

R

D

Taylor Spheroid IIM VOF (a)

−1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5 0 0.5 1

−1

−0.5 0 0.5 1 (b)

Figure 5: (a) Comparison of the steady state deformation (D) versus conductivity ratio R between Taylor’s approximation (Eq. (4.3), solid), the spheroidal model (dashed), the IIM, and VOF results [41]. Other parameters are: Q=10, µr=1, and Ca_E=0.18. (b) Steady state contour plots of the electric potential (left) and circulation plots (right) with(Q,R) = (10,1/1.81)and Ca_E=0.18. The dashed curve represents the initially spherical drop shape.

a larger range of R for agreement with our spheroidal model. Fig. 5(b) shows the electric potential (left) and velocity field (right) for R=0.55 (an oblate drop with a circulation from pole to equator), where our IIM result is in good agreement with VOF result and the spheroidal model.

4.3 Equilibrium deformation as a function of electric capillary numberCa_E We simulate drop deformation for various values of the permittivity and conductivity ratios, with parameter values from Fig. 16, Fig. 18 and Fig. 20 of [17]: (Q, R) = (0.1,0.1) for the prolate A drop, (Q, R) = (2,10) for the oblate drop, and (Q, R) = (50,0.04) for the prolate B drop. The spatial resolution is varied between N=128−512, with higher resolution needed at electric capillary number CaE≥0.25. The step size is h=6/N, and the time step ∆t=h/12.5 (for prolate) and ∆t=h/25 (for oblate). Starting with a spherical drop shape, the simulations are continued until a steady equilibrium shape is obtained.

Fig. 6 shows the equilibrium deformation number versus the electric Capillary num- ber CaEfor the prolate A and B, and oblate drops in panels (a), (b) and (c), respectively. In addition to results from our IIM simulations, we also superimpose the predictions from the spheroidal model [31], the boundary integral approach [17], and two other analytical models by Taylor [39] and Ajayi [1].

In Fig. 6(a), we observe that for small electric capillary number (CaE≤0.05) all the models and numerical results are in good agreement. This is expected, as all models

0 0.05 0.1 0.15 0.2 0.25 0.3

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35

CaE

D

IIM Spheroid Ajayi Taylor BIM (a)

0 0.05 0.1 0.15 0.2 0.25 0.3

0 0.05 0.1 0.15 0.2

CaE

D

IIM Spheroid Ajayi Taylor BIM (b)

0 0.05 0.1 0.15 0.2 0.25 0.3

−0.4

−0.35

−0.3

−0.25

−0.2

−0.15

−0.1

−0.05 0

CaE

D

IIM Spheroid Ajayi Taylor BIM (c)

Figure 6: Equilibrium drop deformation versus Ca_E from our IIM (solid), spheroidal model (dashed), Ajayi’s second-order approximation (dotted-dashed), Taylor’s results (dotted-solid) and the BIM (symbol). (a) Prolate A drop with (Q,R) =(0.1,0.1). (b) Prolate B drop with(Q, R)=(50,0.04). (c) Oblate drop with(Q, R)=(2,10).

−1 −0.5 0 0.5 1 0

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

0 0.5 1 (a) 1.5

−1 −0.5 0 0.5 1

0 0.5 1 1.5

0 0.5 1 (b) 1.5

−1.50 −1 −0.5 0 0.5 1 1.5

0.2 0.4 0.6 0.8 1

0 0.5 1 (c) 1.5

Figure 7: Steady state contour plots of the electric potential (left) and circulation plots (right) for (a) the prolate A drop with(Q, R) = (0.1,0.1) and CaE=0.325, (b) the prolate B drop with (Q, R) = (50,0.04) and Ca_E=0.305, and (c) the oblate drop with (Q, R) = (2,10) and Ca_E=0.304. The dashed curve represents the initially spherical drop shape.

should converge to Taylor’s result in the limit of small deformation. For CaE>_{0.05, Tay-} lor’s result start to deviate significantly from the other approximations. Ajayi’s second- order model is in good agreement with the spheroidal model and numerical simulations up to CaE≈0.1. The spheroidal model shows excellent agreement with the IIM and BIM simulations for the full range of CaE presented in this paper. We make similar observa- tions for the prolate B shape. In Fig. 6(b), Taylor’s result is in good agreement with all the other approximations for a larger range of CaE values, up to≈0.1. Ajayi’s second- order model agrees with the spheroidal and numerical simulations up to CaE≈0.2. For the oblate drop in Fig. 6(c), Taylor’s result is in agreements up to CaE≈0.1. Ajayi’s re- sult agrees for CaEslightly below 0.2, while the spheroidal model matches the numerical results up to CaE≈0.225.

For a ‘clean’ drop (no surface-active agents such as surfactant on the drop surface) in a DC electric field, Taylor predicted and observed that the circulatory motion is from the equator to the pole (counterclockwise in the first quadrant) for the prolate A drop [39].

On the other hand, a circulation from the pole to the equator (clockwise circulation in the first quadrant) is found inside the prolate B and oblate drops. Fig. 7 shows the steady state drop shape, electric potential, and flow patterns for prolate A and B, and oblate drops. For each panel, the electric potential contour is on the left, while the flow field is depicted on the right.

We also investigate large drop deformation and its dynamics leading to break-up un- der a strong electric field (large electric capillary number). For large CaE the axisym- metric drop shape deviates significantly from spheroids as illustrated in [17]. In partic-

−1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5 0 0.5

1 (a)

(b) (c) (d) (e) (a)

0 5 10 15 20 25 30

1 1.5 2

T

Position of equator

0 5 10 15 20 25 30

0 0.5 1

T

Position of pole

(b)

T = 32.3438

−1.5 −1 −0.5 0 0.5 1 1.5

0 0.5 1

0 0.5 (c) 1

Figure 8: Transient dynamics for the oblate drop in Fig. 6 with(Q, R) = (2,10) and Ca_E=0.315. Panel (a):

The dashed curve represents the initially spherical drop shape. Time T=1.4, 4.7, 18.7, 30, and 32.3 for curve (a), (b), (c), (d), and (e), respectively. Panel (b): Position of the equator (top) and pole (bottom) versus time.

Panel (c): Contour plots of the electric potential (left) and circulation plots (right).

ular, there may exist a critical capillary number above which the drop breaks up into smaller droplets. To illustrate that our IIM code captures non-spheroidal drop shape and dynamics leading to drop breakup for large CaE, we use the parameters for Fig. 7(c) ((Q, R) = (2,10)). Fixing CaE=0.315 and N=256, Fig. 8(a) shows the drop shape at times prior to the pinch-off. This sequence of drop shapes is in good agreement with Fig. 18 in [17], where the boundary integral results show that drop will break up for CaE>_0.297.

In addition, our IIM code (with N=256) is able to obtain steady spheroidal shapes up to CaE=0.304, which is slightly higher (by 2.4%) than the threshold 0.297 for break-up reported in [17]. Fig. 8(b) shows the location of the equator (pole) versus time in the top (bottom) panel. The late-time behavior of the pole indicates a pinch-off at the poles around T∼32. Fig. 8(c) depicts the electric potential contour (left) and flow field (right) at the time before pinch-off.

We also conducted simulations for a prolate A drop with CaE=0.37 and the same (Q, R) = (0.1,0.1)as those for the prolate A drop in Figs. 6(a) and 7(a). For this set of Q and R, the prolate A drop is predicted to pinch off for CaE≥0.34 [17]. Fig. 9(a) shows the drop shape at times prior to the pinch-off. The electric potential contour and flow field are depicted in Fig. 9(b). We note that, unlike the oblate case in Fig. 8, the prolate A drop takes longer to elongate and get close to pinch-off.

5 Conclusion

In this work we implemented the IIM to simulate the dynamics of a three-dimensional axisymmetric leaky dielectric drop under a DC electric field. Within the leaky dielectric

−2 0 2

−3

−2

−1 0 1 2 3

(a) (c) (b) (d)

(a) (e) T = 93.5156

−1 0 1

0 0.5 1 1.5 2 2.5 3 3.5 4

0 0.5 1 1.5 2 2.5 3 3.5 4 (b) 4.5

Figure 9: Transient dynamics for the prolate A drop in Fig. 6 with(Q, R) = (0.1,0.1)and Ca_E=0.37. Panel (a):

The dashed curve represents the initially spherical drop shape. Time T=3, 20.3, 50.8, 77.2, and 93.4 for curve (a), (b), (c), (d), and (e), respectively. Panel (b): Contour plots of the electric potential (left) and circulation plots (right).

framework with piecewise constant electrical properties, the electric force is treated as an interfacial force (calculated from the jump in Maxwell stress) along with the capillary force. In this work an IIM solver for the axisymmetric fluid flow [24] is combined with the AIIM solver for the electric field potential [15]. Following the approach in [15], the inter- facial electric force is computed from the interfacial Maxwell stress using the one-sided least square interpolation (see Section 3.1). We first performed numerical checks to con- firm the second-order convergence of the scheme for the electric potential (at grid points and on the interface), the electric force, and the fluid flow variables. This is a significant improvement from [15], which exhibits first-order convergence due to the immersed- boundary fluid solver. The second-order convergence in our IIM can be improved to higher-order by systematically calculating the higher-order corrections across the inter- face.

We then ran simulations with different combinations of permittivity and conductiv- ity ratios, and compared against analytical models and boundary integral simulation re- sults. Our results show excellent agreement for a wide range of electric capillary number.

Quantitative comparison against the boundary integral simulations from [17] with ex- actly the same physical parameters shows excellent agreement not only in equilibrium drop deformation and flow pattern, but also the non-equilibrium electrohydrodynamic of a viscous drop.

An obvious advantage of our approach over the boundary integral code is that our IIM code can be easily extended to electro-diffusion in an axisymmetric configuration. In addition, the interface can also be extended to an elastic interface (such as a membrane) for a capsule immersed in a fluid.

Another advantage of our fluid solver is that it can very easily be extended to com- pute solutions for the incompressible Navier-Stokes equations. In this case, the flow field