A hybrid immersed boundary and immersed interface method for electrohydrodynamic simulations

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

A hybrid immersed boundary and immersed interface method for electrohydrodynamic simulations

Wei-Fan Hu

^a

, Ming-Chih Lai

^b,∗

, Yuan-Nan Young

^c

aDepartmentofAppliedMathematics,NationalChiaoTungUniversity,1001,TaHsuehRoad,Hsinchu300,Taiwan

bCenterofMathematicalModelingandScientificComputing&DepartmentofAppliedMathematics,NationalChiaoTungUniversity, 1001,TaHsuehRoad,Hsinchu300,Taiwan

cDepartmentofMathematicalSciences,NewJerseyInstituteofTechnology,Newark,NJ07102,USA

a r t i c l e i n f o a b s t r a c t

Articlehistory:

Received28February2014

Receivedinrevisedform15August2014 Accepted5November2014

Availableonline11November2014

Keywords:

Immersedboundarymethod Immersedinterfacemethod Ellipticinterfaceproblem Electrohydrodynamics Leakydielectricmodel Navier–Stokesequations

Inthispaper,wedevelopahybridimmersedboundary(IB)andimmersedinterfacemethod (IIM)tosimulatethe dynamicsofadropunder anelectricfieldinNavier–Stokes flows.

Withintheleakydielectricframeworkwithpiecewiseconstantelectricpropertiesineach fluid,theelectricstresscanbetreatedasaninterfacialforceonthedropinterface.Thus, boththeelectricandcapillaryforcescanbeformulatedinaunifiedimmersedboundary framework.TheelectricpotentialsatisfiesaLaplaceequationwhichissolvednumerically by an augmented immersed interface method whichincorporates the jump conditions naturallyalongthenormaldirection.TheincompressibleNavier–Stokesequationsforthe fluidsare solvedusing aprojectionmethodonastaggeredMAC gridand the potential is solvedatthecell center.The interfaceistrackedin aLagrangian mannerwithmesh control byaddinganartificialtangentialvelocitytotransportthe Lagrangianmarkersto ensurethatthespacingbetweenmarkersis uniformthroughoutthecomputations.Aseries ofnumericaltestsforthepresentschemehavebeenconductedtoillustratetheaccuracy andapplicabilityofthemethod.Wefirstcomputethepotentialanditsgradient(electric field)toperformtheaccuracycheck forthepresentaugmentedIIM. Wethencheckthe convergenceoftheinterfacialelectricforceandthefluidvariables.Wefurtherrunaseries ofsimulations withdifferentpermittivityand conductivityratios and comparewiththe resultsobtainedbythesmalldeformationtheoryandothernumericalresultsinliterature.

Inaddition,wealsostudytheelectriceffectforadropundershearflow.

©^{2014ElsevierInc.}All rights reserved.

1. Introduction

Thestudyofhydrodynamicsdrivenbyanelectricfield(electrohydrodynamics)hasmanyindustrialapplicationsinmicro- fluidicsystems[22,33].Inparticular,aweaklyconducting(leakydielectric)dropsuspendedinanotherleakydielectricfluid underanelectricfield(asdepictedinFig. 1)hasbeenextensivelystudiedfromdifferentperspectives.G.I.Taylor[27] con- cluded that the equilibriumdrop shapecan be explained by balancing the viscous stress withthe electricstress on the dropinterfaceaslongashetookintoaccounttheinducedsurfacechargesduetothemismatchinelectricconductivityand permittivitybetweenthetwofluids.Furthermore,itisfoundthat dropcanbedeformedintoeitheraprolateoranoblate

*

Correspondingauthor.

E-mailaddresses:weifanhu.am95g@g2.nctu.edu.tw(W.-F. Hu),mclai@math.nctu.edu.tw(M.-C. Lai).

http://dx.doi.org/10.1016/j.jcp.2014.11.005 0021-9991/©^{2014ElsevierInc.}All rights reserved.

equilibriumshape withcirculatory flowsinside thedrop[20].Thedropdeformation andflowpatternsmainlydependon the electricalpropertiesofthefluidsystem. Acomprehensivereview offurther theoreticaldevelopmentscanbe foundin [24].Mosttheoretical models arelimitedto smalldeformation froma sphericaldropunder moderateelectricfields,and such limitation leads to quantitative discrepancies reportedin [29].Recently a spheroidal model was usedto study the large electrodeformationofa leakydielectricdrop(see[21] andreferencestherein).Resultsfromthesemodelsagreewell withexperimentsuptomoderateelectriccapillary number.However,thesespheroidalmodelsstillcannotpredictthedrop deformationunderlargeelectricfieldstrengths.

Thus,therehasbeenacontinuinginterestinusingnumericalmethodstosimulatetheelectrohydrodynamicsofaviscous drop under an electricfield. Numerous works on thenumerical simulations of dropdeformations underleaky dielectric theory arereportedintheliterature.Basedonhowtheinterface istreated,theseworkscanbe categorizedintothefront trackingmethod[11,30],levelsetmethod[4,31],phasefieldmethod[18,35],andthevolume-of-fluidmethod[9,32].Other numericalapproachesincludelatticeBoltzmannmethod[36]andboundaryintegralmethod[13,25],justtonameafew.

In this paper,a hybrid immersed boundary (IB)and immersedinterface method (IIM) are developed to simulate the dynamics ofaleaky dielectric dropunderan electricfield inNavier–Stokesfluids. The immersedboundarymethod used tosolvethefluidequationsissimilarinspirittothe front-trackingmethodasin[11,30].Otherrecentimmersedboundary methods to study particle–particleinteractions in electrohydrodynamicscan be found in [3,10], inparticular, the former one usedtheavailableimmersedinterfacemethodsolver[17]tosolvethepotentialequationbuttheapplicationsandthe electricforcecalculationsaredifferentfromthepresentstudy.

The majorcontributionsandsignificantdifferencesofourworkfrompreviousworks[4,3,11,10,18,30–32,35]areasfol- lows. Firstly, theMaxwell stress tensorarising fromtheelectriceffect iscastasan interfacial electricforce(the jump of Maxwell tensoracrosstheinterface)ratherthanthevolumeforce intheequationssothatthecapillaryandelectricinter- facialforcescanbeformulatedinaunifiedimmersedboundaryframework.Ascanbeseeninthenext section(Section2), thepresentapproachavoidsapplyingthedivergenceoperatortotheMaxwelltensor.Thisisparticularlyimportantsincethe tensorisdiscontinuousacrosstheinterfaceduetodifferentpermittivityandconductivity.

Secondly,asharpimmersedinterfacemethodisusedtocomputethepotentialanditsgradient.One-sidedinterpolations to the interface from either side are used to obtain the Maxwell stress tensor (thus the interfacial electric force) more accurately. In most aforementioned works in literature, the electric potential equation Eq. (6) is solved numerically by smoothing thepiecewise conductivity

σ

overa narrowregion usingtheharmonicmeans withoutanyspecialtreatments neartheinterface.However,accordingtothejumpconditionEq.(7),one canimmediatelyseethat∇φ ^{hasnonzerojump} acrosstheinterface.Thus,bytheabovesmoothingnumericalmethod,theelectricfieldE= −∇φ^{willhave O}(1)errornear theinterface.Furthermore,thesmoothingofpermittivity

ε

alsocauses O(1)errorneartheinterface.Thecomputationsof

∇

ε

and∇ · (

ε

E)bothhaveO(1/h)errorneartheinterface.Asaresult,thedirectcalculationofvolumeelectricforceEq.(9) isnotaccurate neartheinterface. Here,instead ofcomputingthevolume force,we compute theinterfacial electricforce in a moreaccurate manner.We are alsoable todemonstratethe numericalconvergenceof thisinterfacial force.A more rigoroustheoreticalanalysisontheaccuracyoffinitedifferenceschemesforellipticinterfaceproblemscanbefoundin[2].

The paper is organized as follows. The governing equations for the electrohydrodynamicswithin the leaky dielectric framework arepresented in Section 2.A simpleversion ofaugmented immersedinterface method forsolving piecewise ellipticinterface problemisintroducedandtestedinSection 3.Ahybridnumericalalgorithmofimmersedboundaryand immersedinterfacemethodsforsolvingtheelectrohydrodynamicsequationsisoutlinedinSection4.Thenumericalresults consistingof convergencecheckanddrop deformationsunderDCelectricfieldswithdifferentpermittivitiesandconduc- tivitiesare studiedindetail inSection5.Some concludingremarksandbriefdiscussion onfuturedirectionsare givenin Section6.

2. Governingequationsofelectrohydrodynamics

In this paper, we consider a leaky dielectric drop (Ω⁻) suspended in another immiscibleleaky dielectric fluid (Ω⁺) under a DC (direct current) electricfield E_∞ as depicted in Fig. 1. The governing equationsconsist oftwo-dimensional incompressibleNavier–Stokesequationswithsurfacetensionandelectricforces.We furtherassumethat boththedensity andviscosityareidenticalforthedropandsuspendedfluidasweconcentratemoreontheelectro-deformationofthedrop duetodifferentelectricconductivityandpermittivity.The dropisneutrallybuoyantinthefluiddomainΩ andthegravi- tational forceisneglected. Thefluidinterface Σ isrepresentedbyaLagrangian parametricformX(s,t)= (^X(s,t),Y(s,t)), 0≤^s≤²

π

,wheres istheparameteroftheinitialconfigurationoftheinterface.Undertheimmersedboundaryformulation, thistwo-fluidsystemiscastasasinglefluidwithvariablephysicalpropertiesinasingledomainΩ= Ω⁻∪ Ω⁺

ρ ( ∂

u

∂

t

+ (

^u

· ∇)

^u

) + ∇

^p

= μ

u

+

^fC

+

^fE in

Ω,

(1)

∇ ·

u

=

0 in

Ω,

(2)

∂

X

∂

t

(

s

,

t

) =

^U

(

s

,

t

) =



Ω

u

(

x

,

t

)δ 

x

−

^X

(

s

,

t

) 

dx on

Σ.

(3)

Fig. 1. A leaky dielectric drop suspended in another leaky dielectric fluid under DC electric field.

Eqs. (1)–(2) are the familiar incompressibleNavier–Stokes equations withthe fluid density

ρ

, the velocity u, the pres- sure p,thedynamicviscosity

μ

,thecapillaryforcearisingfromtheinterfacial tensionfC,andtheelectricforcefE.Eq.(3) simply statesthatthe interface movesalongwiththe localfluid velocity (theinterfacial velocity) whichislinked by the two-dimensionalDiracdeltafunctionδ(x)= δ(^x)δ(y).

Thecapillaryforceduetotensionontheinterfaceisimmersedas fC

(

x

,

t

) =



Σ

∂

∂

s

( γ τ ) δ 

x

−

X

(

s

,

t

) 

ds

,

(4)

where

γ

isthesurfacetension,and

τ

_=|^X_X^s

s| istheunittangentvectoralongtheinterface.Here,thesubscriptinXsdenotes the partialderivative of X withrespectto s. Incaseof asurfactant-covereddrop considered,the surface tension

γ

may not bea constantandit dependsonthesurfactant concentration viaLangmuirequation ofstate. Forsuchcase, an extra surfactantconvection–diffusionequationmustbeincorporatedintothefluidsystem.Oneoftheauthors(M.-C.Lai)andhis grouphavedevelopedanumericallyconservativeschemeforsolvingtheinsoluble[15]andsolublesurfactantequations[5].

TheelectricforcefE fortheleakydielectricmodelwillbegiveninthenextsection.

2.1. Leakydielectricmodel

Inelectrohydrodynamics, the dynamic current(withdimension I)is usually smallsothe inducedmagnetic effectcan be neglected.Thus, theelectricfieldintensityE(M LT⁻³I⁻¹) isirrotationalwhich resultsin∇ ×^E=^{0. ByGausslawina} dielectricmaterialwithpermittivity

ε

(M⁻¹L⁻³T⁴I²),thevolumechargedensityq_v(L⁻³T I)canbewrittenasq_v= ∇ · (

ε

E). Theconservationoffreechargedensityinelectrohydrodynamicscanbeexpressedas

Dqv

Dt

+ ∇ · ( σ

E

) =

⁰

,

(5)

where _Dt^D = _∂^∂_t +^u· ∇ ^{is the material derivative and}

σ

(M⁻¹L⁻³T³I²) is the electric conductivity of the medium. In a homogeneousincompressiblefluidwherethepermittivity

ε

andconductivity

σ

are bothconstant,onecanshow thatthe free chargedensitydecaysfromthe initialcharge densitywiththerelaxationtime scaletE=

ε

/

σ

[20].Theviscous time scaleofthefluidmotionisdefinedbyt_F =

ρ

L²/

μ

,where

ρ

and

μ

arefluiddensityandviscosity,andL isthecharacteristic lengthscale.

Inaleakydielectricmodel[20,24]foratwo-fluidsystemwithdifferentelectricpropertiesdepictedinFig. 1,thecharge accumulatesatthefluidinterfacealmostinstantaneouslyascomparedtothetimescaleofthefluidmotion(t_E^tF).Thus, thebulkchargedensityequation(5)canbesimplifiedinto

∇ · ( σ

E

) =

0

,

in

Ω \ Σ.

(6)

Withpiecewiseconstantoftheelectricconductivity

σ

⁻(inΩ⁻)and

σ

⁺(inΩ⁺)respectively,oneconcludesthat∇ ·^E=⁰ inbothdomainsΩ⁻andΩ⁺.SincetheelectricfieldE isirrotational(∇ ×^E=^{0)asmentionedearlier,theelectricfieldcan} beexpressedintermsofelectricpotentialφbyE= −∇φ^{.Therefore,Eq.(6)canbefurtherreducedtoLaplaceequationfor} thepotentialφineachfluiddomainΩ⁻andΩ⁺.Inthefarfield,aconstantelectricfieldE_∞ isappliedtothefluidsystem.

The boundary conditions along the interface Σ separating Ω⁻ and Ω⁺ are based on the continuity of the electric potentialandthenormalcomponentoftheelectricfluxdensityacrosstheinterface:

[φ] =

^{0 and}

[ σ _{∇φ ·}

n

] =

^{0 on}

Σ,

(7)

wherethejump[·]^{indicatesthequantityfromthe}Ω⁺sideminustheoneofΩ⁻side,andthenormalvectorn ispointing outwardfromΩ⁻toΩ⁺side.Thepermittivityofdropandthesuspendedfluidarepiecewiseconstantanddenotedby

ε

⁻

and

ε

⁺,respectively.AsketchoftheproblemsetupisdepictedinFig. 1.

Asdescribedin[20,24],thestressinducedinadielectricmediumunderan electricfieldisgivenbytheMaxwellstress tensoroftheform

M_E

= ε



EE

−

¹

2

(

E

·

^E

)

I



.

(8)

The equivalentvolumeforce densityoftheMaxwellstress canbeobtainedbytakingthedivergenceofaboveequationso wehave

f_E

= ∇ ·

^ME

= −

¹

2

(

E

·

^E

)∇ ε + ∇ · ( ε

E

)

E

.

(9)

The first term on the right-hand side of the above equation is due to polarization stress and it acts along the normal direction of the interface. The second term is due to the interaction of the electriccharges under the direction of the electric field.Since both thepermittivity andconductivityare piecewise constant,one immediatelysees that theelectric force inEq.(9)isnon-zeroonlyinthevicinity oftheinterface duetothefactof∇ ·^E=^{0 inbothdomains} Ω⁻ andΩ⁺. Thereforeforourmodelingpurpose,weregardtheelectriceffectasaninterfacialforcefromthejumpofMaxwellstressin thenormaldirectioninsteadofapplyingthevolumeforceasinEq.(9)intothefluidequations.Theinterfacialelectricforce isdefinedas

F_E

= [

^ME

·

ⁿ

] = 

M⁺_E

−

^M−_E



·

ⁿ

,

(10)

where M⁻_E andM⁺_E stand forMaxwell stress tensor inthe dropand thesuspended fluid, respectively. Thus, theelectric volumeforcefE canbealternativelyrepresentedbytheelectricinterfacialforceusingtheDiracdeltafunctionas

f_E

(

x

,

t

) =



Σ

F_E

(

s

,

t

)δ 

x

−

^X

(

s

,

t

) 

|

^Xs

|

^ds

.

(11)

ItisquiteinterestingtoconcludethattheinterfacialelectricforceinEq.(10)hasexactlythesameformasthecontinuous electricsurfaceforcederivedbyTomaretal.[32]inwhichtheauthorstakethelimitofvanishingtransitionregionthickness tobezero.TheexplicitformofthepresentinterfacialelectricforceinEq.(10)willbederivedindetailinAppendix A.Note also that, this kind of interfacial force approach is especially necessary to numerical simulations for boundary integral methodsuchasin[13]andhasbeenadoptedbyMcConnellandMiksisin[19]forvesicleelectrohydrodynamicsimulations.

3. Asimpleversionofimmersedinterfacemethodforsolvingellipticinterfaceproblem

In thissection, we propose asimple version ofimmersedinterface method forsolving theelectricpotential equation Eq.(6) ina domainΩ withjumpconditionsacross theinterface Σ (Eq.(7)).We considerthe problemona rectangular computational domain Ω= [^a,b]× [^c,d] ^{withan immersedinterface} Σ= {^X(s)= (^X(s),Y(s)),0≤^s<2

π

}^{,where s isa} Lagrangian parameter and X(0)=^X(2

π

). The interface Σ divides the domain Ω into two regions; namely inside (Ω⁻) and outside (Ω⁺) of the interface. Without loss of generality, we consider the followingelliptic interface problemwith inhomogeneousjumpsas

∇ · ( σ ∇φ) =

f in

Ω \ Σ, [φ] =

v

, [ σ ∇φ ·

n

] =

w

,

on

Σ.

(12)

Ofcourse, theabove equation shouldbe accompaniedwithsome suitable boundarycondition (say DirichletorNeumann condition)alongthecomputationaldomain∂Ω.Otherboundaryconditionswillnotchangethemainingredientspresented here. Notethat,thejumpsv(s)andw(s)arebothfunctionsdefinedontheinterfaceΣ.Since

σ

ispiecewiseconstant,we canrewritetheaboveequationintheformofPoissonequationφ= ^f/

σ

indifferentdomainsΩ⁺andΩ⁻withthesame jump conditions.In thefollowing,we use theshorthand φn torepresentthe normalderivative ∇φ ·^{n.Inorder toapply} theschemedevelopedin[23,14] directly,thejumpconditions[φ]^and[φn]^{mustbeknownandgiven. However,thejump} innormalderivative [φn] ^{isunknown.Wethus introduceanaugmented unknownvariable}[φn]=^{g intothe systemasin} [17] suchthatthesolutionsatisfiesthepiecewisejumpcondition [

σ

φn]=^{w.Insummary,Eq.(12)becomesthefollowing} ellipticinterfaceproblemas

φ =



f⁻

/ σ

⁻ in

Ω

⁻

f⁺

/ σ

⁺ in

Ω

⁺

[φ] =

^v

, [φ

n

] =

^g

, [ σ φ

n

] =

^w

,

on

Σ.

(13)

Aftersomesimplecalculations,therelationbetweenthejumpsof[φn]^and[

σ

φn]^{canbeexpressedbyeitheroneof}

φ

_n⁺

+ σ

⁻

[ σ ] [φ

n

] = [ σ φ

n

]

[ σ ] ,

if

σ

⁻

> σ

⁺

,

or (14)

φ

_n⁻

+ σ

⁺

[ σ _] ^[φ

ⁿ

^{] =}

[ σ φ

_n

]

[ σ _] ^,

^if

σ

⁺

> σ

⁻

.

(15)

TheabovechoiceisproposedbyXu[34]whichwillbecomeclearinournumericalschemediscussedlater.

Fig. 2. (a)Thegridpointsusedinfive-pointLaplacianattheirregularpointxi jmightcrosstheinterface.(b)Agridlayoutindicatingtheneighboringgrid pointsusedinleastsquarescubicpolynomialapproximationtotheone-sidednormalderivativeatX^∗_k.

3.1. AnaugmentedIIMtoincorporatethejumpsinthenormaldirection

Inthissubsection,weintroduceanaugmentedapproachforellipticinterfaceproblemsdescribedinprevioussubsection.

Ourtechnique issimilar to the one in [17] butthe computational details of the correction termsare quite different. To proceed,letusfirstlayoutauniformCartesiangridinΩ withmeshwidthh= ^x= ^{y.Thegrid pointsx}i j aredefined at thecell centerwhere the discrete solutions φi j are located.As in[17],we classify the grid point aseither a regular or irregular point. Fora regular grid point, we mean that the standard five-point Laplacian at that point doesnot cut through theinterface. Onthe other hand,ifthefive-point Laplacianofa grid pointinvolves usingthe gridpoints inside andoutsidetheinterface simultaneously, then wecall itan irregular point.Since thesolutionis smooth eitherinsideor outsidethedrop,thefive-pointLaplacianofaregularpointdoesnotneedtobemodifiedinordertohavethesecond-order accuracy.However,amodificationisneededatanirregularpointwherethesolutionisnotsmoothacrosstheinterface.The modificationdependsonthejumpconditions[φ] ^and[φn]^.Thus,thediscretizationofEq.(13)atthegridpointx_{i j} canbe generallywrittenintheformof

_h

φ

i j

+

^{Ci j}

h²

=

fi j

,

(16)

whereCi j isthecorrectiontermwhichisnonzeroonlyifthegridpointisirregular.

Letus briefly describe how the correction termis derived at thisparticular irregular point asillustrated in Fig. 2(a).

Whenweapplythefive-pointLaplaciantox_{i j},thegridpointx_i−1,jfallsintothedifferentsideoftheinterfacesothemain sourceofcorrectioncomesonlyfromthepointx_i−1,j.Toderivethecorrectiontermasin[14,23],wefirstneedtofindthe orthogonalprojectionofxi−¹,jontheinterface(sayX_k^∗=^X(s^∗_k)asshowninFig. 2(a)),andthenapplythetruncatedTaylor’s seriesexpansionalongthenormaldirectionatX^∗_k.Thecorrectiontermthusbecomes

C_{i j}

= [φ]

X^∗_k

+

^d

[φ

n

]

X_k^∗

+

^d2 2



[

^f

] − κ [φ

n

] −

¹

|

^Xs

|

∂

∂

s



1

|

^Xs

|

∂ [φ]

∂

s



X_k^∗

,

(17)

wherethevalued isthesigneddistancebetweenthegridpointx_i−1,janditsorthogonalprojectionX^∗_k,and

κ

isthelocal curvature ofthe interface. Notethat, fora fixed interface, the curvatureis known atthose orthogonal projection points.

However,whentheinterfaceisevolvingaspresentapplications,thecurvatureatthoseprojectionpointscanbe computed bythecubicsplineinterpolationoftheonesattheneighboringLagrangianmarkers.Weleavethedetailsonhowtocompute thelocalcurvatureatLagrangianmarkersinSection4.LetΦandΨ bethesolutionvectorsformedbyφi j and[φn]X^∗_k,then thedifferenceequationhasthematrixform AΦ+^EΨ=^{F .}

Notethat,inthecalculationofthecorrectiontermC_{i j},allthetermsarebasicallyavailableexceptthenormalderivative jump [φn]X^∗_k which can be linked to theknown jump [

σ

φn]X^∗_k through Eqs.(14) or(15). However, an approximation to φ_n⁺(X^∗_k) or φ_n⁻(X^∗_k) isneeded. Here, we use the leastsquarescubic polynomial approach to approximatethose values as follows.First,we choosetheneighboringgridpointsofX^∗_k= (^X∗_k,Y_k^∗)suchastheonesdepictedinFig. 2(b),andthen use the valuesofφi j at thosegrid points to constructthe leastsquarescubic polynomial P(x,y).Thus, the value ofφ_n⁺(X^∗_k) canbeapproximatedbycomputing∇^P(X_k^∗)·ⁿ(X^∗_k)directly.Thisleastsquarescubicpolynomialapproachhasthird-orderof accuracytotheapproximation ofthefunctionitselfwhilethe derivativehassecond-orderofaccuracy. Noticethat,inour presentnumericalexperiments,theaboveapproachtoapproximateφ_n⁺(X^∗_k)hasbetteraccuracythantheone-sidednormal extrapolationformulaproposedbyXuin[34].

Theapproximationforφ_n⁻canbedoneinasimilarmannersoweomitithere.Withthisapproximation,wecanrewrite Eqs.(14)–(15)inamatrixformasB^±Φ+^σ_[σ^∓_]Ψ=^{G,where B+andB−denotethematricesresultingfromtheaboveleast} squares cubic polynomial approximationalong the outward andinward normaldirections, respectively.Note that, we do not formthematrix B^± explicitlyinpracticesincetheGMRESiterationisappliedtosolvetheabovelinearsysteminthe following solution procedures. The vector G is the knownjumps at those orthogonal projectionsX^∗_k. Thus, theresultant linearsystembecomes



A E

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

Φ

Ψ 

=

F

G



.

(18)

Although theabove augmented techniqueis similar totheone developed in[34], thelinearsystemin Eq.(18)issolved differently.In[34],theauthorfirsteliminatedtheunknownΨ toobtaintheequationforΦ;andthensolvedtheresultant linearsystemby pseudotimeiterationto steadystate. Here,wefirst eliminateΦ fromEq.(18)toobtainalinearsystem ofΨ



B^±A⁻¹E

− σ

^∓

[ σ ]

^I



Ψ =

^B±



A⁻¹F



−

^G

.

(19)

Thisisan N_b×^Nb systemforΨ whereN_b isthenumberoftheprojectionpointsusedontheinterface.Theorderofabove linear systemisone-dimension lowerthan theone forsolvingΦ.Wethen usetheGMRES iterativemethodtosolve the abovelinearsystem.SincetheGMRESmethodonlyrequiresthematrix-vectormultiplication,itisnotnecessarytoconstruct thematrices A and B^± explicitly.IneachGMRESiteration,theinversionofA canbeperformedefficientlybyapplyingthe fastPoissonsolverprovidedbyFishpack[1]publicsoftwarepackage.Thedetailednumericalalgorithmcanbesplitintothe followingthreesteps.

Step 1. ApplyonefastPoissonsolvertosolveΦ^∗ in A

Φ

^∗

=

^F

.

Step 2. ApplyGMRESiterativemethodtosolveΨ untilitconvergesin



B^±A⁻¹E

− σ

^∓

[ σ ]

^I



Ψ =

^B±

Φ

^∗

−

^G

.

Step 3. ApplyonefastPoissonsolvertosolveΦ in A

Φ =

^F

−

^E

Ψ.

As mentioned before, the usage ofEq. (14)or Eq.(15)depends on theratio of

σ

r=

σ

⁻/

σ

⁺. If

σ

r>1, then we use Eq. (14)tocompute [φn]^{;otherwise, we useEq.(15).Onecan easily compute that}

σ

⁻/[

σ

_]=

σ

r/(1−

σ

r) and

σ

⁺/[

σ

_]=

1/(1−

σ

r).Thereasonoftheabove choice istoensure thatthecoefficient ofdiagonalmatrix σ^∓

[σ]I haslargermagnitude so thatthematrixtendstobemorediagonallydominant. Insuchway,theGMRESiterativemethodcanconvergequickly.

Here, we setthe stoppingcriteriafor GMRESmethod ash² since thecorrection termis nowwithlocal truncationerror O(h).Thus,theoverallcomputationalcostforSteps1–3inourpresentschemecanbeevaluatedintermsofthenumberof fastPoissonsolversbeingapplied.

3.2. Numericaltests

Example1(Convergenceandefficiencytest).Inthistest,wedemonstratetheaccuracyandefficiencyofthepresentaugmented IIMbychoosingagivenanalyticalsolutionforEq.(12)inΩ= [−¹,1]× [−¹,1]^as

φ (

x

,

y

) =



exp

(

x

+

^y

)

x

∈ Ω

⁻

,

sin x sin y x

∈ Ω

⁺

.

Theinterfaceischosenasanellipse(₀^x_.2)²+ (₀^y_.5)²=^{1 sothatthejump}discontinuities v and w inEq.(12)canbeexactly computedaswell.Forsimplicity,heretheDirichletboundaryconditionsareusedbecauseanexactsolutionisavailable.In Table 1,weshow themeshrefinementresultsforthesolutionanditsderivativeswithconductivityratio

σ

r=^{10 and 0}.1, respectively.Themeshwidthisdefinedash=²/N inbothx and y directionswiththegridsizeN.Notethat,thederivatives φx and φy are definedat thecell edges andcomputed by thestandard centered difference scheme usingthe computed solutions oftwo adjacent grid points. At the irregular point, a correction term has to be added asin [14]. Thiscan be donestraightforwardlyheresincethejump[φn]^{isalsosolvedinournumericalalgorithm.AsshowninTable 1,therateof} convergenceforthesolutionφisaroundsecond-orderwhilethederivativesφx andφy areslightlylessthansecond-order.

As forthe computational complexity,one can seefromthe last columnofTable 1 that thenumberof GMRESiterations plateaus to10 as thegridsize doubles.Allthe numericalcomputationswere carriedouton aPC withIntel Corei7CPU

Table 1

Theaccuracyandefficiencyresultsforthenumericalsolutionφhanditsderivatives(φh)xand(φh)yasthemeshisrefined.Here,φe istheexactsolution soitsderivatives(φe)xand(φe)ycanbecomputeddirectly.

N φh− φe∞ Rate (φh)x− (φe)x∞ Rate (φh)y− (φe)y∞ Rate Iter.

σr=¹⁰

32 3.31E−^{04 – 2.51E}−^{03 – 1.75E}−^{03 – 7}

64 8.46E−^{05 1.97 3.42E}−^{04 2.88 3.84E}−^{04 2.19 9}

128 2.43E−05 1.80 9.85E−05 1.80 1.13E−04 1.76 10

256 5.37E−06 2.18 3.84E−05 1.36 3.50E−05 1.69 10

σr=⁰.1

32 7.99E−^{04 – 5.07E}−^{03 – 5.16E}−^{03 – 6}

64 1.27E−^{04 2.65 5.21E}−^{04 3.28 6.10E}−^{04 3.08 6}

128 3.08E−^{05 2.04 1.91E}−^{04 1.45 2.52E}−^{04 1.28 7}

256 5.37E−06 2.52 2.85E−05 2.74 3.12E−05 3.01 7

Table 2

ThemeshrefinementresultsbetweenthepresentIIM(denotedbyφh)andthesmoothingmethod(denotedbyφ_h^S)withratioσr=⁰.1.

N φh− φe∞ φh^S− φe∞ (φh)x− (φe)x∞ (φh^S)x− (φe)x∞

32 3.07E−^{03 2.84E}−^{02 1.45E}−^{03 6.34E}−⁰¹

64 7.72E−^{04 8.57E}−^{03 4.44E}−^{04 5.91E}−⁰¹

128 1.93E−^{04 4.28E}−^{03 1.10E}−^{04 6.32E}−⁰¹

256 4.84E−^{05 2.45E}−^{03 2.85E}−^{05 6.80E}−⁰¹

and16GRAM.ThecomputationalcostforthepresentfastPoissonsolveris O(N²log N)withthegridsizeN inbothx and y directions.ThetotalCPUrunningtime foreach numericalresultlistedinTable 1iswithinasecond evenwiththegrid sizeN=^256.

Example2(Comparisonwiththesmoothingmethod).Here we consideran examplewithzero jumpconditions [φ]=^{0 and} [

σ

φn]=^{0,andcomparethenumericalresultswiththosefromthesmoothingmethod.Bythesmoothingmethod,wemeant} thatEq.(12)issolvedbyfirstsmoothing

σ

throughharmonicaveragingas

1

σ (

x

) =

^H

(

x

) σ

⁻

⁺

1

−

H

(

x

) σ

⁺

^,

andthenusingtheregular centereddifferenceschemeto discretizetheequation.Here, H(x)istheindicator function(or regularizedHeavisidefunction)whichcanbeobtainedbysolving

H

(

x

) = −∇ ·



Σ n

δ 

x

−

^X

(

s

) 

|

^Xs

|

^ds

.

Inthiscase, theinterface ischosen asa circlewithradius0.5immersedincomputational domain[−¹,1]× [−¹,1]^.The exactsolutionisgivenby

φ (

x

,

y

) =



a2

(

x²

+

y²

) +

a1 x

∈ Ω

⁻

,

(

x²

+

^y2

)

² x

∈ Ω

⁺

.

Theratioistakenby

σ

r=⁰.1.Theparametersa1 anda2 canbedeterminedbysettingthejumpconditions[φ]^and[

σ

φn] to be zeroacross the interface. From Table 2, one can see that the numerical solution fromthe smoothing method φ_h^S hasfirst-orderconvergencewhileitsderivativesdonotseemtoconvergeatall. Thisisnot surprisingsince thederivative is discontinuous across the interface inthe smoothing method.Both the solution (φh) andits derivative (φh)x from our presentIIMaremoreaccuratethanthosefromthesmoothingmethod.

4. Numericalalgorithm

Tosolvethewholeelectrohydrodynamicequationsnumerically,letusfirstperformthenon-dimensionalizationonthose governingequationsasin[11].Toproceed,wescaleallphysicalvariablesbytheassociatedcharacteristicscalesasfollows.

x^∗

=

^x

R

,

t^∗

=

 γ

ρ

R³t

,

u^∗

=

ρ

R

γ

^u

^,

^p

∗

=

^R

γ

^p

^, ε

^∗

= ε

ε

⁺

^, σ

^∗

= σ σ

⁺

^,

^E

∗

=

^E

|

E_∞

| .

Here, R is the initial drop radius. Aftersome simple calculations, the dimensionlessgoverning equations become (after dropping^∗innotations)

∂

u

∂

t

+ (

u

· ∇)

u

+ ∇

p

=

Oh

u

+



Σ



FC

+

CaEFE

|

Xs

|  δ 

x

−

X

(

s

,

t

) 

ds in

Ω,

(20)

∇ ·

^u

=

^{0 in}

Ω,

(21)

F_C

= ∂

∂

s

( γ τ )

on

Σ,

(22)

φ =

^{0 in}

Ω, [φ] =

⁰

, [ σ φ

n

] =

⁰

,

on

Σ,

(23)

E

= −∇φ,

^ME

= ε



EE

−

¹

2

(

E

·

^E

)

I



,

F_E

= 

M⁺_E

−

^M−_E



·

ⁿ

,

(24)

∂

X

∂

t

(

s

,

t

) =

^U

(

s

,

t

) =



Ω

u

(

x

,

t

)δ 

x

−

^X

(

s

,

t

) 

dx on

Σ.

(25)

There aretwodimensionlessnumbers;namely,theOhnesorgenumber,Oh=

μ

/√

γρ

R istheratiooftheviscousforce totheinertiaandsurfacetensionforce;andtheelectriccapillarynumber,Ca_E=

ε

⁺R|^E∞|²/

γ

measuresthestrengthofthe electricfieldrelativetothesurfacetensionforce.

Next,wedescribeahybridnumericalschemetosolvetheelectrohydrodynamicequationsEqs.(20)–(25).Theunderlying idea is that the fluid equations are solved by usual immersed boundary methodwhile the electricpotential equation is solved by theproposed immersedinterface methodinSection 3.Forsimplicity,letusassume thecomputational domain tobe arectangulardomainΩ= [^a,b]× [^c,d]^{.Withinthisdomain,auniformlatticegridwithmeshwidthh isemployed,} andthefluidvariablesaredefinedinastaggeredMACmanner[8].Thatis,thevelocitycomponentsu and v aredefinedat cell-normaledges(x_i−1/2,y_j)= (^a+(ⁱ−¹)h,c+(^j−¹/2)h)and(x_i,y_j−1/2)= (^a+(ⁱ−¹/2)h,c+(^j−¹)h),respectively,while thepressure p andtheelectricpotentialφarebothdefinedatthecellcenterx= (^xi,yj)= (^a+ (ⁱ−¹/2)h,c+ (^j−¹/2)h). SincethedropinterfaceX(s)= (^X(s),Y(s)),s∈ [⁰,2

π

_]isclosed(soperiodic),wecanusethespectralFourierdiscretiza- tiontorepresentX(s)andY(s).Wefirstchooseanevennumberofcollocationpointssuchthattheinterfaceisdiscretized by X(s_k)= (^X(s_k),Y(s_k)),wheres_k=^ks,k=⁰,1,. . .M with s=²

π

/M. The interfacecan berepresentedintruncated Fourierseriesas

X

(

s

) =

M



/2−1 l=−^M/2



X_le^ils

,

and Y

(

s

) =

M



/2−1 l=−^M/2

Y



_le^ils

,

(26)

where X_l and Y_l are the corresponding Fourier coefficients for X(s) and Y(s), respectively, and can be computed very efficiently usingthe FastFourierTransform(FFT). Under theFourierrepresentation,the derivativeswithrespect to s can be computedquiteeasily byusingthepseudospectralmethod[28].TheFouriercoefficientsof p-thderivativeof X(s) are simply computed as (il)^pX_l. Similar procedures can be applied forthe derivatives of Y(s). Inthis fashion, all interfacial quantities such asthetangent vector

τ

(thus, thenormalvector n) andthe curvature

κ

canbe computedwithspectral accuracy. As demonstrated in our recent work [12] for a three-dimensional axisymmetric interface, the presentspectral method forcomputingthecurvatureisindeedmore accuratethan thefinitedifference method(traditionally usedinIB).

In addition, this spectral representation ofthe interface is moreadvantageous when the surfactant is presentalong the interfacesincethesurfactantconvection-diffusionequationcanbesolvedmoreaccuratelyaswell.

In thefollowing,we describe howto advanceone time stepforthe solutionvariables. Atthe beginningofeach time stepn,theinterfacepositionXⁿ andthefluidvelocityuⁿmustbegiven.Thenumericalalgorithmisasfollows.

1. Computetheelectricpotential φⁿ inEq.(23)bytheaugmentedIIMdevelopedinSubsection3.1.Usethevaluesofφⁿ tocompute theelectricfieldEⁿ= (−φxⁿ,−φⁿy)onthegrid.(Attheirregularpoint,acorrectiontermmustbeaddedto achieve thedesiredaccuracy.)Thenweperformone-sidedinterpolationtocompute theMaxwell stressesM⁺_E andM⁻_E toobtaintheinterfacialelectricforceFⁿ_E atLagrangianmarkers.

2. GiventheinterfacemarkersXⁿ,wefirstcompute theunittangentvector

τ

ⁿ andthencompute theinterfacialtension forceFⁿ_C inEq.(22).Theabovebothtermsinvolvingspatialderivativeswithrespecttos canbecomputedwithspectral accuracy.

3. Distributetheinterfacial tensionandelectricforcefromthe Lagrangianmarkers tothefluid gridpoints by usingthe discretedeltafunctionasintraditionalIBmethod.

4. Solve the Navier–Stokes equations by the pressure-increment projection method to obtain new velocity uⁿ⁺¹. This procedure involves solvingtwo Helmholtz-type equationsfor theintermediate velocity andone Poissonequation for thepressure-incrementwhichagaincanbeefficientlydonebyapplyingthefastdirectsolverprovidedbyFishpack.

5. Interpolatethenewvelocityonthefluidgridpointstothemarkerpointsandthenmovethemarkerstonewpositions Xⁿ⁺¹asinEq.(25).

Here,thedetailednumericalimplementationofStep 1isdescribedinSubsection3.1whileSteps 2–5arequitestandardin immersedboundarymethodandcanbefoundinanyrelatedliteraturesuchasin[16,12].

Fig. 3. ElectricforceFE atdifferentnumberofmarkersM=⁶⁴,128,256 and 512.Leftpanelisthex componentandrightpanelisthe y component.

M=^{64 for}dash-dottedlines,M=^{128 fordottedlines,M}=^{256 fordashedlines,andM}=^{512 forsolidlines.}

Remark.Inthecontextofimmersedboundarysimulations,theinterfaceisusuallyevolvedbytrackingtheinterfacemarkers inaLagrangian manner.Oncetheinitialpositions oftheLagrangianmarkers arechosen,the movementofthosemarkers is determinedbytheinterpolatedlocalfluid velocity.Veryoften,astime evolves,theseLagrangianmarkersmaybeeither clusteredordispersedinsome regions.Forexample,ournumericalresultsinthenext sectionshow that ahomogeneous DCelectric field induces acirculatory flow inside theviscous leaky dielectric drop,moving fromthe equator towardthe tipsandthe Lagrangian markersare densely clusteredaround thetip regions. As aresult, the overall numericalstability andaccuracy maysuffer atlater times.To avoidsuch clustering of Lagrangian markers, we haveintroduced an artificial tangentialvelocityU^A(s,t)intoEq.(25)as

∂

X

(

s

,

t

)

∂

t

=

^U

(

s

,

t

) +

^UA

(

s

,

t

) τ ,

(27)

sothattheLagrangianmarkerscanbeuniformlydistributed.Here,thetermU^A(s,t)hastheformof

U^A

(

s

,

t

) =

^s 2

π

2π



0

∂

U

∂

s

· τ

ds

−



s

0

∂

U

∂

s

· τ

ds

.

(28)

Thedetailedderivationfortheabovetangentialvelocitycanbefoundin[16].

5. Numericalresults

Inthis section,we perform a seriesofnumerical testsforthepresentscheme developedin theprevious section. We firstverifytheconvergenceoftheinterfacialelectricforceandthefluidvariables.Then,werunaseriesofsimulationswith differentpermittivityratio

ε

r=

ε

⁻/

ε

⁺andconductivityratio

σ

r=

σ

⁻/

σ

⁺,andcomparetheresultswiththeonesobtained bythesmalldeformationtheory.Wealsostudytheelectriceffectforadropundershearflow.Finally,westudysomedrop deformationsbyslightlyincreasingtheelectriccapillarynumberCaE toseehowthedropbehavesandcomparewithother numericalresultsinliterature.Throughoutthispaper,weputacirculardropwiththeradiusR=^{1 insidethefluiddomain} Ω= [−⁴,4]× [−⁴,4] ^{andsettheinitial velocityto bezero}everywhere. Theflow issimplydriven bytheelectrical field.

Theboundary conditionforthevelocity field iszero(excepttheone fortheshearflow whichisindicated inFig. 7).The farelectricfieldischosenasE_∞= (⁰,−¹)sotheboundaryconditionsforthepotentialφare φ=y (Dirichlet)at y= ±^4, and ^∂φ_∂x =0 (Neumann)atx= ±^{4.TheCartesiangridmeshwidthh forthefluidvariablesandthepotentialisnowdefined} ash=⁸/N in both x and y directionswiththe grid size N. There are M Lagrangian markers distributedinitially along theinterface withtheLagrangianmesh s=²

π

/M<h.Here,wesimplychoose M=^{N inall}simulations. Thetimestep is chosen ast=^h/4.The Ohnesorge numberis chosen asOh=^{1, thesurface tension}

γ

=^{1 and the electriccapillary} numberCaE=⁰.5 unlessotherwisestated.

5.1. ConvergencetestfortheinterfacialelectricforceFE

As mentioned in Section 2.1, the electricfield effects are incorporated as an electric interfacial force FE in Eq. (10) insteadofthe moreconventional approachin[11,31]. Inthissimulation,we set

σ

r=^{3 and}

ε

r=^{2, andplotF}E inFig. 3 withdifferentnumberofLagrangianmarkersM=⁶⁴,128,256 and 512.Onecanseethattheinterfacialelectricforcetends toconvergeasthemarkersizeM increases.Inaddition,weplotthesuccessiveerrorsdefinedby(^FE)2M− (^FE)M∞versus themarkersizeM inFig. 4.Again,roughlysecond-orderconvergencehasbeenobservedinthistest.

Fig. 4. The successive errors of the electric interfacial force FEin both x and y directions (see legends). The rate of convergence is around second-order.

Table 3

Themeshrefinementresultsforthevelocitycomponentsu andv,theinterfaceconfigurationX,andtherelativevolumelossofthedropatT=^2.

N u2N−uN∞ Rate v2N−vN∞ Rate X2N−XN∞ Rate ^|AN_A^−A0|

0

64 1.19E−02 – 1.62E−02 – 3.01E−02 – 8.74E−05

128 2.84E−03 2.07 5.36E−03 1.60 8.80E−03 1.78 5.85E−05

256 7.77E−^{04 1.87 8.85E}−^{04 2.60 1.41E}−^{03 2.64 3.16E}−⁰⁵

5.2. Convergencetestforthefluidvariables

In thissubsection,we performtheconvergencetest forthepresentnumericalalgorithmofthewhole electrohydrody- namicequations.Weagainchoosetheconductivityratio

σ

r=^{3 and}permittivityratio

ε

r=^{2 withdifferentnumberofgrid} points N=⁶⁴,128,256,512.Allthenumericalsolutionsarecomputedup totime T=^{2.Sincetheanalyticalsolutionsare} notavailableinthistest,wecomputetheerrorsbetweentwosuccessivegridsdenotedby^u2N−^uN^{sotherateofconver-} genceiscalculatedbyrate=^log2^uN−^uN/2

u2N−uN.Theerrorsandratiosofv andtheinterfaceconfigurationX areallcomputedin thesamemanner.Table 3showsthemaximumerrorsforthefluidvelocityandtheinterfaceconfiguration.Notethat,the fluidvariablesaredefinedonthestaggeredgrid.Thus,whenwerefinethemesh,thenumericalsolutionswillnotcoincide withthesamegridlocationsandasimplebi-linearinterpolationmustbeimplemented.Duetothefactthattheimmersed boundaryformulationhasthesingular forcingtermintheequations,regularizingthesingulartermby smoothingdiscrete deltafunctionsmightcausethemethodtobefirst-orderaccurate.However,thepresentnumericalresultsshowbetterthan first-orderconvergence.

Meanwhile, we alsocheck the volume (area in2D) loss ofthe drop inour presentcomputations. From Table 3 (last column),oneseesthattherelativevolumelossisoforderofmagnitude10⁻⁵ whichisnegligiblysmallinpresentsimula- tions.Weattributesuchgoodvolume-conservationbehaviortothechoiceoffluidsolverandtheenforcementofuniformly distributedLagrangianmarkers.Veryrecently,BoyceGriffith[7]hasperformedratherdetailednumericaltestsonvolume- conservationfordifferentIBimplementations.Inourworkthepressure-incrementprojectionmethodunderstaggeredgrid discretizationisadoptedforourfluidsolverbecauseofitssuperiorperformanceforvolume-conservation[7].

5.3. Comparisonwithsmall-deformationtheory

In this subsection,we simulatethe drop deformation undera DC electricfield, andcompare withasymptotic results from the small-deformation theory.G.I. Taylor[27] developed a linear modelto compute steadyequilibrium shape ofa leakydielectricdropunderaDCelectricfield.Hefoundthattheviscousdropcanturnintoeitheroblateorprolateshape.

BasedonTaylor’sasymptoticresults,thedropdeformationdependsontheelectriccapillarynumberCaE,theratioofelectric conductivity,permittivity,andfluidviscosity.Theequilibriumdropdeformationcanbequantifiedbythedropdeformation number D as

D

=

^L

−

B L

+

B

,

where L and B arethedropdeformationsalongthemajorandminoraxes,respectively.Thefirst-ordersmall-deformation modelforatwo-dimensionalviscousdropisdevelopedby Feng[6](alsolistedin[36]) andisanalogoustoTaylor’slinear modelforasphericaldrop[27].Inthislinearmodeltheequilibriumdropdeformationisapproximatedby

D

=

^fd

( σ

r

, ε

r

)

3

(

1

+ σ

r

)

^2CaE

,

(29)