Least squares moving prticle semi-implicit

Comp.Part.Mech.(2014)1:277–305DOI10.1007/s40571-014-0027-2

Leastsquaresmovingparticlesemi-implicitmethod

AnarbitraryhighorderaccuratemeshfreeLagrangianapproachforincompressible?owwithfreesurfaces

TasukuTamai·SeiichiKoshizuka

Received:27February2014/Revised:8May2014/Accepted:26May2014/Publishedonline:11July2014©OWZ2014

AbstractInthispaper,aconsistentmeshfreeLagrangianapproachfornumericalanalysisofincompressible?owwithfreesurfaces,namedleastsquaresmovingparticlesemi-implicit(LSMPS)method,isdeveloped.Thepresentmethodologyincludesarbitraryhigh-orderaccuratemesh-freespatialdiscretizationschemes,consistenttimeintegra-tionschemes,andgeneralizedtreatmentofboundarycondi-tions.LSMPSmethodcanresolvetheexistingmajorissuesofwidelyusedstrong-formparticlemethodforincompressible?ow—particularly,thelackofconsistencyconditionforspa-tialdiscretizationschemes,dif?cultyinenforcingconsistentNeumannboundaryconditions,andseriousinstabilitylikeunphysicalpressureoscillation.Applicationsofthepresentproposaldemonstrateremarkableenhancementsofstabilityandaccuracy.

KeywordsLeastsquaresmovingparticlesemi-implicitmethod·LSMPSmethod·Highorderscheme·Meshfreecompactscheme·Movingparticlesemi-implicitmethod

1Introduction

Todaythe?niteelementmethod(FEM)[101],the?nitevol-umemethod(FVM)[98],andthe?nitedifferencemethod(FDM)[87]basedcomputationalmechanicsplayaconspic-uousroleintechnologyadvancement.Anoteworthyfeatureofthemisthattheydivideacontinuumdomainintodiscretesubdivisionusuallycalledmesh/grid,whichrequiresconnec-T.Tamai(B)·S.Koshizuka

GraduateSchoolofEngineering,TheUniversityofTokyo,7-3-1,Hongo,Bunkyo-ku,Tokyo113-8656,Japane-mail:tasuku@mps.q.t.u-tokyo.ac.jpS.Koshizuka

e-mail:koshizuka@sys.t.u-tokyo.ac.jp

tivitybasedonatopologicalmap;however,thecharacteris-ticofthemisnotalwayssuitable.Forinstance,inordertoadapttopologicalandgeometricchangesundergonebytherealmaterial,suchassimulationsof?uid?oworlargestraincontinuumdeformation,Lagrangiandescription(i.e.movingmesh/grid)couldbeapplied;however,onewouldfacedis-tortionofmesh/gridwhichresultsineitherterminationofthecomputationorseveredeteriorationinnumericalaccuracy.ArbitraryLagrangianEulerian(ALE)formulations[26,36]havebeendevelopedtoovercomethedif?cultycausedbydistortionofmesh/grid,ofwhichobjectiveistomovemesh/gridindependentlyfromactualmotionofmaterialsothatdis-tortioncouldbeminimized;nevertheless,distortionofmesh/gridstillremainsandcausesoverwhelmingerrorsinnumer-icalsolutions.

Underthesecircumstances,meshfreemethodsand/orpar-ticlemethods,whichdiscretizeacontinuumbyonlyasetofnodalpointsorparticles,havebeensoughtinorderto?ndbet-terdiscretizationprocedureswithoutmesh/gridconstraints.Sincetheconnectivityamongnodescanbegeneratedany-timedesiredandcanchangewithtime,meshfree/particlemethodscaneasilyhandlesimulationsofverylargedefor-mations,evenwiththechangesofthetopologicalstructureandfragmentation–coalescenceofcontinuum.

Ingeneral,accordingtocomputationalmodelingsandfor-mulations,meshfree/particlemethodscanbecategorizedintotwodifferentclassi?cationsaswell:theweakformformula-tionsofPartialDifferentialEquations(PDEs);andthestrongformformulationsofPDEs.The?rstclassofmeshfree/particlemethodsisusedwithvariousweakformulationssuchasGalerkinmethods,forexample,diffuseelementmethod(DEM)[73],elementfreeGelerkin(EFG)method[11,13–15,67],reproducingKernelparticlemethod(RKPM)[21,63–65],h-pcloudmethod[27,28,58,74],partitionofunitymethod(PUM)[7–9,69],meshlesslocalPetrov–Gelerkin

123

278(MLPG)method[2–6],?nitepointmethod(FPM)[75–77,80],particle?niteelementmethod(P-FEM)[38,78,79],reproducingKernelelementmethod(RKEM)[56,62,66,86].Thesecondcategoryofmeshfree/particlemethodsisusedtoapproximatethestrongformPDEsdiscretizedbyspe-ci?ccollocationtechniques,forinstance,smoothedparti-clehydrodynamics(SPH)method[31,60,68,71,72],movingparticlesemi-implicit(MPS)method[48,49],meshfree?nitedifferencemethod[57,59],?nitepointsetmethod(FPM)[95–97].

Invariousweakformmeshfree/particlemethods,shapefunctions,ormoregeneralmeshfreeinterpolants,arecon-structedbasedontheso-called“partitionofunity”,thencon-sistencyconditions,polynomialcompleteness,orreproduc-ingconditionsaresatis?ed.Ontheotherhand,althoughthestrongformparticlemethodssuchastheSPHmethodandtheMPSmethodhavebeenshowntobeusefulwidelyinengi-neeringapplicationsespeciallyin?uiddynamics,theirstan-dardformulaeofspatialdiscretizationschemesarenotcon-sistentexceptunderverylimitedconditions,anddonotholdpolynomialcompleteness/reproducingconditionsordiffer-entialcompleteness/reproducingconditions.AccordingtotheLax’sequivalencetheorem[51],thismajorissuemustneverbeoverlooked.Moreover,thismatteryieldsadverseeffectsforbothcomputationalaccuracyandstability.Somecorrectionmethodsforresolvingthelackofpolynomialcom-pletenessorreproducingconditionsonthespatialdiscretiza-tionschemes(whichwillbediscussedinChap.3)havebeenproposed;however,theyarefarfromadequateintermsofthecompatibilitywithsatisfactionofhigherorderconsis-tencyconditionsandnumericalstability.Strongformmesh-free/particlemethodsstillhaveadif?cultyrelatingtopro-ceduresofenforcingboundaryconditions,especiallyNeu-mannboundaryconditions.Hence,prevailingstrongformparticlemethods,whoseadvantageisthattheycanhand-ilyrunnumericalanalysisofcontinuumwithlargedefor-mation,evenwiththechangesoftopologicalstructureandfragmentation–coalescence,areinadequatelystudiedasanaccuratemathematicalcomputation.

Withtakingparticularnotetocontroversiesdescribedabove,wedevelopanewconsistentfullyLagrangianmesh-freeparticlemethod,namedleastsquaresmovingparticlesemi-implicit(LSMPS)method,fornumericalanalysisofincompressible?uid?owwithfreesurfaces.Asitsnamesug-gests,LSMPSmethodisbasedonthemethodofweighted“LeastSquares”procedure,andfollowsfundamentalsoftheMPSmethod[48,49]:“MovingParticle”meansmeshfreefullyLagrangianapproach,and“Semi-implicit”representsthetypeoftimeintegrationalgorithmforincompressible?owwhichiswell-knownastheprojectionmethod[33].LSMPSmethodsucceedsthenameoftheexistingMPSmethod;how-ever,allofLSMPSformulaearedifferentfromofthecurrentMPSmethod.

123

Comp.Part.Mech.(2014)1:277–305

Inthispaper,weintroduceanewmethodologyandfor-mulaeasameshfreeLagrangianapproach(particlemethod),includingarbitraryhighorderaccuratemeshfreespatialdis-cretizationschemes,consistenttimeintegrationschemes,andgeneralizedtreatmentofboundaryconditions.Additionally,somenumericaldemonstrationscomparedwiththeconven-tionalMPSsolutionsshowdrasticimprovementofaccuracyandstability.

2Preliminary2.1Notation

Inordertoexpeditethepresentation,weintroducesomepre-liminariesfornotations.Throughoutthispaper,theletterdisapositiveintegeranddenotesthespatialdimension.Ω?Rdisanonempty,open,bounded,andconnectedset.?ΩdenotestheboundaryofΩ,and?ΩisassumedtobeLipschitzcon-tinuousorsmoother,asthecasemaybe.

N0denotesthesetofnonnegativeintegers.Ifα:=(α1,α2,...,αd)∈Nd0

(1)

isand-tupleofnonnegativeintegers,wecallαamulti-index.Then,thequantity|α|:=

??dαi,(2)

i=1

isde?nedtobethelengthofα.Wealsousethefollowing

conventions:α!:=α1!...αd!,

?α∈Nd0.

(3)

Ifα,β∈Nd0,wesayβ≤αprovided1≤?i≤d,

βi≤αi.

(4)

Bythesametoken,

??α??

??????β

:=

α!

α1=β1...αd??βd.(5)

Ifx:=(x1,...,xd)T∈Rdandα∈Nd0,thenxαis

de?nedasfollows:xα:=xα1

1

···xαdd.(6)

Iff(x)isarealvaluedfunctiononanopensubsetofRd,α∈

Nd0andsmoothnessoffisassumedenough,thenDαthorderFréchetderivativeoffasfollows:

xf(x)

denotestheαDαx

f(x):=?|α|f(x)

11d.(7)

d

Comp.Part.Mech.(2014)1:277–305Ifx=(x??1,...,xd)TisanarbitrarypointofRd,then

????x??????

d2=??

xi2(8)

i=1

denotestheEuclidiannormofx,andweusuallyuse??x??as

abridgednotation.

2.2Consistency,completeness,reproducingcondition,and

thepartitionofunityInordertoadvanceconcretediscussionsonthestudyofcon-vergence,accordingtoBelytschkoetal.[12],threetermsarefrequentlyused:(i)Consistency[90]:apropertyofthedis-cretizationschemesforpartialdifferentialequations,whichisusuallyutilizedinFiniteDifferenceapproximations,(ii)Reproducingconditions[64,65]:theabilityoftheapproxi-mationtoreproducespeci?edfunctions,whichareusuallypolynomials,(iii)Completeness[37]:polynomialcomplete-nessorcompleteness,whichisusuallydiscussedintheFEM.2.2.1Consistency

Strikwerda[90]de?nestheconsistencyconditionofpartic-ularoperatorsforapproximatingderivativesasfollows:Theorem2.1(Consistency[12,90])AschemeLhu=fthatisconsistentwiththedifferentialequationLu=fisaccurate(consistent)oforderpifforanysuf?cientlysmoothfunctionv

Lv?Lhv=O(hp)

(9)

Intheabovede?nition,aparameterhdenotesthere?nementofmesh/gridandpreferstotheorderofconsistency.Obvi-ously,itisnecessaryforconvergencethatp>0,andwerequirethatp≥1fortheef?cientnumericalcalculation.Accordingtothewell-knownLax-Richtmyerequivalencetheorem[51],aconsistent?nitedifferenceschemeforawell-posedpartialdifferentialequationisconvergentifandonlyifitisstable.Consequently,anydiscretizationschememustsatisfythepthorder(p>0)consistencyconditiontoobtainconvergence.Theconsistencyisstraightforwardtoverifyandstabilityistypicallymucheasiertoshowthanconvergence;therefore,theconvergenceisusuallystudiedviatheLax-Richtmyerequivalencetheorem.

2.2.2Completeness,reproducingcondition,andthe

partitionofunity–nullitySinceascertainingwhethermeshfreeinterpolantsareconsis-tentforirregularlydistributednodesissigni?cantlymoredif-?cultthanexaminingwhether?nitedifferenceschemesforauniformstructuredgrid,thecompletenessorreproducing

279

conditionswhichplaythesamerole[12]astheconsistencyconditions1areexplored,instead.

Thereproducingconditionsorthecompletenessaretheabilityoftheapproximationtoreproducespeci?edfunctionswhichareusuallypolynomials.Onecansaythatanapproxi-mationfh(x)iscompletetoorderpifanygivenpolynomialuptoorderpcanbereproducedexactly.Ifanapproximationfh(x)isgivenfh(x)=??

by

Φi(x)f(xi),(10)

i

where{Φi(x)}1≤i≤Naretheinterpolantfunctionsand{f(xi)}1≤i≤Naregivennodalvaluesforthesetofnodes{xi}1≤i≤N,thenthecompletenessorthereproducingcondi-tionscanbede?nedasfollows:

Theorem2.2(pthordercompleteness/reproducingcondi-tion)Foramultiindexα∈Nd(x)∈Rholdsthe0:pth0≤order|α|≤polynomialp,aninterplantfunctionΦicom-pleteness/reproducing??

conditionifitsatis?es:(x?xi)αΦi(x)=δα0,(11)

i

or??

equivalently,xiαΦi(x)=xα.

(12)

i

Somemeshfreediscretizationschemesareformulatedtosatisfyalternatives,thedifferentialcompletenessorthedif-ferentialreproducingconditionswhicharerequirementthatthederivativesofapolynomial?eldbereproducedcorrectly.Theycanbede?neddirectlyfromTheorem2.2bytakingderivativesofEqs.(11)and(12),i.e.

Theorem2.3(pthorderdifferentialcompleteness/reproducingcondition)LetaninterplantfunctionΦki|β|(Rd≤)k.(For≤pmulti),anindiciesinterplantα,functionβ∈Nd(x)∈C0:Φ0≤|α|≤p,0≤i(x)holdsthepthorderdifferentialcompleteness/reproducingconditionifitsatis?es:??(x?xi)αDβ

xΦi(x)=(?1)|β|α!δαβ,(13)

i

orequivalently,??xiαDβ

xΦi(x)=

α!

i

xα?β.

(14)

ItshouldbenotedthattheTheorem2.2and2.3aresimi-lartopthorderconsistencyconditionandpthorderdifferen-tialconsistencyconditionforRKPMshapefunction[55,65],respectively.

1

pthorderpolynomialcompletenessconditionsorreproducingcondi-tionsissuf?cientconditionofpthorderconsistency.

123

280Ifα=0,Eqs.(11)and(12)become??

Φi(x)=1,

(15)

i

andifα=0,Eq.(11)does??

(x?xi)αΦi(x)=0.

(16)

i

Thesearetheoriginofthename“ThePartitionofUnity”,and“ThePartitionofNullity”.Thesepropertiesarecloselyrelatedwithnotonlymeshfreeinterpolant[65]butalsomesh-basedone.Obviously,sinceFEMinterpolantfunc-tions,calledshapefunctionssatisfytheKroneckerdeltaprop-erty,theyareconstructedbasedonthepartitionofunity–nullity.Moreover,weightingcoef?cientsoflinearcombina-tiononthe?nitedifferenceschemesarebuiltonthepartitionofunity–nulltytoobtaincertainorderofconsistency.

Interestingly,Liuetal.[65]showedthatleastsquaresbasedinterpolantsprovideachievementofpolynomialcom-pleteness/reproducingconditions,orso-calledthepartitionofunity–nullity—viceversa,correctedinterpolantformulaetoful?llthepolynomialcompleteness/reproducingcondi-tionsarequiteidenticaltotheproductionderivedfromleastsquaresprocedures.Chakravarthy[17]showedthefunda-mentalconceptofthe?nitedifferenceschemesthattodeter-mineunknownsemergedfromTaylorexpansionorpoly-nomialapproximationisgeneralizationof?nitedifferenceoperator.Also,ityieldsnormalequationswhichistheveryideaofthelinearleastsquaresapproaches.

Withfocusingattentiontoacloserelationshipofconsis-tency,completeness/reproducingconditions,andtheparti-tionofunity–nullity,andleveragingthefactthattheleastsquaresschemescontributeachievementofarbitraryhigh-orderconsistencyconditionsformeshfreespatialdiscretiza-tionschemes,wedevelopnewformulaeinthenextsection.

3Meshfreespatialdiscretizationschemes3.1Anoverviewoftheexistingmeshfreespatial

discretizationschemes

Asmentioneditintheintroduction,variousmeshfreeand/orparticlemethodshavebeensought,inordertodiscretizeadomainwithoutmesh/gridconstraints.Therearealotofdis-cretizationschemesformeshfreeinterpolantsandmeshfree?nitedifference.Oneofthemostprevalentparticlemeth-odsistheSPHmethodwhichdiscretizespartialdifferentialequationsbyaintegralrepresentationcollocationtechnique.EventhoughtheSPHmethodhasachievedalotofsuccessincomputationalmechanics,ithasnotbeenviewedasanaccuratemathematicalcomputationwhichstemsfromthe

123

Comp.Part.Mech.(2014)1:277–305

factthatitlacksarigorousconvergencetheoryaswellasasuccessivere?nementprocedure[81].

TheearlySPHinterpolantsdonotsatisfythediscretepar-titionofunityandnullity[63]exceptunderverylimitedcon-ditions.ThismeansthedefectionofSPHinterpolants0thorhigherordercompleteness/reproducingconditionsingeneralparticledistribution,whichresultsinincapabilityofrepre-sentingrigidbodymotioncorrectly,eventhoughitisGalileaninvariant(rigidbodytranslationonly).Also,theprimalSPHgradientoperatorandLaplaceoperatorgenerallyhaveanalo-gousproblem.Inordertosolvethismatter,severalcorrectionschemeshavebeenproposed,forexample,Monaghan’ssym-metrizationonderivativeapproximation[70,71],Johnson–Beisselcorrection[41],Randles–Liberskycorrection[82],Krongauz–Belytschkocorrection[12],Chen–Berauncor-rection[18–20],Bonet-Kulasegaramintegrationcorrection[16],Aluru’scollocationRKPM[1],Zhang–Batracorrection[99,100].TheycorrectSPHkernelinterpolanttosatisfycom-pleteness/reproducingconditionintheinterpolation?eld,orequivalently,tomodifySPHderivativeapproximationoper-atorsdirectlytomeetderivativecompleteness/reproducingconditioninthederivativeoftheinterpolants.Itmustbemen-tionedthatalmostall1storhigherorderconsistentcorrectedformulaeaslistedabovearebasedontheleastsquaresmeth-ods,forinstance,themostwidelyused1storderconsistentgradientapproximationoperatorfortheSPHmethoddevel-opedbyRandles–Libersky[82]isoneoftheleastsquaresbaseddiscretizations.

Leastsquaresproceduresareexcellentwithmeshfreeinterpolantsormeshfree?nitedifferenceschemes.Liuetal.[65]demonstratethatmovingleastsquares(MLS)[50]approximationisequivalenttoreproducedKernel(RK)interpolant[65]whicharecorrectedkernelapproximationbasedonthereproducingconditions.MLSRKinterpolantcanprovidearbitraryhighorderconsistencycondition,andisthecontemporaryversionoftheclassicalMLSonesincethebasicconceptofreproducingkernel,ormoregener-allyspeaking,thefundamentalsofthepartitionofunity–nullityincubatesvariousderivationofmeshfreeinterpolants[54,61,62].WithtakingparticularnoteoftheSPHkernelinterpolantinconsistency,Dilts[24,25]utilizesMLSinter-polant[50]toimprovetheaccuracyoftheSPHkernelapprox-imation,so-calledmovingleastsquaresparticlehydrody-namics(MLSPH).SinceMLScanconstructsuf?cientlysmoothedinterpolationgloballywithoptionalorderofpoly-nomialcompleteness/reproducingconditions,itiswidelyusedinmeshfree/particlemethods,suchasvariousGalerkinmeshfreeapproachlistedinIntroduction(e.g.DEM,EFGM,RKPM,etc.).

SincetheearlyMPSgradientoperatorandLaplaceoper-ator[49]generallylack1storhigherorderdifferentialcom-pleteness/reproducingconditionexceptunderverylimitedconditions(e.g.regularparticledistributionisassumed),

Comp.Part.Mech.(2014)1:277–305variedcorrectiontechniqueshavebeenproposedinordertoenhancetheaccuracyoftheMPSmethod.Forinstance,Khayyer–Gotohgradientoperatoranti-symmetrization[42],Khayyer–Gotohdivergenceoperatorcorrection[43],Khayyer–GotohLaplaceoperatorcorrection[44],Khayyer–Gotohgradientoperatorcorrection[45],andSuzukigradi-entoperatorcorrection[39,91]areproposed.OnlySuzukimethodbasedonweightedleastsquarestechniquecanachieve1storderconsistencyforgradientoperator;however,othersdonotholddifferentialcompleteness/reproducingconditionsingeneralcase.Inotherwords,theyarefarfromsatisfactionofhighorderaccuracy.Ofcourse,sim-ilartoadoptingleastsquaresapproachintoSPHmethod,leastsquaresbasedformulaecanbeintroducedtotheMPSmethod.Kohetal.[46]utilizetwodimensionalsecondordergeneralized?nitedifferenceschemes[57,59]basedonweightedleastsquarestotheMPSmethod.Tamaietal.[93]formulategeneralized?nitedifferenceschemesbasedontheweightedleastsquaresfortheMPSmethod,whichprovidesarbitraryhighorderconsistencyandcanbeappliedforarbitrarydimension.

Consistentleastsquaresbasedspatialdiscretization?nitedifferenceschemes[46,57,59,91,93]resolvetherackofpolynomialcompleteness/reproducingconditionsontheMPSmethodsothathighorderconsistencyconditionsareful?lledandabsoluteenhancementofaccuracywouldbegiven;however,utilizingleastsquaresbasedschemeraiseanewproblem—normalequationsderivedfromtheleastsquaresprocedureswillbeill-conditionedproblemswhichresultsineitherseriousdeteriorationinnumericalaccuracyandstabilityorterminationofthecomputation.Selectingneighborhoodstenciles[59]tocircumventthisill-conditionedproblemsisproposed;however,thistechniquecannotbethefundamentalsolutionsincetheconditionnum-bersofcoef?cientmatricesderivedfromtheleastsquares,so-calledthe“momentmatrices”,arenotindependentfromcharacteristiclengthofcalculationpointsspacing.

Inordertoovercometheweaknessoftheleastsquaresbasedspatialdiscretizationschemesthatnormalequationsshallbeill-conditionedandtoobtainhighorderconsistencyconditionsfromthem,wereformulatetheleastsquaresbasedschemesinthenextsection.

3.2Anewmeshfreespatialdiscretizationschemesbasedon

theweightedleastsquaresprocedure3.2.1Stone–WeierstrasstheoremoflocallycompactversionLetf:Rd→Rbeasuf?cientlysmoothfunction2thatisde?nedonasimplyconnectedopensetΩ?Rd.Accord-ingtotheStone-Weirestrasstheoremoflocallycompactver-2

Atleast,f(x)∈C0(Ω)

¯.281

sion[84,88,89],fora?xedpointx

¯∈Ω¯,oneshouldalwaysbeabletoapproximatef(x)byapolynomialserieslocally.Thus,wecande?nealocalfunction

??

fl(x,x

¯):=f(x),?x∈B(x

¯),0,?x∈B(x¯),(17)where

B(x

¯):=??x????????x?x¯??<re,x∈Rd

??

.(18)

Ifthefunctionf(x)issmoothenoughasassumed,there

existsalocaloperatorLx¯:C0(B(x¯))→Cp(B(x¯))s.t.fl(x,x¯)≈Lx¯f(x):=??pT(x)a(x¯),(19)

where

??p(x):=??xα??????0≤|α|≤p??

,(20)

ispthordercompletepolynomialbasis,anda(x

¯)iscoef-?cientvector.UtilizingTaylorexpansionofapproximated

polynomialfunctionaroundxiwithnearbypointxj∈B(x

¯)yields

??p??1??

xj?xi αα??h|α|=1

Dxf(xi)?{f(xj)?f(xi)}=Rp+1

ij,(21)

whereRp+1

ij

:=Lx¯f(x)?fl

(x),

(22)

istheresidualoflocalpolynomialapproximation.Equation(21)isawell-knownformofTaylorexpansion,andweuseitseveraltimesinthispaperwithoutspecialnoteagain.3.2.2Weightfunction

Weusetheweightedleastsquaresproceduresfornewspatialdiscretizationschemesformulae,thentheweight(window)functionwhichsatisfythefollowingconditionsisde?ned.

w(x,re)∈Ck

(Rd),1≤k,(23)?x∈Rd,

0≤w(x,re)≤Cw<∞,

(24)??x??≥re??w(x,re)=0,

(25)????

x??<??y??≤re???w(x,re)>w(y,re),

(26)w(x,re)dx=C=(Const.),(Typically,C=1),

(27)

Rd

wherereisthedilationparameterandtheradiusofcom-pactsupportoftheweightfunction.IntheLSMPSmethod,

singularweightfunctionlike

123