An immersed boundary method for simulating Newtonian vesicles in viscoelastic fluid

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

An immersed boundary method for simulating Newtonian vesicles in viscoelastic fluid

Yunchang Seol

^a,∗

, Yu-Hau Tseng

^b

, Yongsam Kim

^c

, Ming-Chih Lai

^a

aDepartmentofAppliedMathematics,NationalChiaoTungUniversity,1001TaHsuehRoad,Hsinchu30010,Taiwan bDepartmentofAppliedMathematics,NationalUniversityofKaohsiung,Kaohsiung81148,Taiwan

cDepartmentofMathematics,Chung-AngUniversity,Dongjak-Gu,Heukseok-Dong221,Seoul156-175,SouthKorea

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

Articlehistory:

Received27March2018

Receivedinrevisedform27August2018 Accepted18October2018

Availableonline22October2018

Keywords:

Viscoelasticity Oldroyd-Bfluid Vesicle Tank-treading Tumbling

Immersedboundarymethod

In thispaper, a two-dimensionalimmersed boundary methodis developed to simulate the dynamicsofNewtonian vesicleinviscoelasticOldroyd-Bfluidundershearflow.The viscoelasticity effect of extra stress is well incorporated into the immersed boundary formulation using the indicator function. Our numerical methodology is first validated in comparison with theoretical results in purely Newtonian fluid, and then a series of numerical experiments is conducted to study the effects of different dimensionless parametersonthevesiclemotions.Althoughthetank-treading(TT)motionofNewtonian vesicleinOldroyd-BfluidundershearflowcanbeobservedjustlikeinNewtonianfluid, it is surprising to find that the stationary inclination angle can be negative without the transition to tumbling (TB) motion. Moreover, the inertia effect playsa significant role that isableto turn the vesicleback topositive inclination anglethroughTT-TB-TT transitionastheReynoldsnumberincreases.Tothebestofourknowledge,thisisthefirst numericalworkforthedetailedinvestigationsofNewtonianvesicledynamicssuspendedin viscoelasticOldroyd-Bfluid.Webelievethatournumericalresultscanbeusedtomotivate furtherstudiesintheoryandexperimentsforsuchcouplingvesicleproblems.

©^{2018ElsevierInc.Allrightsreserved.}

1. Introduction

Avesicle isa closedphospholipidmembrane,separatingthe inner fluidfromthe outer medium.Its typical size isan orderof10 μm,sothedimensionlessReynoldsnumberwithtypicalshearratelessthan10 Hzisabout10⁻³ inwater.So far,theexperimentsofvesiclehavebeenconductedinpurelyNewtonianfluid,whereNewtonianvesiclesaresuspendedin surroundingNewtonianfluids.Ingeneral,vesiclesmimictheredbloodcells(RBCs),andthebloodcontainingRBCsexhibits non-Newtonianpropertiessuchasviscoelasticity,shear-thinning,thixotropy,andyieldstress[9,39,40].Motivatedfromthis, weherestudytheviscoelasticfluideffectontwo-dimensionalvesicledynamicsundershearflow.AnOldroyd-Bfluidmodel isadopted toaccount the viscoelasticity,andan immersedboundary (IB)method isdevelopedto simulateprimarily the viscoelastic effectona Newtonian vesicle interactingwithsurrounding Oldroyd-Bfluid in two-phase flows,although our methodisnotonlyrestrictedtothistypeoffluid.

For the past two decades, three types of vesicle motions under shear flow have been widely studied, namely tank- treading,tumbling,andvacillating-breathing(ortrembling,swinging).Thosemotionscanbecharacterizedbytheso-called

*

Correspondingauthor.

E-mailaddresses:[email protected](Y. Seol),[email protected](Y.-H. Tseng),[email protected](Y. Kim),[email protected](M.-C. Lai).

https://doi.org/10.1016/j.jcp.2018.10.027

0021-9991/©^{2018ElsevierInc.Allrightsreserved.}

Fig. 1. SimulationsetupofaNewtonianvesicleinOldryod-BfluidundershearflowdenotedbyN/Ofluid.ThevesicleisfilledwithaNewtonianfluid (viscosityμn)whiletheOldryod-Bfluidisamixtureofpolymer(viscosityμp)andNewtoniansolvent(viscosityμs).Theangleθinradiansdenotesthe inclinationangleofthevesicle.

reducedareaofvesicleconfiguration,theviscositycontrastbetweeninnerandouterfluids,andthecapillarynumberasso- ciatedwithmembranebendingrigidity. So,bycarefullycontrollingthesedimensionlessparameters,extensivestudieswere presented intheories[21,29,25], experiments[18,20,10], andsimulations [24,47,4,11,5],justtonamea few.Refer[1] and references thereinfora review ofvesiclerheology. Recently,the inertia effectinfiniteReynolds numberflowbecomes a popularresearchsubjectaswell[23,36,27].

Among manynumericalmethods,theIB(orfront-tracking)methodisquiteusefulforbiologicalproblemstosimulatea flexible structure interactingwithsurroundingbulk fluid.There arevarious extensionsofIB methodtomodelvesicles or RBCs dynamicsinpurelyNewtonianfluid [22,26,14,45,15,37,35].Anotherextensionincorporatesa richerfluid suchasthe Oldroyd-B fluidto studytheviscoelastic effecton dropletdynamics[38,8,30,17].Thereexists somenumericalmethods to studydropletdynamicsinmultiphaseviscoelasticflow.In[7],thevolume-of-fluid(VOF)methodisemployedtoinvestigate adropletdeformationundershearflow.Thelevel-setmethodin[33] andthefront-trackingmethodin[17] areemployedto studyadropletinbuoyancy- andpressure-drivenflows.Inthesestudies,themultiphaseviscoelasticflowisaccomplished by varying the parameter ofpolymer relaxation time (or the dimensionlessWeissenberg number) inthe polymer stress evolution equation. Inthis framework,a smoothed indicator functionmultiplied tothe parameteris closeto zerointhe regionofNewtonianphase,whichisproblematicwhensolvingtheequation.Tocircumventthisdifficulty,thediffuseinter- face(orphase-field)methodproposedin[46] tostudyNewtoniandropletsuspendedinanOldroyd-Bfluidundercreeping shearflowadoptsadifferentapproach;thatis,thetimeevolutionofpolymerstressequationissolvedontheentiredomain firstandthentheextrapolymerstressisaddedbacktotheviscoelasticregionviathephasefieldfunction.Here,following the similar idea ofhandlingthe extra stress partin [46], we develop an immersedboundary methodforsimulating the dynamicsofNewtonianvesicleinanOldroyd-Bfluidwhichalsotakestheviscositycontrastandinertiaeffectintoaccount.

Tothebestofourknowledge,itseemstobethefirstnumericalworkonthiscouplingflowproblem.Twomajorinteresting findingsarebriefedhere.Firstly,thetank-treadingmotionofNewtonianvesicleinOldroyd-Bfluidundershearflowcanbe observedjustlikeinNewtonianfluid.However,asweincreasetheWeissenbergnumber,thestationaryinclinationanglecan benegativewithoutthetransitiontotumblingmotion.Secondly,byincreasingtheinertiaeffectontheunmatchedviscosity contrast case, thetank-treading motionwithnegative inclinationangle can turninto thesame motionbutwithpositive inclinationangle.Theinterplaybetweentheviscoelasticityandtheinertiaisinterestingandworthsfurtherinvestigations.

The rest ofthis paperis organized asfollows.In next section, we presentthe governing equationsbased onthe im- mersed boundaryframework whichdescribe thetwo-phase flowconsistingofa Newtonianvesicle inan Oldroyd-Bfluid.

InSection 3,wetransformtheequationsintodimensionlessformandoutlinethetime-steppingalgorithm forsolving the whole fluidsystem.Toverifyourcodeimplementation,numericalaccuracytestforanOldroyd-Bfluidsolverisconducted aswell.InSection5,aseriesofnumericalexperimentsisconductedtostudytheeffectsofdifferentdimensionlessparam- etersontheNewtonianvesicledynamicsundershearflowinanOldroyd-Bfluid.Conclusionsandfutureworkaregivenin Section6.

2. Mathematicalmodel

Inthispaper,weconsideranimmiscibletwo-phaseflowconsistingofaNewtonian vesicleinanOldroyd-B fluidunder shear flow ina two-dimensional domain  asdepicted inFig. 1.The vesicle boundary  separatingthose two fluidsis assumed to be a closedsmooth curve immersed in .The modelis formulated by the immersed boundary method in which thefluid partsindifferentregions usedifferentmodels(Newtonian orOldroyd-B)andthe vesicleeffectistermed astheelastictensionandbendingforce alongthevesicle boundary.Certainly,thoseforcesare derivedby theenergylaw definedonthevesicleboundarydescribedinliterature.Meanwhile,therearerelatedenergylawstodescribetheviscoelastic effectoffluidssuchasphase-fieldmethoddevelopedin[46].Basedonthephase-fieldvariableanditsgradient,themixing energy can be represented by Ginzburg–Landau form which produces the interfacial tension force between two phases.

The extrastress forpolymercontributionin anOldroyd-B fluidcan be derivedfromthedumbbellelastic energythrough a variationalprocedure.Byaddingtheviscous stressofsurroundingfluid,thetotalstresstensoristhesumofthosethree

terms.Thefluidvelocityu(x,t)andpressure p(x,t)aredescribedinEulerianmannerwhilethevesicleboundaryX(

α

,t)is describedinLagrangianmanner,sothegoverningequationsofthistwo-phaseflowcanbewritteninasinglefluidsystem asfollows.

ρ

 ∂

u

∂

t

+ (

^u

· ∇)

^u



= −∇

^p

+ ∇ ·

^[

((

1

−

^H

) μ

n

+

^H

μ

s

) D(

^u

) +

^H

σ

]

+

^{f in}

,

(1)

∇ ·

^u

=

^{0 in}

,

(2)

f

(

x

,

t

) =





(

F_γ

+

Fb

)( α ,

t

) δ(

x

−

X

( α ,

t

))|

X_α

|

d

α

in

,

(3)

∂

X

∂

t

( α ,

t

) =

^U

( α ,

t

) =





u

(

x

,

t

) δ(

x

−

^X

( α ,

t

))

dx on

,

(4)

γ ( α ,

t

) = γ

0

 |

^Xα

( α ,

t

)|

|

^Xα

( α ,

0

)| −

¹



on

,

(5)

F_γ

( α ,

t

) =

¹

|

X_α

|

∂( γ _τ)

∂ α ^,

^Fb

⁽ α ,

t

) =

cb



κ

ss

+ κ

³

2



n on

,

(6)

ξ σ

^

+ σ = μ

p

D(

^u

)

in

,

(7)

where

σ



₌ ^∂ σ

∂

t

+

^u

· ∇ σ ₋ ^ _∇

u

σ ₊ σ ( ∇

^u

)

^T



,

(8)

D(

^u

) = ∇

^u

+ (∇

^u

)

^T

, ∇

^u

=



ux uy

v_x v_y



.

(9)

Eqs. (1) and(2) aretheNavier–StokesequationsinwhichH istheindicator(orHeaviside)functiondefinedby H

(

x

,

t

) =

1 if x is in Oldroyd-B fluid,

0 if x is in Newtonian vesicle

.

⁽¹⁰⁾

Herewe assumebothfluidshavethesameconstantdensitydenotedby

ρ

,whereastheNewtonian fluidviscosity

μ

n and theOldroyd-Bfluidviscositymaybedifferent.TheNewtoniansolventviscosity

μ

sandthepolymerviscosity

μ

pdetermine the total viscosity of Oldroyd-B fluid, see Fig. 1. The extra stress

σ

contributes the viscoelasticity to the Navier–Stokes equationsonlyinthe Oldroyd-Bfluid partwhere H(x,t)=^{1.InEq. (3),theEulerianforce f arisesfromthespreading of} interfacial force F via theDiracdeltafunction δ(x)= δ(^x)δ(y).Similarly, asshowninEq. (4), theinterfacial velocity U is interpolatedfromthelocalEulerianfluidvelocityvia thedeltafunction.Heretheinextensiblevesicle boundaryisrelaxed bythenearlyinextensibleapproachinwhichaspring-liketensiongiveninEq. (5) isintroduced.Thisapproachismotivated fromtheequality∂|Xα|/∂^t= |Xα|∇s·^{U andhasbeenvalidatedtobeeffectivefor}applicationsundervariousflowsin2D[3], axisymmetriccase[15],and3D[35].Theinterfacialforce inEq. (6) consistsoftwocomponents;namely,thetensionforce Fγ andthebendingforce F_b.Inthebendingforce inEq. (6),c_b isthe bendingrigidityconstant,

κ

isthe localcurvature, andthesurfaceLaplacianof

κ

isrepresentedby

κ

ss

=

¹

|

^Xα

|

∂

∂ α



1

|

^Xα

|

∂ κ

∂ α



= −

^Xα

·

^Xαα

|

^Xα

|

⁴

∂ κ

∂ α ⁺

1

|

^Xα

|

²

∂

²

κ

∂ α

²

^.

⁽¹¹⁾

The latter expansion formula will be used to enhance the numerical accuracy, refer to Appendix D in [44]. The vesicle boundary(t) betweentwofluidsisa one-dimensionalcurve representedby (t)= {^X(

α

,t)= (^X(

α

,t),Y(

α

,t))|⁰≤

α

<

2

π

_}, where

α

is the Lagrangian coordinate. The subscript

α

of a function denotes its partial derivative of the function withrespectto

α

.We then definethe unit tangentvector asτ = (τ1,τ2)=Xα/|Xα| ^{andtheoutward unit normalvector} as n= (τ2,−τ1), where |· | ^{stands for the usual Euclidean norm. So the signed curvature ofthe curve can be written} as

κ

= (Xα Yαα−Yα Xαα)/

Xα²+^Y2_α3/2

.Theextra stress

σ

(x,t) inEq. (1) representsthe viscoelasticcontributionofthe Oldroyd-BfluidandisgovernedbyEq. (7),wheretheparameterξ isthepolymerrelaxationtime,

σ

^istheupperconvected time derivative of

σ

defined in Eq. (8), and D(^u) is the rate of deformation tensor as in Eq. (9). Of course, the above governingequations(1)–(7) shouldaccompanywithsomesuitableinitialandboundaryconditionswhichwillbedescribed later.

3. Numericalmethodandaccuracy 3.1. Dimensionlessgoverningequations

The governingequations(1)–(7) canbe writteninthedimensionlessformasfollows.Foravesicleboundaryenclosing an area A withitsperimeter L,we definethe effectiveradius ofthevesicle R=√

A/

π

asthe characteristiclength scale.

Thecharacteristictimescaleisdefinedbyt_c=¹/

γ

˙,where

γ

˙ istheshearrateofflow.Thenallphysicalvariablesarescaled bytheassociatedcharacteristicquantitiesas

x^∗

=

^x

R

,

t^∗

=

^t tc

,

u^∗

=

^u R

/

tc

,

p^∗

=

^tc

μ

s

p

, σ

^∗

₌

^tc

μ

s

σ .

(12)

After performing substitutions and a few calculations, the dimensionless governing equations for Newtonian vesicle in Oldroyd-Bfluidsystem(N/Ofluid)canbewrittenas(droppingtheasterisk)

Re

 ∂

u

∂

t

+ (

u

· ∇)

u



= −∇

p

+ ∇ ·

[

((

1

−

H

)λ +

H

) D(

u

) +

H

σ

]

+

f in

,

(13)

∇ ·

^u

=

^{0 in}

,

(14)

f

(

x

,

t

) =





F_γ

+

^Fb

 ( α ,

t

) δ(

x

−

^X

( α ,

t

)) |

^Xα

|

^d

α

in

,

(15)

∂

X

∂

t

( α ,

t

) =

^U

( α ,

t

) =





u

(

x

,

t

) δ(

x

−

^X

( α ,

t

))

dx on

,

(16)

γ ( α ,

t

) = γ

0

 |

^Xα

( α ,

t

)|

|

^Xα

( α ,

0

)| −

¹



on

,

(17)

F_γ

( α ,

t

) =

¹

|

^Xα

|

∂( γ τ)

∂ α ^,

^Fb

⁽ α ,

t

) =

¹ Ca





s

κ + κ

³

2



n on

,

(18)

W i

σ

^

+ σ = β D(

u

)

in

,

(19)

where

Re

= ρ

R²

μ

stc

,

Ca

= μ

_sR³ c_btc

,

W i

= ξ

tc

, λ = μ

_n

μ

s

, β = μ

_p

μ

s

.

(20)

Here, Re isthefluidReynoldsnumber,thecapillarynumberCa measurestheratioofviscousforceandbendingforce,and theWeissenbergnumberW i measurestheratioofelasticforceandviscousforceinOldroyd-Bfluid.Wealsodenoteλ by the viscosity contrastbetweenthe interiorNewtonian fluid andthe Newtonian solventof Oldroyd-B fluid,and β by the viscositycontrastbetweenthepolymerandtheNewtoniansolvent.Forinstance,inasingle-phaseNewtonianfluid,weset theparametersbyβ=^{W i}=^{0 and}λ=^{1.Noticethat,thereisno}dimensionlessparameterintherighthandsideofFγ term sinceitcanbeabsorbedintothestiffnessparameter

γ

0 whichwillbechosenlater.

3.2. Numericalscheme

Wenextpresentthenumericalschemeandsomerelatedimplementationdetailsforsolvingthedimensionlessgoverning equations(13)–(19). Throughoutthispaper, theperiodicboundarycondition isimposed forall physicalvariablesinthe x direction, whilethe Dirichlet andzero Neumannboundary conditions are imposed for thevelocity and the extra stress, respectively, in the y direction. To solve the Navier–Stokes equations (13)–(14) in a computational domain ⊆ R^{2, we} layout auniformCartesian gridwithmeshwidthh= ^x= ^{y,andallocatethefluid velocityu}= (^u,v)andthepressure p inastaggeredmarker-and-cell(MAC)gridmanner[16].Expressingtheextrapolymerstressby

σ =

 σ

^a

σ

^b

σ

^b

σ

^c

 ,

we allocatethediagonalcomponents,

σ

^a and

σ

^c,tothecellcenterandtheoff-diagonalcomponent

σ

^b tothecellcorner, so that the divergenceform ∇ · (^H

σ

) iscompatible withthe staggeredgrid associated withfluid velocitiesasshown in Fig.2.Thosestressgridallocationscanbefoundinthereference[8].Besides,theindicatorfunction H isdefinedatthecell center,sowhenavariableisnecessarilyshiftedbyahalf-meshwidthforrelevantcomputations,thelinearinterpolationis simplyemployed.Forexample,wedefine

Fig. 2. Grid allocations for fluid variables u,v,p, the extra stressσ, and the indicator function H .

σ

^b

i+¹₂,j

= σ

^b

i+¹₂,j−¹₂

+ σ

^b

i+¹₂,j+¹₂

2 and

σ

_i^b_,j

= σ

^b

i−¹₂,j−¹₂

+ σ

^b

i−¹₂,j+¹₂

+ σ

^b

i+¹₂,j−¹₂

+ σ

^b

i+¹₂,j+¹₂

4

.

Notethatwedefinea variableatthecellcenterby

σ

_i^b_,j=

σ

^b(xi,yi)=

σ

^b(x0+ (ⁱ−¹/2)h, y0+ (^j−¹/2)h)duetotheuse ofstaggeredgrid.Inthesimilarfashion,othervariablescanbelinearlyinterpolatedfromthevaluesofneighboringpoints.

FortheclosedvesicleboundaryX,weuseFourierrepresentationtodiscretizeitas

X

( α ,

t

) =

N



/2−¹ k=−N/2

X

(

k

,

t

)

e^ikα

.

(21)

SothevesicleboundaryisrepresentedbyasetofLagrangianmarkersX_=^X(

α

),where

α

= 

α

and

α

=²

π

/N.The associatedgeometricquantitiesusedintheaboveformulationssuchasinterfacetangentandnormalvector,thecurvature, andtheirderivativesareallcomputedattheLagrangianmarkersXusingtheFourierspectraldifferentiations[42]. These spectralderivativeswithrespectto

α

canbecomputedefficientlyusingFastFourierTransform(FFT)inaspectralaccuracy.

Toremovethealiasingerrorincomputations,we adoptthede-aliasing2/3-rulefilter.Duetotheinextensibilityofvesicle boundary, the Lagrangian markers clustering during the simulations isnot that severe comparing withthe droplet case.

Therefore,wedonotemploythemarkersredistributiontechniqueinpresentsimulations.

WenowpresenthowtomarchtheLagrangianmarkersXⁿ=^X(nt)fromtimeleveln to Xⁿ⁺¹=^X(nt+ ^t)attime leveln+^{1 withthetimestepsize} t.Inthefollowing,thepolymerstress

σ

ⁿ,thefluidvelocity uⁿ,thepressure pⁿ,and theLagrangianmarkersXⁿ areallgiveninadvance,andfromthesevariablesweaimtoupdate

σ

ⁿ⁺¹,uⁿ⁺¹,pⁿ⁺¹,andXⁿ⁺¹. Thestep-by-stepnumericalprocedurecanbedoneasfollows.

Step 1. AttheLagrangianmarkersXⁿ_,wefirstcompute allthegeometricquantitiesusingthespectral differentiationsand thespring-liketensionby

γ

_ⁿ

₌ γ (

Xⁿ_

) = γ

0

 |(

^Xα

)

ⁿ_

|

|(

^Xα

)

⁰_

| −

¹

 ,

andthenincorporatethosetocomputetheinterfacialforces F_γ

(

Xⁿ_

) =

¹

|(

X_α

)

ⁿ_

|

 ∂ γ _τ

∂ α



n



,

F_b

(

Xⁿ_

) =

¹ Ca



− (

X_α

)

ⁿ_

· (

^Xαα

)

ⁿ_

( κ

α

)

ⁿ_

|(

X_α

)

ⁿ_

|

⁴

+ ( κ

αα

)

ⁿ_

|(

X_α

)

ⁿ_

|

²

+ ( κ

_ⁿ

)

³ 2



nⁿ_

.

Step 2. DistributethetensionandbendingforcesactingonLagrangianmarkersintotheEuleriangridbyusingthesmoothed Diracdeltafunctionδhas

fⁿ

(

x

) =

N−¹



=0

(

F_γ

(

Xⁿ_

) +

Fb

(

Xⁿ_

))δ

h

(

x

−

Xⁿ_

)

ds

(

Xⁿ_

),

wherex= (^x,y)istheEuleriangridpointandthearc-lengthelementisobtainedbyds(Xⁿ_)= |(Xα)ⁿ_|

α

.Forthe discrete deltafunction δh(x)= _h¹2φ_x

h

φy h

, we employ the 6-point supported C³ function φ newly developed in [2].

Step 3. Solve theNavier–Stokes equationsbythe second-orderincremental pressure-correction projection methodin[12]

asfollows.

Re



3u^

−

4uⁿ

+

uⁿ⁻¹

2



t

+

2

uⁿ

· ∇

h



uⁿ

−

uⁿ⁻¹

· ∇

h



uⁿ⁻¹



= − ∇

hpⁿ

+ λ

hu^

+ ∇

h

· 

(

1

− λ)

^Hn



∇

huⁿ

+ (∇

huⁿ

)

^T

 + ∇

h

· (

^Hn

σ

ⁿ

) +

^fn

,



_hp^

=

^3Re

2



t

∇

h

·

^u

, ∂

p^

∂

n

=

^{0 on}

∂ 

D

,

u^

=

^un+1on

∂ 

D

,

uⁿ⁺¹

=

^u

−

²



t

3Re

∇

hp^

, ∇

hpⁿ⁺¹

= ∇

hp^

+ ∇

hpⁿ

−

²

λ

t

3Re



h

(∇

hp^

),

wherethediscreteoperators∇h and∇h·approximatethegradientanddivergenceoperators,respectively,usingthe second-order finitedifferenceinstaggeredgrid.Forthenonlinearterms,theskew-symmetricformisemployedas (u· ∇h)u=¹₂(u· ∇h)u+¹₂∇h(uu).TheDirichletboundaryisdenotedby∂D.ThesmoothedindicatorfunctionHⁿ canbeobtainedbysolvingthediscretizedPoissonequation,asdescribedin[43],



hHⁿ

(

x

) = ∇

h

·

N



−¹

=⁰

nⁿ_

δ(

x

−

Xⁿ_

)

ds

(

Xⁿ_

).

The numerical schemeforviscous termwritten above issuggestedwhen λ≥^{1, so thatthe more viscous partis} treatedimplicitlywiththe constant coefficientλ whilethelessviscous partis treatedexplicitly.Therefore,a fast solver using FFT can be implemented directly. Notice that, as the traditional IB method, the interfacial force fⁿ is computed explicitly inStep 2 so herethe extra stress termis also computedexplicitly for thesimplicity and compatibility.

Step 4. Once uⁿ⁺¹ is determined in the Eulerian grid points, we interpolate the fluid interfacial velocity Uⁿ_⁺¹ =

xuⁿ⁺¹(x)δh(x−^Xn_)h².Therefore,theLagrangianmarkerscanbeupdatedbyXⁿ_⁺¹=^Xn_+ ^tUn_^+1. Step 5. Asthelaststep,wesolvetheextrastressequation

∂ σ

∂

t

= F( σ ,

u

),

where

F( σ ,

u

) = −∇

h

· (

^u

σ ) + 

∇

hu

σ ₊ σ ( ∇

hu

)

^T

 − σ

W i

+ β

W i

 ∇

hu

+ (∇

hu

)

^T

 .

Following in[8],we use thesecond-order Runge–Kutta(RK2) methodwithzero-Neumannboundary condition at Dirichletboundaryoffluidvelocitytoupdate

σ

ⁿ⁺¹.Introducingtwointermediatevariables

σ

1 and

σ

2 by

σ

1

= σ

ⁿ

_{+ }

t

F( σ

ⁿ

,

uⁿ

), σ

2

= σ

1

+ 

^t

F( σ

1

,

uⁿ⁺¹

),

the updatedextrastress isobtainedby

σ

ⁿ⁺¹= (

σ

ⁿ₊

σ

2)/2.InthiscomputationusingRK2method,we notethat all the componentsof

σ

arelocated atcell centersby linearinterpolation,andthe convectiontermis writtenin divergenceformduetotheidentityu· ∇h

σ

= ∇h· (^u

σ

).

4. Codeverificationandvalidation

Wetesttheconvergenceratesofsingle-phaseOldroyd-Bfluidsolverandimmersedboundarymethodforvesicledynam- icsintwo-phaseN/Ofluid.ThecodevalidationisalsoperformedbycomparisonofdropletdeformationinN/Ofluid.

4.1. Convergencestudyforsingle-phaseOldroyd-Bfluid

Here,wetestthesingle-phaseOldroyd-Bfluidsolverpresentedintheprevioussectionbycheckingthenumericalaccu- racyofthefollowinganalyticalsolutioninacomputationaldomain= [⁰,2

π

_{]× [}

π

/2,5

π

/2]^.Thesingle-phaseOldroyd-B fluidequationsarewrittenas

Re

 ∂

u

∂

t

+ (

^u

· ∇)

^u



= −∇

^p

+ 

^u

+ ∇ · σ

in

,

(22)

∇ ·

^u

=

^{0 in}

,

(23)

W i

σ

^

+ σ = β D(

^u

) + φ

ⁱⁿ

,

(24)

whereanauxiliaryfunctionφisaddedinEq. (24) anddefinedby

Table 1

TemporalconvergenceanalysisinL_∞normforthefluidvelocityu= (^u,v),pressurep,andthepolymerstressσattwodifferenttimest=⁰.6 andt=¹.2.

M ^uM−^ue ∞ Rate ^vM−^ve ∞ Rate ^pM−^pe ∞ Rate

(t=⁰.6)

64 1.287E-02 – 1.173E-02 – 4.270E-02 –

128 6.291E-03 1.03 5.735E-03 1.03 2.181E-02 0.96

256 3.108E-03 1.01 2.833E-03 1.01 1.100E-02 0.98

512 1.544E-03 1.00 1.408E-03 1.00 5.522E-03 0.99

(t=¹.2)

64 9.512E-03 – 8.290E-03 – 3.775E-02 –

128 4.697E-03 1.01 4.091E-03 1.01 1.935E-02 0.96

256 2.334E-03 1.00 2.033E-03 1.00 9.780E-03 0.98

512 1.163E-03 1.00 1.013E-03 1.00 4.914E-03 0.99

M σ_M^a−σe^a ∞ Rate σ_M^b−σe^b ∞ Rate σ_M^c−σe^c ∞ Rate

(t=⁰.6)

64 8.391E-02 – 1.035E-02 – 8.508E-02 –

128 4.697E-02 0.99 5.286E-03 0.96 4.234E-02 1.00

256 2.099E-02 1.00 2.661E-03 0.99 2.108E-02 1.00

512 1.049E-02 1.00 1.334E-03 0.99 1.052E-02 1.00

(t=¹.2)

64 4.583E-02 – 1.765E-02 – 4.660E-02 –

128 2.329E-02 0.97 8.937E-03 0.98 2.351E-02 0.98

256 1.173E-02 0.98 4.486E-03 0.99 1.179E-02 0.99

512 5.887E-03 0.99 2.247E-03 0.99 5.903E-03 0.99

φ =

^{W i}

β

 φ

^a

φ

^b

φ

^b

φ

^c



,

G

(

t

) =

^{e−2(1+β)t/Re}

,

φ

^a

= −

^2Gsin x sin y

−

^{2G2cos2x sin2y}

+

^{2G2sin2x cos2y}

−

^{4G2sin2x sin2y}

, φ

^b

= −

^4G2sin x cos x sin y cos y

,

φ

^c

=

^2Gsin x sin y

+

^{2G2cos2x sin2y}

−

^{2G2sin2x cos2y}

−

^{4G2sin2x sin2y}

.

Theprime symbolof G denotesthe derivativewithrespectto timet. Aftersome carefulcalculationsreferringto Taylor–

Greenvortex,onecaneasilycheckthatthefollowinganalyticalsolutionsatisfiestheaboveequations(22)–(24).

u_e

= (

G cos x sin y

, −

G sin x cos y

),

p_e

= −

^Re

4

(

cos 2x

+

^{cos 2 y}

)

G²

,

(25)

σ

_e

_{= β D(}

u

) = β

 −

2G sin x sin y 0 0 2G sin x sin y

 .

Table1showstheconvergenceanalysisforallsolutionvariablesattwodifferenttimest=⁰.6 andt=¹.2.Wesetthe parametersbyRe=^{W i}= β =^{1 and}t=^h/8,whereh istheEulerianmeshwidthdefinedbyh=²

π

/M withgridsizeM.

Theconvergencerateofthefluidvelocityu isobtainedby Rate

=

^log2

(

^uM

−

^ue

∞

/

^u2M

−

^ue

∞

),

andso are theother variables. The numericalresults for all variables in Table1 show clean first-order accuracy dueto theexplicittreatmentofthepolymerstressterm∇ ·

σ

ofEq. (22) whichisdescribedinStep3.Ingeneral, theimmersed boundarymethodisfirst-orderaccurate,unlessspecialalgorithmisemployedforhigher-ordersuchasinimmersedinterface method.Thus,thisaccuracybehaviorisacceptableandconsistentwiththepresentIBmethod.Onemightwanttofurther extendour methodto achieve thesecond-order accuracy insolving theOldroyd-B fluid flow by treating theextra stress terminthesamewayasthenonlinearadvectionterminprojection methodofStep3inSubsection3.2.However,we do notadoptthistreatmenthereduetothepresentimmersedboundarymethodisonlyfirst-orderaccuratewhenthevesicle boundaryforceistakenintoaccount.Wefurtherverifythespatialconvergenceofourfluidsolver.Todothis,smallertime stepsizeischosenbyt=^h/128 sothatthedominanterrormainlyarisesfromthespace,whileotherparametersremain thesame.Table2confirmsthatoursolverhasthesecondorderaccuracyinL2 norm.Thesolutionvariablesarealldefined onstaggered gridso whenthemesh isrefined,thesolution locationsatcoarseandfinegrids donot coincidewitheach other.Hence,theusageofL2 errorinsteadofL_∞ isanaturalchoice.

4.2. ConvergencestudyforNewtonianvesicledynamicsinOldroyd-Bfluid

Inthissubsection,we performtheconvergencestudytoverifythepresentnumericalalgorithmforvesicledynamicsin two-phase flow. We puta Newtonian vesicle into Oldroyd-B fluidunder shearflow. Wechoose the vesiclewith

ν

=⁰.8

Table 2

SpatialconvergenceanalysisinL2normforthefluidvelocityu= (^u,v),pressurep,andthepolymerstressσ attwodifferenttimest=⁰.6 andt=¹.2.

M ^uM−^ue 2 Rate ^vM−^ve 2 Rate ^pM−^pe 2 Rate

(t=⁰.6)

64 2.958E-04 – 2.749E-04 – 5.352E-04 –

128 6.366E-05 2.22 5.981E-05 2.20 1.253E-04 2.09

256 1.466E-05 2.12 1.384E-05 2.11 3.042E-05 2.04

512 3.510E-06 2.06 3.324E-06 2.06 7.501E-06 2.02

(t=¹.2)

64 1.628E-04 – 1.483E-04 – 4.084E-04 –

128 3.982E-05 2.03 3.680E-05 2.01 1.028E-04 1.99

256 9.854E-06 2.01 9.175E-06 2.00 2.587E-05 1.99

512 2.451E-06 2.01 2.291E-06 2.00 6.490E-06 2.00

M σ_M^a−σe^a 2 Rate σ_M^b−σe^b 2 Rate σ_M^c−σe^c 2 Rate

(t=⁰.6)

64 9.637E-04 – 1.625E-04 – 9.563E-04 –

128 2.295E-04 2.07 3.554E-05 2.19 2.283E-04 2.07

256 5.597E-05 2.04 8.275E-06 2.10 5.575E-05 2.03

512 1.381E-05 2.02 1.994E-06 2.05 1.377E-05 2.02

(t=¹.2)

64 4.004E-04 – 1.621E-04 – 3.898E-04 –

128 1.012E-04 1.98 3.778E-05 2.10 9.920E-05 1.97

256 2.545E-05 1.99 9.097E-06 2.05 2.504E-05 1.99

512 6.383E-06 2.00 2.230E-06 2.03 6.291E-06 1.99

Fig. 3. The convergence rates for the fluid velocity(u,v)and the vesicle configuration X. Left: M=64; Right: M=^128.

and the dimensionless parameters are fixed by Ca=^{1 and W i}=^{1. We shall choose four different Cartesian grid sizes} M=⁶⁴,128,256,512 toexamineconvergencerates.Undershearflow(u,v)= (^y,0)inthecomputationaldomain[−⁸,8]^2, we settheparameters by Re=^10−3,λ=^2,β=^1,h=¹⁶/M,andt=^h/128.Forthevesicle interfacediscretization,as we doublethegridsizeM, weneedtohalve thelocalarclengthsothat thenumberofLagrangianmarkersusedmustbe increasedby afactoroftwoaccordingly.Tobe consistent,oneshould alsobe carefultoadjustthestiffnessparameter

γ

0 used inEq. (17) for interfacial inextensibility. Forthis,we choose

γ

0=⁵×^103M/512,soit becomes

γ

0=⁵×^{103 when} M=^{512 whichisusedthroughoutSection5.}

Fig. 3showstheconvergenceratesforthevelocity componentsu= (^u,v) andthe vesicleconfigurationX atdifferent timesuptot=^{50.Unlikeintheprevious} subsection,theanalyticsolutioninthistestisnotavailable,sowecomputethe convergenceratebyestimatingthetwosuccessivemaximumerrorsas

Rate

=

^log2

(

^uM

−

^u2M

_∞

/

^u2M

−

^u4M

_∞

).

Theratesforv andX aredefinedinasimilarmanner.Wecanseefromthefigurethattheratesarelargerthanonewhen M=⁶⁴,128,256 andcloseto onewhen M=¹²⁸,256,512,so theaveragerateisaround one forall variables indicating that ournumerical scheme isroughly first-orderaccurate. This resultis consistent withthe IB method beingfirst-order accurategenerally.

Fig. 4. DeformationparameterD ofdropletinN/OfluidasafunctionoftheWeissenbergnumberW i fortwoCapillarynumbers:(a)Ca=⁰.1,(b)Ca=⁰.2.

OurresultiscomparedwithnumericalresultofYueetal.(2005),andalsowithexperimentaldataofGuidoetal.(2003)andofMaffettoneandGreco (2004).

4.3. ComparisonofdropletdeformationinN/Ofluid

Intheprevioussubsection,wehaveshowntheconvergenceratesoffluidvariablesforaNewtonianvesicleinOldroyd-B fluid under shear flow. In order to verify that the present numericalscheme is physically plausible, we hereperform a comparisonwithexperimental andnumericalstudies inliterature fordropletdeformation inN/Ofluid undershear flow.

Inthistest,we denoteN/Ofluid byaNewtonian dropletsuspendedinan Oldroyd-Bfluid.Todothis, inEq. (3),we omit thebendingforceandreplacethetensionforceby Fγ(

α

,t)= −_Ca²_d

κ

n,whereCad isthecapillarynumberassociatedwith dropletinterfacialtensionandthefactor2 comes fromournon-dimensionalizationprocess.Thedeformationparameteris definedbyD= (^L−^B)/(L+^B),whereL andB arethelongestandshortestlengthsofthedropletmeasuredfromitscenter to theinterface, respectively. In thistest, we initially configurea circulardroplet withunity radius whichis centered in thecomputationaldomain[−⁸,8]× [−⁴,4]^{withthegridsize2048}×^{1024,asshowninFig.1.ThenumberofLagrangian} markersis chosenas N=^{4096,theuniformEulerianmeshwidthish}=¹/128,andthetime stepsizeis t=^h/32.We alsofixtheReynoldsnumberby Re=^{10−4andtheviscosityratiobetweenpolymerandNewtoniansolventby}β=^1,thus λ=^2.

Fig.4 showsour numericalresultson the dropletdeformation asa function ofthe Weissenbergnumber W i fortwo different capillary numbers Ca=⁰.1 and 0.2. It is also compared with experimental results [13,28] and the numerical results basedon phase-fieldmethod [46]. Notice that, the experimental results are available only forsmallWeissenberg number.In thefigure,ournumericalresultis quitecomparabletothe experimentbyGuido etal. [13] for smallW i and Ca=0.1,while the result ofYue etal. [46] is comparable to the one of Maffettoneand Greco[28]. Forsmall W i and Ca=⁰.2, ourresult becomes quite comparableto bothexperiments, whilethe resultof Yue etal.is still comparableto theoneofMaffettoneandGreco.IncomparisonwithnumericalresultsofYueetal.,ourresultsshowslightdifferencebut bothplotsofdeformation parameterversusWeissenbergnumbershowconcaveupward behaviorasWeissenbergnumber increases. So we can conclude that our present numerical method is validated by simulatingthe Newtonian droplet in Oldroyd-B fluid.Inthenext section,westudyonthevesicle dynamicsmainlyintwo-phaseN/Ofluidundersteadyshear flow.

5. Numericalresultsanddiscussion

Inthissection,we performnumericaltestsfortwo-dimensional vesicledynamicsunderviscoelasticshearflow. Toex- ploretheeffectsofviscoelasticityonvesiclemotions,theinnerfluidisalwayssettobeNewtonian,whiletheouterfluidis givenbyeitherNewtonianornon-Newtonian(Oldroyd-B).Wedenotethesetwo-phasefluidsbyN/NandN/O,respectively.

Inthefollowingsubsections,wefirstvalidate ourIBmethodbycomparingwiththeKeller–Skalak(KS)theoryinN/Nfluid [21]. We then intensively investigatethe viscoelastic effecton vesicle motions inN/O settingby varying the parameters suchas W i (theWeissenbergnumber),β (theviscosityratiobetweenpolymer andNewtonian solvent), Re (theReynolds number),andλ(theviscositycontrastbetweeninnerandouterfluids).

Throughoutthispaper,weinitiallyconfigureanellipticalvesiclewiththedesiredreducedarea

ν

=⁴

π

A/L² byfixingthe enclosedarea A=

π

andchangingthetotalvesicleperimeterL.Weputthevesicleundershearflow(u,v)= (^y,0)inthe computationaldomain[−⁸,8]^{2 withthegridsize5122,seeFig.1.TheuniformEulerianmeshwidthish}=¹⁶/512 andthe timestep sizeist=^h2/4.ThenumberofLagrangianmarkersrepresentingthevesicleboundaryischosen asN=^1024.

In mostcases,thenumerical resultsare obtainedup to time t=^{50 whichtakesabout30 hours tofinish a run usinga} single-core computerequipped withaCPU3.7 GHz.Unless otherwisestated, we fix theReynolds numberby Re=¹⁰⁻³

Fig. 5. Comparison ofN/N (leftcolumn,(a,c))and N/O(rightcolumn, (b,d))cases. (a)–(b):thetimeevolutionaryplotsofinclinationangleθ/πfor variousν.(c)–(d)thetimeevolutionaryplotsoftank-treadingfrequencyωforvariousν.Hereνdenotesthereducedarea.

to be inStokes regime,the capillarynumberby Ca=^{1,andtheviscosity ratiobetweenpolymerandNewtonian solvent} by β=^{1. Although mostly omitted here, we tested for different Ca ranged from1 to 100 and initial} tilt-angles 0 and

π

/4.Sincethequantitativeresultsassociatedwithshear-inducedvesicleTTmotionsarealmostidenticalfordifferentinitial configurationofvesicle,wetiltthevesicleinitiallyonlybyanangleθ=

π

/4.Inordertosupportthisstatement,inthenext subsection,we show thetime-dependent variationof inclinationangleandangularfrequencyoftank-treadingvesicle for differentconfigurationsinitiallytilted by0,

π

/8,and

π

/4.Thestiffnessparameterusedforthenearlyinextensiblevesicle boundary ischosenas

γ

0=⁵×^{103 withwhichtherelative changesinvesiclearea A anditstotallength L arelessthan} 0.5% duringallsimulationswehaveperformed.

For shear-induced vesicle motions in N/N fluid, two major types of motion have been well investigated; namely (1) tank-treading(TT) motion,and(2)tumbling (TB)motion.A vesiclewithTTmotiontakesastationaryshapeandconstant inclination angle θ independentof time,while a vesiclewith TBmotionperiodically rotateslike a rigidbody.Thesebe- haviors havebeenfound inasymptotictheory [21] andalsonumericalsimulationssuch asin[24].Anothermotioncalled avacillating-breathing(ortrembling,swinging)isfoundbetweenTTandTBstatesinthetheory[29],theexperiment[19], andnumericalsimulation[31].However,hereonlyTTandTBmotionsareobservedsincethevacillating-breathingismainly attributedtothepresenceofvorticityaxisinthreedimensions.

5.1. ComparisonofmatchedviscosityN/NandN/Ofluids

As afirsttest,weexaminetheeffectofviscoelasticitybycomparingthevesicletank-treading(TT)motionundershear flow forN/NandN/Ofluids. Throughoutthispaper,we callthematched viscosityN/Nfluid iftheparameters λ=^{1 and} β=^{0 areused inN/Nfluid setup;while wecall thematchedviscosityN/Ofluid if}λ=¹+ β ^andβ >0 in N/Ofluid. In thelattercase,theviscosityofNewtonianfluidmatcheswiththetotalviscosityofOldroyd-Bfluidby

μ

n=

μ

s+

μ

p.Here wechooseβ=^{1 andthe}WeissenbergnumberW i=^{1 for}convenienceinthisstudy.Notethat,theusageofβ=^{1 isquite} commoninother numericalstudies[7,46,8]. Meanwhile,asreportedin[41],the valueof W i canbemeasured closeto1 orevenhighercloseto10 inbloodflowforhardenedcellsandthroughmicrotubes.