Vesicle electrohydrodynamic simulations by coupling immersed boundary and immersed interface method

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

Vesicle electrohydrodynamic simulations by coupling immersed boundary and immersed interface method

Wei-Fan Hu

^a,∗

, Ming-Chih Lai

^b

, Yunchang Seol

^c

, Yuan-Nan Young

^d

aDepartmentofAppliedMathematics,NationalChungHsingUniversity,145,XingdaRoad,TaichungCity402,Taiwan

bCenterofMathematicalModelingandScientificComputing&DepartmentofAppliedMathematics,NationalChiaoTungUniversity, 1001, Ta HsuehRoad,Hsinchu 300,Taiwan

cNationalCenterforTheoreticalSciences,No.1,Sec.4,Road.Roosevelt,NationalTaiwanUniversity,Taipei 10617,Taiwan dDepartmentofMathematicalSciences,NewJerseyInstituteofTechnology,Newark,NJ 07102,USA

a rt i c l e i n f o a b s t ra c t

Articlehistory:

Received21August2015

Receivedinrevisedform29February2016 Accepted16April2016

Availableonline19April2016

Keywords:

Immersedboundarymethod Immersedinterfacemethod Vesicleelectrohydrodynamics Leakydielectricmodel Navier–Stokesequations

In this paper, we develop a coupledimmersed boundary (IB) and immersed interface method(IIM)to simulatethe electrodeformationand electrohydrodynamics ofavesicle inNavier–StokesleakydielectricfluidsunderaDCelectricfield.Thevesiclemembraneis modeledas aninextensible elasticinterfacewithanelectriccapacitanceand anelectric conductance.Within theleaky dielectric framework and the piecewiseconstantelectric properties in each fluid, the electric stress can be treated as an interfacial force so thatboth themembrane electricand mechanical forces can beformulated in aunified immersed boundary method. The electric potential and transmembrane potential are solved simultaneously via an efficient immersed interface method. The fluid variables in Navier–Stokes equations are solved using a projection method on astaggered MAC gridwhilethe electricpotential issolved atthecell center. Aseriesofnumerical tests havebeencarefullyconductedto illustratetheaccuracy and applicabilityofthepresent methodtosimulatevesicleelectrohydrodynamics.Inparticular,weinvestigatetheprolate–

oblate–prolate(POP)transition andthe effectofelectricfieldand shearflow onvesicle electrohydrodynamics.Ournumericalresultsare ingoodagreementwiththoseobtained inpreviousworkusingdifferentnumericalalgorithms.

©^{2016ElsevierInc.}All rights reserved.

1. Introduction

Vesicleelectrohydrodynamics(dynamicsofthefluidflowaroundavesicleunderanexternalelectricfield)hasreceived much attentionduetoits wideapplicationsinbothbiological andbiomedicalsystemsthat involvevesiclemanipulations, such aselectroporation[16,34],creatingcrack[18],andelectrofusion[13].Acomprehensivereview formoreapplications canbefoundin[3,4].Inthispaper,ouraimistodevelopanumericalschemeforsimulationsofvesicleelectrodeformation;

inparticular,we focusonaleakydielectricvesiclesuspendedinanotherleakydielectric fluidunderauniformDC(direct current)electricfield.

Varioustheoreticalmodelshavebeendevelopedfortheelectrohydrodynamicsandelectrodeformationofanon-porated vesicle underan external electricfield. Fora DC electricfield, theoreticalmodels havebeen developedfora slightlyde-

*

Correspondingauthor.

E-mailaddresses:wfhu@nchu.edu.tw(W.-F. Hu),mclai@math.nctu.edu.tw(M.-C. Lai),ycseol@ntu.edu.tw(Y. Seol),yyoung@njit.edu(Y.-N. Young).

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

Fig. 1. A vesicle suspended in a fluid subject to a uniform DC electric field E_∞= (⁰,E_∞).

formedvesicle[10,11,16,37]andaspheroidalvesicle[35,46,47].TheextensiontoavesicleunderanAC(alternativecurrent) electricfield canbe found in [6,12,32,43]. These theoretical resultsare limitedby themodel assumptions such assmall deformationfromanearlysphericalvesicleandthevesiclecanonlytakeaspheroidalshape.

Inexperimentsthedependenceofthevesicleshapeanddynamicsonthephysicalparametersprovidesameanstomea- surethemembranebendingmodulus,membraneelectricalcapacitance,andmembraneconductance[2–4,39].Forexample, whenthevesicleinteriorislessconductivethantheexteriorfluid,thevesiclemightundergoaprolate–oblate–prolate(POP) transition[39].

Numericalmethods forsimulatingvesicle electrohydrodynamicsarefarandfew between.Focusing onStokesflow dy- namics,McConnellandMiksisdevelopedaboundary-integralnumericalschemeforsimulatingthetwo-dimensionalvesicle electrohydrodynamics[30].Theyfurtherextendedtheir worktostudythevesiclebreathingdynamicsunderanACelectric field [31]. Using a level-set/immersed interface hybrid approach, Kolahdouz andSalac proposed a numerical schemefor three-dimensionalvesicleelectrohydrodynamicsinNavier–Stokesflow[22]andusedittostudytheeffectsofdifferentfluid andmembranepropertiesonthePOPtransition[23].Theimmersedinterfacemethodisusedtosolvetheelectricpotential [21]andthelevelsetmethodisusedtosolveforthefluidflow[22].Inthepresentwork,ourhybridnumericalschemeuses animmersedinterfacemethodfortheelectricpotentialandanimmersedboundarymethodfortheNavier–Stokesflow.

The new numerical implementations and techniquesreported here for both the immersed interface method (for the electricpart)andthecouplingbetweenelectricandfluidflowviaimmersedboundarymethodarequitedifferentfromthe schemesin[22].Firstly,basedonourprevious workforsimulatingdropelectrohydrodynamics[15],we castthe Maxwell stress tensorasan electricinterfacial force (the jumpofMaxwell tensoracross theinterface)ratherthan avolume force intheequationssothatthevesiclemembraneforceandinterfacialelectricstresscanbeformulatedinaunifiedimmersed boundaryframework.Thisstepissignificantlyadvantageousbecauseitallowsustocomputetheelectricforce inan accu- ratemanner,andalsogreatly simplifiesthe numericalalgorithmforthe fullsystem. Secondly,we developan augmented immersedinterfacemethodforsolvingtheelectricpotential.Differentfromdropelectrohydrodynamics,theellipticinterface equationfortheelectricpotentialforthevesicleelectrohydrodynamic couplesatime-varyingtransmembranepotential.As aresult,thevesicletransmembranepotentialissolvedsimultaneouslywiththeelectricpotential.

Thirdly,wetreatthevesiclemembraneasanearlyinextensibleinterfacesothatsolvinganunknowntension(whichen- forcestheinextensibilityconstraint)canbeavoided.Furthermore,theFourierspectralmethodisusedtoevaluateinterfacial derivativessothatwe cancomputethegeometricalquantitiesontheinterface moreaccurately.Theconvergencestudyof inextensibleapproachcanbefoundinourpreviousworkforaxisymmetricvesicles[14] andanotherrecentworkforfully three-dimensionalvesicles [40].Lastly, asaforementioned,weconsider thevesicleelectrohydrodynamicsinNavier–Stokes flow asopposedtomostofexisting worksforStokes flow. Thisenablesusto investigatepossibleinertia effectsatfinite Reynoldsnumbers.

The paper is organized as follows. The formulation for the vesicle electrohydrodynamics within the leaky dielectric framework ispresentedinSection 2.An efficientimmersedinterface methodforsolvingthe electricandtransmembrane potentials is introduced andvalidated in Section 3. The new numerical algorithm that couples immersed boundary and immersed interface method for solving the vesicle electrohydrodynamicsis outlined in Section 4. A series of numerical simulationstoinvestigatethedynamicsofprolate–oblate–prolatetransitionandthestudyofhowtheelectricfieldaffects the vesicledynamics ina shearflow are presentedinSection 5. Some concludingremarksandfuture work aregivenin Section6.

2. Governingequationsofvesicleelectrohydrodynamics

Inthispaper,weconsideravesiclecontainingaleakydielectricfluid suspendedinanotherleakydielectricfluidunder a uniform DC electric field E_∞= (⁰,E_∞) asillustrated in Fig. 1. Separated by the elastic vesicle membrane  with a membranecapacitance C_m andconductanceG_m,thefluidsinside(domain⁻)andoutside(domain⁺)ofthevesicleare characterizedby electricconductivity

σ

andpermittivity

ε

.Here,we denote

σ

⁻ and

ε

⁻ asthe fluidelectricconductivity andpermittivity insidethe vesicle,and

σ

⁺ and

ε

⁺ forthe suspendingfluid.The conductivityandpermittivity ratiosare denotedas

σ

r=

σ

⁻/

σ

⁺and

ε

r=

ε

⁻/

ε

⁺,respectively.ThedynamicsofthevesicleunderaDCelectricfieldisdeterminedby theinteractionbetweenthehydrodynamics,elasticbending,tension,andelectricstressesonthemembrane.Forsimplicity, weassumethatboththefluiddensityandviscosityareidenticalfortheinteriorandsuspendingfluids.

The vesicle electrohydrodynamics is governed by the two-dimensional incompressible Navier–Stokes equations with an external DC electric field andforces on the vesicle membrane. The immersed boundary formulation is an Eulerian–

Lagrangianframeworkinwhichthefluidvariables(suchasvelocityfield,pressure,andelectricpotential)aredefinedinthe Eulerianmanner,whilethevesiclemembranerelatedvariables(suchasmembranepositionandtheelectricMaxwellstress) aredefinedintheLagrangianmanner.Forexample,thevesiclemembraneisrepresentedinaLagrangianparametricform asX(s,t)= (^X(s,t),Y(s,t)),0≤^s≤²

π

,wheres isaparameteroftheinitialconfigurationoftheimmersedboundary(vesi- clemembrane).Underthisrepresentation,thearc-lengthelement isdefinedas|^Xs|^ds=

X_s²+^Y2sds, wherethesubscript s denotesthepartialderivativewithrespectto s.Thus,thegoverningequationscanbewrittenasasinglefluidsystemin thedomain= ⁻∪ ⁺

ρ

 ∂

u

∂

t

+ (

^u

· ∇)

^u



= −∇

^p

+ μ

u

+

^{f in}

,

(1)

∇ ·

^u

=

^{0 in}

,

(2)

f

(

x

,

t

) =





(

F_γ

+

^Fb

+

^FE

)(

s

,

t

) δ(

x

−

^X

(

s

,

t

)) |

^Xs

|

^{ds in}

,

(3)

∂

X

∂

t

(

s

,

t

) =

U

(

s

,

t

) =





u

(

x

,

t

)δ(

x

−

X

(

s

,

t

))

dx on

,

(4)

∇

s

·

^U

=

^{0 on}

,

(5)

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

ρ

,the velocity u, the pressure p,andthedynamic viscosity

μ

.The externalforcing termf dueto immersedboundaryforceinEq.(3)consists of the membraneelastic tensionforce Fγ , thebending force Fb,and theelectric force FE.Eq.(4) simplystates that the membrane movesalong with thelocal fluid velocity (the interfacial velocity). The Eulerianfluid andthe Lagrangian im- mersed boundary variables are linked by thetwo-dimensional Diracdelta functionδ(x)= δ(^x)δ(y).Eq.(5) describesthe inextensibilityconstraintofthevesiclemembranewhichisequivalentlywrittenaszerosurfacedivergence(definedbelow) ofthevelocityalongthemembrane.Next,wedescribetheimmersedboundaryforcesinEq.(3)indetail.

A vesicle lipid bilayer membraneis inextensible and resistant to bending. Thus the total membrane elastic energy E consists of two parts; namely, a surface energy Eγ thatinvolves the membrane tension arising from the inextensibility constraint,andaHelfrichtypeenergyEb [9]toresistthebendingofmembrane.Thetotalmembraneenergyishence

E

(

t

) =

^Eγ

(

t

) +

^Eb

(

t

) =





 γ +

^cb

2

κ

²

^ |

^Xs

|

^ds

,

(6)

where

γ

is the membrane tension (which acts as a Lagrange multiplier) determined from the inextensibility constraint Eq.(5),c_b isthebendingmodulus,and

κ

isthemembranecurvature.Bytakingthevariationalderivativetothemembrane energyinEq.(6),thevesicletensionforceFγ andbendingforceF_bareobtainedas

F_γ

=

¹

|

^Xs

| γ

s

τ − γ κ

n

,

F_b

=

^cb



s

κ + κ

³

2



n

,

(7)

where

τ

= (

τ

1,

τ

2)=^Xs/|^Xs|^andn= (

τ

2,−

τ

1)aretheunittangentandnormalvectoralongthevesiclemembranerespec- tively.The surfacedivergenceoperator ∇s·^{inEq.(5)andthesurface Laplace(or}Laplace–Beltrami) operator s inEq.(7) aredefined,respectively,as

∇

s

· =

¹

|

^Xs

|

∂

∂

s

· τ

and

s

=

¹

|

^Xs

|

∂

∂

s



1

|

^Xs

|

∂

∂

s



.

(8)

ThedetailedderivationofEq.(7)canbefoundin[17].Noticethat,thesignofcurvature

κ

ispositiveforacircleunderthe currentparametricrepresentation.

The electricstressinadielectricmedium underan electricfield isgivenbytheMaxwell stresstensoroftheform[29, 36]

ME

= ε



EE

−

¹

2

(

E

·

^E

)

I



.

(9)

Since the permittivity

ε

and electricfield E are discontinuous across the membrane, we regard the electric effect as a membraneforcearisingfromthejumpoftheMaxwelltensorasinourpreviousworkforthedropletelectrohydrodynamics [15].Thatis,

F_E

= [

^ME

·

ⁿ

] = (

^M+_E

−

^M−_E

) ·

ⁿ

,

(10)

whereM⁻_E andM⁺_E standforMaxwellstresstensorinsidethevesicleandsurroundingfluid,respectively.

Wealsonotethat thesamemembraneelectricforce hasbeenusedintheboundaryintegralmethod[30,31] andlevel setmethod[22] forvesicleelectrohydrodynamicsimulations.TheremainingissueishowtofindtheelectricfieldE inthe domain,whichweaddressasfollows.

2.1. Vesicleleakydielectricmodel

Inthispaperweadopttheleakydielectricmodelwhereelectricchargesinthebulkrelaxonamuchshortertimescale thanthe characteristicfluid dynamictimescale [36,41],while thechargesare allowed toaccumulate onan interface that separatestwomediawithmismatchedelectricalproperties.Accordingly,theelectricfieldE isirrotational andcanberep- resentedbyE= −∇φ^{,wheretheelectricpotential}φsatisfiestheLaplaceequation

φ =

^{0 in}

 \ .

⁽¹¹⁾

Underan electricfield the lipidbilayermembranecharges asacapacitor.Opposite chargesaccumulate on eitherside ofthevesiclemembraneandthusformatransmembranepotential Vm (unlikethedropletinterfacebehavior, see[15] for example).Thatis,

[φ] = φ

⁺

− φ

⁻

=

^Vm

(

s

,

t

)

on

,

(12)

wherethe bracket[·]^{standsforthe jumpofthequantity} approaching fromthe⁺ side tothe⁻ side.Alsounderthe leakydielectricassumption,thenormalcomponentofohmic currentJ=

σ

E iscontinuous acrossthemembrane.Thus,we have

[ σ φ

n

] = σ

⁺

φ

_n⁺

− σ

⁻

φ

_n⁻

= −(

J⁺

−

J⁻

) ·

n

=

0 on

,

(13)

where the shorthand φn denotes the normal derivative ∇φ ·^n.The transmembranepotential Vm is calculated fromthe conservationlawofcurrentdensityacrossthemembrane[38]

Cm

∂

Vm

∂

t

+

GmVm

= σ

⁺

φ

_n⁺

= σ

⁻

φ

_n⁻ on

,

(14)

whereCmandGmarethemembranecapacitanceandconductance,respectively.Itisworthmentioningthatthetransmem- branepotentialisdefinedbythepotentialjumpfromthe⁻ sidetothe⁺sidesothereisasigndifferencebetweenthe presentV_m andtheoneusedin[38].Asaresult,therighthandsideofEq.(14)in[38]willbe

σ

⁺E⁺·ⁿ=

σ

⁻E⁻·^n.

3. AnimmersedinterfacemethodforthebulkelectricpotentialφandtransmembranepotentialVm

Theimmersedinterfacemethod(IIM)isajump-capturingmethodonaCartesiangridforsolvingellipticequationswith discontinuities across the interfaces,see Liand Ito[26] formore details. Inthis section, we will develop an augmented immersed interface methodfor solving the electric potential φ and the transmembrane potential V_m. That is, we need to solveEqs. (11)and(14)simultaneouslywith thejumpconditions Eqs.(12)–(13).Ourpresentnumericalschemeis an extensionofourpreviouselectricpotentialsolverforaviscousdrop[15],wherethereisnojumpinelectricpotentialacross thefluidinterface(noneedtosolveEq.(14)there).

Toproceed,letusfirstlayoutauniformCartesiangridinthecomputationaldomainwithmeshwidthh= ^x= ^{y for} simplicity.Thegridpointx_{i j} islocatedatthegridcenterwherethediscretepotentialφi j= φ(^xi j)isdefined.Theinterface

 embeddedincutsthroughsomegridcellssothepotentialsolutionisnotsmoothacrosstheinterface aswecansee fromthejumpconditionsEqs.(12)–(13).Wethusclassifythegridpointaseitheraregularorirregularpoint:Foraregular point,itmeansthatthestandardfive-pointLaplaciandiscretization(denotedby h) forEq.(11)doesnotcutthroughthe interface so the second-order local truncation erroris achieved. On the other hand, atan irregularpoint, the five-point Laplaciancutsthroughtheinterfacesothegridpointsusedinvolvebothinsideandoutsidetheinterface.Sincethesolution anditsderivativeshavejumps acrosstheinterface,a correctionfortheLaplaciandiscretizationisneededattheirregular pointtomaintainthedesiredaccuracy.Thus,thediscretizationofEq.(11)atthegridpointx_{i j} canbegenerallywrittenin theformof

_h

φ

_{i j}

= −

^{Ci j}

h²

,

(15)

whereC_{i j} isthecorrectiontermwhichisnonzeroonlyifthegridpointisirregular.

Letus describe whatthe correction term isat theparticular irregular point asdepictedin Fig. 2(a). When we apply thefive-pointLaplaciantoxi j,thegridpointxi−¹,jfallsindifferentsideoftheinterface sothecorrectionofdiscretization comes onlyfromthepoint xi−¹,j.Toderive thecorrection term, oneneeds tofind theorthogonalprojection of xi−¹,j at the interface (say X_∗=^X(s_∗) in Fig. 2(a)), andthen apply the Taylor’s expansion along the normal directionat X_∗.The correctiontermthusbecomes

Fig. 2. (a) The five-point Laplacian of the irregular point xi j. (b) A diagram showing least squares approximation for one-sided normal derivative at Xk.

C_{i j}

= [φ]

X_∗

+

^d

[φ

n

]

X_∗

+

^d2 2

 − κ

X_∗

[φ

n

]

X_∗

−

s

[φ]

X_∗

,

(16)

wherethevalued isthesigneddistancebetweenthegridpointxi−¹,jandtheorthogonalprojectionX_∗,and

κ

isthelocal curvatureoftheinterface.Themoredetailedderivationoftheabovecorrectiontermcanbefoundin[25,33].

Since theabovecorrection terminvolvestheusageofnormalderivative jump[φn]^{,itisquite naturalto usethejump} conditionEq.(13)toderivethefollowingequationasin[15,44]

φ

_n⁻

+ σ

⁺

[ σ ] [φ

n

] = [ σ φ

n

]

[ σ ] =

0

.

(17)

Here, weassumethat thevesicleislessconductivethanthesurroundingfluid(

σ

⁻<

σ

⁺). Otherformulacanbeusedfor the case(

σ

⁻>

σ

⁺) as shown in [15,44]. The terms in the righthand side of transmembrane potential Eq. (14) involve one-sidednormalderivatives

σ

⁺φ_n⁺and

σ

⁻φ_n⁻.Onecanfurtherrewritetheequationseparatelyby

C_m

σ

⁺

∂

V_m

∂

t

+

^Gm

σ

^+Vm

^{= φ}

n+

,

^Cm

σ

⁻

∂

V_m

∂

t

+

^Gm

σ

^−Vm

^{= φ}

n−

,

andthensubtracttheabovefirstequationbythesecondonetoobtain



1

−

¹

σ

r

 

C_m

σ

⁺

∂

V_m

∂

t

+

^Gm

σ

^+Vm



= [φ

n

]

^on

,

(18)

where the conductivityratio is defined by

σ

r=

σ

⁻/

σ

⁺. The advantage of using the above Eq. (18) instead ofEq. (14) is to establishthe relationshipbetween thetwo unknown interfacial functions Vm= [φ]^and [φn]^{.This new}formulation simplifiesournumericalalgorithmaswecanseelater.McConnelletal. exploitedthesimilarformulationfortheboundary integralmethodtosimulatevesicleelectrohydrodynamicsinStokesregime[30].

3.1. Implementationdetails

We now write down the discretization details for Eqs. (15)–(18), and solve the resultant matrix equation. Since the vesicle membrane is self-enclosing, we can use the Fourier spectral discretization to represent the configuration X. We firstchooseacollectionofdiscreteLagrangianmarkersXk=^X(sk)withequallydistributedinterfacialcoordinatessk=^k s, k=⁰,1,· · · ,^{M with} s=²

π

/M.ThemembranethuscanberepresentedintruncatedFourierseries(invectorform)as

X

(

s

) =

M

/2−¹ l=−M/2

X



_le^{i l s}

,

(19)

whereXl aretheFouriercoefficientsforX(s),andcanbe computedveryefficientlyusingtheFastFourierTransform(FFT).

Under thisFourierrepresentation,the derivativeswithrespect tos canalsobe computedquiteeasily by usingthepseu- dospectralmethod[42].Inthisfashion,allgeometricalquantitiessuchasthetangentvector

τ

(thus,thenormalvector n) andthecurvature

κ

canbecomputedwithspectralaccuracy.As notedbefore,all thevesiclemembranevariablesarede- finedinLagrangianmannersothejumpsof[φ]=^Vm and[φn]^{arealldefinedatX}k.ThankstotheFourierrepresentation, thejumps[φ]^and[φn]^{attheorthogonalprojectionpointsX}∗ inthecorrectionterminEq.(16)canbeeasilyinterpolated through the valuesat the markers X_k. Meanwhile, the curvature

κ

X_∗ and s[φ]X_∗ inEq. (16) can be also computed in a spectral manner.One shouldmentionthat the presentrepresentation fortheunknown derivative jump[φn] ^definedat the markerpoints X_k differs fromour previous work onthe dropletcase [15] inwhich [φn] ^{is defineddirectlyat those} projectionpointsX_∗.

Let,1and2bethesolutionvectorsformedbyφi j,[φ]Xkand[φn]Xk,thenthedifferenceEq.(15)canbealternatively writteninthematrixformof

A

 +

^C1



1

+

^C2



2

=

^F

,

(20)

where C1 and C2 are theformal matrix operators involvingthe Fourierinterpolation ofjumps inthe correction termas described in the previous subsection.The righthand side vector F simply results fromthe boundary conditions forthe potentialφ= −^E∞y.Tobecomplete,weneedtwomorematrixequationsfor1 and2,respectively.

As discussedearlier, thevalue of thenormal derivative jump[φn]Xk (denoted by the vector 2) can be linked tothe knownjump[

σ

φn]X_k throughEq.(17).However,aswecanseefromEq.(17),aone-sidednormalderivativeapproximation φ_n⁻(X_k)is needed. Toapproximatethat, we firstconstructa leastsquares cubicpolynomial by usingthe valuesof φi j at the grid points satisfying ^xi j−^Xk≤^{4h inthe domain} ⁻ as depictedin Fig. 2(b).In general, thenumber of chosen gridpoints ismorethanthe numberofunknown coefficientsintheleastsquaresapproximationsothe cubicpolynomial P(x,y) canbe obtained. Fortheinterface point Xk ina highcurvatureregion,we haveto increasethe spatialresolution to makesure enough gridpoints fall into thesame region sothe leastsquares cubicpolynomial is achievable.Once we have the polynomial P(x,y),the one-sided normal derivative φ_n⁻(X_k) can be computed by ∇^P(X_k)·ⁿ(X_k) directly. This least squarescubic polynomial approach roughly has third-order ofaccuracy to the approximation of the function itself while the derivative has second-order of accuracy. With this approximation, we can rewrite Eq. (17) in a matrix form as

B⁻

 + σ

⁺

[ σ ] 

2

=

⁰

,

(21)

whereB⁻denotestheresultantmatrixfromtheaboveleastsquarescubicpolynomialapproximationformally.

Nowweneedtoconstructthematrixequationfortheunknownvector1whichthetransmembranepotentialEq.(18) shouldbeused.Here,we providetheBackwardEuler(BE)andtheCrank–Nicholson(CN)methodforthetimeintegration asfollows.

• ^BE:



1

−

¹

σ

r

 

C_m

σ

⁺



1

− 

₁ⁿ

t

+

^Gm

σ

⁺

^

¹



= 

2

.

(22)

• ^CN:



1

−

¹

σ

r

 

Cm

σ

⁺



1

− 

₁ⁿ

t

+

^Gm

σ

⁺



1

+ 

₁ⁿ 2



= 

2

+ 

ⁿ₂

2

.

(23)

Here,₁ⁿandⁿ₂ aregivenintheprevioustimestepwiththetimestepsize t.

WecancouplethematrixEqs.(20)–(21)witheitherEq.(22)(BE)orEq.(23)(CN)intoonelinearsystemas

⎡

⎣

B^{A− C}0^{1 σ}_[^C_σ⁺_]²I 0

α

I

β

I

⎤

⎦

⎡

⎣ 



1



₂

⎤

⎦ =

⎡

⎣

^F0 G

⎤

⎦ ,

(24)

where

α

andβ are nonzero(if

σ

r=^{1)constantcoefficientsdepending ontheusageofEq.(22)or(23),andthevector G} collectstheknownterms.Inpractice,wedonotformthematricesC1,C2andB⁻explicitlyaswecanseefromtheiterative procedureofourmatrixsolverbelow.Tosolvetheabovelinearsystem,wefirsteliminateand2 fromEq.(24)bySchur complementtechnique(aftersomecarefulmatrixcalculations)toobtainanewlinearsystemof1 as



B⁻A⁻¹

(

C1

− α β

^C2

) + α

β σ

⁺

[ σ _]

^I





1

=

^B−A−1

(

F

−

¹

β

^C2G

) +

¹

β

σ

⁺

[ σ _]

^G

^.

⁽²⁵⁾

ThisisaM×^{M systemfor}1 whichismuchsmallerthantheoriginal one.WethenusetheGMRESiterativemethodto solvetheabovelinearsystem.Since theGMRES methodonlyrequiresthematrix-vectormultiplication,itisnotnecessary to construct the matrices A⁻¹, B⁻, C₁ and C₂ explicitly. Note that, the matrix A comes from the five-point Laplacian discretizationsothe inversionof A can beperformedefficientlyby applyingthefastPoissonsolverprovidedby Fishpack publicsoftwarepackage[1].

Insummary,thedetailednumericalalgorithmforsolvingthelinearsystemEq.(24)tofind,1 and2 canbesplit intothefollowingfoursteps.

Table 1

MeshrefinementresultsforaccuracyandefficiencyforExample 1.Thenumericalsolutionisde- notedbyφhandexactsolutionbyφe.

N BE CN

φh− φe∞ Rate Iter. φh− φe∞ Rate Iter.

80 8.046E–04 – 3 2.128E–04 – 3

160 4.958E–04 0.69 3 3.168E–05 2.74 4

320 2.587E–04 0.93 3 7.959E–06 1.99 4

640 1.313E–04 0.97 3 1.994E–06 2.00 4

Step 1. ApplyonefastPoissonsolvertoobtain^∗ in A



^∗

=

F

−

¹

β

^C2G

.

Step 2. UseGMRESiterationtosolve1 in



B⁻A⁻¹

(

C₁

− α β

^C2

) + α

β σ

⁺

[ σ ]

^I





1

=

^B−



^∗

+

¹

β

σ

⁺

[ σ ]

^G

.

Step 3. Update2 explicitlybyEq.(22)or(23).

Step 4. ApplyonefastPoissonsolvertosolvein A

 =

F

−

C1



1

−

C1



2

.

During the GMRES iteration in Step 2,we set the stopping criterion as h² which has the same magnitudeas the local truncation errorofour scheme.Theoverall computational cost forSteps1–4 in ourpresentscheme canbe evaluated in termsofthenumberoffastPoissonsolverbeingapplied.

Let usconclude thissubsection by remarking the comparison betweenthe presentnumerical algorithm and theone proposedbyKolahdouzandSalac[21]forsolvingtheaboveelectricandtransmembranepotentialequations.Althoughboth algorithms sharethesameIIMspirit developedearlierin[25,33],theauthorsin[21] usetheone-sidednormalderivative φ_n⁻astheunknownsinGMRESiterationswhilehereweusethenormalderivativejumpof[φn]^{astheunknownsdirectly.}

Meanwhile, in[21],thesurface Laplacianjump s[φ]^{attheprojectionpoint neededinthecorrection termofEq.(16)is} obtainedbyfirstextendingthejump[φ]=^Vmalongthenormaldirectiontothegridpoints viaclosestpointmethod[28]

withlevelsetrepresentationoftheinterfaceandthen evaluatingthestandardCartesian Laplacian [φ] ^(equalsto s[φ]⁾ onthegridpoints.However,thepresentFourierspectralrepresentationwithLagrangianmarkersisabletocalculate s[φ]

directly.Inthenextsubsection,weshallprovideamoresystematicnumericalcheckontheaccuracyandefficiencyforour presentnumericalalgorithm.

3.2. Numericaltest:accuracyandefficiencystudy

Inthissubsection,wedemonstratetheaccuracyandefficiencytestforthenumericalschemeoftheIIMforsolvingthe electricandtransmembranepotentialsdevelopedintheprevioussubsection.Allthenumericalcomputationswere carried outonadesktopPCwithIntelCorei7CPUand16 GRAM.

Example1.WeconstructananalyticalsolutionforEqs.(11)and(14)withthejumpconditionsEqs.(12)–(13)onacircular interface X(s)= (⁰.5cos s,0.5sin s)as

φ (

x

,

y

,

t

) =



_e−t

σ⁻

(

x²

−

^y2

)

x

∈ 

⁻

,

e^−t

σ⁺

(

x²

−

^y2

)

x

∈ 

⁺

,

Vm

(

s

,

t

) =



1

−

¹

σ

r



e^−t

σ

⁺

⁽

^X

⁽

^s

⁾

2

−

^Y

(

s

)

²

).

Thecomputationaldomainischosenas= [−¹,1]×[−¹,1]^.Theconductivitiesarechosenas

σ

⁻=⁰.1,

σ

⁺=^{1 (}

σ

r=⁰.1).

ThemembranecapacitanceandconductanceareCm=^{1 andG}m=₁₋⁴_σr+^Cm,respectively.WesetN tobethegridsizeand thus themeshsizeish=²/N. ThetimestepforthediscretizationoftransmembranepotentialEq.(22)or(23)ischosen inthesameorderofmeshwidthas t=h/4 andthesimulationisranuptoT=0.2.Wepresenttherateofconvergence forthemaximumerrorbetweenthenumericalsolutionφh andtheexactsolutionφe inTable 1.Asexpected,thefirst- and second-order convergenceareachievedforBEandCNscheme,respectively.Asforthecomputationalcomplexity,one can seethattheaveragenumberofGMRESiterationsinthetimeinterval[⁰,0.2]^{islessthan4andbecomessteadyevenasthe} gridnumberdoubles.

Fig. 3. Membrane shape X(s)= (^r(s)cos s,r(s)sin s)with the radius function r(s)=⁰.1(4+^{cos 3s}).

Table 2

MeshrefinementresultsforaccuracyandefficiencyforExample 2.Thenumericalsolutionisde- notedbyφhandexactsolutionisbyφe.

N BE CN

φh− φe∞ Rate Iter. φh− φe∞ Rate Iter.

80 8.299E–03 – 8 7.249E–04 – 9

160 4.033E–03 1.04 9 6.313E–05 3.52 10

320 2.017E–03 0.99 9 1.229E–05 2.36 9

640 1.010E–03 0.99 8 3.204E–06 1.93 8

Example2.Inthistest,we constructanotheranalytic solutionforEqs.(11)–(13)onamore complicatedinterface X(s)= (r(s)cos s,r(s)sin s)withtheradiusfunctionr(s)=⁰.1(4+^{cos 3s})(seeFig. 3)as

φ (

x

,

y

,

t

) =



_e−t

σ⁻e^xcos y x

∈ 

⁻

e^−t

σ⁺e^xcos y x

∈ 

⁺

,

Vm

(

s

,

t

) =



1

−

¹

σ

r



e^−t

σ

^+e

X(s)cos Y

(

s

).

Againthecomputationaldomainis= [−¹,1]× [−¹,1]^{.Weusethesame}conductivitiesasExample1butwithacapaci- tanceofCm=^{1 and}conductanceofGm=⁰.1.However,thetransmembranepotential Vm doesnotsatisfy Eq.(18)exactly so an extratermmust beadded inthe righthand sideof theequation. The gridsize andthetime step are thesameas theonesin Example 1andthesolution iscomputedup to T=⁰.2. Table 2showsthe meshrefinementanalysisforour numericalresults.Again,we obtainthe correspondingfirst andsecond-orderconvergenceforBEandCNscheme andthe averagenumberofGMRES iterationsstillstays steadyaswedoublethegridsize N.Thisnumericaltestdemonstratesthe robustcapabilityofourpresentalgorithmtohandleaninterfacewithmorecomplexgeometry.

4. Numericalalgorithm

Thegoverningequations(1)–(5),(7),(10),(11)and(14)arerenderednon-dimensionalasfollows.Weusetheeffective radiusofvesicle R=√

A0/

π

asascalinglength,where A0 istheenclosedareaofvesicle.Weusethemembranecharging time ofaquasi-spherical non-conductingmembrane[16]t_mm= ^RC_σ+^m

1 2+_σ¹_r

asour characteristictimescale.All physical variablesarescaledbytheassociatedcharacteristicscalesinthefollowing,

x^∗

=

^x

R

,

t^∗

=

^t

t_mm

,

p^∗

=

^t2mm

ρ

R²p

, γ

^∗

=

^R2

c_b

γ , ε

^∗

= ε

ε

⁺

^, σ

^∗

= σ σ

⁺

^,

^E

∗

=

^E E_∞

.

Thedimensionlessgoverningequationsforthevesicleelectrohydrodynamics(afterdropping*innotations)canbesumma- rizedas

∂

u

∂

t

+ (

u

· ∇)

u

= −∇

p

+

¹ Re

u

+







F_γ

+

¹

ReCaF_b

+

^Mn ReF_E



δ(

x

−

^X

(

s

,

t

)) |

^Xs

|

^{ds in}

,

⁽²⁶⁾

∇ ·

^u

=

^{0 in}

,

(27)

∂

X

∂

t

(

s

,

t

) =

^U

(

s

,

t

) =





u

(

x

,

t

)δ(

x

−

^X

(

s

,

t

))

dx on

,

(28)

where

γ = γ

0

(|

^Xs

| − |

^Xs

|

⁰

)

on

,

(29)

F_γ

= |

X¹s

| γ

s

τ ₋ γ κ

n

,

F_b

=



s

κ ₊ κ

³

2



n on

,

(30)

φ =

^{0 in}

 \ , [φ] =

^Vm

, [ σ φ

n

] =

^{0 on}

,

(31)



1

−

¹

σ

r

 

C

∂

Vm

∂

t

+

^{G V}m



= [φ

n

]

^on

,

(32)

E

= −∇φ,

^ME

= ε



EE

−

¹

2

(

E

·

^E

)

I



,

F_E

= (

^M+_E

−

^M−_E

) ·

ⁿ

.

(33)

Notice that, as our previous work on vesicle dynamics [14,40], we relax the inextensibility constraint Eq. (5) of vesicle membraneandreplacetheunknowntensionbyaspring-liketensionasshowninEq.(29)withalargeelasticstiffness

γ

0. Here |^Xs|^{0 isthelocalstretchingfactorofthevesicle initial}configuration.Thespring-liketensionisintendedtokeep the stretchingfactor|^Xs|^{closetoitsinitial}configurationsothecurrentvesiclemembraneisnearlyinextensible.Thenumerical validation andconvergenceofthisnearly inextensibleapproachhavebeenperformedinourrecentworkson3Daxisym- metric[14]andfully3Dvesicledynamics[40].

Thereareseveraldimensionlessnumbers;namely,theReynoldsnumber,Re=_μ^ρ_t^R_mm^{2 ;thecapillarynumber,Ca}=^μR_tmm^3/cb, measures the ratioof restoringbending andmembrane chargingtimescale; the Masonnumber, Mn= _μ/^t_ε^mm+E²_∞ measures thestrengthoftheelectricfield;thedimensionlessmembranecapacitance,C= ^RC_t^m_mm^{/σ+ measures themembranecharging} rate;thedimensionlessmembraneconductivity,G=^{R G}m/

σ

⁺ measurestheratioofchargingandmembraneconductivity rate. For 2D vesicle, another dimensionless parameter to measure the deflation of its shape is called the reduced area

ν

=⁴

π

A₀/L²₀,whereL₀ istheperimeterofvesicle.Thedimensionlessnumberisnothingbutthearearatioofthevesicle to a circle with the same perimeter. The typical values observed fromthe experimental data [31,37,39] are as follows:

the vesicle radius R≈^{10−5 m, the strength of electric field E}∞≈^{105 V/m, the fluid viscosity}

μ

_≈10⁻³ Pa s, the fluid density

ρ

≈^{103 kg/m3,thebulk fluid}conductivity

σ

⁺≈^{10−4 S/m,the} permittivity

ε

⁺≈^{10−10F/m, andthemembrane} conductivity Cm≈^{10−2 F/m2.Referring theabove physicalparameters andtaking}

σ

r=⁰.1,we obtainthe dimensionless numbers Re=⁰.02,Mn=^20,andC=⁰.1. Throughoutthispaper,we setthecapillarynumber Ca=10 correspondingto cb=¹⁰⁻¹⁷J (typicallycb=^{10−19J)sincelargerbendingforcecanstabilizethevesicleshapeinpresent}simulationsduring POPtransitionwithoutwrinklingasweshallexplainlater.

Next we describe the numerical scheme to simulate the full vesicle electrohydrodynamic systems Eqs. (26)–(33). As our dropletwork in[15],the ideaisto solvethe fluidequationsby theusual immersedboundary methodandtheelec- tric potential by the immersedinterface method proposed in Subsection 3.1. We consider the computational domain as a rectangle = [^a,b]× [^c,d]^{. Withinthis domain,a uniform Cartesian grid withmesh width h inboth x and y direc-} tions isemployed. The fluid variables are defined onthe standard staggered marker-and-cell(MAC) manner [8]. Thatis, the velocity component u and v are defined at the cell-normal edges (xi−¹/2,yj)= (^a+ (ⁱ−¹)h,c+ (^j−¹/2)h) and (xi,yj−¹/2)= (a+ (i−1/2)h,c+ (j−1)h)respectively,whilethepressurep andtheelectricpotentialφarebothdefinedat thecellcenter(xi,yj)= (^a+ (ⁱ−¹/2)h,c+ (^j−¹/2)h).Asmentionedbefore,forthevesiclemembraneinterface,weusea spectralcollectionpointssk=^k s,k=⁰,1,. . .M with s=²

π

/M torepresenttheLagrangianmarkers Xk=^X(sk)sothat any spatial derivativescan be performedspectrally accurate by using FastFourierTransform(FFT). In ourprevious work [14],it hasbeenconfirmedthat usingtheFourierspectral methodto compute interfacial derivativesindeedoutperforms thefinitedifferencemethodusedintraditionalIBmethodintermsofaccuracy.Thisaccuracyeffectisespeciallyprofound in the calculation of s

κ

inthe bending force F_b since the fourth-order derivative is involved. Here thepseudospectral methodisusedtoincrease theaccuracyofthecomputationofgeometricalquantitiesonvesicle membrane.However, the overallconvergencerateofthefluidvariableswouldremainonlyfirst-orderaccurateduetotheappliedimmersedboundary methodforthefluidsolver.

Thetime-steppingfortheoverallvesicleelectrohydrodynamicsystemcanbedescribedinthefollowing.Atthebeginning ofeachtimestepn,thefluidvelocityuⁿ,thevesiclemembraneconfigurationXⁿ,andthetransmembranepotential V_mⁿ must begiven.Thedetailednumericalalgorithmisgivenasfollows.

1. Solvetheelectricpotentialφⁿ andthetransmembranepotential(usingtheCrank–Nicholsonscheme)bytheimmersed interface method proposed inSubsection 3.1. Then we performthe one-sideddifference by leastsquarespolynomial approachintroducedinSubsection3.1tocomputeEⁿ= (−φⁿx,−φⁿy)attheLagrangianmarkersandusethemtocompute theMaxwellstresstensorM⁺_E andM⁻_E toobtaintheinterfacialelectricforceFⁿ_E (alsoatLagrangianmarkers).

2. For the given vesicle membrane configurationXⁿ, we compute the tensionforce Fⁿγ associated withthe spring-like tension

γ

ⁿ=

γ

0(|^Xs|ⁿ− |^Xs|⁰)andthebendingforceFⁿ_b inEq.(30).

3. Distributetheinterfacial forcetermsFⁿ_E,Fⁿγ ,andFⁿ_b fromtheLagrangian markerstothefluidgridpoints byusingthe discretedeltafunctionasintraditionalIBmethod.Here,thediscretedeltafunctionischosenasthesmoothed4-point piecewisefunctionin[45]whichhasbettersuppressionofnon-physicaloscillationsinthevolumeforcecomputations.