A Momentum-Conserving Implicit Material Point Method for Surface Energies with Spatial Gradients
Jingyu Chen, Victoria Kala, Alan Marquez-Razon, Elias Gueidon, David A. B. Hyde, Joseph Teran
AA Momentum-Conserving Implicit Material Point Method for SurfaceEnergies with Spatial Gradients
JINGYU CHEN,
UCLA, USA
VICTORIA KALA,
UCLA, USA
ALAN MARQUEZ-RAZON,
UCLA, USA
ELIAS GUEIDON,
UCLA, USA
DAVID A. B. HYDE,
UCLA, USA
JOSEPH TERAN,
UC Davis, USA
Fig. 1. Our method enables the simulation of a wide variety of thermomechanical and surface-tension-driven effects. (Top)
Letter-shaped candles melt andinteract. (Bottom)
A large melting candle; soap spreading on a water surface; water droplets falling and streaking on ramps; partial rebound of a water dropletimpact; wine settling in a glass; a water droplet settling on a hydrophobic surface.
We present a novel Material Point Method (MPM) discretization of surfacetension forces that arise from spatially varying surface energies. These varia-tions typically arise from surface energy dependence on temperature and/orconcentration. Furthermore, since the surface energy is an interfacial prop-erty depending on the types of materials on either side of an interface, spatialvariation is required for modeling the contact angle at the triple junctionbetween a liquid, solid and surrounding air. Our discretization is based onthe surface energy itself, rather than on the associated traction conditionmost commonly used for discretization with particle methods. Our energybased approach automatically captures surface gradients without the explicitneed to resolve them as in traction condition based approaches. We includean implicit discretization of thermomechanical material coupling with anovel particle-based enforcement of Robin boundary conditions associatedwith convective heating. Lastly, we design a particle resampling approachneeded to achieve perfect conservation of linear and angular momentumwith Affine-Particle-In-Cell (APIC) [Jiang et al. 2015]. We show that ourapproach enables implicit time stepping for complex behaviors like the
Authors’ addresses: Jingyu Chen, UCLA, USA; Victoria Kala, UCLA, USA; Alan Marquez-Razon, UCLA, USA; Elias Gueidon, UCLA, USA; David A. B. Hyde, UCLA, USA; JosephTeran, UC Davis, USA.Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than ACMmust be honored. Abstracting with credit is permitted. To copy otherwise, or republish,to post on servers or to redistribute to lists, requires prior specific permission and/or afee. Request permissions from [email protected].© 2021 Association for Computing Machinery.0730-0301/2021/2-ART $15.00https://doi.org/10.1145/nnnnnnn.nnnnnnn
Marangoni effect and hydrophobicity/hydrophilicity. We demonstrate therobustness and utility of our method by simulating materials that exhibithighly diverse degrees of surface tension and thermomechanical effects,such as water, wine and wax.CCS Concepts: •
Mathematics of computing → Discretization ; Partialdifferential equations ; Solvers ; •
Applied computing → Physics .Additional Key Words and Phrases: Surface Tension, Momentum Conserva-tion, Melting, Marangoni Effect, Material Point Methods, Particle-In-Cell
ACM Reference Format:
Jingyu Chen, Victoria Kala, Alan Marquez-Razon, Elias Gueidon, David A. B.Hyde, and Joseph Teran. 2021. A Momentum-Conserving Implicit MaterialPoint Method for Surface Energies with Spatial Gradients.
ACM Trans. Graph.
1, 1 (February 2021), 15 pages. https://doi.org/10.1145/nnnnnnn.nnnnnnn
Surface tension driven flows like those in milk crowns [Zheng et al.2015], droplet coalescence [Da et al. 2016; Li et al. 2020; Thüreyet al. 2010; Wojtan et al. 2010; Yang et al. 2016a] and bubble forma-tion [Da et al. 2015; Huang et al. 2020; Zhu et al. 2014] comprisesome of the most visually compelling fluid motions. Although theseeffects are most dominant at small scales, increasing demand forrealism in computer graphics applications requires modern solverscapable of resolving them. Indeed surface tension effects have beenwell examined in the computer graphics and broader computationalphysics literature. We design a novel approach for simulating sur-face tension driven phenomena that arise from spatial variations in
ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021. a r X i v : . [ c s . G R ] J a n • Jingyu Chen, Victoria Kala, Alan Marquez-Razon, Elias Gueidon, David A. B. Hyde, and Joseph Teran cohesion and adhesion forces at the interface between two liquids.This is often called the Marangoni effect [Scriven and Sternling 1960;Venerus and Simavilla 2015] and perhaps the most famous exampleis the tears of wine phenomenon [Thomson 1855]. Other notableexamples of the Marangoni effect include repulsive flows inducedby a soap droplet on a water surface as well as the dynamics ofmolten waxes and metals [Farahi et al. 2004; Langbein 2002].The spatial variation in the surface forces can be characterizedin terms of the potential energy Ψ 𝑠 associated with surface tension: Ψ 𝑠 = ∫ Γ 𝑘 𝜎 ( x ) 𝑑𝑠 ( x ) . (1)Here the surface tension coefficient 𝑘 𝜎 is proportionate to the rela-tive cohesion and adhesion at the interface between the two fluids.Typically this coefficient is constant across the multi-material in-terface Γ , however with the Marangoni effect the coefficient varieswith x ∈ Γ . These variations are typically driven by temperature orconcentration gradients and give rise to many subtle, but importantvisual behaviors where the variation typically causes fluid to flowaway from low surface energy regions towards high surface energyregions. For example with tears of wine, the spatial variations in thesurface energy arise from inhomogeneity in the mixture of alcoholand water caused by the comparatively rapid evaporation of alcoholand high surface tension of water.There are many existing techniques in the computational physicsliterature that resolve spatial variations in the surface tension. Forparticle-based methods like Smoothed Particle Hydrodynamics (SPH)[Monaghan 1992] and Particle-In-Cell (PIC) [Harlow and Welch1965], most of these approaches are based on the Continuum Sur-face Force (CSF) model of Brackbill et al. [1992]. Marangoni effectshave not been addressed in computer graphics, other than by Huanget al. [2020] where it was examined for material flows in soap films.To capture the Marangoni effect, most approaches do not workwith the potential energy Ψ 𝑠 in Equation (1) but instead base theirdiscretization on its first variation. This variation results in theinterfacial traction condition t = 𝑘 𝜎 𝜅 n + ∇ 𝑆 𝑘 𝜎 . (2)Here t is the force per unit area due to surface tension at the inter-face Γ , ∇ 𝑆 is the surface gradient operator at the interface and 𝜅 and n are the interfacial mean curvature and normal, respectively.The original CSF technique of Brackbill et al. [1992] resolves themean curvature term in Equation (2), but not the surface gradientterm. Tong and Browne [2014] show that the CSF approach can bemodified to resolve the surface gradient. However, while this andother other existing approaches in the SPH and PIC literature arecapable of resolving the spatial variation, none support implicit timestepping for the surface tension forces.We build on the work of Hyde et al. [2020] and show that effi-cient implicit time stepping with Marangoni effects is achievablewith PIC. As in Hyde et al. [2020] we observe that similarities withhyperelasticity suggest that the Material Point Method (MPM) [Sul-sky et al. 1994] is the appropriate version of PIC. We show that by building our discretization from the energy in Equation (1) ratherthan the more commonly adopted traction condition in Equation (2),we can naturally compute the first and second variations of thepotential needed when setting up and solving the nonlinear systemsof equations associated with fully implicit temporal discretization.Interestingly, by basing our discretization on the energy in Equa-tion (1), we also show that no special treatment is required for theinterfacial spatial gradient operator ∇ 𝑠 as was done in e.g. [Tongand Browne 2014]. Furthermore, we show that our approach todiscretizing the Marangoni forces can also be used to impose thecontact angle at liquid/solid/air interfaces [Young 1805]. We showthat this naturally allows for simulation of droplet streaking effects.While our method is a generalization of the MPM technique inHyde et al. [2020], we also improve on its core functionality. TheHyde et al. [2020] approach is characterized by the introductionof additional surface tension particles at each time step which aredesigned to represent the liquid interface Γ and its area weightedboundary normals. These surface tension particles are temporaryand are deleted at the end of the time step to prevent excessivegrowth in particle count or macroscopic particle resampling. How-ever, the surface tension particles are massless to prevent violationof conservation of particle mass and momentum. We show that thisbreaks perfect conservation of grid linear and angular momentumwhen particles introduce grid nodes with no mass that the surfacetension forces will act upon. This is an infrequent occurrence, butbreaks the otherwise perfect conservation of grid linear and angularmomentum expected with conservative MPM forces. We design anovel mass and momentum resampling technique that, with theintroduction of two new types of temporary particles, can restoreperfect conservation of grid linear and angular momentum. We callthese additional temporary particles balance particles. We showthat our novel resampling provides improved behavior over theoriginal approach of Hyde et al. [2020], even in the case of standard,non-Marangoni surface tension effects. Furthermore, although otherresampling techniques exists for PIC methods [Edwards and Bridson2012; Gao et al. 2017b; Yue et al. 2015], we note that ours is the firstto guarantee perfect conservation when using generalized particlevelocities associated with the Affine Particle-in-Cell (APIC) method[Fu et al. 2017; Jiang et al. 2015, 2016].Since variations in surface energy are typically based on temperatureand/or concentration gradients, we couple our surface tension coef-ficients with thermodynamically driven quantities. Furthermore, weresolve solid to liquid and liquid to solid phase changes as a functionof temperature since many Marangoni effects arise from melting andcooling. Notably, we show that our novel conservative resamplingnaturally improves discretization of Robin and Neumann boundaryconditions on the interface Γ needed for convection/diffusion of tem-perature and concentration. In summary, our primary contributionsare: • A novel implicit MPM discretization of spatially varying sur-face tension forces. • A momentum-conserving particle resampling technique forparticles near the surface tension liquid interface.
ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
Momentum-Conserving Implicit Material Point Method for Surface Energies with Spatial Gradients • 3
P2GSplitBoundary Sampling Active GridsParticle Group x p A nq A np m p m q x np x nq v np v nq x np m p A np v np x nq A nq m q v nq m p m p m p A np A np v np v np v np s nr b nr s nr A np b nr m q v nq A nq x nq Fig. 2. Splitting. After surface particles (yellow) are created, the mass andmomentum of the interior MPM particles (blue) that are closest to thesurface particles are immediately distributed. Particles in each particlegroup are assigned equal mass. MPM particles (black) that are not pairedwith any surface particles remain intact for the splitting process. Surfaceparticles (yellow) and balance particles (red) are assigned the same linearvelocity and affine velocity of the original particle (blue). • An implicit MPM discretization of the convection/diffusionevolution of temperature/concentration coupled to the sur-face tension coefficient including a novel particle-based Robinboundary condition.
We discuss relevant particle-based techniques for simulating Marangoniand surface tension effects, contact angle imposition, thermody-namic evolution of temperature and/or concentration, as well asresampling techniques in particle-based methods.
Particle Methods:
Particle-based methods are very effective forcomputer graphics applications requiring discretization of surfacetension forces. Hyde et al. [2020] provide a thorough discussion ofthe state of the art. Our approach utilizes the particle-based MPM[de Vaucorbeil et al. 2020; Sulsky et al. 1994] PIC technique, largelydue to its natural ability to handle self collision [Fei et al. 2018,2017; Guo et al. 2018; Jiang et al. 2017], topology change [Wanget al. 2019; Wolper et al. 2020, 2019], diverse materials [Daviet andBertails-Descoubes 2016; Klár et al. 2016; Ram et al. 2015; Schreckand Wojtan 2020; Stomakhin et al. 2013; Wang et al. 2020c; Yue et al.2015] as well as implicit time stepping with elasticity [Fang et al.2019; Fei et al. 2018; Stomakhin et al. 2013; Wang et al. 2020b]. Weadditionally use the APIC method [Fu et al. 2017; Jiang et al. 2015,2016] for its conservation properties and beneficial suppression ofnoise. Note that our mass and momentum remapping techniqueis designed to work in the context of the APIC techniques whereparticles store generalized velocity information.SPH is very effective for resolving Marangoni effects. The approachesof Tong and Browne [2014] and Hopp-Hirschler et al. [2018] areindicative of the state of the art. Most SPH works rely on the CSFsurface tension model of Brackbill et al. [1992], which transformssurface tension traction into a volumetric force that is only non-zeroalong (numerically smeared) material interfaces. CSF approachesgenerally derive surface normal and curvature estimates as gradientsof color functions, which can be very sensitive to particle distribu-tion. Also, CSF forces are not exactly conservative [Hyde et al. 2020].SPH can also be used to simulate the convection and diffusion oftemperature/concentration that gives rise to the spatial variation
G2PMergeParticle Contributionsto the Grid Momentum m p m q x np x nq A n +1 p v n +1 p v n +1 q A n +1 q p s i αr = ˜ m p N i ( s nr )ˆ v n +1 i α p i αp = ˜ m p N i ( x np )ˆ v n +1 i α p b i αr = ˜ m p N i ( b nr )ˆ v n +1 i α G2P (cid:31) r ∈ Π p (cid:31) i Q i αβγ p s i αr + (cid:31) i Q i αβγ p i αp + (cid:31) r ∈ Π p (cid:31) i Q i αβγ p b i αr = (cid:31) i Q i αβγ m p N i ( x np ) Q i αδ(cid:27) G n +1 δ(cid:27)p x np x nq b nr s nr Fig. 3. Merging. The merging process is a modified version of G2P. For theparticles that are not associated with surface particles (black), a regularG2P is performed. Among each particle group, we calculate each particle’scontribution to the grid momentum and the generalized affine momentsof their summed momenta about their center of mass. Then, we restorethe mass of the original particle associated with the group prior to thesplit and compute its generalized affine inertia tensor from its grid massdistribution. Using the affine inertia tensor of the original particle, wecompute generalized velocity of the particle after the merging from thegeneralized moments of the group. in surface energy in the Marangoni effect. Hu and Eberhard [2017]simulate Marangoni convection in a melt pool during laser weldingand Russell [2018] does so in laser fusion additive manufacturingprocesses. Both approaches use SPH with Tong and Browne [2014]for discretization of Equation (2). Although SPH is very effective forresolving Marangoni effects, all existing approaches utilize implicittreatment of Marangoni forces.
Marangoni effect and contact angle:
The Marangoni effect is vi-sually subtle and has not been resolved with particle-based meth-ods in computer graphics applications. Perhaps the most visuallycompelling example of the Marangoni effect is the tears of winephenomenon on the walls of a wine glass [Scriven and Sternling1960; Venerus and Simavilla 2015]. Tears of wine were simulated inAzencot et al. [2015], however the authors modeled the fluid usingthin film equations under the lubrication approximation and didnot model surface tension gradients. The fingering instabilities theyobserved are stated to occur due to the asymmetric nature of theirinitial conditions. The Marangoni effect was resolved by Huang etal. [2020] recently with thin soap films to generate compelling dy-namics of the characteristic rainbow patterns in bubbles. Relatedly,Ishida et al. [2020] simulated the evolution of soap films includ-ing effects of thin-film turbulence, draining, capillary waves, andevaporation. Outside of computer graphics, the Marangoni effecthas been recently studied in works like Dukler et al. [2020], whichmodels undercompressive shocks in the Marangoni effect, and deLangavant et al. [2017], which presents a spatially-adaptive levelset approach for simulating surfactant-driven flows. Also, Nas andTryggvason [2003] simulated thermocapillary motion of bubblesand drops in flows with finite Marangoni numbers (flows with sig-nificant transport due to Marangoni effects).Our approach for the Marangoni effect also allows for impositionof contact angles at air/liquid/solid interfaces. This effect is impor-tant for visual realism when simulating droplets of water in contactwith solid objects like the ground. Contact angles are influencedby the hydrophilicity/hydrophobicity of the surface on which these
ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021. • Jingyu Chen, Victoria Kala, Alan Marquez-Razon, Elias Gueidon, David A. B. Hyde, and Joseph Teran droplets move or rest [Cassie and Baxter 1944; Johnson Jr. and Dettre1964]. Wang et al. [2007] solve General Shallow Wave Equations,including surface tension boundary conditions and the virtual sur-face method of Wang et al. [2005], in order to model contact anglesand hydrophilicity. Yang et al. [2016b] use a pairwise force model[Tartakovsky and Meakin 2005] for handling contact angles in theirSPH treatment of fluid-fluid and solid-fluid interfaces. Clausen etal. [2013] also consider the relation between surface tension andcontact angles in their Lagrangian finite element approach.
Particle Resampling:
Reseeding or resampling particles is a com-mon concern in various simulation methods; generally speaking,particles need to be distributed with sufficient density near dynamicareas of flow or deformation in order to accurately resolve the dy-namics of the system [Ando et al. 2012; Losasso et al. 2008; Narainet al. 2010]. Edwards and Bridson [2012] use a non-conservativerandom sampling scheme to seed and reseed particles with PIC.Pauly et al. [2005] resample to preserve detail near cracks/fracture,but their resampling does not attempt to conserve momentum. Yueet al. [2015] use Poisson disk sampling to insert new points in low-density regions and merge points that are too close to one anotherwith MPM. However, their resampling method is not demonstratedto be momentum-conserving. A conservative variant of this split-and-merge approach is applied in Gao et al. [2017b] where massand linear momentum are conserved during particle splitting andmerging. However, angular momentum conservation is not con-served. Furthermore, these techniques use PIC not APIC and neitheris designed to guarantee conservation with the generalized velocitystate in APIC techniques [Fu et al. 2017; Jiang et al. 2015, 2016].
Thermomechanical Effects:
Thermodynamic effects in visual sim-ulation date back to at least Terzopoulos et al. [1991]. More recently,melting and resolidification for objects like melting candles havebeen simulated using various methods, including Lattice Boltzmann[Wang et al. 2012] and SPH [Lenaerts and Dutré 2009; Paiva et al.2009], though these results leave room for improved visual and phys-ical fidelity. FLIP methods have also been used for thermodynamicproblems, such as Gao et al. [2017a], which adapts the latent heatmodel from Stomakhin et al. [2014]. Condensation and evaporationof water were considered in several works based on SPH [Hochstet-ter and Kolb 2017; Zhang et al. 2017]. SPH was also applied to theproblem of simulating boiling bubbles in Gu and Yang [2016], whichmodels heat conduction, convection and mass transfer. Recently,particle-based thermodynamics models were incorporated into anSPH snow solver [Gissler et al. 2020a], with temperature-dependentmaterial properties such as the Young’s modulus. In another vein,Pirk et al. [2017] used position-based dynamics and Cosserat physicsto simulate combustion of tree branches, including models for mois-ture and charring. Yang et al. [2017] used the Cahn-Hilliard andAllen-Cahn equations to evolve a continuous phase variable formaterials treated with their phase-field method. Maeshima et al.[2020] considered particle-scale explicit MPM modeling for addi-tive manufacturing (selective laser sintering) that included a latentheat model for phase transition. For a detailed review of thermody-namical effects in graphics, we refer the reader to Stomakhin et al.[2014].
Surface Tension:
Many methods for simulating non-Marangonisurface tension effects have been developed for computer graph-ics applications. We refer the reader to Hyde et al. [2020] for adetailed survey. More recently, Chen et al. [2020] incorporated sub-cell-accurate surface tension forces in an Eulerian fluid frameworkbased on integrating the mean curvature flow of the liquid inter-face (following Sussman and Ohta [2009]). With an eye towardsresolving codimensional flow features, such as thin sheets and fila-ments, Wang et al. [2020a] and Zhu et al. [2014] simulated surfacetension forces using moving-least-squares particles and simplicialcomplexes, respectively. Most related to the present work, Hyde etal. [2020] proposed an implicit material point method for simulatingliquids with large surface energy, such as liquid metals. Their sur-face tension formulation follows [Adamson and Gast 1967; Brackbillet al. 1992; Buscaglia and Ausas 2011] and incorporates a potentialenergy associated with surface tension into the MPM framework.Material boundaries are sampled using massless MPM particles.
We first define the governing equations for thermomechanicallydriven phase change of hyperelastic solids and liquids with variablesurface energy. As in Hyde et al. [2020] we also cover the updatedLagrangian kinematics. Lastly, we provide the variational form ofthe governing equations for use in MPM discretization. We notethat throughout the document Greek subscripts are assumed to runfrom , , . . . , 𝑑 − for the dimension 𝑑 = , of the problem. Re-peated Greek subscripts imply summation, while sums are explicitlyindicated for Latin subscripts. Also, Latin subscripts in bold are usedfor multi-indices. We adopt the continuum assumption [Gonzalez and Stuart 2008]and updated Lagrangian kinematics [Belytschko et al. 2013] used inHyde et al. [2020]. At time 𝑡 we associate our material with subsets Ω 𝑡 ⊂ R 𝑑 , 𝑑 = , . We use Ω to denote the initial configuration ofmaterial with X ∈ Ω used to denote particles of material at time 𝑡 = . A flow map 𝝓 : Ω × [ ,𝑇 ] → R 𝑑 defines the material motionof particles X ∈ Ω to their time 𝑡 locations x ∈ Ω 𝑡 as 𝝓 ( X , 𝑡 ) = x .The Lagrangian velocity is defined by differentiating the flow mapin time V ( X , 𝑡 ) = 𝜕 𝝓 𝜕𝑡 ( X , 𝑡 ) . The La-grangian velocity can be difficult to work with in practice sincereal world observations of material are made in Ω 𝑡 not Ω . TheEulerian velocity v : Ω 𝑡 → R 𝑑 is what we observe in practice.The Eulerian velocity is defined in terms of the inverse flow map 𝝓 − (· , 𝑡 ) : Ω 𝑡 → Ω as v ( x , 𝑡 ) = V ( 𝝓 − ( x , 𝑡 ) , 𝑡 ) where 𝝓 − ( x , 𝑡 ) = X . In general, we can use the flow map and its inverse to pull backquantities defined over Ω 𝑡 and push forward quantities definedover Ω , respectively. For example, given 𝐺 : Ω → R , its pushforward 𝑔 : Ω 𝑡 → R is defined as 𝑔 ( x ) = 𝐺 ( 𝝓 − ( x , 𝑡 )) . Thisprocess is related to the material derivative operator 𝐷𝐷𝑡 where
𝐷𝑔𝐷𝑡 ( x , 𝑡 ) = 𝜕𝐺𝜕𝑡 ( 𝝓 − ( x , 𝑡 )) = 𝜕𝑔𝜕𝑡 ( x , 𝑡 ) + (cid:205) 𝑑 − 𝛼 = 𝜕𝑔𝜕𝑥 𝛼 ( x , 𝑡 ) 𝑣 𝛼 ( x , 𝑡 ) (seee.g. [Gonzalez and Stuart 2008] for more detail). ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
Momentum-Conserving Implicit Material Point Method for Surface Energies with Spatial Gradients • 5
In the updated Lagrangian formalism [Belytschko et al. 2013] wewrite quantities over an intermediate configuration of material Ω 𝑠 with ≤ 𝑠 < 𝑡 . For example, we can define ^ 𝑔 (· , 𝑠 ) : Ω 𝑠 → R as ^ 𝑔 ( ˜x , 𝑠 ) = 𝐺 ( 𝝓 − ( ˜x , 𝑠 )) for ˜x ∈ Ω 𝑠 . As shown in Hyde et al. [2020],this is particularly useful when discretizing momentum balanceusing its variational form. The key observation is that the updatedLagrangian velocity can be written as ^v ( ˜x , 𝑠, 𝑡 ) = V ( 𝝓 − ( ˜x , 𝑠 ) , 𝑡 ) = v ( ^ 𝝓 ( ˜x , 𝑠, 𝑡 ) , 𝑡 ) with ^ 𝝓 ( ˜x , 𝑠, 𝑡 ) = 𝝓 ( 𝝓 − ( ˜x , 𝑠 ) , 𝑡 )) for ˜x ∈ Ω 𝑠 . Intu-itively, ^ 𝝓 (· , 𝑠, 𝑡 ) : Ω 𝑠 → Ω 𝑡 is the mapping from the time 𝑠 con-figuration to the time 𝑡 configuration induced by the flow map.This has a simple relation to the material derivative as 𝜕 ^v 𝜕𝑡 ( ˜x , 𝑠, 𝑡 ) = 𝜕 V 𝜕𝑡 ( 𝝓 − ( ˜x , 𝑠 ) , 𝑡 ) = 𝐷 v 𝐷𝑡 ( ^ 𝝓 ( ˜x , 𝑠, 𝑡 ) , 𝑡 ) . As in Hyde et al. [2020] we willgenerally use upper case for Lagrangian quantities, lower case forEulerian quantities and hat superscripts for updated Lagrangianquantities. The deformation gradient F = 𝜕 𝝓 𝜕 X isdefined by differentiating the flow map in space and can be used toquantify the amount of deformation local to a material point. Weuse 𝐽 = det ( F ) to denote the deformation gradient determinant. 𝐽 represents the amount of volumetric dilation at a material point.Furthermore, it is used when changing variables with integration.We also make use of similar notation for the ^ 𝝓 mapping from Ω 𝑠 to Ω 𝑡 , i.e. ^F = 𝜕 ^ 𝝓 𝜕 ˜x , ^ 𝐽 = det (cid:0) ^F (cid:1) . Our governing equations primarily consist of conservation of massand momentum which can be expressed as 𝜌 𝐷 v 𝐷𝑡 = ∇ · 𝝈 + 𝜌 g , 𝐷𝜌𝐷𝑡 = − 𝜌 ∇ · v , x ∈ Ω 𝑡 (3)where 𝜌 is the Eulerian mass density, v is the Eulerian materialvelocity, 𝝈 is the Cauchy stress and g is gravitational acceleration.Boundary conditions for these equations are associated with a freesurface for solid material, surface tension for liquids and/or pre-scribed velocity conditions. We use 𝜕 Ω 𝑡𝑁 to denote the portion ofthe time 𝑡 boundary subject to free surface or surface tension condi-tions and 𝜕 Ω 𝑡𝐷 to denote the portion of the boundary with Dirichletvelocity boundary conditions. Free surface conditions and surfacetension boundary conditions are expressed as 𝝈 n = t , x ∈ 𝜕 Ω 𝑡𝑁 (4)where t = for free surface conditions and t = 𝑘 𝜎 𝜅 n + ∇ 𝑆 𝑘 𝜎 from Equation (2) for surface tension conditions. Velocity boundaryconditions may be written as v · n = 𝑣 𝑛 bc , x ∈ 𝜕 Ω 𝑡𝐷 . (5) Each material point is either a solid orliquid depending on the thermomechanical evolution. For liquids,the Cauchy stress 𝝈 is defined in terms of pressure and viscousstress: 𝝈 = − 𝑝 I + 𝜇 (cid:18) 𝜕 v 𝜕 x + 𝜕 v 𝜕 x 𝑇 (cid:19) , 𝑝 = − 𝜕 Ψ 𝑝 𝜕𝐽 , with Ψ 𝑝 ( 𝐽 ) = 𝜆 𝑙 ( 𝐽 − ) . Here 𝜆 𝑙 is the bulk modulus of the liquidand 𝜇 𝑙 is its viscosity. For solids, the Cauchy stress is defined in terms of a hyperelastic potential energy density Ψ 𝑠 as 𝝈 = 𝑗 𝜕𝜓 ℎ 𝜕 F f 𝑇 where f ( x , 𝑡 ) = F ( 𝝓 − ( x , 𝑡 ) , 𝑡 ) and 𝑗 ( x , 𝑡 ) = 𝐽 ( 𝝓 − ( x , 𝑡 ) , 𝑡 ) are theEulerian deformation gradient and its determinant, respectively. Weuse the fixed-corotated constitutive model from [Stomakhin et al.2012] for 𝜓 ℎ . This model is defined in terms of the polar SVD [Irvinget al. 2004] of the deformation gradient F = U 𝚺 V 𝑇 with 𝜓 ℎ ( F ) = 𝜇 ℎ 𝑑 − ∑︁ 𝛼 = ( 𝜎 𝛼 − ) + 𝜆 ℎ ( 𝐽 − ) , where the 𝜎 𝛼 are the diagonal entries of 𝚺 and 𝜇 ℎ , 𝜆 ℎ are the hyper-elastic Lamé coefficients. We assume the internal energy of our materials consists of potentialenergy associated with surface tension, liquid pressure and hypere-lasticity and thermal energy associated with material temperature.Conservation of energy together with thermodynamic consider-ations requires convection/diffusion of the material temperature[Gonzalez and Stuart 2008] subject to Robin boundary conditionsassociated with convective heating by ambient material: 𝜌𝑐 𝑝 𝐷𝑇𝐷𝑡 = 𝐾 Δ 𝑇 + 𝐻𝐾 ∇ 𝑇 · n = − ℎ ( 𝑇 − ¯ 𝑇 ) + 𝑏. (6)Here 𝑐 𝑝 is specific heat capacity, 𝑇 is temperature, 𝐾 is thermaldiffusivity, 𝐻 is a source function, ¯ 𝑇 is the temperature of ambientmaterial and n is the surface boundary normal. ℎ controls the rateof convective heating to the ambient temperature and 𝑏 representsthe rate of boundary heating independent of the ambient materialtemperature ¯ 𝑇 .The total potential energy Ψ in our material is as in Hyde et al.[2020], however we include the spatial variation of the surface en-ergy density in Equation (1) and the hyperelastic potential for solidregions: Ψ ( 𝝓 (· , 𝑡 )) = Ψ 𝜎 ( 𝝓 (· , 𝑡 )) + Ψ 𝑙 ( 𝝓 (· , 𝑡 )) + Ψ ℎ ( 𝝓 (· , 𝑡 )) + Ψ 𝑔 ( 𝝓 (· , 𝑡 )) . Here Ψ 𝜎 is the potential from surface tension, Ψ 𝑙 is the potentialfrom liquid pressure, Ψ ℎ is the potential from solid hyperelasticityand Ψ 𝑔 is the potential from gravity. As in typical MPM discretiza-tions, our approach is designed in terms of these energies: Ψ 𝑔 ( 𝝓 (· , 𝑡 ) = ∫ Ω 𝑅 g · 𝝓 𝐽𝑑 X , Ψ 𝜎 ( 𝝓 (· , 𝑡 )) = ∫ 𝜕 Ω 𝑡 𝑘 𝜎 ( x , 𝑡 ) 𝑑𝑠 ( x ) Ψ 𝑙 ( 𝝓 (· , 𝑡 )) = ∫ Ω 𝜆 𝑙 ( 𝐽 − ) 𝑑 X , Ψ ℎ ( 𝝓 (· , 𝑡 )) = ∫ Ω 𝜓 ℎ ( F ) 𝑑 X . Here 𝑅 is the pull back (see Section 3.1) of the mass density 𝜌 . Notethat in the expression for the surface tension potential it is usefulto change variables using the updated Lagrangian view as in Hydeet al. [2020]: Ψ 𝜎 ( 𝝓 (· , 𝑡 )) = ∫ 𝜕 Ω 𝑡 𝑘 𝜎 ( x , 𝑡 ) 𝑑𝑠 ( x ) = ∫ 𝜕 Ω 𝑠 𝑘 𝜎 ( ^ 𝝓 ( ˜x , 𝑠, 𝑡 ) , 𝑡 ) | ^ 𝐽 ^F − 𝑇 ˜n | 𝑑𝑠 ( ˜x ) . (7) ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021. • Jingyu Chen, Victoria Kala, Alan Marquez-Razon, Elias Gueidon, David A. B. Hyde, and Joseph Teran
Here ˜n is the outward unit normal at a point on the boundary of Ω 𝑠 and the expression 𝑑𝑠 ( x ) = | ^ 𝐽 ^F − 𝑇 ˜n | 𝑑𝑠 ( ˜x ) arises by a changeof variables from an integral over Ω 𝑡 to one over Ω 𝑠 . Notably, thespatial variation in 𝑘 𝜎 does not require a major modification of theHyde et al. [2020] approach. The strong form of momentum balance in Equation (3), togetherwith the traction (Equation (4)) and Dirichlet velocity boundaryconditions (Equation (5)), is equivalent to a variational form that isuseful when discretizing our governing equations using MPM. Toderive the variational form, we take the dot product of Equation (3)with an arbitrary function w : Ω 𝑡 → R 𝑑 satisfying w · n = for x ∈ 𝜕 Ω 𝑡𝐷 and integrate over the domain Ω 𝑡 , applying integrationby parts where appropriate. Requiring that the Dirichlet velocityconditions in Equation (4) hold together with the following inte-gral equations for all functions w is equivalent to the strong form,assuming sufficient solution regularity: ∫ Ω 𝑡 𝜌 𝐷𝑣 𝛼 𝐷𝑡 𝑤 𝛼 𝑑 x = − 𝑑𝑑𝜖 𝑃𝐸 (
0; w ) − 𝜇 𝑙 ∫ Ω 𝑡 𝜖 𝑣𝛼𝛽 𝜖 𝑤𝛼𝛽 𝑑 x . (8)Here 𝜖 𝑤 = (cid:16) 𝜕𝑤 𝛼 𝜕𝑥 𝛽 + 𝜕𝑤 𝛽 𝜕𝑥 𝛼 (cid:17) and 𝜖 𝑣 = (cid:16) 𝜕𝑣 𝛼 𝜕𝑥 𝛽 + 𝜕𝑣 𝛽 𝜕𝑥 𝛼 (cid:17) andPE ( 𝜖 ; w ) = Ψ ( 𝝓 (· , 𝑡 ) + 𝜖 W ) . Here W is the pull back of w (see Section 3.1). Note that this no-tation is rather subtle for the surface tension potential energy. Forclarification, Ψ 𝜎 ( 𝝓 (· , 𝑡 ) + 𝜖 W ) = ∫ 𝜕 Ω 𝑠 𝑘 𝜎 ( ^ 𝝓 + 𝜖 ^w )| ^ 𝐽 𝜖, ^w ^F − 𝑇 𝜖, ^w ˜n | 𝑑𝑠 ( ˜x ) , where ^w ( ˜x , 𝑠, 𝑡 ) = w ( ^ 𝝓 ( ˜x , 𝑠, 𝑡 )) is the updated Lagrangian pull backof w , ^F 𝜖, ^w = 𝜕 ^ 𝝓 + 𝜖 ^w 𝜕 ˜x is the deformation gradient of the mapping ^ 𝝓 + 𝜖 ^w and ^ 𝐽 𝜖, ^w = det (cid:0) ^F 𝜖, ^w (cid:1) is its determinant. Lastly, for discretizationpurposes, in practice we change variables in the viscosity term 𝜇 𝑙 ∫ Ω 𝑡 𝜖 𝑣𝛼𝛽 𝜖 𝑤𝛼𝛽 𝑑 x = 𝜇 𝑙 ∫ Ω 𝑠 ^ 𝜖 𝑣𝛼𝛽 ^ 𝜖 𝑤𝛼𝛽 ^ 𝐽𝑑 ˜x . (9)We can similarly derive a variational form of the temperature evo-lution in Equation (6) by requiring ∫ Ω 𝑡 𝜌𝑐 𝑝 𝐷𝑇𝐷𝑡 𝑞𝑑 x = − ∫ 𝜕 Ω 𝑡 𝑞ℎ𝑇𝑑 S ( x ) + ∫ 𝜕 Ω 𝑡 𝑞 (cid:0) ℎ ¯ 𝑇 + 𝑏 (cid:1) 𝑑 S ( x )+ ∫ Ω 𝑡 𝐻𝑞𝑑 x − ∫ Ω 𝑡 ∇ 𝑞 · 𝐾 ∇ 𝑇𝑑 x (10)for all functions 𝑞 : Ω 𝑡 → R . The thermomechanical material dependence is modeled by allowingthe surface tension coefficient 𝑘 𝜎 , the liquid bulk modulus 𝜆 𝑙 , theliquid viscosity 𝜇 𝑙 and the hyperelastic Lamé coefficients 𝜇 ℎ , 𝜆 ℎ tovary with temperature. When 𝑇 exceeds a user-specified meltingpoint 𝑇 melt , the solid phase is changed to liquid and the deformationgradient determinant 𝐽 is set to . Similarly, if liquid temperature drops below 𝑇 melt , The phase is updated to be hyperelastic solidand the deformation gradient F is set to the identity matrix. Inpractice, resetting the deformation gradient and its determinanthelps to prevent nonphysical popping when the material changesfrom solid to liquid and vice versa. We remark that incorporating amore sophisticated phase change model, such as a latent heat buffer,is potentially useful in future work [Stomakhin et al. 2014]. The contact angle between a liquid, a solid boundary and the ambientair is governed by the Young equation [Young 1805]. This expressionrelates the resting angle 𝜃 (measured through the liquid) of a liquidin contact with a solid surface to the surface tension coefficientsbetween the liquid, solid and air phases: 𝑘 𝜎𝑆𝐺 = 𝑘 𝜎𝑆𝐿 + 𝑘 𝜎𝐿𝐺 cos ( 𝜃 ) . (11)The surface tension coefficients are between the solid and gas phases,solid and liquid phases, and liquid and gas phases, respectively. Asin Clausen et al. [2013], we assume 𝑘 𝜎𝑆𝐺 is negligible since we areusing a free surface assumption and do not explicitly model the air.Under this assumption, the solid-liquid contact angle is determinedby the surface tension ratio − 𝑘 𝜎𝑆𝐿 / 𝑘 𝜎𝐿𝐺 . We note that, while onewould expect surface tension coefficients/energies to be positive,this ratio can be negative under the assumption of zero solid-gassurface tension. Furthermore, we note that utilizing this expressionrequires piecewise constant surface tension coefficients where thevariation along the liquid boundary is based on which portion isin contact with the air and which is in contact with the solid. Thedistinct surface tension coefficients on different interfaces providecontrollability of the spreading behavior of the liquid on the solidsurface. As in Hyde et al. [2020], we use MPM [Sulsky et al. 1994] and APIC[Jiang et al. 2015] to discretize the governing equations. The domain Ω 𝑡 𝑛 at time 𝑡 𝑛 is sampled using material points x 𝑛𝑝 . These pointsalso store approximations of the deformation gradient determinant 𝐽 𝑛𝑝 , constant velocity v 𝑛𝑝 , affine velocity A 𝑛𝑝 , volume 𝑉 𝑝 , mass 𝑚 𝑝 = 𝜌 ( x 𝑝 , 𝑡 ) 𝑉 𝑝 , temperature 𝑇 𝑛𝑝 , and temperature gradient ∇ 𝑇 𝑛𝑝 . Wealso make use of a uniform background grid with spacing Δ 𝑥 whendiscretizing momentum updates. To advance our state to time 𝑡 𝑛 + ,we use the following steps:(1) Resample particle boundary for surface tension and Robinboundary temperature conditions.(2) P2G: Conservative transfer of momentum and temperaturefrom particles to grid.(3) Update of grid momentum and temperature.(4) G2P: Conservative transfer of momentum and temperaturefrom grid to particles. The integrals associated with the surface tension energy in Equa-tion (7) and the Robin temperature condition in Equation (10) aredone over the boundary of the domain. We follow Hyde et al. [2020]and introduce special particles to cover the boundary in order to
ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
Momentum-Conserving Implicit Material Point Method for Surface Energies with Spatial Gradients • 7
Particle Group Π p b nr x np s nr Fig. 4. A portion of an MPM fluid in the simulation domain. Surface par-ticles (yellow) are sampled on faces of the zero isocontour of the level setformed by unioning spherical level sets around each MPM particle. Eachsurface particle generates an associated balance particle (red) such that theclosest MPM particle (blue) to a boundary particle lies on the midpoint of aline segment between the surface particle and balance particle. A single blueparticle at 𝒙 𝑝 may be paired with multiple surface particles and balanceparticles, and they are considered to be in a particle group 𝚷 𝑝 . MPM parti-cles that are not associated with any surface tension particles are markedas black. serve as quadrature points for these integrals. As in Hyde et al.[2020] these particles are temporary and are removed at the end ofthe time step. However, while Hyde et al. [2020] used massless sur-face particles, we design a novel conservative mass and momentumresampling for surface particles. Massless particles easily allow formomentum conserving transfers from particle to grid and vice versa;however, they can lead to loss of conservation in the grid momen-tum update step. This occurs when there is a grid cell containingonly massless particles. In this case, there are grid nodes with nomass that receive surface tension forces. These force componentsare then effectively thrown out since only grid nodes with mass willaffect the end of time step particle momentum state (see Section 6.1).We resolve this issue by assigning mass to each of the surface parti-cles. However, to conserve total mass, some mass must be subtractedfrom interior MPM particles. Furthermore, changing the mass ofexisting particles also changes their momentum, which may lead toviolation of conservation. In order to conserve mass, linear momen-tum and angular momentum, we introduce a new particle for eachsurface particle. We call these balance particles, and like surfaceparticles they are temporary and will be removed at the end of thetime step. We show that the introduction of these balance particlesnaturally allows for conservation both when they are created at thebeginning of the time step and when they are removed at the endof the time step. We first introduce surface particlesusing the approach in Hyde et al. [2020]. A level set enclosing theinterior MPM particles is defined as the union of spherical level setsdefined around each interior MPM particle. Unlike Hyde et al. [2020],we do not smooth or shift the unioned level set. We compute thezero isocontour of the level set using marching cubes [Chernyaev1995] and randomly sample surface particles along this explicitrepresentation. In Hyde et al. [2020], three-dimensional boundarieswere sampled using a number of sample points proportional to the surface area of each triangle. Sample points were computedusing uniform random barycentric weights, which leads to a non-uniform distribution of points in each triangle. Instead, we employa strategy of per-triangle Monte Carlo sampling using a robustPoisson distribution, as described in Corsini et al. [2012], whereuniform triangle sample points are generated according to Osada etal. [2002]. We found that this gave better coverage of the boundarywithout generating particles that are too close together (see Figure5). We note that radii for the particle level sets are taken to be . Δ 𝑥 (slightly larger than √ Δ 𝑥 ) in 2D and . Δ 𝑥 (slightly larger than √ Δ 𝑥 ) in 3D. This guarantees that even a single particle in isolationwill always generate a level set zero isocontour that intersects thegrid and will therefore always generate boundary sample points.Note also that as in Hyde et al. [2020], we use the explicit marchingcubes mesh of the zero isocontour to easily and accurately generatesamples of area weighted normals 𝑑 A 𝑟 where (cid:205) | 𝑑 A 𝑟 | ≈ ∫ 𝑡 𝑛 Ω 𝑑 x are chosen with direction from the triangle normal and magnitudebased on the number of samples in a given triangle and the trianglearea. For each surface particle s 𝑛𝑟 , weadditionally generate a balance particle b 𝑛𝑟 . First, we compute theclosest interior MPM particle for each surface particle x 𝑛𝑝 ( s 𝑛𝑟 ) . Thenwe introduce the corresponding balance particle as b 𝑛𝑟 = s 𝑛𝑟 + (cid:16) x 𝑛𝑝 ( s 𝑛𝑟 ) − s 𝑛𝑟 (cid:17) . (12) After introducing the surface s 𝑛𝑟 and balance b 𝑛𝑟 particles, we assign them mass and momentum(see Figure 2). To achieve this in a conservative manner, we firstpartition the surface particles into particle groups Π 𝑝 defined as theset of surface particle indices 𝑟 such that x 𝑛𝑝 is the closest interiorMPM particle to s 𝑛𝑟 (see Figure (4)). We assign the mass 𝑚 𝑝 of theparticle x 𝑛𝑝 to the collection of x 𝑛𝑝 , s 𝑛𝑟 and b 𝑛𝑟 for 𝑟 ∈ Π 𝑝 uniformlyby defining a mass of ˜ 𝑚 𝑝 = 𝑚 𝑝 | Π 𝑝 |+ to each surface and balancepoint as well as to x 𝑛𝑝 . Here | Π 𝑝 | is the number of elements in the set.This operation is effectively a split of the original particle x 𝑛𝑝 withmass 𝑚 𝑝 into a new collection of particles x 𝑛𝑝 , s 𝑛𝑟 , b 𝑛𝑟 , 𝑟 ∈ Π 𝑝 withmasses ˜ 𝑚 𝑝 . This split trivially conserves the mass. Importantly, byconstruction of the balance particles (Equation (12)) we ensure thatthe center of mass of the collection is equal to the original particle x 𝑛𝑝 : 𝑚 𝑝 (cid:169)(cid:173)(cid:171) ˜ 𝑚 𝑝 x 𝑛𝑝 + ∑︁ 𝑟 ∈ Π 𝑝 ˜ 𝑚 𝑝 s 𝑛𝑟 + ˜ 𝑚 𝑝 b 𝑛𝑟 (cid:170)(cid:174)(cid:172) = x 𝑛𝑝 . (13)With this particle distribution, conservation of linear and angularmomentum can be achieved by simply assigning each new particlein the collection the velocity v 𝑛𝑝 and affine velocity A 𝑛𝑝 of the originalparticle x 𝑛𝑝 . We note that the conservation of the center of mass(Equation (13)) is essential for this simple constant velocity split toconserve linear and angular momentum (see [Annonymous 2021]). ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021. • Jingyu Chen, Victoria Kala, Alan Marquez-Razon, Elias Gueidon, David A. B. Hyde, and Joseph Teran
Fig. 5. Isocontour and sampled boundary particles for an ellipsoid. (Left)
Using the method of Hyde et al. [2020]. Note how low-quality triangles areundersampled and how sample points often clump near triangle centers. (Right)
The present method, which does not suffer from similar issues.
After the addition of the surface and balance particles, we transfermass and momentum to the grid in the standard APIC [Jiang et al.2015] way using their conservatively remapped mass and velocitystate 𝑚 𝑛 i = ∑︁ 𝑝 ˜ 𝑚 𝑝 (cid:169)(cid:173)(cid:171) 𝑁 i ( x 𝑛𝑝 ) + ∑︁ 𝑟 ∈ Π 𝑝 𝑁 i ( s 𝑛𝑟 ) + 𝑁 i ( b 𝑛𝑟 ) (cid:170)(cid:174)(cid:172) ,𝑚 𝑛 i v 𝑛 i = ∑︁ 𝑝 ˜ 𝑚 𝑝 𝑁 i ( x 𝑛𝑝 ) (cid:16) v 𝑛𝑝 + A 𝑛𝑝 ( x i − x 𝑛𝑝 ) (cid:17) + ∑︁ 𝑝 ˜ 𝑚 𝑝 ∑︁ 𝑟 ∈ Π 𝑝 𝑁 i ( s 𝑛𝑟 ) (cid:16) v 𝑛𝑝 + A 𝑛𝑝 ( x i − s 𝑛𝑟 ) (cid:17) + ∑︁ 𝑝 ˜ 𝑚 𝑝 ∑︁ 𝑟 ∈ Π 𝑝 𝑁 i ( b 𝑛𝑟 ) (cid:16) v 𝑛𝑝 + A 𝑛𝑝 ( x i − b 𝑛𝑟 ) (cid:17) . Here N i ( x ) = N ( x − x i ) are quadratic B-splines defined over theuniform grid with x i living at cell centers [Stomakhin et al. 2013].Note that for interior MPM particles far enough from the boundarythat Π 𝑝 = ∅ . This reduces to the standard APIC [Jiang et al. 2015]splat since ˜ 𝑚 𝑝 = 𝑚 𝑛𝑝 . We also transfer temperature from particlesto grid using 𝑇 𝑛 i 𝛼 = ∑︁ 𝑝 𝑚 𝑝 𝑁 i ( x 𝑛𝑝 )( 𝑇 𝑛𝑝 + ( 𝑥 i 𝛼 − 𝑥 𝑛𝑝𝛼 )∇ 𝑇 𝑛𝑝𝛼 ) . Note that for the temperature transfer, we only use surface particlesto properly apply the thermal boundary conditions, and we do notuse these particles to transfer mass-weighted temperature to thegrid.
We discretize the governing equations in the standard MPM mannerby using the particles as quadrature points in the variational forms.The interior MPM particles x 𝑛𝑝 are used for volume integrals and the surface particles s 𝑛𝑟 are used for surface integrals. By choosing 𝑠 = 𝑡 𝑛 , 𝑡 = 𝑡 𝑛 + and by using grid discretized versions of ^w ( ˜x ) = (cid:205) j w j 𝑁 j ( ˜x ) , ^v ( ˜x , 𝑡 𝑛 , 𝑡 𝑛 + ) = (cid:205) i ^v 𝑛 + 𝑁 i ( ˜x ) , ^ 𝑞 ( ˜x ) = (cid:205) j 𝑞 j 𝑁 j ( ˜x ) and ^ 𝑇 ( ˜x ) = (cid:205) i ^ 𝑇 i 𝑁 i ( ˜x ) . As in Hyde et al. [2020], the grid momen-tum update is derived from Equation (8): 𝑚 𝑛 i ^v 𝑛 + − v 𝑛 i Δ 𝑡 = f i ( x + Δ 𝑡 ^q ) + 𝑚 𝑛 i g , (14) f i ( ^x ) = − 𝜕𝑒𝜕 ^x i ( ^x ) − 𝜇 𝑙 ∑︁ 𝑝 𝝐 𝑣 ( ^x; x 𝑛𝑝 ) (cid:18) 𝜕𝑁 i 𝜕 x ( x 𝑛𝑝 ) (cid:19) 𝑇 𝑉 𝑛𝑝 , (15)where f i is the force on grid node i from potential energy and vis-cosity, 𝜖 𝑣 ( ^x; x 𝑛𝑝 ) = (cid:18)(cid:205) j ^x j 𝜕𝑁 j 𝜕 x ( x 𝑛𝑝 ) + (cid:16) ^x j 𝜕𝑁 j 𝜕 x ( x 𝑛𝑝 ) (cid:17) 𝑇 (cid:19) is the strainrate at x 𝑛𝑝 , g is gravity, and ^q is either (for explicit time integra-tion) or ^v 𝑛 + (for backward Euler time integration). x represents thevector of all unmoved grid node positions x i . We use 𝑒 ( y ) to denotethe discrete potential energy Ψ where MPM and surface particlesare used as quadrature points: 𝑒 ( y ) = ∑︁ 𝑝 (cid:32) 𝜓 ℎ ( F 𝑝 ( ^y )) + 𝜆 𝑙 ( 𝐽 𝑝 ( ^y ) − ) (cid:33) 𝑉 𝑝 + ∑︁ 𝑟 𝑘 𝜎 ( s 𝑛𝑟 )| ^ 𝐽 𝑟 ( ^y ) ^F − 𝑇𝑟 ( ^y ) 𝑑 A 𝑛𝑟 | , where, as in [Stomakhin et al. 2013], F 𝑝 ( ^y ) = (cid:205) i y i 𝜕𝑁 i 𝜕 x ( x 𝑛𝑝 ) F 𝑛𝑝 and as in Hyde et al. [2020], 𝐽 𝑝 ( ^y ) = (cid:16) − 𝑑 + 𝑦 𝛼 𝜕𝑁 i 𝜕𝑥 𝛼 ( x 𝑛𝑝 ) (cid:17) 𝐽 𝑛𝑝 and ^F 𝑝 ( y ) = (cid:205) i y i 𝜕𝑁 i 𝜕 x ( x 𝑛𝑝 ) . With these conventions, the 𝛼 componentof the energy-based force on grid node i is of the form − 𝜕𝑒𝜕𝑥 i 𝛼 ( y ) = − ∑︁ 𝑝 𝜕𝜓 ℎ 𝜕𝐹 𝛼𝛿 ( F 𝑝 ( ^y )) 𝐹 𝑛𝑝𝛾𝛿 𝜕𝑁 i 𝜕𝑥 𝛾 ( x 𝑛𝑝 ) 𝑉 𝑝 − ∑︁ 𝑝 𝜆 𝑙 ( 𝐽 𝑝 ( y ) − ) 𝜕𝑁 i 𝜕𝑥 𝛼 ( x 𝑛𝑝 ) 𝐽 𝑛𝑝 𝑉 𝑝 − ∑︁ 𝑟 𝑘 𝜎 ( s 𝑛𝑟 ) 𝜕 | det (cid:16) ^G 𝑟 (cid:17) ^G − 𝑇𝑟 𝑑 A 𝑛𝑟 | 𝜕 ^ 𝐺 𝛼𝛿 ( ^F 𝑟 ( ^y )) 𝜕𝑁 i 𝜕𝑥 𝛿 ( x 𝑛𝑝 ) . (16) We note that the viscous contribution to the force in Equation (15)is the same as in Ram et al. [2015]. We would expect 𝑉 𝑛 + 𝑝 inthis term when deriving from Equation (9), however we approxi-mate it as 𝑉 𝑛𝑝 . This is advantageous since it makes the term lin-ear; and since ^ 𝐽 𝑝 ( ^x ) ≈ from the liquid pressure and hypere-lastic stress, it is not a poor approximation. Lastly, we note thatthe surface tension coefficient 𝑘 𝜎 ( s 𝑛𝑟 ) will typically get its spa-tial dependence from composition with a function of temperature 𝑘 𝜎 ( s 𝑛𝑟 ) = ˜ 𝑘 𝜎 ( ^ 𝑇 ( s 𝑛𝑟 )) = ˜ 𝑘 𝜎 ( T 𝑠,𝑛𝑟 ) .In the case of implicit time stepping with backward Euler ( ^q = ^v 𝑛 + ),we use Netwon’s method to solve the nonlinear systems of equations.This requires linearization of the grid forces associated with poten-tial energy in Equation (16). We refer the reader to Stomakhin et al.[2013] and Hyde et al. [2020] for the expressions for these terms, aswell as the definiteness fix used for surface tension contributions. ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
Momentum-Conserving Implicit Material Point Method for Surface Energies with Spatial Gradients • 9
We discretize Equation 10 in a simi-lar manner which results in the following equations for the gridtemperatures 𝑇 i : 𝑐 𝑝 𝑚 i ^ 𝑇 𝑛 + − 𝑇 𝑛 i Δ 𝑡 = − ∑︁ 𝑝 𝐾 𝜕𝑁 i 𝜕𝑥 𝛼 ( x 𝑛𝑝 ) ^ 𝑇 𝑛 + 𝜕𝑁 j 𝜕𝑥 𝛼 ( x 𝑛𝑝 ) 𝑉 𝑛𝑝 − ∑︁ 𝑟 ℎ𝑁 i ( s 𝑛𝑟 ) ^ 𝑇 𝑛 + 𝑁 j ( s 𝑛𝑟 )| 𝑑 A 𝑛𝑟 |+ ∑︁ 𝑟 𝑁 i ( s 𝑛𝑟 ) (cid:2) ℎ ¯ 𝑇 ( s 𝑛𝑟 ) + 𝑏 ( s 𝑛𝑟 ) (cid:3) | 𝑑 A 𝑛𝑟 |+ ∑︁ 𝑝 ℎ𝑁 i ( x 𝑛𝑝 ) 𝐻 ( x 𝑛𝑝 ) 𝑉 𝑛𝑝 . Note that by using the surface particles s 𝑛𝑟 as quadrature pointsin the variational form, the Robin boundary condition can be dis-cretized naturally with minimal modification to the Laplacian andtime derivative terms. Also note that inclusion of this term modifiesboth the matrix and the right side in the linear system for ^ 𝑇 𝑛 + . Wefound that performing constant extrapolation of interior particletemperatures to the surface particles provided better initial guessesfor the linear solver. Once grid momentum and temperature have been updated, wetransfer velocity and temperature back to the particles. For interiorMPM particles with no associated surface or balance particles ( Π 𝑝 = ∅ ), we transfer velocity, affine velocity and temperature from gridto particles in the standard APIC [Jiang et al. 2015] way: v 𝑛 + 𝑝 = ∑︁ i 𝑁 i ( x 𝑛𝑝 ) ^v 𝑛 + , A 𝑛 + 𝑝 = Δ 𝑥 ∑︁ i 𝑁 i ( x 𝑛𝑝 ) ^v 𝑛 + ( x i − x 𝑛𝑝 ) 𝑇 . For interior MPM particles that were split with a collection of surfaceand balance particles ( Π 𝑝 ≠ ∅ ), more care must be taken sincesurface and balance particles will be deleted at the end of the timestep. First, the particle is reassigned its initial mass 𝑚 𝑝 . Then wecompute the portion of the grid momentum associated with eachsurface and balance particle associated with 𝑝 as p 𝑠 i 𝑟 = ˜ 𝑚 𝑝 𝑁 i ( s 𝑛𝑟 ) ^v 𝑛 + , p 𝑏 i 𝑟 = ˜ 𝑚 𝑝 𝑁 i ( b 𝑛𝑟 ) ^v 𝑛 + , 𝑟 ∈ Π 𝑝 . We then sum this with the split particle’s share of the grid momen-tum to define the merged particle’s share of the grid momentum p i 𝑝 = ˜ 𝑚 𝑝 𝑁 i ( x 𝑛𝑝 ) ^v 𝑛 + + ∑︁ 𝑟 ∈ Π 𝑝 p 𝑠 i 𝑟 + p 𝑏 i 𝑟 . Note that the p i 𝑝 may be nonzero for more grid nodes than theparticle would normally splat to (see Figure 3). We define the par-ticle velocity from the total momentum by dividing by the mass v 𝑛 + 𝑝 = 𝑚 𝑝 (cid:205) i p i 𝑝 . To define the affine particle velocity, we use ageneralization of Fu et al. [2017] and first compute the generalizedaffine moments 𝑡 𝑝𝛽𝛾 = (cid:205) i 𝑄 i 𝛼𝛽𝛾 𝑝 i 𝑝𝛼 of the momentum distribution 𝑝 i 𝑝𝛼 where 𝑄 i 𝛼𝛽𝛾 = 𝑟 i 𝑝𝛾 𝛿 𝛼𝛽 is the 𝛼 component of the 𝛽𝛾 linearmode at grid node i . Here r i 𝑝 = x i − x 𝑛𝑝 is the displacement from thecenter of mass of the distribution to the grid node x i . We note thatthese moments are the generalizations of angular momentum toaffine motion, as was observed in [Jiang et al. 2015], however in ourcase we compute the moments from a potentially wider distribution of momenta 𝑝 i 𝑝𝛼 . Lastly, to conserve angular momenta (see [An-nonymous 2021] for details), we define the affine velocity by invert-ing the generalized affine inertia tensor (cid:205) i 𝑄 i 𝛼𝛾𝛿 𝑚 𝑝 𝑁 i ( x 𝑛𝑝 ) 𝑄 i 𝛼𝜖𝜏 of the point x 𝑛𝑝 using its merged mass distribution 𝑚 𝑝 𝑁 i ( x 𝑛𝑝 ) . How-ever, as noted in [Jiang et al. 2015], the generalized inertia tensor 𝑚 𝑝 Δ 𝑥 I is constant diagonal when using quadratic B-splines for 𝑁 i ( x 𝑛𝑝 ) and therefore the final affine velocity is A 𝑛 + 𝑝 = 𝑚 𝑝 Δ 𝑥 t 𝑝 .Temperature and temperature gradients are transferred in the sameway whether or not a MPM particle was split or not: 𝑇 𝑛 + 𝑝 = ∑︁ i ^ 𝑇 𝑛 + 𝑁 ( x i ) , ∇ 𝑇 𝑛 + 𝑝 = ∑︁ i ^ 𝑇 𝑛 + ∇ 𝑁 ( x i ) . To demonstrate our method’s ability to fully conserve momentumand center of mass, we simulate a two-dimensional ellipse thatoscillates under zero gravity due to surface tension forces, andwe compare these results to those obtained using the method ofHyde et al. [2020]. As seen in Figure 6, both methods are dissipativein terms of kinetic energy due to the APIC transfer [Jiang et al.2015]. However, the present technique perfectly conserves totallinear momentum, total angular momentum, and the center of massof the ellipse, unlike Hyde et al. [2020]. In this example, explicittime stepping was used, and the CFL time step was restricted to bebetween × − and × − seconds. A grid resolution of Δ 𝑥 = / was used, along with a bulk modulus of . and a surfacetension coefficient of 𝑘 𝜎 = . . 8 particles per cell were randomlysampled, then filtered to select only the ones inside the ellipse;this resulted in a center of mass for the ellipse at approximately ( . , . ) . Frame T o t a l K i n e t i c E n e r g y P i m i v i Hyde Ours
Frame T o t a l L i n e a r M o m . P i m i v i -3 Hyde(x) Hyde(y) Ours(x) Ours(y)
Frame -2-1012 T o t a l A n g u l a r M o m . P i x i m i v i -3 Hyde Ours
Frame C e n t e r o f M a ss P o s i t i o n Hyde(x) Hyde(y) Ours(x) Ours(y)
Fig. 6. The present method (blue) conserves total mass, total linear andangular momentum, and center of mass, unlike Hyde et al. [2020] (red).Both methods are energy-dissipative.
ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
We demonstrate our method’s ability to handle highly dynamicsimulations with a wide range of surface tension strengths. We sim-ulated several spherical droplets with the same material parameters,each of which free falls from a fixed height and impacts a dry, fric-tionless, hydrophobic surface. The bulk modulus is 83333.33 andthe gravitational acceleration is 9.8. The same grid resolution of Δ 𝑥 = / was used for all the simulations, and the time step wasrestricted between − to × − seconds by the CFL condition.With different surface tension coefficient 𝑘 𝜎 , the droplet showeddistinct behaviors upon impact, as shown in Figure 7. Fig. 7. (Top)
Spherical droplets with different surface tension coefficient freefall from the same height. In the top figure, from left to right, the surfacetension coefficients are 𝑘 𝜎 = , , , . , . . (Middle Row) full rebound ofthe droplet (initial height: . and 𝑘 𝜎 = ). (Bottom Row) partial reboundof the droplet (initial height: . and 𝑘 𝜎 = ). We also captured the partial rebound and the full rebound behaviorsof the droplet after the impact. The middle and bottom rows ofFigure 7 show the footage of a droplet with 𝑘 𝜎 = dropped from aheight of . and a droplet with 𝑘 𝜎 = dropped from a height of . ,respectively. With a higher surface tension coefficient and a higherimpact speed, the droplet is able to completely leave the surfaceafter the impact. Our results qualitatively match the experimentoutcomes from Rioboo et al. [2001]. As discussed in Section 3.6, our method allows for distinct 𝑘 𝜎 valuesat solid-liquid and liquid-air interfaces. Tuning the ratio between 𝑘 𝜎 at these interfaces allows simulating different levels of hydrophilic-ity/hydrophobicity. Figure 8 shows an example of several liquiddrops with different 𝑘 𝜎 ratios falling on ramps of . ◦ angle. Thegrid resolution was Δ 𝑥 = / , and the time step was restrictedbetween − and − seconds by the CFL condition. Coulombfriction with a friction coefficient of . was used for the rampsurface. When there is a larger difference between solid-liquid andliquid-air surface tension coefficients (i.e., a smaller 𝑘 𝜎 ratio), theliquid tends to drag more on the surface and undergo more sepa-ration and sticking. The leftmost example, with a 𝑘 𝜎 ratio of . ,exhibits hydrophobic behavior. Fig. 8. Liquid drops fall on a ramp with varying ratios between the solid-liquid and liquid-air surface tension coefficients. From left to right: ratios of . , . , . , . . (Top) Frame 60. (Bottom)
Frame 100.
The two-dimensional lid-driven cavity is a classic example in theengineering literature of the Marangoni effect [Francois et al. 2006;Hopp-Hirschler et al. 2018]. Inspired by works like these, we simu-late a square unit domain and fill the domain with particles up toheight − Δ 𝑥 ( Δ 𝑥 = / ), which results in a free surface nearthe top of the domain. A linear temperature gradient from 1 onthe left to 0 on the right is initialized on the particles. To achievethe Marangoni effect, the surface tension coefficient 𝑘 𝜎 is set todepend linearly on temperature: 𝑘 𝜎 = − 𝑇 𝑝 . 𝑘 𝜎 is clamped to be ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
Momentum-Conserving Implicit Material Point Method for Surface Energies with Spatial Gradients • 11 in [ , ] to avoid artifacts due to numerical precision. Gravity is setto zero, dynamic viscosity is set to × − and implicit MPM isused with a maximum Δ 𝑡 of . . Results are shown in Figure 9.We note that the center of the circulation drifts to the right overthe course of the simulation due to uneven particle distributionresulting from the circulation of the particles; investigating particlereseeding strategies to stabilize the flow is interesting future work. Fig. 9. Frame 500 of a two-dimensional lid-driven cavity simulation. Thesimulation is initially stationary, but velocity streamlines (red) show theflow pattern characteristic of Marangoni convection that develops due to atemperature-dependent surface tension coefficient. The contour plot showsthe evolving temperature field (initially a linear horizontal distribution).
Figure 10 shows that our method enables simulation of variouscontact angles, emulating various degrees of hydrophobic or hy-drophilic behavior as a droplet settles on a surface. We adjust thecontact angles by assigning one surface tension coefficient, 𝑘 𝜎𝐿𝐺 , tothe surface particles on the liquid-gas interface, and another one, 𝑘 𝜎𝑆𝐿 , to those on the solid-liquid interface. Following the Youngequation (Equation (11)) and our assumption that 𝑘 𝜎𝑆𝐺 is negligible,the contact angle is given by 𝜃 = arccos (cid:0) − 𝑘 𝜎𝑆𝐿 / 𝑘 𝜎𝐿𝐺 (cid:1) . Note thatthe effect of gravity will result in contact angles slightly smallerthan targeted. The grid resolution was set to Δ 𝑥 = / , and thetime step restricted between . and × − . Each droplet isdiscretized using k interior particles and k surface particles.The bulk modulus of the fluid is . , 𝑘 𝜎𝐿𝐺 is set to , and thedynamic viscosity is . . We demonstrate a surface tension driven flow by simulating the soapreducing the surface tension of the water. We initialize a × . × rectangular water pool and identify the particles in the middle of thetop surface as liquid soap. Boundary particles associated with thewater particles have higher surface tension than those associatedwith the soap particles. In order to visualize the effect of the surfacetension driven flow, we randomly selected marker particles on thetop surface of the pool. Due to the presence of the soap, the centerof the pool has lower surface tension than the area near the edge Fig. 10. As our droplets settle, we are able to obtain contact angles of ap-proximately 45, 90, 135 and 180 degrees, using a 𝑘 𝜎𝑆𝐿 / 𝑘 𝜎𝐿𝐺 ratio of −√ / , , √ / and , respectively. of the container. Figure 11 shows footage of this process. A gridresolution of Δ 𝑥 = / was used, and the time step was limitedbetween − and × − by the CFL number. The bulk modulusof water is set to be . . The surface tension 𝑘 𝜎 for water isset to be . and 𝑘 𝜎 for soap is . . Fig. 11. The soap in the center of the pool surface reduces the surfacetension. The surface tension gradient drives the markers towards the wallsof the container. Frames 0, 10, 20, 40 are shown in this footage.
We consider an example of wine flowing on the surface of a pre-wetted glass. The glass is an ellipsoid centered at ( . , . , . ) withcharacteristic dimensions 𝑎 = . , 𝑏 = . , and 𝑐 = . . We initial-ize a thin band of particles with thickness of Δ 𝑥 on the surfaceof the wine glass and observe the formation of ridges and fingersas the particles settle toward the bulk fluid in the glass. The gridresolution is Δ 𝑥 = / and a maximum allowable time step is − constrained by the CFL number. We set the surface tensioncoefficient on the liquid-gas interface is 𝑘 𝜎𝐿𝐺 = . and the one on ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021. the solid-liquid interface is 𝑘 𝜎𝑆𝐿 = . . The piecewise constantsurface tension gradient leads to a more prominent streaking behav-iors of the liquid on the glass wall. The results are shown in Figure12. Fig. 12. Wine is initialized in a glass with part of the interior pre-wetted.The falling wine forms tears and ridges, and the tears eventually connectwith the bulk fluid. Frames 30 and 90 are shown.
We simulate several scenarios with wax candles. In these examples,wax melts due to a heat source (candle flame) and resolidifies whenit flows away from the flame. Ambient temperature ^ 𝑇 is taken to be K, and the melting point is
K. Thermal diffusivity 𝐾 is takento be . , and specific heat capacity 𝑐 𝑝 is set to . No internal heatsource is used ( 𝐻 = ); instead, heating and cooling are applied onlyvia the boundary conditions.To simulate the candle wicks, we manually construct and samplepoints on cubic splines. As the simulation progresses, we delete par-ticles from the wick that are too far above the highest ( 𝑦 -direction)liquid particle within a neighborhood of the wick. The flames arecreated by running a separate FLIP simulation as a postprocessand anchoring the result to the exposed portion of each wick. Werendered these scenes using Arnold [Georgiev et al. 2018] and post-processed the renders using the NVIDIA OptiX denoiser (based onChaitanya et al. [2017]).We consider the effect of varying 𝑘 𝜎 on the overall behavior ofthe flow. Figure 13 compares 𝑘 𝜎 values of . , . , . , and . . Inthese examples, a grid resolution of Δ 𝑥 = / was used, alongwith boundary condition parameters ℎ = . and 𝑏 = . The fig-ure demonstrates that as surface tension increases, the molten waxspreads significantly less. As the wax cools and resolidifies, visuallyinteresting layering behavior is observed.Figure 14 shows an example of several candle letters melting in acontainer. Wicks follow generally curved paths inside the letters.Melt pools from the different letters seamlessly interact. This simu-lation used a surface tension coefficient 𝑘 𝜎 = . , ℎ = . , 𝑏 = ,dynamic viscosity of . , Δ 𝑥 = / , and a bulk modulus of . for liquid and solid phases. Fig. 13. Various 𝑘 𝜎 values ( . , . , . , . ) are simulated in the case ofa melting candle. Frame 1202 is shown.Fig. 14. Letter-shaped candles melt inside a container. (Top) Frame 1, beforeflames are lit. (Middle)
Frame 60, in the middle of melting. (Bottom)
Frame200, as flames are extinguished and wax pools resolidify.
We simulated a liquid metal droplet that moves under the Marangonieffect, i.e. due to a temperature-induced surface tension gradient.The droplet first falls to the ground and spreads on the dry surface.We then turn on the heating while the droplet is still spreading;only one side of the droplet is subjected to heating. At its originaltemperature, the surface tension coefficient 𝑘 𝜎 of the droplet is ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
Momentum-Conserving Implicit Material Point Method for Surface Energies with Spatial Gradients • 13
Fig. 15. A liquid metal droplet subjected to heating on one side. The surfacetension coefficient increases as the temperature increases. (Top) the liquidmetal at frame 45 and frame 130. (Bottom) the particle view of temperaturedistribution at frame 45 and frame 130. The red color indicates highertemperature.Table 1. Performance measurements for one time step of several of ourexamples, broken down by (1) generating surface and balance particles, (2)transferring MPM particle quantities to the background grid and implicitgrid update, and (3) transferring background grid quantities to MPM par-ticles, merging surface, balance and MPM particles, and updating MPMparticle states. Note that the merging time does not exceed . of (3). Alltimes are in milliseconds. Example → Grid Grid → Part.Droplet Impact ( 𝑘 𝜎 = ) 2M 794K 100K 2224 7705 2360Droplets on Ramps ( 𝑘 𝜎𝑆𝐿 / 𝑘 𝜎𝐿𝐺 = . ) 1.5M 70K 100K 258 1147 287Contact Angles ( 𝑘 𝜎𝑆𝐿 / 𝑘 𝜎𝐿𝐺 = ) 256K 230K 250K 492 3715 571Soap Droplet in Water 1M 4M 200K 2166 32099 3205Wine glass 2M 1.6M 500K 1549 11268 1172Candle ( 𝑘 𝜎 = . ) 2M 618K 50K 1420 27500 1662Candle Letters 256K 3.1M 100K 4601 170048 2739Droplet with Marangoni Effect 4.1M 235K 200K 29991 8812 3244 . . As the temperature increase, 𝑘 𝜎 increases linearly with thetemperature. When the change in temperature is greater than 𝐾 ,the surface tension coefficient reaches its maximum value of .Since the hotter side of the droplet has higher surface tension, thestronger surface tension drives the particles to flow to the colderside as shown in Figure 15. This surface tension gradient results inan interesting self-propelled behavior of the liquid metal droplet.A domain of size × × with a grid resolution of Δ 𝑥 = / is used for this example. The time step is restricted between − to × − seconds by the CFL condition. The bulk modulus ofthe liquid is . . ℎ = . and 𝑏 = are used for the thermalboundary conditions. Table 1 shows average per-timestep runtime details for several of ourexamples. For this table, all experiments were run on a workstationequipped with 128GB RAM and with dual Intel ® Xeon ® E5-2687Wv4 CPUs at 3.00Ghz.
Our method allows for simulation of surface tension energies withspatial gradients, including those driven by variation in temperature. Our MPM approach to the problem resolves many interesting char-acteristic phenomena associated with these variations. However,while we provide for perfect conservation of linear and angularmomentum, our approach to the thermal transfers is not perfectlyconservative. Developing a thermally conservative transfer strategyis interesting future work. Also, although we simulate tears of wineon the walls of a glass, we did not simulate the effect of alcoholevaporation on the surface energy variation. Adding in a mixturemodel as in [Ding et al. 2019] would be interesting future work.Lastly, although our approach was designed for MPM, SPH is morecommonly used for simulation of liquids. However, SPH and MPMhave many similarities as recently shown by the work of Gissler etal. [Gissler et al. 2020b] and it would be interesting future work togeneralize our approach to SPH.
REFERENCES
A. Adamson and A. Gast. 1967.
Physical chemistry of surfaces . Vol. 150. IntersciencePublishers New York.R. Ando, N. Thurey, and R. Tsuruno. 2012. Preserving Fluid Sheets with AdaptivelySampled Anisotropic Particles.
IEEE Trans Vis Comp Graph
18, 8 (Aug. 2012), 1202–1214.Annonymous. 2021.
Supplementary Technical Document . Technical Report.O. Azencot, O. Vantzos, M. Wardetzky, M. Rumpf, and M. Ben-Chen. 2015. Functionalthin films on surfaces. In
Proc 14th ACM SIGGRAPH/Eurograph Symp Comp Anim .137–146.T. Belytschko, W. Liu, B. Moran, and K. Elkhodary. 2013.
Nonlinear finite elements forcontinua and structures . John Wiley and sons.J. Brackbill, D. Kothe, and C. Zemach. 1992. A continuum method for modeling surfacetension.
J Comp Phys
Comp Meth App Mech Eng
Transactions of theFaraday society
40 (1944), 546–551.C. R. A. Chaitanya, A. S. Kaplanyan, C. Schied, M. Salvi, A. Lefohn, D. Nowrouzezahrai,and T. Aila. 2017. Interactive reconstruction of Monte Carlo image sequences usinga recurrent denoising autoencoder.
ACM Trans Graph
36, 4 (2017), 1–12.Y.-L. Chen, J. Meier, B. Solenthaler, and V.C. Azevedo. 2020. An Extended Cut-CellMethod for Sub-Grid Liquids Tracking with Surface Tension.
ACM Trans Graph
Marching cubes 33: Construction of topologically correct isosurfaces .Technical Report.P. Clausen, M. Wicke, J. R. Shewchuk, and J. F. O’brien. 2013. Simulating liquids andsolid-liquid interactions with Lagrangian meshes.
ACM Transactions on Graphics(TOG)
32, 2 (2013), 17.M. Corsini, P. Cignoni, and R. Scopigno. 2012. Efficient and Flexible Sampling withBlue Noise Properties of Triangular Meshes.
IEEE Trans Vis Comp Graph
18, 6 (2012),914–924. https://doi.org/10.1109/TVCG.2012.34F. Da, C. Batty, C. Wojtan, and E. Grinspun. 2015. Double bubbles sans toil and trouble:discrete circulation-preserving vortex sheets for soap films and foams.
ACM TransGraph (SIGGRAPH 2015) (2015).F. Da, D. Hahn, C. Batty, C. Wojtan, and E. Grinspun. 2016. Surface-only liquids.
ACMTrans Graph (TOG)
35, 4 (2016), 1–12.G. Daviet and F. Bertails-Descoubes. 2016. A Semi-implicit Material Point Method forthe Continuum Simulation of Granular Materials.
ACM Trans Graph
35, 4 (2016),102:1–102:13.C. C. de Langavant, A. Guittet, M. Theillard, F. Temprano-Coleto, and F. Gibou. 2017.Level-set simulations of soluble surfactant driven flows.
J Comp Phys
348 (2017),271–297.A. de Vaucorbeil, V. P. Nguyen, S. Sinaie, and J. Y. Wu. 2020. Chapter Two - Materialpoint method after 25 years: Theory, implementation, and applications. Advancesin Applied Mechanics, Vol. 53. Elsevier, 185 – 398. https://doi.org/10.1016/bs.aams.2019.11.001M. Ding, X. Han, S. Wang, T. Gast, and J. Teran. 2019. A thermomechanical materialpoint method for baking and cooking.
ACM Trans Graph
38, 6 (2019), 192.Y. Dukler, H. Ji, C. Falcon, and A. L Bertozzi. 2020. Theory for undercompressive shocksin tears of wine.
Phys Rev Fluids
5, 3 (2020), 034002.E. Edwards and R. Bridson. 2012. A high-order accurate Particle-In-Cell method.
Int JNumer Meth Eng
90 (2012), 1073–1088.Y. Fang, M. Li, Ming Gao, and Chenfanfu Jiang. 2019. Silly rubber: an implicit materialpoint method for simulating non-equilibrated viscoelastic and elastoplastic solids.ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
ACM Trans Graph
38, 4 (2019), 1–13.R. Farahi, A. Passian, T. Ferrell, and T. Thundat. 2004. Microfluidic manipulation viaMarangoni forces.
Applied Phys Let
85, 18 (2004), 4237–4239.Y. Fei, C. Batty, E. Grinspun, and C. Zheng. 2018. A multi-scale model for simulatingliquid-fabric interactions.
ACM Trans Graph
37, 4 (2018), 51:1–51:16. https://doi.org/10.1145/3197517.3201392Y. Fei, H. Maia, C. Batty, C. Zheng, and E. Grinspun. 2017. A multi-scale model forsimulating liquid-hair interactions.
ACM Trans. Graph.
36, 4 (2017), 56:1–56:17.https://doi.org/10.1145/3072959.3073630M. M. Francois, J. M. Sicilian, and D. B. Kothe. 2006. Modeling of thermocapillary forceswithin a volume tracking algorithm. In
Modeling of Casting, Welding and AdvancedSolidification Processes–XI (Opio, France). 935–942.C. Fu, Q. Guo, T. Gast, C. Jiang, and J. Teran. 2017. A Polynomial Particle-in-cell Method.
ACM Trans Graph
36, 6 (Nov. 2017), 222:1–222:12.M. Gao, A. Tampubolon, C. Jiang, and E. Sifakis. 2017b. An adaptive generalizedinterpolation material point method for simulating elastoplastic materials.
ACMTrans Graph
36, 6 (2017), 223:1–223:12. https://doi.org/10.1145/3130800.3130879Y. Gao, S. Li, L. Yang, H. Qin, and A. Hao. 2017a. An efficient heat-based model forsolid-liquid-gas phase transition and dynamic interaction.
Graphical Models
ACM Transactions onGraphics (TOG)
37, 3 (2018), 1–12.C. Gissler, A. Henne, S. Band, A. Peer, and M. Teschner. 2020a. An Implicit CompressibleSPH Solver for Snow Simulation.
ACM Trans Graph
39, 4, Article 36 (July 2020),16 pages. https://doi.org/10.1145/3386569.3392431C. Gissler, A. Henne, S. Band, A. Peer, and M. Teschner. 2020b. An implicit compressibleSPH solver for snow simulation.
ACM Trans Graph (TOG)
39, 4 (2020), 36–1.O. Gonzalez and A. Stuart. 2008.
A first course in continuum mechanics . CambridgeUniversity Press.Y. Gu and Y.-H. Yang. 2016. Physics Based Boiling Bubble Simulation. In
SIGGRAPHASIA 2016 Technical Briefs (Macau) (SA ’16) . Association for Computing Machinery,New York, NY, USA, Article 5, 4 pages. https://doi.org/10.1145/3005358.3005385Q. Guo, X. Han, C. Fu, T. Gast, R. Tamstorf, and J. Teran. 2018. A material pointmethod for thin shells with frictional contact.
ACM Trans Graph
37, 4 (2018), 147.https://doi.org/10.1145/3197517.3201346F. Harlow and E. Welch. 1965. Numerical calculation of time-dependent viscous incom-pressible flow of fluid with free surface.
Phys Fl
8, 12 (1965), 2182–2189.H. Hochstetter and A. Kolb. 2017. Evaporation and Condensation of SPH-Based Fluids.In
Proc ACM SIGGRAPH/Eurographics Symp Comp Anim (Los Angeles, California) (SCA ’17) . Association for Computing Machinery, New York, NY, USA, Article 3,9 pages. https://doi.org/10.1145/3099564.3099580M. Hopp-Hirschler, M. S. Shadloo, and U. Nieken. 2018. A Smoothed Particle Hydro-dynamics approach for thermo-capillary flows.
Comp Fluids
176 (2018), 1 – 19.https://doi.org/10.1016/j.compfluid.2018.09.010H. Hu and P. Eberhard. 2017. Thermomechanically coupled conduction mode laserwelding simulations using smoothed particle hydrodynamics.
Comp Part Mech
ACM Trans Graph
39, 4, Article 41 (July 2020), 14 pages. https://doi.org/10.1145/3386569.3392094D.A.B. Hyde, S.W. Gagniere, A. Marquez-Razon, and J. Teran. 2020. An Implicit UpdatedLagrangian Formulation for Liquids with Large Surface Energy.
ACM Trans Graph
39, 6, Article 183 (Nov. 2020), 13 pages. https://doi.org/10.1145/3414685.3417845G. Irving, J. Teran, and R. Fedkiw. 2004. Invertible Finite Elements for Robust Simulationof Large Deformation. In
Proc ACM SIGGRAPH/Eurograph Symp Comp Anim . 131–140.S. Ishida, P. Synak, F. Narita, T. Hachisuka, and C. Wojtan. 2020. A Model for Soap FilmDynamics with Evolving Thickness.
ACM Trans Graph
39, 4, Article 31 (July 2020),11 pages. https://doi.org/10.1145/3386569.3392405C. Jiang, T. Gast, and J. Teran. 2017. Anisotropic elastoplasticity for cloth, knit and hairfrictional contact.
ACM Trans Graph
36, 4 (2017), 152.C. Jiang, C. Schroeder, A. Selle, J. Teran, and A. Stomakhin. 2015. The Affine Particle-In-Cell Method.
ACM Trans Graph
34, 4 (2015), 51:1–51:10.C. Jiang, C. Schroeder, J. Teran, A. Stomakhin, and A. Selle. 2016. The Material PointMethod for Simulating Continuum Materials. In
ACM SIGGRAPH 2016 Course . 24:1–24:52.R. E. Johnson Jr. and R. H. Dettre. 1964. Contact angle hysteresis. III. Study of anidealized heterogeneous surface.
J Phys Chem
68, 7 (1964), 1744–1750.G. Klár, T. Gast, A. Pradhana, C. Fu, C. Schroeder, C. Jiang, and J. Teran. 2016. Drucker-prager Elastoplasticity for Sand Animation.
ACM Trans Graph
35, 4 (2016), 103:1–103:12.D. Langbein. 2002.
Capillary surfaces: shape – stability – dynamics, in particular underweightlessness . Vol. 178. Springer Science & Business Media. T. Lenaerts and P. Dutré. 2009. An architecture for unified SPH simulations.
CW Reports (2009).W. Li, D. Liu, M. Desbrun, J. Huang, and X. Liu. 2020. Kinetic-based Multiphase FlowSimulation.
IEEE Trans Vis Comp Graph (2020).F. Losasso, J. Talton, N. Kwatra, and R. Fedkiw. 2008. Two-Way Coupled SPH and ParticleLevel Set Fluid Simulation.
IEEE Trans Visu Comp Graph
14, 4 (2008), 797–804.T. Maeshima, Y. Kim, and T. I. Zohdi. 2020. Particle-scale numerical modeling of thermo-mechanical phenomena for additive manufacturing using the material point method.
Computational Particle Mechanics (2020). https://doi.org/10.1007/s40571-020-00358-xJ. Monaghan. 1992. Smoothed particle hydrodynamics.
Ann Rev Astron Astroph
30, 1(1992), 543–574.R. Narain, A. Golas, and M. Lin. 2010. Free-flowing granular materials with two-waysolid coupling.
ACM Trans Graph
29, 6 (2010), 173:1–173:10.S. Nas and G. Tryggvason. 2003. Thermocapillary interaction of two bubbles or drops.
Int J Multiphase Flow
29, 7 (2003), 1117–1135. https://doi.org/10.1016/S0301-9322(03)00084-3R. Osada, T. Funkhouser, B. Chazelle, and D. Dobkin. 2002. Shape Distributions.
ACMTrans. Graph.
21, 4 (Oct. 2002), 807–832. https://doi.org/10.1145/571647.571648A. Paiva, F. Petronetto, T. Lewiner, and G. Tavares. 2009. Particle-based viscoplasticfluid/solid simulation.
Computer-Aided Design
41, 4 (2009), 306–314.M. Pauly, R. Keiser, B. Adams, P. Dutré, M. Gross, and L. Guibas. 2005. Meshlessanimation of fracturing solids.
ACM Trans Graph
24, 3 (2005), 957–964. https://doi.org/10.1145/1073204.1073296S. Pirk, M. Jarząbek, T. Hädrich, D.L. Michels, and W. Palubicki. 2017. Interactive WoodCombustion for Botanical Tree Models. 36, 6, Article 197 (Nov. 2017), 12 pages.https://doi.org/10.1145/3130800.3130814D. Ram, T. Gast, C. Jiang, C. Schroeder, A. Stomakhin, J. Teran, and P. Kavehpour. 2015.A material point method for viscoelastic fluids, foams and sponges. In
Proc ACMSIGGRAPH/Eurograph Symp Comp Anim . 157–163.R. Rioboo, C. Tropea, and M. Marengo. 2001. Outcomes from a Drop Impact on SolidSurfaces.
Atomization and Sprays
11, 2 (2001).M. A. Russell. 2018.
A Smoothed Particle Hydrodynamics Model for the Simulation ofLaser Fusion Additive Manufacturing Processes . Ph.D. Dissertation. UC Berkeley.C. Schreck and C. Wojtan. 2020. A practical method for animating anisotropic elasto-plastic materials.
Computer Graphics Forum - Eurographics 2020
39, 2 (2020).L. Scriven and C. Sternling. 1960. The marangoni effects.
Nature
Proc Symp Comp Anim . 25–32.A. Stomakhin, C. Schroeder, L. Chai, J. Teran, and A. Selle. 2013. A Material PointMethod for snow simulation.
ACM Trans Graph
32, 4 (2013), 102:1–102:10.A. Stomakhin, C. Schroeder, C. Jiang, L. Chai, J. Teran, and A. Selle. 2014. AugmentedMPM for phase-change and varied materials.
ACM Trans Graph
33, 4 (2014), 138:1–138:11.D. Sulsky, Z. Chen, and H. Schreyer. 1994. A particle method for history-dependentmaterials.
Comp Meth App Mech Eng
SIAM J Sci Comp
31, 4 (2009), 2447–2471.A. Tartakovsky and P. Meakin. 2005. Modeling of surface tension and contact angleswith smoothed particle hydrodynamics.
Phys Rev E
72 (Aug 2005), 026301. Issue 2.https://doi.org/10.1103/PhysRevE.72.026301D. Terzopoulos, J. Platt, and K. Fleischer. 1991. Heating and melting deformable models.
The Journal of Visualization and Computer Animation
2, 2 (1991), 68–73.J. Thomson. 1855. XLII. On certain curious motions observable at the surfaces of wineand other alcoholic liquors.
London, Edinburgh, and Dublin Phil Mag J Sci
10, 67(1855), 330–333.N. Thürey, C. Wojtan, M. Gross, and G. Turk. 2010. A multiscale approach to mesh-basedsurface tension flows.
ACM Trans Graph (TOG)
29, 4 (2010), 1–10.M. Tong and D. J. Browne. 2014. An incompressible multi-phase smoothed particlehydrodynamics (SPH) method for modelling thermocapillary flow.
Int J Heat MassTransfer
73 (2014), 284 – 292. https://doi.org/10.1016/j.ijheatmasstransfer.2014.01.064D. C. Venerus and D. N. Simavilla. 2015. Tears of wine: New insights on an oldphenomenon.
Scientific reports
Comp Anim Virtual Worlds
23, 3-4 (2012), 279–289. https://doi.org/10.1002/cav.1457 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/cav.1457H. Wang, Y. Jin, A. Luo, X. Yang, and B. Zhu. 2020a. Codimensional Surface TensionFlow Using Moving-Least-Squares Particles.
ACM Trans Graph
39, 4, Article 42 (July2020), 16 pages. https://doi.org/10.1145/3386569.3392487H. Wang, G. Miller, and G. Turk. 2007. Solving General Shallow Wave Equations onSurfaces (SCA ’07) . Eurographics Association, Goslar, DEU, 229–238.H. Wang, P. J. Mucha, and G. Turk. 2005. Water Drops on Surfaces.
ACM Tran. Graph
24, 3 (July 2005), 921–929. https://doi.org/10.1145/1073204.1073284ACM Trans. Graph., Vol. 1, No. 1, Article . Publication date: February 2021.
Momentum-Conserving Implicit Material Point Method for Surface Energies with Spatial Gradients • 15
S. Wang, M. Ding, T. Gast, L. Zhu, S. Gagniere, C. Jiang, and J. Teran. 2019. Simulationand Visualization of Ductile Fracture with the Material Point Method.
Proceedingsof the ACM on Computer Graphics and Interactive Techniques
2, 2, 18.X. Wang, M. Li, Y. Fang, X. Zhang, M. Gao, M. Tang, D. Kaufman, and C. Jiang. 2020b.Hierarchical optimization time integration for CFL-rate MPM stepping.
ACM TransGraph (TOG)
39, 3 (2020), 1–16.X. Wang, Y. Qiu, S.R. Slattery, Y. Fang, M. Li, S.-C. Zhu, Y. Zhu, M. Tang, D. Manocha,and C. Jiang. 2020c. A Massively Parallel and Scalable Multi-GPU Material PointMethod.
ACM Trans Graph
39, 4, Article 30 (July 2020), 15 pages. https://doi.org/10.1145/3386569.3392442C. Wojtan, N. Thürey, M. Gross, and G. Turk. 2010. Physics-inspired topology changesfor thin fluid features.
ACM Trans Graph
29, 4 (2010), 50:1–50:8. https://doi.org/10.1145/1778765.1778787J. Wolper, Y. Chen, M. Li, Y. Fang, Z. Qu, J. Lu, M. Cheng, and C. Jiang. 2020. AnisoMPM:animating anisotropic damage mechanics.
ACM Trans. Graph.
39, 4, Article 37(2020).J. Wolper, Y. Fang, M. Li, J. Lu, M. Gao, and C. Jiang. 2019. CD-MPM: continuum damagematerial point methods for dynamic fracture animation.
ACM Trans. Graph.
38, 4,Article 119 (2019).S. Yang, X. He, H. Wang, S. Li, G. Wang, E. Wu, and K. Zhou. 2016a. Enriching SPHsimulation by approximate capillary waves. In
Symp Comp Anim . 29–36. T. Yang, J. Chang, M. C. Lin, R. R. Martin, J. J. Zhang, and S.-M. Hu. 2017. A unifiedparticle system framework for multi-phase, multi-material visual simulations.
ACMTransactions on Graphics (TOG)
36, 6 (2017), 224.T. Yang, M. C. Lin, R. R. Martin, J. Chang, and S.-M. Hu. 2016b. Versatile Interactions atInterfaces for SPH-Based Simulations (SCA ’16) . Eurographics Association, Goslar,DEU, 57–66.T. Young. 1805. III. An essay on the cohesion of fluids.
Phil TransRoyal Soc London
95 (1805), 65–87. https://doi.org/10.1098/rstl.1805.0005arXiv:https://royalsocietypublishing.org/doi/pdf/10.1098/rstl.1805.0005Y. Yue, B. Smith, C. Batty, C. Zheng, and E. Grinspun. 2015. Continuum foam: amaterial point method for shear-dependent flows.
ACM Trans Graph
34, 5 (2015),160:1–160:20.T. Zhang, J. Shi, C. Wang, H. Qin, and C. Li. 2017. Robust Gas Condensation Simulationwith SPH based on Heat Transfer. In
Pacific Graphics Short Papers , Jernej Barbic,Wen-Chieh Lin, and Olga Sorkine-Hornung (Eds.). The Eurographics Association,27–32. https://doi.org/10.2312/pg.20171321W. Zheng, B. Zhu, B. Kim, and R. Fedkiw. 2015. A new incompressibility discretizationfor a hybrid particle MAC grid representation with surface tension.
J Comp Phys
280 (2015), 96–142.B. Zhu, E. Quigley, M. Cong, J. Solomon, and R. Fedkiw. 2014. Codimensional surfacetension flow on simplicial complexes.