A fast moving least squares approximation with adaptive Lagrangian mesh refinement for large scale immersed boundary simulations

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

A

fast

moving

least

squares

approximation

with

adaptive

Lagrangian

mesh

refinement

for

large

scale

immersed

boundary

simulations

Vamsi Spandan

a

,

∗

,

Detlef Lohse

a

,

b

,

Marco

D. de Tullio

c

,

Roberto Verzicco

a

,

d

,

∗

a_{PhysicsofFluidsandMaxPlanck-UniversityofTwenteCenterforComplexFluidDynamics,UniversityofTwente,Enschede,POBox217,} 7500 AE,Netherlands

b_{MaxPlanckInstituteforDynamicsandSelf-Organization,37077,Göttingen,Germany} c_{PolitecnicodiBari,ViaReDavid200,70125,Bari,Italy}

d_{UniversityofRome‘TorVergata’,ViadelPolitecnico,Rome00133,Italy}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received1December2017

Receivedinrevisedform2August2018 Accepted23August2018

Availableonline28August2018

Keywords:

Immersedboundarymethod Movingleastsquares Multiphaseflows

In this paper we propose and test the validity of simple and easy-to-implement algorithms within the immersed boundary framework geared towards large scale simulations involving thousands of deformable bodies in highly turbulent flows. First, we introduce a fast moving least squares (fast-MLS) approximation technique with which we speed up the process of building transfer functions during the simulations which leads to considerable reductions in computational time. We compare the accuracy of the fast-MLS against the exact moving least squares (MLS) for the standard problem of uniform flow over a sphere. In order to overcome the restrictions set by the resolution coupling of the Lagrangian and Eulerian meshes in this particular immersed boundary method, we present an adaptive Lagrangian mesh refinement procedure that is capable of drastically reducing the number of required nodes of the basic Lagrangian mesh when the immersed boundaries can move and deform. Finally, a coarse-grained collision detection algorithm is presented which can detect collision events between several Lagrangian markers residing on separate complex geometries with minimal computational overhead.

©2018 Elsevier Inc. All rights reserved.

1. Introduction

Numericalsimulationsofflowsaroundmovinganddeformablebodieswithcomplexgeometrieshavebeenanimmense challenge to thefluid dynamics communityfora longtime. Computationalstudies usingbody-fitted methodswhere the

mesh needs to be adjusted according to the instantaneous shape andposition of the immersed body become not only

mathematically complex but also terminally expensive [1]. This makes non body-fitted methods such asthe immersed

boundary method (IBM) very attractive to deal with complex flow systems especially when there is coupling between

a carrierfluid anda deformable interface, membrane or a structure [2]. IBM uses a single time-invariant mesh for the solutionoftheNavier–Stokes(NS)equationsanda separatesetofnodesdistributedonthesurface ofanimmersedbody. The influenceofan immersedrigidordeformablebodyonthesurroundingfluidisachievedthroughaspatially andtime varyingforcingfunctionintheNSequations.ThismethodwasfirstusedbyPeskin[3] tostudyflowaroundheartvalvesand

*

Correspondingauthors.

E-mailaddress:vamsispandan@gmail.com(V. Spandan).

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

lationsinvolving IBM(the termlargescale herearbitrarily refers tonumericalsimulations involving O

(10

)

fullyresolved deformablebodieswithO(103

)

LagrangianmarkerseachinaturbulentflowonaEulerianmeshofO

(10

9

)

ormorenodes). Inthispaper,we tacklesome ofthesebottlenecks andintroducesimple andeffectiveoptimisation algorithmswithin the immersedboundaryframeworkwhicharegearedtowardsscalingupsimulationsofhighly turbulentflowswiththousands offullyresolveddeformableparticles.

As mentioned before,IBM requirestwo different meshes;an Eulerianone for thesolution ofNS equationsanda La-grangianmeshonthesurfaceoftheimmersedbodywhichisprimarilyusedtoenforcetheinterfacialboundarycondition. When theimmersedbodiescan deform,theLagrangian meshalsoserves asthe computationalgrid forthe discretisation ofthestructureequations;suchastheCauchy–Navierequationsinelasticdeformationproblems.Whiletherearemultiple variantsofIBMwhichvary inhowtheinfluenceoftheimmersedbodyis imposedontothe carrierfluid, themost pop-ularapproach inthecaseofmoving bodiesisthe oneintroduced by [12].In thisapproach,theforcing functionon each Lagrangiannode iscomputedfromacloud ofsurroundingnodesoftheEulerianmesh.SincetheLagrangiannodesdonot coincidewiththeEulerianmeshateverytimeinstant,interpolationandextrapolationalgorithmsareimplementedintothe flowsolver;thisbringsintoplayacoupleofissues.Firstly,itisimportanttoensurethatthetransferfunctionsusedto

ex-changeinformationbetweentheEulerianandLagrangianmeshesconservemomentumandenergybetweentheimmersed

body and thesurrounding fluid in addition to maintaininga specific order of accuracycompatible withthe flow solver. Secondly,thecalculationofthetransferfunctionshould notbecomecomputationallyintensiveasthishamperstheability oftheflowsolvertohandlelargescale systems.AnothermajorissueinnumericalsimulationsinvolvingIBMisthecoupling ordependenceoftheresolutionoftheLagrangianmeshwiththeEulerianmesh.Accordingtothestandardimplementation theremustbe,atleast,one-to-onecorrespondenceofEulerianandLagrangiannodesifwehavetopreventholes orflux leak-age attheimmersedsurface.Inotherwords,havingaLagrangianresolutioncoarserthantheEuleriancounterpartresultsin adiscontinuousrepresentationoftheimmersedboundarywhichleadstoartificialeffectsintheflow.Typically,the accept-ablechoiceisthatthelengthscaleofanindividualLagrangianelementissetto0.7–0.9timesthatofthelocalEuleriangrid size.However,itisimportanttonotethatone-to-onecorrespondenceofEulerianandLagrangiannodesbecomesterminally expensiveinwall-boundedflowswherethemeshisrefinedinthenearwallregionandthisforcestheLagrangianmeshto beunnecessarilyrefinedwhenthemovingimmersedbodyevolvesinthebulkregion.

Inrecentworks, themovingleastsquares(MLS)approximationhasbeenusedtobuild thetransferfunctionsbetween theEulerianandLagrangianmeshesandinparticularithasbeenfavouredforthesimulationofdeformableinterfacesand membranes[8,9,13].MLSallowsforinterpolatingboththevalueandthespatialderivativesofaflowvariableatanygiven pointbasedona setoffieldnodesdefinedwithin alocalsupportdomaininadditiontoensuring thetransferfunctionis smoothinpresenceofmovinganddeformingboundaries.MLSwasfirstintroducedbymathematicianstosolvetheproblem ofdatafitting andsurface reconstruction andgivenits highly genericnature ithasbeen usedinseveralfieldsthereafter [8,9,13–23].WhileMLSprovidesseveraladvantagesinhandlingmovinganddeformablebodiesininhomogeneousunsteady flows, there are certain inherent issues in the algorithm which make it difficult to use. In comparison to conventional interpolation techniques,MLSis a computationallyintensive algorithmas itrequires theconstruction ofmultiple weight matricesanda matrixinversionforeveryLagrangian markerwheretheapproximationisrequired.Thisseverelyhampers thewall-clock time ofthesimulationwhenthereare O

(10

6–107

)

Lagrangianelements immersedintheflow. Inthefirst partofthispaper,weintroduce asimplealgorithm(referredtoasfast-MLS)wherewecompute theMLSshapefunctions only once at the beginning of the simulation at pre-defined nodes and then use a relatively inexpensive interpolation technique (Shepard’sinterpolationhere) toestimate theshape function atanyposition.Withthiswe aimtoachieve the smoothnessandaccuracyofMLSshapefunctionswhileeliminatingtheoperationcostsofconstructingandinvertingweight matricesasintheexact-MLSapproach.Wecomparetheresultsandthetotalcomputationaltimefrombothapproachesfor thewell-studiedproblemofuniformflowoverastationarysphere.

Asmentionedearlier,duetothecouplednatureoftheEulerianandLagrangianmeshesinIBMsimulations,anyincrease intheEulerianmeshforcesanincreaseintheLagrangianmesh.Thiscanseverelyhampertheperformanceofthecodewhen theimmersedbodies aremoving andtheLagrangian meshhasto beevolvedindifferentregions ofthe flowwhichmay haverelativelycoarseorrefinedEulerianelements.Inordertotacklethisproblem,weproposesegregatingtheLagrangian

mesh into two parts (i) base mesh which can be used to deform the immersed body (ii) refined mesh forcomputing

onlyonceduringthestartofthesimulationwhiletherefinedmeshisdefineddynamicallyaccordingtothelocalEulerian resolution. This approach not only makes the flow solver more flexible since the deformation governing equations can now be discretisedindependentlyof theEulerian mesh,butitalso reducesthe memoryloadofthe codeaswell asthe computationalcosts.Totestthevalidityoftheproposedapproach,wecompareresultsofadeformingdropinacross-flow by bothfullyresolving theLagrangianmeshaccordingtotheEulerianresolutionandalsousingtherefinementtechnique whichrequiresrelativelylessoperationsandmemory.

Finally,whensimulatingthousandsofdeformablebodiesdispersedinaturbulentflow,detectingandmodellingcollision betweenthembecomesachallengingcomputationaltask.Tothisextent,inthelastpartofthepaper,wediscussasimple andeasy-to-implementcoarsegrainedcollisiondetectionalgorithmwhich canbeusedto detectcollisioneventsbetween markersresidingondifferentdeformableimmersedbodieswithminimalcomputationaloverhead.Itisworthtostressthat, ourprimary focusisoncollisiondetectionbetweendeformablebodieswithin theimmersedboundaryframeworkandnot onthespecificphysicalcollisionmodellingthatcanbeparametrisedbyseparatemodels.

The paper is organised as follows. In the next section we describe the fast-MLS approach and compare the results from thisapproach withthe conventionalMLS-IBM simulations. Insection 3, we describe theadaptive Lagrangian mesh refinementtechniqueandcomparetheresultsfromthistechniquewithpreviousworks.Insection4,wediscussthedetails andresultsfromthecoarse-grainedcollisiondetectionalgorithmandthenprovideasummaryandoutlookinsection5. 2. Fastmovingleastsquaresfortheimmersedboundarymethod

First,wedescribebrieflytheMLSapproximationusedinthecontextofIBM[8,9,13].InMLS,afieldvariableq givenon apre-definedsetofnodesisapproximatedatanydesiredpointlyingwithinthesenodesandpositionedatxl asfollows:

Q

(

xl

)

=

pT

(

xl

)

a

(

xl

)

(1)

Q istheapproximatedquantitywhilepT

(

xl

)

isabasisfunctionvectorwithdimensionm.xl isthepositionvectorofthe

La-grangianmarkerl.pT

₍

_x

l

)

= [

p1

(

xl

),

...,

pm

(

xl

)

]

isbuiltusingthePascal’striangleandpyramidintwoandthree-dimensions,

respectively. Inthiswork, we considera linearbasis functionwithpT

₍

_x

₎

_{= [}

_1,x,y,z

_]

_,i.e.m

₌

_{4 which hasalreadybeen}

showntobeacost-efficientandreliablechoiceforflowsimulationsincorporatingtheimmersedboundarymethod[8,9,13]. Itistoberememberedthatthefast-MLSapproachdescribedlaterisindependentofthechoiceofm.

Thevectora

(

xl

)

iscomputedbyminimisingaweightedL2norm J whichisdefinedas: J

=

Ne



k=1

W

(

xl

−

xk

)

[

pT

(

xk

)

a

(

xl

)

−

qk

]

2 (2)

Ne is the number of nodes in the Eulerian grid which are in support of the Lagrangian marker l and the selection of

these nodesis described later. W

(

xl

−

xk

)

isa weight function forpoint xk andcan be computedindifferent ways;for

exampleusingexponentialfunctions,cubicsplines,quadraticsplinesetc.[24].Thechoiceoftheweightfunctionisarbitrary and primarily depends on the desired levels of accuracy provided the conditions for continuity and smoothness of the weightfunctionwithinthesupportdomainaresatisfied.Asplineweightfunctionwithadesiredorderofcontinuitycanbe constructedusingthemethoddescribedin[25].Minimising J withrespecttoa

(

xl

)

leadstotherelationA

(

x

)

a

(

x

)

=

B

(

x

)

q,

where q is a vector containing the valuesof the field variable q onthe pre-defined setof nodes; A

(

x

),

B

(

x

)

andq are definedas: A

(

xl

)

=

Ne



k=1 W

(

xl

−

xk

)

p

(

xk

)

pT

(

xk

)

(3) B

(

xl

)

= [

W

(

xl

−

x1

)

pT

(

x1

)...W

(

xl

−

xNe

)

p T

₍

_x Ne

)

]

(4) q

= [

q1

...q

Ne

]

T ₍₅₎

Minimisingequation(2) withrespecttoa

(

x

)

givesthecoefficientvectora

(

x

)

=

A−1

(

x

)

B

(

x

)

q andtheapproximatedvalueof thefieldvariableq atxl canbecomputedasfollows:

Q

=

pT

(

x

)

A−1

(

x

)

B

(

x

)

q (6)

Theapproximationinequation(6) canalsobewrittenasQ

= 

T

₍

_x

₎

_q,where

_

T

₍

_x

₎

_{isavectorofshapefunctioncoefficients}

ofsize1

×

Ne andiscomputedas



T

(

x

)

=

pT

(

x

)

A−1

(

x

)

B

(

x

).

ThisshapefunctioniscomputedforeveryLagrangianmarker

intheflowandisusedtoapproximatethevalue oftherequiredEulerianquantity(typicallyfluidvelocity)attheposition ofthemarker.Approximationofthefluidvelocitygivesinformationonthemagnitudeoftheforcingthatneedstobe trans-ferred totheEuleriangridwhichensures thattheinterfacial boundaryconditionissatisfied.Ifthevectorq

= [

u1

...u

Ne

]

T_,

uk being thefluid velocity definedatthe Eulerian nodesandits corresponding approximatedquantity Q

=

U ,the total

IBM force corresponding tothe Lagrangian markerl is written as F_li

= (

V_bi

−

Ui

_)/t,

_{where F}i l, V

i

Fig. 1. (a)Schematicofatwo-dimensionalsupportdomainoftheLagrangianmarker(shownbyredtriangle);dottedlinesshowalltheEuleriancellcentres insupportoftheLagrangianmarker.(b)HomeEuleriancelloftheLagrangianmarkerdiscretisedbyNi=9 and(c)Ni=25 sub-Euleriannodes,respectively.

(Forinterpretationofthecoloursinthefigure(s),thereaderisreferredtothewebversionofthisarticle.)

force,desiredvelocityofthefluidattheLagrangianmarkerandtheinterpolatedfluid velocityintheith direction.Inthis way,the actualvelocity ofthe fluid atthemarkeris forcedto beequal tothe desiredvelocity Vb,in orderto imposea

no-slipcondition.Similarstrategiescanbe usedtoimposeafree-slipboundarycondition [26].Oncetheforcinghasbeen determined,thenextstepistodistributeitontotheEuleriangridinsuchawaythatthetotalforce duringthetransferis conservedwhichleadstothefollowingrelation:

Ne



k=1 fk

V

k

=

Nl



l=1 Fl

V

l

,

(7)

V

k isthe volumeofthe Euleriancellk and

V

l

=

Alhl isthe forcing volumeassociatedwiththe Lagrangianmarkerl;

Al isthesurfaceareaoftheLagrangianmarkerl andhl

=

Ne



k=1

φ

_kl

(x

k

+ 

yk

+ 

zk

)/3.

Oncetheshapefunctioniscomputed,

theforcingtobe includedintheEuleriancellk iscomputedasfk

=

Nl



l=1

cl

φ

lkFl.Thescalingfactorcl iscalculatedsuchthat

equation(7) issatisfiedandisgivenasfollows: cl

=

V

l Ne



k=1

φ

l_k

V

k (8)

φ

l_kareindividualelementsofthevector



T

₍

_x

₎

_{andrepresentstheshapefunctionvaluefortheEuleriancellk corresponding}

totheLagrangianmarkerl.TheconstructionoftheauxiliarymatricesA

(

x

)

andB

(

x

)

followedbyaninversionofA

(

x

)

arethe computationallymostintensivepartsoftheMLSapproximation.Forstationarybodies,onlyasingleMLSinterpolationper Lagrangianmarkeratthebeginningofthesimulationisrequired.However,whentheimmersedbodiesmoveordeform,in additiontoperformingtheMLSinterpolationevery timestep,wealsoneedtoevaluatehydrodynamicforces,interpolating pressureandvelocity gradientcomponentsatthe markerposition,whichleads toadditionalMLSoperations. Sincethese stepshavetobeperformedforeveryindividualLagrangianmarker,itresultsinasteepincreaseofthecomputationaltime withincreasingLagrangianresolution.Inthefast-MLSapproachweeliminatetheprocess ofcomputingtheshapefunction

(

x)during thesimulationandconstructitonlyatthestartofeverysimulationforapre-definedsetofnodes.Theshape function for anyarbitrary Lagrangian marker is then computed through a relatively inexpensive interpolation procedure suchastri-linearinterpolationortheShepard’sinterpolation.Theprocedureforselectingthenodesusedforcomputingthe shapefunctionsandinterpolationoftheshapefunctionsisnowdetailed.

InFig.1(a)weshow theschematicofatwo-dimensionalsupportdomainused foran arbitraryLagrangian marker.For thesakeofbrevity,inthefollowing,werestricttoatwo-dimensionalsystem;althoughthenumericalexamplesgiveninthe resultsareallthreedimensionalflows.Thesupportdomainforourexact-MLSapproximationinatwo-dimensionalsystem consistsofNe

=

9 (3

×

3)Euleriancells[8,9].ThesupportdomainisbuiltbyfirstidentifyingthehomecelloftheLagrangian

marker(indicatedbythegreyregioninFig.1(a))andthentakingitsneighbouringcells.Infast-MLSwetakeadvantageof thefact thatthe shapefunctionof Lagrangianmarkersatthe samerelativedistancesto theEuleriancellcentres oftheir respectivesupportdomainswillbeexactlythesame.Whilethecalculationsgivenbelowareforauniformmesh,thesame methodcan alsobeapplied whenthemesh isstretchedin asingle direction,the onlydifference beinga requirementto storeadditionalshapefunctionsatthestartofthesimulation.

Toincorporatethefast-MLSapproach,wefirstdiscretiseasingleEuleriancellusinganimaginarymesh(hereaftercalled thesub-Eulerianmesh)beforethestartofthesimulationasshownbyhollowcirclesinFig.1(b);inthiscase,theresolution

Fig. 2. (a)Maximumrelativedifferenceintheshapefunctioncoefficientscomputedfromthefast-MLS(φf)andtheexact-MLS(φe)techniquesversusthe

numberofgridpointsusedinthesub-Eulerianmesh.(b)DragcoefficientCDofastationaryrigidspherefacingauniformflowstreamversustheReynolds

numberRe.Inboth(a)and(b),themarkersizeisinclusiveoftheerror-bar.Referencedatain(b)isfrom[27].

ofthe sub-Eulerianmeshis Ni

=

9 points.Next,we computetheexact shapefunctionsateachofthesepoints usingthe

algorithms described in equations (3)–(6). This leads to nine different exact shape vector functionswhich contain nine coefficientseach.Duringthesimulation,thelocationoftheLagrangianmarkerisidentifiedwithrespecttothesub-Eulerian meshinitshomecellandtheshapefunctionisestimatedasfollows:



f

(

xl

)

=

4



i=1

α

i



ie

(

xi

)

(9)

α

i

=

1/d2_j 4



j=1 1/d2 j (10)



f

(

xl

)

isthe shapefunction estimatedfromthe fast-MLS approach,



e

(

xi

)

is theexact shapefunction calculated atthe

nodesofthesub-Eulerian meshatthestartofthesimulation,xi is thepositionofthenode i onthesub-Eulerian mesh;

α

iaretheinterpolationcoefficientscomputedthroughShepherd’sprocedurewhichisarelativelyinexpensiveinterpolation

technique anddependsonlyondiwhichisthedistancebetweentheLagrangianmarkerandthenodesofthesub-Eulerian

meshinconsideration.ThenumberofEuleriansub-nodescanbefurtherrefinedanddiscretisedusing25points asshown inFig.1(c)inwhichcasetheshapefunctionson25(5

×

5)differentnodeswith9coefficientseachwouldneedtobestored at thestart of thesimulation.The interpolation ofthe shape function accordingtoequation 9wouldthen be computed basedonadifferentsetofnodes.Thenextlevelofrefinementofthesub-Eulerianmeshwouldconsistof81(9

×

9)nodes. Acommentisnecessaryregardingthechoiceofthenumberofsub-Euleriannodesusedinthecalculationof

α

iinequations

(9) and(10).Whileestimatesfromallthesub-Euleriannodescanbeusedinthecalculationof

α

i(forexample9inFig.1(b)

and25in

1

(c)),wefindthatgiventhe1/d2 dependenceofthecoefficients, includingadditionalsub-Euleriannodesinthe calculationof

α

idoesnotimprovetheestimationsignificantly.

We nowcomparetheshapefunctioncoefficientscalculatedusingboth fast-MLSandexact-MLStechniques.InFig.2(a) weshowthenormalisedmaximumdifferenceintheshapefunctioncoefficientsversusthenumberofgridpointsusedinthe sub-Eulerian mesh.Inorderto ensurethat themaximumrelativedifference isnot affectedby positionoftheLagrangian marker, we run the tests for 1000Lagrangian markers randomly positioned in the Eulerian cell. The maximum relative difference issampledover the1000randomlyplaced Lagrangianmarkers.Weseethatwithincreasingsub-Eulerianmesh resolution there is a rapid decrease in the difference between the shape functions computed from the exact-MLS and thefast-MLSapproaches. Whileasub-Eulerianmeshdiscretisedusing9(3

×

3)or25(5

×

5) pointsleadstounacceptable deviationsintheshapefunctioncoefficients,theerrorwith81(9

×

9)pointsisalreadybelow3%.Wenowuseasub-Eulerian meshdiscretisedinthree-dimensionsusing729(9

×

9

×

9)pointsinthefast-MLStechniquetocomputethedragcoefficient ofastationarysphereinauniformflow.Theerrorintheshapefunctioncoefficientsinathreedimensionalsystemfollows thesametrendasshowninFig.2(a),whiletheonlydifferenceistheincreaseinthenumberofNi,thesub-Euleriangrid

nodes. InFig.2(b),we comparethe spheredrag coefficient versus theReynolds numberwhichis definedas Re

=

U D/

ν

(U is the free stream velocity, D is the diameter of the sphere and

ν

is the kinematic viscosity of the fluid). A good agreementcanbeseenbetweenthedrag coefficientsobtainedfromthesimulationsinvolvingthefast-MLSandexact-MLS approach.Thescalingfactorclofequation(8) whichdependsontheshapefunctioncoefficientsestimatedfromthefast-MLS

approach, volumeofLagrangian markers andthecorresponding Euleriancellensures that thereis nolossof momentum whiletransferringinformationbetweentheEulerianandLagrangianmeshes.

InFig.3(a),weshowacomparisonofthecomputationalperformanceoftheflowsolverwhenusingfast-MLSand exact-MLStechniques.The parallelisationalgorithms anddatastructuresused inthissolver havealreadybeen describedin [9]. Thesesimulations wereperformedonthethinnodesoftheDutchsupercomputing facility‘Cartesius’whereeach node is

Fig. 3. (a)Scalingplotshowingthetimestepforonefulliterationofthesolverversusthenumberofcoresusedforboththeexact-MLSandfast-MLS techniques.(b)ComparisonofthetotaltimetakenforonefulliterationversusthetotalnumberofLagrangianmarkersimmersedintheflow.t isnormalised usingthetimetakenforthefirstdatapointatNf=320.Circlescorrespondtoexact-MLSwhilesquaremarkerscorrespondtofast-MLS.

composed of2

×

12 core2.6 GHzIntel XeonE5-2690v3 CPU’s. Forthistest case, weinitialise atotalof16000 spherical particles eachdiscretisedusing320faces (i.e.atotal of5MillionLagrangianmarkers) ina doubly-periodicCartesianbox (discretisedusing anEulerian resolutionof720

×

720

×

3840).The flowis driveninone ofthe periodicdirections using auniformpressuregradient whiletheimmersedparticlecanmoveanddeformbasedonthetechniquesdescribed in[9]. Itcan beseen thatthefast-MLS approachwhileensuring goodscalabilityalsoreducesthe totalwallclock timeper time step.Here,onefulliterationconsistsofintegratingthefluidandimmersedbodygoverningequationsforonetime-step.In ordertofurtheranalysethecomputationalboostprovidedbythefast-MLSapproachwenowfixtheEulerianresolutionand increasethe totalnumberofLagrangianmarkers (Nf) intheflow.Forthesesimulations theEuleriangridiskeptfixedto

120

×

120

×

720 andatotalof120computingprocessorswereused.InFig.3(b)weplotthenon-dimensionaltimetakenfor onefulliterationwithincreasingnumberoffaceswherethetimeisnormalisedusingthetimetakenforthefirstdatapoint i.e.Nf

=

320.Here itisimportanttopoint outthatforthesesimulationstheimmersedbodiescanbothmoveanddeform

accordingtothe deformationalgorithm discussedby [9]. Itimpliesthat withoutthe fast-MLSapproach asingle iteration wouldconsistofmultipleexact-MLSoperations.Thegain incomputationaltime achievedwiththefast-MLSapproachcan beseenclearlyfromFig.3(b).

Acommentisnecessaryonthememoryrequirementsofthesub-Euleriancellapproachwithinthefast-MLSformulation. Thememoryloadismanageablefortworeasons.First,inthecaseofauniformmesh,thesub-Euleriancellapproachtakes advantage of the fact that the shape functions depend entirely on the relative position of the Lagrangian marker with the neighbouring Eulerian cells. Thus, only the shape coefficients of the sub-Eulerian cells of a single Eulerian cell will need to be storedby all processors. Second, inthe caseofnon-uniform grid stretching inone direction, the coefficients ofsub-Eulerian cells inthe directionofthe stretched grid isstored. Although thisincreasesthe memoryrequirementin comparisontoacasewithuniformgridspacing,forasimulationwithmultipleprocessors,eachprocessorwillonlyneedto storetheinformationoftheEuleriancellsallocatedtoitsmemoryandnotofalltheEuleriancellsinthestretcheddirection. ThememoryrequirementforthelargestcasesvariesfromO(1) MBforuniformgridspacingtoO(10) MBfornon-uniform gridspacing,whichiseasilymanageableonmoderncomputingarchitectures.

3. DecouplingtheLagrangianmeshes

WenowmoveontodescribingthetechniqueofadaptiveLagrangianmeshrefinementandhowitcanbe usedin sim-ulations involvingIBM.Thepurposeofthe refinementisto tackletwo challengesinIBMsimulations. Firstly,inMLS-IBM simulations,theresolutionoftheLagrangianmeshisinherentlycoupledtotheEulerianmeshduetothecontinuous trans-ferofinformationbetweenbothmeshes.IftheLagrangianmeshistoocoarse,theEulerianflowfieldisnotforceduniformly whichleadstodiscontinuities(oftencalled‘holes’)intheinterfacialboundarycondition.ThisisillustratedinFig.4,where wesimulateflowoverastationaryrigidsphereintwodifferentcases(i) coarseLagrangianmeshi.e.l/h

∼

2.0,and(ii) re-fined Lagrangian mesh, l/h

∼

0.8, where l is the mean edge length of the triangular elements andh is the mean grid spacingintheEulerianmesh.Itcan beclearly seenfromtherightpanel ofFig.4(a)that whentheLagrangian discretisa-tioniscoarser thanthelocalEulerianone,theflowdevelopsdiscontinuities attheinterfacebetweentheimmersedbody andthefluid.Incontrast,whentheLagrangianmeshissufficientlyfine(asseeninFig.4(b)),i.e.themeanedgelength of thetriangularmarkersis smallerthanthat ofthelocalEuleriangridspacing, theflow issmoothandtheintegrity ofthe interfacebetweentheimmersedbodyandthefluidismaintained.

Tackling the problemof discontinuityin the flow withmoving immersed bodiesusing a brute-force approachwould meanresolvingtheLagrangian meshaccordingtothesmallestEuleriangrid cell,sincethepositionofthesebodies isnot knowninadvanceandthey couldendupinthefinestmeshregions. ALagrangian resolutionasfineasthesmallestlocal Eulerianmesh solvesthis problembutthis leadsto inclusion ofadditionalLagrangian markerswhich becomeredundant

Fig. 4. FlowoverastationaryrigidsphereatRe

=

100 for(a)CoarseLagrangianmesh(b)RefinedLagrangianmesh.Theflowisin

ˆ

ezdirection.Theright

panelsshowazoomoftheflowaroundthespherewhere‘holes’or‘discontinuities’intheflowcanbeeasilydetectedinthecaseofacoarseLagrangian mesh.Thecolormaprepresentsthevelocityinthe

ˆ

eydirectionwhileU istheambientflowstreamin

ˆ

ezdirection.

Fig. 5. TriangulationsfortheadaptiveLagrangianmeshrefinement.(a)initialbasetriangularelement.(b),(c),(d)showthefirst,secondandthirdlevelof refinements,respectively.(e)Aspherewheretherighthalfofthesphereisrefinedusingthethirdlevelwhilethelefthalfstillhasitsinitialbasemesh.

when theimmersedbody evolvesinregions wheretheEulerianmeshiscoarse. Animmediateconsequenceisa massive increase of the total computational time. Secondly, if the immersed body can also deform, the same Lagrangian mesh thatisusedforenforcingtheinterfacialboundarycondition(IBMmesh)isusedforthestructuralsolver.Thisimpliesthata refinementintheIBMmeshautomaticallyimpliesarefinementinthestructuralsolvermeshwhichmayleadtounnecessary computational nodesfor discretisingthe deformation governingequationsandsmall time stepsto trackthedynamics of smallstructuralelements.Here,wefocusondescribingtherefinementprocedureandtheresultsfromthesimulationsofa deformingdropwithandwithoutsuchrefinement.

The basicidea ofthe Lagrangian meshrefinement isshowninFigs. 5(a)–(d). Considera triangularelement asshown in Fig.5(a) whichis composed ofthe threevertices (V1

,

V2

,

V3) anda centroid (represented bya square symbol). Given

the vertices of any triangular element, the position of the centroid can be directly computed as C

= (

V1

+

V2

+

V3

)/3

We canfurtherrefine thetriangularelementinmultiplelevelsasshowninFigs. 5(b)–(d).Inthefirst levelofrefinement (Fig.5(b)),thetriangleissplitintotwosuchthatL1

=

0.5(V1

+

V2

).

Inthenextlevelofrefinement(Fig.5(c))thetriangular

element issplit insuch awaythat L2 is atthecentroidofthe triangleelement i.e.L2

= (

V1

+

V2

+

V3

)/3.

Similarly the

trianglecanbe furtherrefinedasshownintheFig.5(d);L31

=

0.5(V1

+

V2

),

L32

=

0.5(V2

+

V3

),

L33

=

0.5(V3

+

V1

).

Here

it is importantto note that such a refinement procedureis computationally inexpensive andfora given basemesh,the positions oftheverticesandcentroids oftherefined trianglescanbe easilycomputeddynamicallyduringthe simulation. Therefinementproceduredescribedaboveholdswhenthe‘mass’isdistributeduniformlyoverthetriangularnetwork,while

Fig. 6. (a)Snapshots ofadropimmersed inacrossflow attwo timeinstants;topand bottompanelcorrespondstoitsinitialanddeformedstate, respectively.(b)ComparisonoftheinverseofthemeanaspectratioofthedropversustheWebernumberwithandwithouttherefinementtechnique.The referencedataisgivenbyafitfromseveralexperimentalmeasurementsgivenin[28].

incaseon a non-uniformlydistributedmass,appropriate weights haveto be giventothe individual vertices.Thelocally refinedLagrangianmarkerstakepartonlyinenforcingtheinterfacialboundaryconditionbutnotinthedeformationofthe immersedbody,i.e.noadditionalcurvature informationisaddedto thebody.Thus,ifthe initialLagrangian resolutionis alreadysufficienttoaccuratelyrepresentthesurfaceduringdeformation,thedecouplingensuresnoredundantcalculations areperformedforthestructuralelements.

We now test this approach for a deforming drop withits centroidfixed immersedin a cross-flow and compute the

deformation dynamics arising from the resulting flow conditions with and without refinement. The control parameters are theReynolds, Re

=

U deq

/

ν

f andthe Webernumbers, W e

=

ρ

fU2deq

/

σ

,deq isthediameterofthedropin itsinitial

spherical shape. For such a flow, the aspect ratio of the deforming drop depends on the Weber number, which is the ratioofinertia forcesactingonthe dropincomparisonto thesurface tensionforces.The responseofthe systemcanbe measuredthroughthequantificationofthedropmorphology.Thecombinedactionofthedynamicpressureandtheshear stressesonthesurfaceofthedropleadstoitsdeformation. Weadoptainteractionpotential baseddeformationapproach gearedtowards liquid–liquid interfacesashasbeendescribed indetailby [9].The computationaldomain istakenofsize L

= (

10,5,5)deq.Thedetails onthebaseandrefinedresolutionsoftheimmersedsphereisgiveninTable1.Thespherical

dropisplacedatxxx

= (

0.5,0.5,0.5)Lz,Lz istheheightintheverticaldirection(e

ˆ

z)anditisboundedbystationaryfree-slip

walls; e

ˆ

y directionis periodic and a uniformflow of UUU

=

Ue

ˆ

x is imposed in the e

ˆ

x direction. In Fig. 6(a) we show the

morphologyofthedropforW e

=

1.0 attwodifferenttimeinstantsshowingitsinitialsphericalstateanddeformedstate. Toquantifythe shapeoftheimmersed dropwecompute themeanaspect ratioofthedropas X

=

lz

/l

x wherelz andlx

arethelengthsofaboxsurroundingthedropin

ˆ

ez,e

ˆ

x directions,respectively.InFig.6(b),wecomparetheinverseofthe

meanaspectratioversustheWeber numberfordropsusingtherefinementtechniquewiththesimulationswithoutusing refinementandalsowithreferencedatatakenfromexperimentalmeasurements andagoodagreementbetweenallthree data-sets is observed.At high Weber numbers (W e

∼

5), there isa finite discrepancy inthe aspect ratios calculated by theapproach whichisattributedtorelatively highdeformationsofthe dropsurface inthisregimewhich theinteraction potentialapproachisunabletocaptureaccurately.Overall,thisshowsthattheadaptiverefinementtechniquecanbeused successfullytodecouplethedeformationmeshfromtheIBMmesh.Inordertoshowtheusefulnessofthisapproachwealso simulatetwostationarybodies,inaTaylor–Couette(TC)cellwherethefluidisconfinedanddrivenusingtwoindependently rotatingcylinders.In Fig.7(a)we showtwo spheres,one closetothe wallandtheother away fromthewall,interacting witha flow being driven inthe e

ˆ

θ direction. As evidenced inFig. 7(b),the grid is stretched in theradial (wall-normal)

direction.One caneasilyobservethat theflowprofiles nearthefluid-solidinterfaceare smoothandnodiscontinuities as seeninFig.4(a)developintheflow.

4. Coarse-grainedcollisiondetection

Asmentionedinsection1,whenmorethanone(deformable)bodyisevolvedinatime-dependentchaoticflow,collision amongtheimmersed bodiesisunavoidable.Since theenhancements proposed inthepaperaims atsimulatingturbulent flows with a dispersed immisciblephase, an efficient procedure to detect collisions becomes mandatory. The algorithms

Fig. 7. (a)SnapshotofradialvelocitycontoursinaTaylor–Couettecellwithtwosphericalparticles.(b)Gridspacingversusthenormalisedradial (wall-normal)position.

andmodelsrequiredtosatisfactorilydetectcollisions havereceivedconsiderableattentioninthepastgivenitswiderange ofapplicationssuchasvideogames,physicsbasedanimations,eventdrivensimulations,multiphaseflowsetc.There exist a numberof problemspecifictechniques intheliterature which canbe used todetect collisionbetweenanytwo points residing ondifferentobjects;see[29] foradetailedreview.The bruteforce approachistocomputethedistancebetween amarkerononebodyagainsteveryothermarkeronotherbodiestodetermineacollisionevent.Usingsuchasimpleand straight-forwardtechnique comesatacost,asthenumberofoperationsrequiredtodetectacollisionbecomes terminally expensive with increasing number of markers. Here ‘terminally expensive’ means an O(N2

)

algorithm, where N is the total numberof Lagrangiannodes usedforall theimmersed boundaries. Significantprogress hasbeenmadein thefield of moleculardynamicswherethe typicalapproachis tomaintain alistofneighbouring particleswithin apredetermined supportradiuswithwhomcollisionislikelyandthenupdatethislistperiodicallythroughoutthesimulation[30].Whilethis a moreefficient andcomputationallyinexpensive procedure wherethetotal operationcount is O(Nlog(N)) ascompared to O(N2

)

operationsinthebruteforcemethod,itstilladdsadditionaloverheadtothetotalsimulationtimeper iteration. WeproposeacomputationallyinexpensivedetectionalgorithmwhichmakesuseoftheEulerianmeshrequiredfortheflow solutionwithminimaloverhead.Here,itisimportanttonotethatinthissection,wefocusonlyoncollisiondetectionand notonmodellingthecollisionforceswhichisanextensivefieldonitsown.

The basic idea of the algorithm is as follows. A collision event occurs when any two arbitrary Lagrangian markers belonging to two separate immersed bodies occupy the same Eulerian cell. This assumption is reasonable since within the immersedboundaryframework thethinfilm fluiddynamicsbetweentwo bodieswhich arejustabouttocollidecan never beresolved [31]. Thisproblemcanbe solvedby modellingthelubrication forceswhentheimmersed bodiescome closer thanthe widthofa single Euleriancell. Onceit isdetermined that twomarkers belongingto separate bodiesare occupying thesame Euleriancell,a collisionforce is activatedonthe respectiveLagrangian markersalong thenormal to thesurfaceoftheimmersedbody.Hereitisimportanttonotethatthealgorithminparticularlyusefulandcomputationally efficientwhenthereareseveraldispersedbodiesintheflowwhichcaninteractwitheachotherandthesurroundingfluid simultaneously.

First,we initialiseathreedimensionalintegerarray

coll[1:N1,1:N2,1:N3]=0

,where

N1,N2,N3

arethenumber

ofEuleriangridpointsineachdirection.Thisarraycanhavethreestates;

coll[i,j,k]=0

impliesthereisnoLagrangian marker occupying the cell,

coll[i,j,k]>0

implies there are one or moreLagrangian markers belongingto the same

immersed body in the cell,

coll[i,j,k]=-1

implies markers from different immersed bodies occupy the cell anda

collision eventmustoccur inthiscell. InFig.8(a) weshow a schematic ofan Euleriangridwiththreeimmersed bodies occupyingthesameEuleriancelli.e.thevalueofthearray

coll

inthiscellissetto

−

1.Thevaluesofthearray

coll

is shownnext totherespectiveEuleriancell centresinFig.8(a).Insucha situationacollisionforce isactivatedoneach of thesemarkers.Forthesimulationsshownhere,weassumethatthemagnitudeofthecollisionforceisinverselyproportional to the distance betweenthe Lagrangian marker andthe Eulerian centre.The formulation of theforce can be of various typeswherelubricationforcescanbeincludedbasedontheapproachvelocityofthemarkerwithrespecttotheEulerian cell centre; however, the focus hereis on the more general problem of detecting a collision event. We now presenta pseudo-codeforthealgorithmwhichshowsitseaseofimplementationintoanygeneralNavier–Stokessolver.

coll[1:N1,1:N2,1:N3] = 0

! Routine for enforcing interfacial boundary condition

do pid = 1,N_particle

do mid = 1,N_f

ii = index_of_marker(1,mid,pid)

jj = index_of_marker(2,mid,pid)

kk = index_of_marker(3,mid,pid)

Fig. 8. (a)Aschematicofthegridshowingthevaluesofthearraycollwiththreedifferentdeformablebodiesapproachingeachotherinthedarkened Euleriancell.(b)Snapshotsatdifferenttimeinstantsofadeformabledrop(showninred)drivenbyacrossflowcollidingwithastationaryrigidsphere (showninblue).(c)Instantaneoussnapshotof120deformabledropsinteractingwiththewallandamongthemselvesinalaminarTaylor–Couetteflow.

if (coll[ii,jj,kk].eq.0)

coll[ii,jj,kk] = pid

if (coll[ii,jj,kk].ne.pid) coll[ii,jj,kk] = -1

end do

end do

! Routine for computing forces from fluid and internal forces

do pid = 1,N_particle

do mid = 1,N_f

ii = index_of_marker(1,mid,pid)

jj = index_of_marker(2,mid,pid)

kk = index_of_marker(3,mid,pid)

if (coll[ii,jj,kk].eq.-1) -> Activate and compute collision force

end do

In Fig.8(b),weshow snapshotsatdifferenttime instantsofsimulationwherea freelymoving dropdriven bya cross flow collideswithastationaryrigidsphere.Ascanbe seen,thecollisioneventoccurssmoothly withoutbothbodies pen-etrating eachother.Inthe supplementarymaterial wehaveattachedananimation whichshowscollisioneventsbetween 120 deforming dropswith2562 markerseach immersedin alaminar Taylor–Couetteflow. An instantaneous snapshotof thisanimationisshowninFig.8(c).Hereitisimportanttoemphasizethattheoverheadintheoverallsimulationtime on implementing thecollision-detectionalgorithmdescribed above islessthan 2%in comparisonto other morefine-grained techniqueswhich haveaminimumoverheadofatleast10% forsuchlarge scalesimulations.Withtheimprovements de-scribed above, we have been able to simulate turbulent Taylor–Couette flow up to Reynolds numbers of approximately 2

×

104withhundredsoffinite-sizebubblestounderstandthegoverningphysicalmechanismsbehinddragreduction[32]. The next challenge in numericalsimulations of finite-size bubblesand drops inturbulent flows is themodelling and implementationofcoalescenceandbreakupphenomena.Withintheimmersedboundaryformulation,bothcoalescenceand breakupeventshavetobemodelledastheinterfacedynamicsarenotsolvedusingfirstprinciples.Thisbringsintoplayan arrayofnewchallengesasgivenbelow,onwhichwearecurrentlyworkingon.

•

Astrategy tomerge/breakuptwo ormoreunstructuredmesheswhileensuringthetriangularelements donotdistort beyondacriticallimit.Here,a choicecanbe madebetween(i)stitchingthemesheswithinthe numericalsimulation whichcanbenumericallyexpensiveforthousandsofinteractingbubbles/drops(fore.g.[33,34]),or(ii)resorttocoarse grainedapproachessuchasimportingpre-definedgeometriesintothesimulationduringcoalescenceorbreakup.

•

Adeterministic strategytodeterminewhencoalescence orbreakup occurs.Thiswoulddependonseveralparameters

suchasapproachvelocityofthebubbles/drops,theirinstantaneousshapes,surfacetensionsetc.

•

Conservationofmass,momentumandenergyduringthenumericalprocessofcoalescenceorbreakup. 5. Summary

The immersedboundarymethod (IBM)isbecoming apopulartool amongcomputational scientiststo tackleproblems

involving moving anddeformable bodies immersed ina fluid. While IBM provides a simple andefficient computational

frameworkforsuchproblems,therearecertaininherentissueswhichlimitsthescaleofthesystemsaffordable.Forexample, currentstateofartsimulationscanhandleO

(10

2

₎

_{disperseddeformablebodiesinweaklyturbulentflows.Inthiswork,we}

proposeandtestdifferentalgorithmswhichcanbeusedtohandleO

(10

4

₎

_{disperseddeformablebodiesinhighlyturbulent}

flows withintheimmersedboundaryframework. Wedescribeafastmovingleastsquaresapproximationtechnique where thetransferfunctionsusedtoexchangeinformationbetweentheEulerianandLagrangianmeshesarecomputedonlyonce atthebeginningofthesimulationandacomputationallyinexpensiveinterpolationtechniqueisusedduringthesimulation toestimatethetransferfunctionatanyarbitrarypointinthedomain.Wefindthatwithsufficienttransferfunctionsstored atthestartofthesimulation,wearenotonlyabletoestimatethetransferfunctionswithin acceptablelevelsoferrorbut alsoreducedrasticallythetotalcomputationaltimerequiredwithincreasingimmersedLagrangianmarkers.Wealsoshow that through asystematicrefinement ofthe Lagrangianmesh, itcan bedecoupled into twomeshes;a basemesh which can be used for the solution of the deformation governing equations and a refined mesh which is used to enforce the interfacialboundarycondition.Additionally,suchaprocedurecanalsobeusedwhenthedispersedbodiesareimmersedin flows withstretchedEuleriangridsthuseliminating redundantLagrangiannodes.Finally,we presentasimpleandeasyto implementalgorithmtodetectcollisioneventsbetweenLagrangianmarkersresidingondifferentdeformablebodieswitha computationaloverheadoflessthan1%foranynumberofLagrangianmarkersimmersedintheflow.

Acknowledgements

TheauthorsthankValentinaMeschiniforvaluablediscussionsduringthecourseofthiswork.Thisworkwassupported by theNetherlands CenterforMultiscaleCatalyticEnergy Conversion (MCEC),an NWOGravitationprogrammefundedby the Ministry of Education, CultureandScience ofthe government ofthe Netherlands. The simulations were carriedout onthenationale-infrastructureofSURFsara,asubsidiaryofSURFcooperation, thecollaborativeICTorganizationforDutch educationandresearch.We alsoacknowledgePRACEforawarding usaccesstoMarconi supercomputer,based inItalyat

CINECAunderPRACEprojectnumber2016143351.

References

[1]E.Dowell,K.Hall,Modelingoffluid-structureinteraction,Annu.Rev.FluidMech.33 (1)(2001)445–490.

[2]R.Mittal,G.Iaccarino,Immersedboundarymethods,Annu.Rev.FluidMech.37(2005)239–261.

[3]C.S.Peskin,Flowpatternsaroundheartvalves:anumericalmethod,J.Comput.Phys.10 (2)(1972)252–271.

[4]M.Uhlmann,T.Doychev,SedimentationofadilutesuspensionofrigidspheresatintermediateGalileonumbers:theeffectofclusteringuponthe particlemotion,J.FluidMech.752(2014)310–348.

[5]F.Picano,W.-P.Breugem,L.Brandt,Turbulentchannelflowofdensesuspensionsofneutrallybuoyantspheres,J.FluidMech.764(2015)463–487.

[6]A.Prosperetti,Lifeanddeathbyboundaryconditions,J.FluidMech.768(2015)1–4.

[7]W.Fornari,F.Picano,L.Brandt,Sedimentationoffinite-sizespheresinquiescentandturbulentenvironments,J.FluidMech.788(2016)640–669.

[8]M.deTullio,G.Pascazio,Amoving-least-squaresimmersedboundarymethodforsimulatingthefluid-structureinteractionofelasticbodieswith arbitrarythickness,J.Comput.Phys.325(2016)201–225.

[19]S.Schaefer,T.McPhail,J.Warren,Imagedeformationusingmovingleastsquares,in:ACMTrans.onGraphics,vol. 25,ACM,2006,pp. 533–540.

[20]S. Fleishman, D. Cohen-Or,C. Silva, Robustmoving least-squares fittingwith sharpfeatures,in: ACM Trans.on Graphics, vol. 24,ACM, 2005, pp. 544–552.

[21]R.Kolluri,Provablygoodmovingleastsquares,ACMTrans.Algorithms4 (2)(2008)18.

[22]Q.-H.Zeng,D.-T.Lu,Curveandsurfacefittingbasedonmovingleast-squaresmethods,J.Eng.Graph.1 (3)(2004)84–87.

[23]L.Kobbelt,M.Botsch,Asurveyofpoint-basedtechniquesincomputergraphics,Comput.Graph.28 (6)(2004)801–814.

[24]G.Liu,Y.Gu,AnIntroductiontoMeshlessMethodsandTheirProgramming,Springer,Berlin,Germany,2005.

[25]G.-R.Liu,M.B.Liu,SmoothedParticleHydrodynamics:AMeshfreeParticleMethod,WorldScientific,2003.

[26]T.Kempe,M.Lennartz,S.Schwarz,J.Fröhlich,Imposingthefree-slipconditionwithacontinuousforcingimmersedboundarymethod,J.Comput.Phys. 282(2015)183–209.

[27]F.F.Abraham,FunctionaldependenceofdragcoefficientofasphereonReynoldsnumber,Phys.Fluids13 (8)(1970)2194–2195.

[28]E.Loth,Quasi-steadyshapeanddragofdeformablebubblesanddrops,Int.J.Multiph.Flow34 (6)(2008)523–546.

[29]H.Sigurgeirsson,A.Stuart,W.-L.Wan,Algorithmsforparticle-fieldsimulationswithcollisions,J.Comput.Phys.172 (2)(2001)766–807.

[30]M.Allen,D.Tildesley,ComputerSimulationofLiquids,OxfordUniversityPress,1989.

[31]S.Thomas,A.Esmaeeli,G.Tryggvason,Multiscalecomputationsofthinfilmsinmultiphaseflows,Int.J.Multiph.Flow36 (1)(2010)71–77.

[32] V.Spandan,R.Verzicco,D.Lohse,PhysicalmechanismsgoverningdragreductioninturbulentTaylor–Couetteflowwithfinite-sizebubbles,J.Fluid Mech.(2018),https://doi.org/10.1017/jfm.2018.478,inpress.

[33]F.Da,C.Batty,E.Grinspun,Multimaterialmesh-basedsurfacetracking,ACMTrans.Graph.33 (4)(2014)112.

[34]M.Razizadeh,S.Mortazavi,H.Shahin,Dropbreakupanddroppaircoalescenceusingfront-trackingmethodinthreedimensions,ActaMech.229 (3) (2018)1021–1043.