Particle-Based Simulation of Fluids

1 EUROGRAPHICS2003/ P. (GuestEditors)Volume22(2003),Number3 Particle-BasedSimulationofFluidsSimonPre mo e1, TolgaTasdizen2, JamesBigler2, AaronLefohn2andRossT. Whitaker21 ComputerScienceDepartment,UniversityofUt ah2 ScientificComputingandImagingInstitute,U niversityofUtahAbstractDuetoourfamiliari tywithhowfluidsmoveandinteract,aswellast heircomplexity, plausibleanimationoffluidsremainsa presenta mixingfluidswithdifferentphysicalpropert ies,fluidinteractionswithstationaryobjec ts, suitedformediumscaleproblemssincecomputa tionalelementsexistonlywhere ,plausiblesimulationoffluidsremainsa challengingprob-lemdespiteenormousimprov ements7;6.

2 Advancesincom-putationalfluiddynamics(CF D)oftencannotbedirectlyappliedtocomputer graphics,becausethey have ,controlovermotionandappearanceis necessaryforartisticpurposessuchthatphys icallyimpos-siblethingsbecomepossibleand thatfluidmotionis script-ableforuser s , existingmethodsarecomputation-allyveryex pensive andveryslow beenmadeto make fluidsimulationmorecontrollable7, , ,multiphaseflows,multiplefluidsmixing,an dsedi-mentaryflowsarenoteasyto hybridrepresentation:frag-mentationandme rgingoffluids,numericaldiffusionin con-vectioncomputation, approachestogrid-basedmethodsforsimulati ngfluidflows:LargeEddySimulation,vortici tyconfinement,vortex vethods, (LES) methods, largetimestepsareallowedandcomputa-tiona lelementsexistonlywhereinterestingflow addresssomeofthedeficienciesofthegrid-ba sedmethods,wedescribea (theNavier-Stokesequations) fitsbetterintocurrentuser-interactionpar adigmandsettingthesim-ulationandcontroll ingit is (inflow andoutflow bound-aries,obstacles,forces)atcoarseres olution(smallnumberoflargeparticles) ,thesimulationis (fluidparticles)

3 Areonlyusedinpartsofthescenewherethey ,if moredetailis re-quiredinpartofthescene, , beenpresentedin computergraphicsliteraturec ,108 Cowley Road,OxfordOX41JF, UKand350 MainStreet,Malden,MA02148, eet Particle-BasedSimulationofFluidsbefore,n oneofthosemethodsdealtwithincompressible flu-idsandwater-like is difficulttocreatea particlesareneces-sarytoobtaina differentmeth-odscanbeusedtocreatea surfacefromparticles, reconstructedona ,forpreview purposes,blendingofpotentialfieldsaround eachparticlecanbeusedtogive fastfeedbackonhow results,theparticle-basedmethodcouldprov ideanalternative forsimu-lationandanimationofvarietyofflu idswithdifferentphys-icalpropertieswhile allowingusercontrolandfastfeedbackat Earlymethodsweregearedtowardssimplify-in gthecomputationbyusingFouriersynthesis18 orprovid-ingspecializedsolutionfora specificproblem9;25.

4 BrienandHodgins22alsouseda heightfieldrepresentationcoupledwitha Marker-And-Cell(MAC)method11tosolve (SPH) col-lectionofparticles,whichmove relevantengineeringproblems,includinghea tandmasstransfer, moleculardynamics, flexibleLagrangianmethod20;21thatcaneasi lycapturelargeinterfacedeformation,break ingandmerging, modelofdeformablesurfacesandmetamorphosi swithactive varietyofsubstancesusingastateequationto computethedynamicsofthesubstance2. Adaptive sam-plingofparticlesimprovedthecomputati onalefficiency ofthemethodbysubdividingparticleswheresu bstanceunder-goeslargedeformation, thegoverningstateequationforanimatedlava flexible,it canonlysolve com-pressiblefluidflow.

5 Severalextensionshave , anothergridlessparticlemethodcalledtheMo ving-ParticleSemi-Implicit(MPS) capableofsimulatinga widevarietyoffluidflowproblemsincludingp hasetransitions,multiphaseflow, sediment-ladenflowsandelasticstructures3 7;15;12. Thecomputationalalgorithminthispaperis a ,thecomputationalelementsintheMPSmethoda rediscretenumberofparticlesoffluidfollow edin , allproblemvari-ablesareobtainedfromvalue sat thesepointsthroughanin-terpolationfuncti on(kernel). fluidcanbedescribedbytheNavier-c eet Particle-BasedSimulationofFluidsuVelocit yrPositionrDistancebetweenparticlesreInt eractionradiusdNumberofspacedimensionstT imedtTimestepwWeightfunction Viscosity Surfacetension Surfacecurvature Densityn0 FluidparticledensityTable1 a velocityfieldofthefluid,thecon-tinuityeq uationstatesthatthemassmis constant:r u=0;(1) : u t+u ru= 1 rp+ r2u+f.

6 (2)where is densityofthefluid,pis thepressure, is ,surfacetension,conservationofenergyandm any otherrelationshipscouldalsobewrittenfora ,mass, parti-clesthataredistancerapartisw(r) = rerif 0 r re0ifre r(3)whereris thedistancebetweentwo particlesiandj,r=jrj rij:(4)Ifallparticlesareallowedtointerac t,thecomplexityisO(N2). Incontrast,ifinteractionradiusisrestrict ed,thecomplexityisonlyO(NM), whereMisthenumberofin-teractingparticles 24withintheradiusofinteractionre. Theparticlenumberdensityncanbecomputedas hnii= j6=iw(jrj rij):Thenumberofparticlesina unitvolumeis approximatedbytheparticlenumberdensityh nii=hnii w(r)dv:Foranincompressiblefluid,thefluid densitymustbecon- solve theNavier-Stokesequations, is theaverageofscalargradientbetweenparticl eiandneighboringparticlej:hr ii=dn0 j6=i j ijrj rij2(rj ri)w(jrj rij):(5)Similarly, thevectorgradientruistheaverageofvectorg radientbetweenparticleiandneighboringpar ticlej:hr uii=dn0 j6=i(uj ui) (rj ri)jrj rij2w(jrj rij):(6) canbeseenasif a fractionofa quantityatpar-ticleiis distributedtoneighboringparticlej.

7 Hr2 ii=2d n0 j6=i( j i)w(jrj rij);(7)where : = Vw(r)r2dv Vw(r)dv:(8)Notethatthismodeldoesnotrequi reany complex boundariesandobjects(stationaryormoving) problemsina , theforcesandviscos-ityinthemomentumconse rvationequationarecomputedexplicitly. Temporaryparticlelocationsr andvelocitiesu arecomputedfromthepositionsrnandvelociti esrnfromtheprevioustimestepnasfollows:u =un+dt r2un+ ( n)n+ g ;(9)andr =rn+u dt:(10)c eet Particle-BasedSimulationofFluidsThesurfa cetensionmodelandcomputationis ,theincompressibilityofthefluidis isnotequalton0. Theparticledensityn needstobemodifiedbyn0suchthatthecontinui tyequationis relatedtomodificationoftheve-locityu0:1d tn0n0=r u0:(11)Themodificationvelocityu0isobtain edfromtheimplicitpressureterminthemoment umconservationequation:u0= dt rpn+1:(12)Notethatthisis +n0;(13)a Poissonequationforpressureis obtained:hr2pn+1ii= dthn ii n0n0:(14)Therighthandsideofequation14isa nalogoustothedi-vergenceofthevelocityvec tor.

8 Equation14issolvedbyusingtheLaplaciandif ferentialoperatoranddiscretizingitintoa sparseandsymmetric; +1iscomputed,thecorrectionvelocityu0also becomesknown:u0= dt hrpn+1i:(15)Newparticlevelocitiesandposi tionsthatsatisfybothcon-servationofmassa ndmomentumcanthenbeupdated:un+1=u +u0(16)rn+1=rn+un+1dt:(17) simpleconditionhn ii< n0;(18) , weuse =0 , (oratmosphericpressure,if applicable)is representfixedobjectstoensurethatparticl edensitynumberiscomputedaccurately. Notethatthereis noexplicitcollisiondetectionAlgorithm1 TheMovingParticleSemi-Implicit(MPS) :u0, ,r Computeparticlenumberdensityn usingnewparticleloca-tionsr Setupandsolve PoissonpressureequationusingConjugateGra dientmethodComputevelocitycorrectionu0fr omthepressureequationComputenew particlepositionsandvelocities:un+1=u +u0rn+1=r + computedat , isa La-grangianmethod,inflow andoutflow offluidis , it canbeeasilycombinedwiththeEulerianap-pro achtohandleinflow andoutflow.

9 We willdescribea (energy, etc.)is onetrulycaresaboutcon-servation, ,largeaspectratiois impossibleat have particlesthataccuratelyrepresentthefluid motionandotherinterestingquantities,but forrenderingthefluida surfacerepresentationis ,oneofthemainproblemswitha purelyLagrangianapproachis ,localresolutionenhancementis hard, moreparticlesareintroducedtoimprove resolu-tion,they willsoonmove ,a gridlesshybridmethodhasbeendeveloped37. :particlepositionsarereconfig-uredc eet Particle-BasedSimulationofFluidsPSfragre placementsrLihrii rInterpolationregionFigure 1 one-dimensionalgridiscreatedintheparticl e s onlyinterpolatedinthecutdiskarea.

10 (AfterHeoet ) :particleconvectionis computedona one-dimensionalgridTheLagrangianphaseis exactlythesameasdescribedin denotetheparticlepositionsandvelocitieso btainedafterthisphaseasuLandrL. Thereconfigurationphaseandconvection(Eul erian) ,theparticlepositionshave fixedboundary(fixedobjects,walls,etc.)or inletoroutletboundary, shouldgobacktotheiroriginalpositionsrn. ThemovingboundarycanbetracedthroughLagra ngianmotionofpointsdescribingthefreesurf acewithoutcomput-ingtheconvectionterm:rn +1s=rLswheresubscriptsdenotesthatthepoin tis ,it is likelythatthepointsonthesurfacewillclust ertogether. To fixthisprob-lem,wemake ,it is Therefore, +1at thenew timestepis deter-minedarbitrarilyandthevelocityofth ecomputingpointucisuc=rn+1 rndt:(19)Theconvectionvelocityis thengivenbyua=rL rndt rn+1 rndt= rn+1 rLdt:(20)EulerianPhaseAfterarbitraryconv ectionvelocityuais computed,proper-tiesatrn+1arecomputedbyi nterpolationofphysicalprop-ertyf(flow velocity, temperature,etc.)