An immersed boundary method for simulating vesicle dynamics in three dimensions

Yunchang Seol

^a

, Wei-Fan Hu

^b

, Yongsam Kim

^c

, Ming-Chih Lai

^d,∗

aNationalCenterforTheoreticalSciences,No. 1,Sec. 4,RooseveltRoad, NationalTaiwanUniversity,Taipei10617,Taiwan bDepartmentofAppliedMathematics,NationalChungHsingUniversity,145,XingdaRoad,TaichungCity402,Taiwan cDepartmentofMathematics,Chung-AngUniversity,Dongjak-Gu,Heukseok-Dong221,Seoul156-175,SouthKorea

dCenterofMathematicalModelingandScientificComputing&DepartmentofAppliedMathematics,NationalChiaoTungUniversity,1001, Ta HsuehRoad,Hsinchu300,Taiwan

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

Articlehistory:

Received5July2015

Receivedinrevisedform16June2016 Accepted22June2016

Availableonline28June2016

Keywords:

Immersedboundarymethod Incompressiblemembrane Three-dimensional vesicle Navier–Stokesequations

Weextendourpreviousimmersedboundary(IB)methodfor3Daxisymmetricinextensible vesicleinNavier–Stokesflows(Huetal.,2014[17]) togeneralthreedimensions.Despitea similar spirit innumerical algorithms to theaxisymmetric case,the fully 3Dnumerical implementation is much more complicated and is far from straightforward. A vesicle membranesurfaceisknowntobeincompressibleand exhibitsbending resistance.Asin 3Daxisymmetriccase,instead ofkeeping thevesicle locallyincompressible,weadopt a modifiedelastictensionenergytomakethevesiclesurfacepatchnearlyincompressibleso thatsolvingtheunknowntension(Lagrangemultiplierfortheincompressibleconstraint) can be avoided. Nevertheless, the new elasticforce derived from the modified tension energyhasexactlythe samemathematicalformas theoriginalone exceptthedifferent definitions oftension.Thevesicle surfaceisdiscretized onatriangular meshwherethe elastic tension and bending force are calculated oneach vertex (Lagrangian marker in the IBmethod)of thetriangulation. A seriesofnumerical testsonthe present scheme are conducted toillustrate the robustnessand applicability ofthe method.We perform theconvergencestudyfortheimmersedboundaryforcesandthefluidvelocityfield.We then studythe vesicle dynamicsin various flows such as quiescent, simple shear, and gravitationalflows.Ournumerical resultsshowgood agreementswiththoseobtainedin previoustheoretical,experimentalandnumericalstudies.

©^{2016ElsevierInc.Allrightsreserved.}

1. Introduction

Avesicleis aliquiddropletwitharadius ofabout10 μm enclosed byaphospholipidmembranesuspended inanin- compressibleviscous fluid media.Such phospholipidmembraneconsistsoftwo-layer tightlypacked lipidmoleculeswith hydrophilicheadsfacingtheexterior andinteriorfluidswhilethehydrophobictailshideinthemiddle.Thisbilayermem- brane hasthe thicknessabout6 nm andexhibitsresistanceagainst membranedilationandbending.Therefore,itisquite naturaltoregardthismembraneasanincompressiblesurfacewithmechanicalfunctionsdeterminedbysomeenergyfunc- tional [16]. Thus, the dynamicsof vesicle influids can be determined by the membrane incompressibility,bending, and hydrodynamicalforces.Thestudyofvesicledynamicsinfluidflowhasbecomean activeresearchareainthecommunities

*

Correspondingauthor.

E-mailaddresses:ycseol@ntu.edu.tw(Y. Seol),wfhu@nchu.edu.tw(W.-F. Hu),kimy@cau.ac.kr(Y. Kim),mclai@math.nctu.edu.tw(M.-C. Lai).

http://dx.doi.org/10.1016/j.jcp.2016.06.035 0021-9991/©^{2016ElsevierInc.Allrightsreserved.}

ofsoftmatterphysicsandcomputationalfluidmechanicsinthepastyears.Forexample,understandingofvesiclebehaviors influid flowsmightleadtoa betterknowledgeofredbloodcells (RBCs)inbloodsimplybecausetheyboth sharesimilar mechanicalbehaviors[33].Certainly,ithasother practical applicationssuchasa drug-deliveryvehicleforcancer therapy [40]andamicro-reactor[13]forenzymaticmRNAsynthesisinbioengineering.

The interaction between the vesicle andsurrounding fluid makes the dynamics rich fromphysical point ofview. For instance,a vesiclecanundergotank-treading,tumbling,ortrembling motionundershear flow,see[6]andthereferences therein.Forthepasttwodecades,thevesicledynamicsingeneralflows(particularlyinshearflow)havebeenextensively studiedbyexperiments[22,24,6],theories[25,32,28,11],andnumericalsimulations,seethedetailedreferencesbelow.

The numericalsimulationsofthevesicleproblemnot onlyinvolveatwo-phaseincompressibleflowbutalsorequireto enforceanincompressibilityconstraintofthemembranesurface,whichmakestheproblemmorechallenging.Thenumerical methods forsimulating vesicleproblems inliterature can be characterized by howthe membranesurface is represented and how the fluid equations are solved. Based on this characterization, several methods have been developed such as boundary integral method [26,43,44,3,45,50,12], level set method [29,37,27,30,8], phase field method [9,2,30,1], particle collisionmethod[34],immersedinterfacemethod[21,41],andimmersedboundarymethodorfront-trackingmethod[18,20, 48,17];justtonameafewrecentones.Inallofthesenumericalmethods,howtoimposethemembraneincompressibility constraint is an important issue. The surface tension in vesicle problems, which has a different physical meaning from that ingeneraltwo-phaseflowproblems,isunknownaprioriandinfactactslikeLagrangemultipliertoenforcethelocal incompressibilityalongthesurface.Thisisexactlythesameroleplayedbythepressuretoenforcethefluidincompressibility inNavier–Stokesequations.

There aretwodifferentapproachestoenforcethelocalincompressibilityconstraintinliterature.Thefirstoneneeds to discretizethewholeequationsfirst(regardlessofusingboundaryintegral,finiteelement,orfinitedifferencemethod)and then tosolve the discretized equations simultaneously forthe tension andfluid variables. This approach can be explicit orsemi-implicitdepending onhowwe treatthetensionforce computations.Thereusually existsatrade-off betweenthe time-stepstabilityandefficiencyinthosealgorithmssimplybecauseiterativeproceduresareneeded.Mostoftheboundary integral method [45,3,50] or level set method [37,27] fall into this category. Another approach, which was used in our previous3Daxisymmetriccase[17],iscalledapenaltyidea.Insteadofkeepingthevesiclemembranelocallyincompressible, the penalty idea makes thevesicle surface patchnearly incompressibleby introducing a modified elastictension energy.

Thisapproachreplacestheunknowntensionbyaspring-liketensiondependingonthesurfaceconfigurationsothatwecan avoidsolvingthewholesystemtoobtainthevariabletension,whichsignificantlysimplifiesthenumericalalgorithm.Inthis paper,weextendourprevious immersedboundary(IB)methodforsimulatingincompressiblevesicles in3Daxisymmetric Navier–Stokesflows[17] togeneralthreedimensions.Weshallshowthatthenewelasticforce derivedfromthemodified tensionenergyhasexactlythesamemathematicalformastheoriginalelastic forceexceptforthedifferentdefinitionsof thetension.Wevalidatethisapproachbyperformingseveralnumericaltestsinoursimulations.

Therestofthepaperareorganizedasfollows.InSection2,wepresentthegoverningequationsforthevesicleproblem under the immersedboundary formulation. We also providesome notions inclassical differential geometrythat will be usedto computethegeometricalquantitiesofthevesiclesurfacemathematicallyandnumerically.Thenweintroduce our approachforanearlyincompressiblevesiclesurfaceandthemodifiedenergy.Thedetailednumericalalgorithmisdescribed inSection3;wefirstexplainhowtoevaluatethemeancurvaturevectorandbendingforcetermsona triangulatedsurface, then outlinethecomplete time-steppingschemeforthealgorithm,andfinallydiscusshowto maintainthe surfacemesh quality duringthesimulations.AseriesofnumericalteststovalidateourpresentalgorithmisgiveninSection4whichis followedbyconclusionandfutureworkinSection5.

2. Equationsofmotion

We consider a single incompressible vesicle (t) suspended in a three-dimensionaldomain  filled withviscous in- compressible Navier–Stokesfluid. FortheIB formulationinwhichthe fluid-relatedquantities arerepresented inEulerian mannerwhilethevesicle-relatedonesareinLagrangianmanner,thegoverningequationscanbewrittenasfollows.

ρ

 ∂

u

∂

t

+ (

^u

· ∇)

^u



= −∇

^p

+ μ 

u

+

^{f in}

,

(1)

∇ ·

^u

=

^{0 in}

,

(2)

f

(

x

,

t

) =





F

(

r

,

s

,

t

)δ(

x

−

^X

(

r

,

s

,

t

))

d A

,

(3)

∂

X

∂

t

(

r

,

s

,

t

) =

^U

(

r

,

s

,

t

) =





u

(

x

,

t

)δ(

x

−

^X

(

r

,

s

,

t

))

dx

,

(4)

∇

s

·

^U

=

^{0 on}

.

(5)

Here, weassume that thefluidsinsideandoutsidevesiclehavethesamedensity

ρ

andviscosity

μ

.Eqs.(1)and(2)are the incompressibleNavier–Stokes equationswiththe fluid velocity u(x,t) andthepressure p(x,t),where x= (^x,y,z) is

E

=

^Eb

+

^Eσ

=







_cb

2 H²

+ σ ^

d A

,

(6)

wherec_b isthebendingrigidity, H isthesurfacemeancurvature,and

σ

istheunknowntensionwhichactsasaLagrange multipliertoenforcetheincompressibilityconditioninEq.(5).Bytakingthevariationalderivativetothesurfaceenergy(6), onecanderivethevesicleboundaryforceF=^Fb+Fσ consistingofthebendingforceFb andtheelasticforceFσ as

F_b

=

^cb 2





_sH

+

^2H

(

H²

−

^K

)



n

,

F_σ

= ∇

s

σ −

^2H

σ

n

,

(7)

where K istheGaussiancurvatureofthemembranesurface,n istheunitoutward normalvectortothesurface,∇s isthe surfacegradientoperator,andsisthesurfaceLaplacian(orLaplace–Beltrami)operator,see[46]forthedetailedderivation.

Inclassicaldifferentialgeometry[7],thesurfacegradient ∇s

σ

ofthetension

σ

onasurface patchX(r,s)can becom- putedas

∇

s

σ ₌

^GXr

⁻

^{F Xs}

E G

−

^F2

σ

r

+

^EXs

−

^{F X}r

E G

−

^F2

σ

s

,

(8)

where E=^Xr·^Xr, F=^Xr·^Xs, andG=^Xs·^Xs are thewell-known coefficientsof thefirst fundamental form ofX. Note thatthesubscriptsr ands ofafunctiondenotethepartialderivativesofthefunctionwithrespecttor ands,respectively.

Similarly,thesurfacedivergenceofavectorfieldU canbedefinedby

∇

s

·

^U

=

^GXr

−

^{F X}s

E G

−

^F2

·

^Ur

+

^EXs

−

^{F X}r

E G

−

^F2

·

^Us

.

(9)

Notingthat the localsurface dilationfactoris|^Xr×^Xs|=√

E G−^{F2,andusingthevector triple productformulaa}× (b×^c)= (^a·^c)b− (^a·^b)c,onecaneasilyobtainthefollowingusefulrelations

X_s

×

ⁿ

=

^GXr

−

^{F X}s

|

^Xr

×

^Xs

| ,

n

×

^Xr

=

^EXs

−

^{F X}r

|

^Xr

×

^Xs

| .

(10)

Thus,thesurfacegradientofthetension

σ

inEq.(8)canalsobeimmediatelywrittenasasuccinctformwithoutusingthe coefficientsofthefirstfundamentalformas

∇

s

σ ₌ ⁽

^Xs

^×

ⁿ

⁾ σ

r

+ (

ⁿ

×

^Xr

) σ

s

|

^Xr

×

^Xs

| .

(11)

Onecanseelaterthatthepresentcomputationfor∇s

σ

isbasedonthediscretizationofEq.(11).

2.1. Approachofnearlyincompressiblemembrane

Asinour previous3D axisymmetric work[17], weadoptanapproach ofnearlyincompressiblemembraneinorderto avoidsolvingtheunknowntensiontoenforcetheincompressibilityconditiongiveninEq.(5).Toproceed,westartwiththe derivationoftherelationshipbetweentheincompressibilityconditionandthesurface dilationfactor|^Xr×^Xs|^.Therateof changeofsurfacedilationfactorcanbederivedinthefollowing:

∂

∂

t

|

^Xr

×

^Xs

| =

^Xr

×

Xs

|

^Xr

×

^Xs

| · (

^Xrt

×

^Xs

+

^Xr

×

^Xst

)

=

n

· (

Xrt

×

Xs

) +

n

· (

Xr

×

Xst

)



since n

=

^Xr

×

Xs

|

^Xr

×

^Xs

|



= (

^Xs

×

ⁿ

) ·

^Xrt

+ (

ⁿ

×

^Xr

) ·

^Xst

(

using

(

a

×

^b

) ·

^c

= (

^b

×

^c

) ·

^a

)

= (

^Xs

×

ⁿ

) ·

^Ur

+ (

ⁿ

×

^Xr

) ·

^Us

(

since X_t

=

^U

)

=

^GXr

−

^{F X}s

|

^Xr

×

^Xs

| ·

^Ur

+

^EXs

−

^{F X}r

|

^Xr

×

^Xs

| ·

^Us

(

using Eq.(10)

)

= (∇

s

·

^U

) |

^Xr

×

^Xs

| (

^{by Eq.(9)}

) .

(12)

Thisderivationnotonlyrevealstheexpressionofthesurfacedivergence∇s·U explicitlybutalsoshowsthatthecondition

∇s·^U=^{0 isequivalenttotheconstantdilationfactorovertime,whichisexactlywhatthelocal}incompressibilitymeans.

In ourmodel forthevesicle membrane,we relax theincompressibility condition in Eq.(5)andreplace theunknown tension

σ

inEq.(7)byaspring-likeelastictension

σ ₌ σ

0

 |

^Xr

×

^Xs

| − |

^X0r

×

^X0s

| 

,

(13)

where|^Xr×^Xs|^and|^X0r×^X0s|^{arethelocaldilationfactorsatlaterandinitialtimes,}respectively.Bychoosingthestiffness constant

σ

0tobesufficientlylarge,wecankeepthosetwolocaldilationfactorsascloseaswelike.However,sinceweallow thedeviationofthedilationfactoratlatertimesfromtheinitialone,thevesiclemembraneisnownearlyincompressible.

Thelegitimacyofpresentnearlyincompressibleapproachcanbevalidatedinournumericalresultslater. Underthisnearly incompressibleapproach,wenolongerneedtosolvetheunknowntensionnumerically.Bysimplysubstitutingthisexplicit formof

σ

inEq.(13)intothesecondtermofEq.(7),weareabletocomputetheelasticforceFσ directlywhichsimplifies numericalprocedure significantly.The aboveelasticforce Fσ canalsobe derivedby takingthe variationalderivativewith respecttoX tothefollowingelastictensionenergy

E_σ

(

X

) = σ

0

2

  |

Xr

×

Xs

| − |

X⁰_r

×

X⁰_s

| 

2

drds

,

(14)

butweomitthedetailhere.

3. Numericalscheme

Inthissection,wedescribeournumericalschemeandsomeimplementationdetailsforsolvingtheequationsofmotion presented intheprevious section. Beforeto proceed,we performanon-dimensionalization onthose governingequations first.Weusetheeffectiveradiusofvesicle R0=√

A/4

π

_{= (}3V/4

π

)^1/3 asthescalinglengthscale,where A and V arethe surface areaandtheenclosedvolumeofthevesicle,respectively.Thecharacteristictimescalet_c canbechosendepending ondifferentflowsituations.Allphysicalvariablesarescaledbytheassociatedcharacteristicscalesasfollows:

x^∗

=

^x R0

,

t^∗

=

^t tc

,

u^∗

=

^u R0

/

tc

,

p^∗

=

^t2c

ρ

R²₀p

, σ

^∗

₌

^R

20

c_b

σ .

Aftersomesimplecalculations,thedimensionlessgoverningequationsforthevesicledynamics(afterdropping*notation) canbesummarizedby

∂

u

∂

t

+ (

^u

· ∇)

^u

= −∇

^p

+

¹ Re



u

+







F_σ

+

¹ ReCaF_b



δ(

x

−

^X

(

r

,

s

,

t

))

d A

,

(15)

∇ ·

^u

=

⁰

,

(16)

F_σ

= ∇

s

σ −

2H

σ

n

,

Fb

=

¹ 2





sH

+

2H

(

H²

−

K

)



n

,

(17)

∂

X

∂

t

(

r

,

s

,

t

) =





u

(

x

,

t

)δ(

x

−

X

(

r

,

s

,

t

))

dx

.

(18)

There are two dimensionless numbers; namely, the Reynolds number Re=

ρ

R²₀/

μ

t_c and the capillary number Ca=

μ

R³₀/c_btc. As discussed before, we adoptthe nearly incompressible idea in which we relax the membrane incompress- ibilityconstraintEq.(5)andreplacetheunknowntensionbyaspring-likeelastictension

σ

asshowninEq.(13).Sincethe elasticstiffness

σ

0 ischosensufficientlylargeandisadjustable,thedimensionlessnumberinfrontofFσ canbeabsorbed inthenumber

σ

0.WeleavethedetailedstudyondifferentstiffnessnumberinSection4.2.

3.1. GridlayoutsforEulerianandLagrangianvariables

In ordertosolvetheNavier–Stokesequations(1)–(2),welay outauniformCartesian gridwithmeshwidthh= ^x=

y= ^{z inthewhole} computationaldomain  anddefine theEulerianfluid variablesonthe staggeredmarker-and-cell (MAC)gridintroducedbyHarlowandWelsh[15].IntheMACgridlayout,thepressurep isdefinedatthecellcenterwhile

Fig. 1. Fluid variables on a staggered MAC grid in 3D (left). Triangular surface patches that share the vertex Xk(right).

thevelocitycomponentsu= (^u,v,w)aredefinedatthecellfacesassociatedwith(x,y,z)directions,seetheleftofFig. 1 indetail.

Forthevesiclesurface,welayoutasurfacetriangulationanddefineallLagrangianvariablessuchasmeancurvatureH, elasticforceFσ ,andbendingforceFbattheverticesofthetriangleswhicharecalledtheimmersedboundary(Lagrangian) markers X.When we choosea marker(ora vertex, sayX_k), thetriangles thatshare thevertexmake1-ring (say T(X_k)), seetherightofFig. 1.Onaparticular triangle(say the-thtrianglein T(X_k)),we denotetheverticesof-thtrianglein counterclockwise orderviewingfromoutside ofthe surface by Xk=^{X1, X2,and X3,} respectively.Then the unit outward normalvectorofthe-thtrianglecanbecalculatedasn=_(^(X_X²2^{−X1)×(X3−X1)}

−^X1)×(^X3−^X1),andtheareaofthe-thtrianglecanalsobe easilycomputedasd A=(X²−X¹)× (X³−X¹)/2.

3.2. Elasticforcecomputation

As mentioned before, the discrete elastic force Fσ is computedat the vertices of the triangles. In the following, we shalldescribehowwecomputetheelastictension

σ

,itssurfacegradient∇s

σ

,andthemeancurvaturevectorH=^{2Hn at} those vertexpoints. Oncewe havethose values,theelasticforce densityatthe vertexX_k can becomputed byFσ(X_k)=

∇s

σ

(X_k)−

σ

(X_k)H(X_k).

Ina continuous level,thephysicalmeaning ofthespring-liketensioninEq.(13)is tokeep thesurface patchdilation factor(thus, surface patchareaby multiplyingdrds)closeto theinitial one.Therefore,itis quitenaturaltocompute the tensiononeachtriangle(saythe-thtriangle)by

σ

= ˜ σ

0



d A

−

^{d A0}



,

(19)

whereformally

σ

˜0=

σ

0/(drds),andd Aandd A⁰ refertotheareasofthesametriangleattwodifferenttimes.Usingthe valueof

σ

oneachtriangle,thetensionatvertexcanbecalculatedby

σ (

X_k

) = 

∈T(X_k)

σ

/

3

,

where1/3 comesfromthefactthatatriangleconsistsofthreedifferentvertices.Thenwecancomputethesurfacegradient

∇s

σ

oneachtriangleby

∇

s

σ

= (

X³

−

^X1

) ×

ⁿ

2 d A

( σ

²

− σ

¹

) +

ⁿ

× (

^X2

−

^X1

)

2 d A



σ

³

− σ

¹

^ ,

(20)

where

σ

ⁱ(i∈ {¹,2,3}^{)arethevaluesof}

σ

atverticesXⁱ (seetheverticeslayoutintherightofFig. 1).Theaboveformula issimply adirectdiscretizationofsurface gradient inEq.(11)undertriangularmesh.Oncewehavethe surfacegradient approximationonthetriangles,wecanestimatethesamequantityatthevertexXkby

∇

s

σ (

Xk

) = 

∈T(Xk)

ω

∇

s

σ

,

where

ω

=_{d A}^{d A}_(X^/_k³₎ ^isthecorrespondingarea-weight forthesurroundingtriangle.Notethat,theeffectivelocalarea atthe vertex X_k is approximatedby d A(X_k)=

∈^T(X_k)d A/3. Thiskind of area-weightedapproach to obtain some geometric quantitiesatvertexfromthevaluesofsurroundingfacetsisquitepopularinliteraturesuchasin[48,47,3,12].

Asforthemeancurvaturevector,weadoptthesameformulausedinthe3DfoampaperbyKimet al.[19].Thatis,the meancurvaturevectorH=^{2Hn ateachvertexX}kcanbecomputedby

H

(

X_k

) =

¹ 2 d A

(

X_k

)



∈T(Xk) n

× 

X³

−

X²



,

(21)

where X¹=^Xk andXⁱ fori∈ {²,3}^{have the samegeometric locations as}illustrated in the rightof Fig. 1. Surprisingly, themeancurvatureobtainedbyaboveformulaisequivalenttothewell-knowncotangentformulainMeyeret al.[31].We leave thisequivalence derivationfor both formulasin Appendix A. One should alsomentionthat, the above formulafor meancurvaturevectorcanberegardedasadiscreteversionof

S2Hnd A=

∂Sn×^{dX overasmoothsurfacepatchS [42].}

3.3. Bendingforcecomputation

In ordertocomputethebending forceFb(Xk) atthe vertexXk,we followtheideausedin[47].Toproceed,webegin withthediscreteversionofbendingenergyinEq.(6)as

E_b

[

^X

] =

¹ 8

Nv



k=1

|

^H

(

X_k

)|

^{2d A}

(

X_k

),

(22)

whereH=^{2Hn isthemeancurvaturevectorasbeforeandN}v isthetotalnumberofverticesinthesurfacetriangulation.

SubstitutingthemeancurvaturevectorformulaofEq.(21)into(22),Wuetal.[47]derivedthefollowingdiscretebending forceactingatthevertexX_k,

Fb

(

Xk

)

d A

(

Xk

) =

¹ 8



∈T(Xk)



(

H

−

n

·

C

)



1 2n

×

E^k



+

¹

2C

×

E^k

+

n

×

h^k



,

(23)

where H

=

¹

3



p∈^V()

|

^H

(

X_p

)|

²

,

C

=

¹ d A



p∈^V()

E^p

×

^H

(

X_p

),

E^k

=

^X3

−

^X2

,

h^k

=

^H

(

X³

) −

^H

(

X²

).

(24) Notethat,asusual,Xⁱfori∈ {²,3}^{areorderedina}counterclockwisedirection(viewingfromoutside)byfixingX¹=^Xkto calculateE^kandh^kinEq.(24).

3.4. Time-steppingscheme

Oncewe knowhow torepresentthevesicle surfaceandhowtocompute theelastictensionandthebendingforce on thetriangularmesh,thetimeintegrationofourschemecompletelyfollowsthetraditionalIBframework[35].Letusdenote the time stepsize byt andthe superscriptn asthe timestep index. Atthe beginningofeach time stepn,the vesicle configurationrepresentedbytheLagrangianmarkersXⁿ_k (theverticesofsurface triangularmesh)andthefluidvelocityuⁿ on Euleriangridare bothgiven.SothenumericalprocedureoftheIBmethodtoobtainuⁿ⁺¹ andXⁿ_k⁺¹ canbeperformed asfollows.

1. ComputetheLagrangianboundaryforceattheLagrangianmarkersby F

(

Xⁿ_k

)

d A

(

Xⁿ_k

) =

^Fσ

(

Xⁿ_k

)

d A

(

Xⁿ_k

) +

¹

ReCaF_b

(

Xⁿ_k

)

d A

(

Xⁿ_k

)

(25)

wheretheelastictensionforceFσ andthebendingforceF_barecomputedusingtheformulasdescribedabove.

2. SpreadtheLagrangianboundaryforceintotheEuleriangridusingthesmoothedDiracdeltafunctionδh as

fⁿ

(

x

) =

Nv



k=1

F

(

Xⁿ_k

)

d A

(

Xⁿ_k

)δ

h



x

−

Xⁿ_k

,

(26)

wherex= (^x,y,z)isthefluidmeshpoint.Forδ_h(x)=_h¹3φ_x

h

φ_y

h

φ_z

h

,weemploy

φ (

r

) =

⎧ ⎪

⎪ ⎪

⎨

⎪ ⎪

⎪ ⎩

1 8



3

−

²

|

^r

| + 

1

+

⁴

|

^r

| −

^4r2



,

if

|

^r

| ≤

¹

,

1 8



5

−

²

|

^r

| − 

−

⁷

+

¹²

|

^r

| −

^4r2



,

if 1

< |

^r

| ≤

²

,

0

,

otherwise

.

(27)

Fig. 2. Re-meshing triangulation by edge addition (top) and deletion (bottom).

3. SolvetheNavier–Stokesequations.Thiscanbedonebythefollowingbackwardsecond-orderprojectionmethod[14] in whichthenonlineartermisapproximatedbytheAdams–Bashforthscheme.

3u^∗

−

^4un

+

^un−1 2



t

+

²



uⁿ

· ∇

h

uⁿ

− 

uⁿ⁻¹

· ∇

h



uⁿ⁻¹

= −∇

hpⁿ

+

¹

Re



_hu^∗

+

^fn

,

(28)



hp^∗

=

³

2



t

∇

h

·

^u∗

, ∂

p^∗

∂

n

=

^{0 on}

∂,

u^∗

=

^un+1on

∂,

uⁿ⁺¹

=

^u∗

−

²



t

3

∇

hp^∗

,

∇

hpⁿ⁺¹

= ∇

hp^∗

+ ∇

hpⁿ

−

²



t

3Re



h

(∇

hp^∗

),

where∇h,∇h·^,andh arethesecond-ordercentraldifference operatorswhichapproximatethegradient,divergence, andLaplaceoperator,respectively,understaggeredMACgridsetting.Forthenonlinearterms,theskew-symmetricform isappliedas(u· ∇h)u=¹₂(u· ∇h)u+¹₂∇h(uu).Onecanseethatthemosttime-consumingpartinthisstepistosolve onePoissonequationforthepressureandthreemodifiedHelmholtzequationsforthevelocityfield.Fortunately,these ellipticsolverscanbeimplementedefficientlyusingFastFourierTransform(FFT).

4. InterpolatethenewvelocityfromtheEulerianfluidpointstotheLagrangianmarkersandthenadvancethemarkersto newpositionsbyusing

Xⁿ_k⁺¹

=

^Xn_k

+ 

^t



x

uⁿ⁺¹

(

x

)δ

_h

(

x

−

^Xn_k

)

h³

.

3.5. Themaintenanceoftriangularmeshquality

As mentioned before, we use triangular mesh to track the shape of the vesicle. To make an initial triangulation of the surface, we start by producing an icosahedron (20 triangles) on a unit sphere and then divide each triangle into 4 smallertrianglesbyconnectingthemidpointsofalledges.Wethenprojectallthenewvertexpointsontotheunit sphere.

Thisrefinement-projection procedureiscontinueduntil properresolutionis achievedon theunit sphere.Then theinitial configurationofvarioustypesofavesiclecanbeconstructedbyapplyingagivenfunctioninsphericalcoordinatesof(θ,φ), where(θ,φ)∈ [⁰,2

π

_{]× [}0,

π

]^.

Sincethe vesiclesurface mightdeformseverely during thenumericalsimulations indynamic flow,we apply aproper re-meshing technique to the triangularmesh to maintain the grid resolution over the surface. If the triangular meshis too coarse, there could be some volume leakage through the boundary surface, and ifit is too fine, a severe time step constraintfornumericalstabilitycouldbeoccurred.Wemaintainproperresolutionofthetriangularmeshbysimplyadding ordeletingmarkerpoints(andthusthetriangles)asneeded,inthefollowingway.

Ateachtimestep,wescanthelengthofthreeedgesforeachtriangleandmakesurethateachedge-lengthisbetween L_max=⁰.5h and L_min=⁰.24h,where h isthespatial meshwidth. Ifthe lengthofan edge islarger than L_max, weadd a newmarkerinthemiddleofthatedgeandthencreatemorenewtrianglesasillustrated intheupperpanelofFig. 2.We remarkthatthesumofareaofalltrianglesthroughthisrefinementprocessisequaltotheareaoftheoriginalconfiguration.

Meanwhile,ifthelengthofanedgeislessthanLmin,wedeletethisedgeandreplacethetwoendpointsoftheedgebya pointlocatedinhalfwaybetweenthem,seethebottompanelofFig. 2.Thisprocessdeletestwotriangleswhichpreviously sharethesameedge.Thisstrategyforre-meshingwasappliedtothenumericalsimulationsfor3Dfoamdynamicsin[19].

Fig. 3. Left:Atriangulatedbiconcavesurfaceanditscross-sectionalview.Right:Thecomparisonofmeancurvaturebetweennumericalvalues(symbols) andexactvalues(solidlines)forthreedifferentsurfaces.Here,thenumberofverticesisNv=2562 correspondingto5120trianglesused.

If we re-mesh the triangulationof a vesicle surface at some time step through the addition or deletion of an edge, thearea ofsome triangleswouldchangeaccordingly,seeFig. 2.Thiscertainly makestheentirenumericalalgorithmmore complex sinceweadoptthenearly incompressibleapproachofavesiclewhich usesthearea differenceofthetriangleat laterandinitialtimes.Toovercomethisdifficulty,whenare-meshingofthevesiclesurfaceoccursattimet,wealsoapply the same re-meshingprocedure to theinitial surface so that a one-to-one mappingcontinuesto exist betweend A and d A⁰inEq.(19).

4. Numericalresults

We performa seriesofnumericaltestsforthree-dimensionalvesiclesimulations influidflows. Webeginby checking therateofconvergenceofthedevelopedschemeforthecomputationofthemeancurvature,bendingforce,fluidvariables, andvesicleconfiguration.Then,wechoosedifferentstiffnesscoefficient

σ

˜0tostudyitseffectonthenearlyincompressibility bycomparinglocalandglobalsurfaceareas.Asapplications,wesimulatethedynamicsofasuspendedvesicleinquiescent, simple shear,andgravitationalflows. Throughoutthispaper,we choosetheeffectiveradius R0=^{1,theReynoldsnumber} Re=^{1,the capillarynumberCa}=^{50,andthegrid size1283, unlessotherwisestated. The numberoftriangles usedfor} an initialvesiclesurfaceischoseneither81920or327680withthecorrespondingnumberofvertices40962or163842so that the volumeleakage error islessthan 1% forall simulationswe haveperformed. Because ofthe vesicledeformation underdynamicflows,theredistributionoftriangularmeshpoints describedintheprevioussectionisappliedtomaintain a proper resolution of a discretized vesicle. For a given initial vesicle shape, another dimensionless numbernamed the reducedvolumeisdefinedas

ν

=_4π(A^3V_/4π)3/2 where A and V arethesurfaceareaandtheenclosed volumeofthevesicle, respectively.Thisnumber

ν

representsthevolumeratiobetweenthevesicleandthespherewiththesamesurfacearea.

4.1. Accuracycheckformeancurvatureandbendingforce

Asthefirsttest,wechecktheaccuracyofthenumericalschemesinEqs.(21)and(23)forcomputingthemeancurvature andthebendingforce,respectively.Here,weconsiderthefollowingthreedifferentsurfacesdescribedinsphericalparamet- ric coordinates(θ,φ)∈ [⁰,2

π

_{]× [}0,

π

]^{forwhichthemeancurvatureandother} geometricalquantitiesintheformulafor thebendingforcecanbecomputedanalytically.

1. Unitsphere:X(θ,φ)= (^cosθsinφ,sinθsinφ,cosφ); 2. Ellipsoid:X(θ,φ)= (⁰.5 cosθsinφ,0.5 sinθsinφ,cosφ); 3. Biconcavesurface:X(θ,φ)=

cosθsinφ,sinθsinφ, (0.1242+⁰.8012 sin²φ−⁰.4492 sin⁴φ)cosφ

 .

We drawthebiconcave surfacetogether withits triangularmeshintheupper-leftpanel ofFig. 3.Inordertocompare themeancurvatureH attheverticesonthesurface,weprojecttheverticesontothexy-planeandcalculatethenormalized distancer fromtheprojectedverticestotheoriginO .Forthebiconcavesurface,thesymbolsinthelower-leftpanelofFig. 3 indicatetheprojectedvertices.Intherightpanel,weshowthepoint-wisevaluesofthemeancurvaturewithrespecttothe distancer forallthethreedifferentsurfaces.Thesymbolsrepresentthenumericalvalueswhilethesolidlinesrepresentthe exact valuesofthemeancurvatures.Onecanseethatthenumericalvaluesforthemeancurvaturematchwiththeexact valuesquitewell exceptattwopointswithr≈⁰,0.9.Weattributethisdiscrepancytothedefectofcurrenttriangulation, sincethesetwopointshaveonlyfiveneighboringtrianglesratherthansix.

Eventhough thepoint-wisevalues ofthecomputational meancurvaturecan behavenot wellatsome vertices,aswe canseefromtheleftpanelofFig. 4,theL² errorsforthemeancurvaturecomputationbehavelikefirst-orderconvergence

underthecurrenttriangulationasweincrease thenumberoftriangles.(TheL² errorofsome functionψ iscalculatedas

 Nv

k=¹|ψ^e(X_k)− ψ(^Xk)|^{2d A}(X_k),whereψ^e(X_k)istheexact valueandψ(X_k) isthecomputedvalue.Foravector-valued function, the absolutevalue should be replaced by the vector norm.) We also show the L² errors forthe bending force F_bd A intherightpanel ofFig. 4.Itisimportanttomentionthat,insteadofstudyingtheconvergenceofF_b,westudythe convergenceofF_bd A.ThisissimplybecausethelatteroneistheimmersedboundaryforceusedinEq.(25).Again,onecan seethatthebendingforcecomputationstillbehaveslikefirst-orderconvergenceasweincreasethenumberoftriangles.

Weshouldmakearemarkthat,forgeneralunstructuredsurfacetriangulations,theabovecomputationalvaluesformean curvatureorbendingforce maynot convergeatall. Inordertohavethe point-wiseconvergenceforthe meancurvature, oneshouldconstructatleastalocalquadraticsurfaceateachvertextoestimatethosegeometricquantities[48].However, that wouldmake ourcode implementations morecomplicatedandwe leave a furtherimprovement forcomputationsof thosegeometricalquantitiestofuturework.

4.2. Studyondifferentstiffnessparameter

σ

˜0

IntheusageofnearlyincompressibleapproachintroducedinSection 2.1,we shouldchoosethestiffnessparameter

σ

˜0 in Eq.(19)large enough to keep the vesiclenearly incompressible. We heretest three differentvaluesof

σ

˜0=⁶×¹⁰⁴, 6×¹⁰⁵,6×^{106 toseehowthischangeaffectsthelocalandglobal}preservationsofsurfacearea.Weputanoblatevesicle withtheinitialshapeas

X

(θ, φ) = (

cos

θ

sin

φ,

sin

θ

sin

φ,

cos

φ/

3

.

5

)

(29)

inaquiescentflowofwhichthereducedvolumeisaround

ν

=⁰.643.Thecomputationaldomainischosenas= [−²,2]^3, andthegridnumberN=^{128 isusedforallspatialdirectionssothemeshwidthish}=¹/32.Thenumberoftrianglesused intheinitialvesiclesurfaceisabout81920 with40962 numberofvertices.Thetimestepsizeischosenast=^h/16.

Fig. 5showsfourtime evolutionaryplots;namely(a)themaximumrelativeerrorofthelocalsurface area(^{d At}(X)− d A⁰(X))/d A⁰(X)∞,(b)therelativeerroroftheglobalsurfacearea|^At−^A0|/^{A0,(c)therelativeerroroftheglobalvolume}

|^Vt−^V0|/^{V0,and(d)thetotalenergy.OnecanseefromFig. 5(a)and(b)thatthelocalandglobalerrorsofthesurfacearea} decreaseastheelasticstiffness

σ

˜0increases.Moreprecisely,whenthecasesofstiffnessparameter

σ

˜0=⁶×^105and6×^106, themaximumrelativeerrorsofthelocalsurface areaareabout1% whiletherelativeerrorsoftheglobalsurface areaare lessthan0.1% and 0.01%,respectively.It isimportanttomentionthattheerrorsshowninFig. 5(a)areinthemaximum normmeaning thatall thetriangularelements inourdiscretizationsatisfy theareadilationorcompressionwithin 1% as timeup tot=^{40 whichshowsthelegitimacy ofpresentnearly}incompressibleapproach.Itisalsointerestingtomention thatthe globalerrordecreasesfromtheorderof10⁻² to theorderof10⁻⁴ inthefactorof1/10 asweincrease

σ

˜0 bya factorof10.Ofcourse,themaintenance ofgoodmeshquality asdescribedinSection 3.5hastobeadoptedfromtimeto time.

Fig. 5(c)showsthat therelative erroroftheglobalvolumeofthevesiclealsodecreasesaftersometime astheelastic stiffness increases. Since the stiffness parameter

σ

˜0 is adjustable, the readersmight wonder whetherdifferent

σ

˜0 cause different flow behaviors. In Fig. 5 (d), we can see that all three different

σ

˜0, the total energy are almost in the same magnitudeuptotimet=^{40.Inotherwords,evenwe usedifferentstiffness parameter}

σ

˜0,theenergycontributedtothe fluidsystemarealmostthesame.Althoughitis notshownhere,thevesicleshapedynamicsforthosethreecasesarequite similar.

4.3. Studyonvesicleconfigurationandfluidvelocity

Wehereperformtheconvergencestudyforthepresentnumericalalgorithm.Again,weputan oblatevesiclewiththe sameshapeasintheprevioussubsectioninquiescentflowinitially.Thecomputational domainisthesameasbeforeandwe

Fig. 5. Thecomparisonforthreedifferentstiffnessparameters:σ˜0=6×10⁴(),6×10⁵(2),and6×10⁶ ().(a)The maximumrelativeerrorofthe localsurfacearea;(b)therelativeerroroftheglobalsurfacearea;(c)therelativeerroroftheglobalvolume;(d)thetotalenergy.

Fig. 6. The ratios of convergence for the fluid velocity(u,v,w)and the vesicle configuration X. Left: N=32; Right: N=^64.

choosefourdifferentCartesiangridsizesN=³²,64,128,256 tostudytheconvergence.Forthevesiclesurfacetriangulation, aswedoublethegridsizeN,weneedtohalvealltheedgesofthetrianglessothatthenumberoftrianglesusedmustbe increasedbyafactoroffourcorrespondingly.(ForthecaseofN=^32,weuse5120 numberoftriangleswhichcorresponds to 2562 numberofvertices.)Tobe consistent,we alsoneedto increasethestiffness parameterby afactorof fouraswe doublethegridsize,thatis

σ

˜0= (^N/32)²×^{104.Thetimestepsizeis}t=^h/8.

Fig. 6 showstherateofconvergenceforthevelocitycomponentsu= (^u,v,w)andthevesicleconfiguration X atsome different timesup tot=^{10.Since the analytic solution isnot known, we estimate the erroron each grid by using the} numericalsolutioninthenextfinergridasthereferencesolution.Inthiscase,forthenumericalvalueofthex-component ofthe velocityuN withthegrid size N,theconvergent ratiocanbe computedusingthetwosuccessiveerrors asratio= log₂(^uN−^u2N∞/^u2N−^u4N∞).(Theratiosforthevelocitycomponentsv,w andthevesicleconfigurationX aredefined similarly.)OnecanseefromFig. 6thattheaverageratioisaroundone,whichindicatesthatournumericalmethodisroughly first-orderaccurate.ThisconfirmsthattheIBmethodisfirst-orderaccurateingeneral.Itisalsointerestingtomentionthat the velocity components u and v have prettysimilar convergencebehaviors dueto theaxis-symmetric property of the vesiclesurfaceinthepresent3Dsimulations.

4.4. Asuspendedvesicleinquiescentflow

We choose a suspended vesicle withthree differentinitial shapes andstudytheir equilibriumin quiescentflow. The initialshapesofthesuspendedvesiclearethefollowing:

Fig. 7. Suspendedvesiclesinquiescentflow.Oblatevesicle(twoupper-left panels);Prolatevesicle(twoupper-rightpanels);Oscillatory vesicle(lower panels).

• ^{Oblatesurface:X}(θ,φ)= (^cosθsinφ,sinθsinφ,cosφ/3.5)

• ^{Prolatesurface:X}(θ,φ)= (⁰.2cosθsinφ,0.2sinθsinφ,cosφ)

• Oscillatorysurface[44]:X(θ,φ)=

3

20r(φ)cosθsinφ,₂₀³r(φ)sinθsinφ,₂₅³r(φ)cosφ



,wherer(φ)=

cos²φ+^{9 sin2}φ+ cos²(4φ).

Thesevesicleshaveapproximatelythereducedvolume

ν

=⁰.64,0.63,and0.60 correspondingly.Thecomputationalsetup uses= [−²,2]^3,

σ

˜0=⁶×^105,andt=^h/16.InFig. 7,thetwoupper-leftpanels showtheoblatevesicleattwodifferent times(t=^{0 andt}=¹⁹.5) wherethelatteroneis inequilibrium.The equilibriumshapeturnsout tobe abiconcave one which is a typical shape of RBCs. Forthe prolate vesicle in the two upper-rightpanels, the equilibrium shape turns to be a dumbbellone.For theoscillatory vesicle inthelower panel, thetransition ismorecomplicatedso we providefour snapshotsatfourdifferenttimes.One cansee thatthe equilibriumshapealsobecomes a typicalbiconcave one whichis ingoodagreementwiththoseobtainedinliterature[44,17]despitethefactthatthepresentsimulationsarefully3D.The presentresultsshownotonlythenumericalstabilityofourschemebutalsoagoodpreservationofthesymmetricstructure.

4.5. Avesicleinshearflow

Wesimulatethedynamicsofaprolatevesicleinasimpleshearflowu= (^z,0,0),seetheschematicviewoftheproblem setup in Fig. 8. It is well-known that there are three different flow regimes; namely,(1) a stationary tank-treading(TT) regimewherethevesicleexhibitsatank-treadingmotionwithnochangeofits shapeintime, (2)atumbling (TB)regime wherethe vesicletumbleslikea rigidbody motion,and(3) atrembling (TR)regime wherethemajor axisofthevesicle oscillatesabouttheshearflowdirection.Ithasbeenwellstudiedinliterature(see[6])thattheabovethreeregimescanbe fullycharacterizedbythenon-dimensionalshearrate

χ

(equalstothecapillarynumberCa inEq.(15)whenthetimescale tc ischosenastheinverseoftheshearrate),theviscositycontrastλ(theviscosityratiobetweentheinteriorandexterior fluids),andthereducedvolume

ν

.Throughoutthissubsection,weusealargercomputationaldomain= [−³,3]^{3 anda} smallerReynoldsnumberRe=^{10−31 sothatwecancompareourpresentnumericalresultsagainstwiththe}theoretical andother numericalresults inStokes regime [26].The mesh widthandthetime step sizeare chosen ash=⁶/128 and

t=^h/64 forthesesimulations.Here,werestrictourflowregimetobeTTbychoosingtheviscositycontrastλ=^{1 andTB} bychoosingλ=^40.

IntheTTregime,weinvestigatethedependenceofthetank-treadingmotionofthevesicleonthenon-dimensionalshear rate

χ

andthereducedvolume

ν

.Theinitialvesicleischosenasaprolateshapeasshownintheprevioussubsectionwith fivedifferentreducedvolume

ν

₌0.8,0.85,0.9, 0.95,and0.975 bykeepingtheeffectiveradiusR0=^{1 fixedbutadjusting} theaspectratio.Astheshearflowturnson,thevesiclealignswiththeflowdirectionsothatastationarytank-treading(TT) motionoccurswhilethevesicleshapedoesnotchangeintime,seethebottompanelsofFig. 8.Thisfeaturediscriminates incompressiblevesiclefromliquid dropletthat elongatesits shape asthelinearshear rateincreases [36]. Theinclination angleisdefinedastheanglebetweenthe x-axisandthemajoraxisL of avesicleprojectedontothe xz-plane,seeFig. 8.

Thistank-treadingmotioncanbecharacterizedbytheinclinationangleθ,andthescaledtank-treadingfrequency

ω

.Here, thetank-treadingfrequency

ω

canbecomputedusingtheformula

ω

=_N¹_v Nv

i=1|^r×^v|

|^r|² ,wherer andv arethepositionand velocity ofthevertices projectedon thexz-plane, respectively.Fig. 9 showsthetank-treadinginclination angle(left) and thescaledtank-treadingfrequency(right)versusthereducedvolumeforthreedifferentdimensionlessshearrates:

χ

=^1, 10,and100.As in[26],we canconcludethat both theinclinationangle θ andthe scaledtank-treadingfrequency

ω

are