$A$-ULMPM: An Arbitrary Updated Lagrangian Material Point Method for Efficient Simulation of Solids and Fluids

The MPM discretization has been widely adopted in computer graphics applications (Jiang et al., 2016). The seminal work of Zhu and Bridson (Zhu and Bridson, 2005) first introduced the FLIP method for sand simulation. Subsequent works further explored its strength in simulating a broader spectrum of material behaviors including snow (Stomakhin et al., 2013), granular materials (Daviet and Bertails-Descoubes, 2016; Klár et al., 2016; Tampubolon et al., 2017; Gao et al., 2018a), foam (Ram et al., 2015; Yue et al., 2015), complex fluids (Fang et al., 2019; Gao et al., 2017), cloth, hair and fiber collisions (Jiang et al., 2017a; Fei et al., 2017, 2018), fracture (Wolper et al., 2019; Wolper et al., 2020) and phase change (Stomakhin et al., 2014; Gao et al., 2018b; Su et al., 2021). We also note the related works of (McAdams et al., 2009) for hair simulation, (Sifakis et al., 2008) for cloth simulation, (Narain et al., 2010) for sand simulation and (Patkar et al., 2013) for bubble simulation, which bear similarities to MPM due to their hybrid nature. Various works have improved or modified aspects of the standard MPM techniques commonly used in graphics (Fang et al., 2020; Yue et al., 2018; Xue et al., 2020; Ding et al., 2019). Among them, notably, Jiang et al. (2015; 2017b) proposed an Affine Particle-In-Cell (APIC) approach that conserves angular momentum and prevents visual artifacts such as noise, instability, clumping and volume loss/gain existing in both FLIP and PIC methods. Furthermore, APIC was enhanced in (Fu et al., 2017; Hu et al., 2018) to improve the kinetic energy conservation in particle/grid transfers.

In the engineering community, MPM was first introduced in (Sulsky et al., 1994) as an extension of the FLIP method (Brackbill et al., 1988), and substantial improvements and variants have been proposed thereafter, including experimental validation for studying dynamic anticrack propagation in snow avalanches (Gaume et al., 2018) and the use of MPM for designing differentiable physics engines for robotics applications (Hu et al., 2019). Different strategies for updating the stress were compared to investigate the energy conservation error in MPM in (Bardenhagen, 2002). The quadrature error and cell-crossing error of MPM was investigated in (Steffen et al., 2008). A generalized interpolation material point method (GIMPM) (Bardenhagen and Kober, 2004) was proposed to obtain a smoother field representation by combining the shape functions of the grid with the particle characteristic function. The cell-crossing error in the MPM discretization was alleviated by introducing the local tangent affine deformation of particles in the convected particle domain method (CPDI) (Sadeghirad et al., 2011). However, CPDI does not completely remove the numerical fracture issue due to gaps between particles and grid cells. Recently, the standard MPM was reformulated with respect to the initial topology, which provides the total Lagrangian formulation for MPM (TLMPM) (de Vaucorbeil et al., 2020), which has been proven to completely eliminate the numerical fracture issue and cell-crossing instability and has been further extended in (de Vaucorbeil and Nguyen, 2021) to support multi-body contacts.

The Lagrangian $\mathcal{L}$ for holonomic systems is defined as:

where $\boldsymbol{q}$ is the generalized displacement, $\boldsymbol{\dot{q}}$ is the generalized velocity, $\mathcal{K}(\dot{\boldsymbol{q}})$ is the kinetic energy and $\mathcal{U}(\boldsymbol{q})$ is the potential energy. By omitting the energy due to external body forces and traction for simplicity, the kinetic and potential energies can be defined as:

where $\rho_{0}$ is the material density at time $t_{0}$, $\Psi(\boldsymbol{F})$ denotes the Helmholtz free energy per unit mass in homogeneous materials, and $\boldsymbol{F}$ is the deformation gradient tensor. Consequently, the Lagrangian density function can be defined as follows:

Based on the Lagrangian framework (Zienkiewicz et al., 1977), the governing equations of motion at time $t_{n}$ can be described as:

where

Note that the term ${\partial\Psi(\boldsymbol{F})}/{\partial\boldsymbol{F}}$ represents a stress tensor that is determined by the material constitutive model, and the term ${\partial\boldsymbol{F}}/{\partial\boldsymbol{q}}$ can be further expressed in terms of a divergence operator. These two terms together define the internal force as the material deforms.

In this formulation, the stress and strain in the material are measured relative to the original configuration at time $t_{0}$ (see Figure 5), such that the deformation gradient tensor $\boldsymbol{F}$ is the displacement derivative of $\boldsymbol{q}_{n}$ with respect to $\boldsymbol{q}_{0}$. Substituting $\boldsymbol{F}_{0n}$ to $\partial\bar{\mathcal{L}}/\partial\boldsymbol{q}_{n}$ gives:

and have the following governing equation of motion:

where the divergence operator $\nabla_{0}$ is also evaluated with respect to the original configuration at time $t_{0}$. In the total Lagrangian formulation, the stress tensor $\mathbb{P}_{0}$ is the first Piola–Kirchhoff stress.

In this formulation, the stress and strain in the material are measured relative to the current configuration at time $t_{n}$ (see Figure 5). Thus, the expression for ${\partial\bar{\mathcal{L}}}/{\partial\boldsymbol{q}}$ can be expanded as follows:

Besides, $\rho_{0}$ can be mapped to $\rho_{n}$ using the relation $\rho_{0}=J_{0n}\rho_{n}$, where $\rho_{n}$ is the density at time $t_{n}$ and $J_{0n}=\text{det}(\boldsymbol{F}_{0n})$. Consequently, the Eulerian formulation gives the following equation of motion:

where $\mathbb{P}_{n}=\mathbb{P}_{0}\boldsymbol{F}_{0n}^{T}/{J_{0n}}$ provides the definition of the Cauchy stress.

A general formulation for measuring stress and strain with respect to an arbitrary reference configuration at time $t_{s}$ has been derived in (Zienkiewicz et al., 1977; Shabana, 2018). In this formulation, the expression for ${\partial\bar{\mathcal{L}}}/{\partial\boldsymbol{q}}$ can be expanded as follows:

Similar to the Eulerian formulation, we map $\rho_{0}$ to $\rho_{s}$ using the relation $\rho_{s}=J_{0s}\rho_{s}$, where $\rho_{s}$ represents the density at time $t_{s}$ and $J_{0s}=\text{det}(\boldsymbol{F}_{0s})$. Consequently, the arbitrary updated Lagrangian formulation gives the following governing equation of motion:

where $\mathbb{P}_{s}=\mathbb{P}_{0}\boldsymbol{F}_{0s}^{T}/{J_{0s}}$ defines a stress measured at time $t_{s}$. By setting $s=0$ and using the defining properties of the initial configuration $J_{00}=1$ and $\boldsymbol{F}_{00}=\boldsymbol{I}$, where $\boldsymbol{I}$ is the identity matrix, equation (9) recovers the total Lagrangian formulation in equation (6). Likewise, setting $s=n$ yields the Eulerian formulation in equation (7).

We show integration of our arbitrary updated Lagrangian formulation with Kernel-MPM (Stomakhin et al., 2013) in Appendix B and MLS-MPM (Hu et al., 2018) in Section 5. In Kernel-MPM, the interpolation is performed using shape functions and stress divergence is evaluated by derivatives of the shape function, while MLS-MPM utilizes APIC particle-to-grid velocity rasterization and uses MLS shape functions to derive the internal force term. Using these definitions, our terminology for different methods is provided as follows:

1. (1)

   Kernel-MPM: In kernel-MPM, we have kernel-EMPM and kernel-TLMPM represent the kernel-MPM in Eulerian and total Lagrangian frameworks, respectively. Kernel-EMPM was presented in (Stomakhin et al., 2013) and kernel-TLMPM was present in (de Vaucorbeil et al., 2020). These two methods can be recovered in our kernel-$A$-ULMPM framework by setting $s=n$ for EMPM and $s=0$ for TLMPM.

2. (2)

   MLS-MPM: In MLS-MPM, we have MLS-EMPM and MLS-TLMPM represent MLS-MPM in Eulerian and total Lagrangian frameworks, respectively. MLS-EMPM was presented in MLS-MPM (Hu et al., 2018). Besides, MLS-EMPM and MLS-TLMPM can be recovered in MLS-$A$-ULMPM by setting $s=n$ for MLS-EMPM and $s=0$ for MLS-TLMPM.

We use a global grid that covers the entire computational domain. Each object has its own configuration map $\phi$ that maps points $\boldsymbol{x}$ within the object to specific locations $\phi(\boldsymbol{x})$ in the global grid. To reduce the associated memory overhead, we implement the global grid using sparsity aware data structures (Setaluri et al., 2014).

In contrast to EMPM, the velocity gradients $\nabla_{s}\boldsymbol{v}_{p}$ and weights $W_{pi}^{s}$ are all evaluated with respect to the configuration map $\phi$ at time $t_{s}$. We use a first order Taylor expansion to reconstruct particle velocities and rasterize particle momentum to the grid as follows:

where $\boldsymbol{r}_{ip}^{s}=\left(\boldsymbol{q}_{p}^{s}-\boldsymbol{q}_{i}^{s}\right)$ and $\nabla_{s}\boldsymbol{v}_{p}^{n}=\frac{\partial\boldsymbol{v}_{p}^{n}}{\partial\boldsymbol{q}_{s}}$, which is further expanded out in equation (31). The mass/velocity at grid nodes are given by:

We also rasterize particle positions to the grid to simplify the theoretical derivations for the deformation gradient in equation (29) and internal forces in equation (19) as shown below:

The idea of introducing a first-order Taylor expansion to enhance velocity rasterization is not new and can also be found in APIC (Jiang et al., 2015) and MLS-MPM (Hu et al., 2018).

With explicit time stepping, the velocity is updated as $\boldsymbol{v}_{i}^{n+1}=\boldsymbol{\hat{v}}_{i}^{n+1}$, where $\boldsymbol{\hat{v}}_{i}^{n+1}$ is given by:

For a semi-implicit update, we follow (Stomakhin et al., 2013; Hu et al., 2018) and take an implicit step on the velocity update by utilizing the Hessian of $\mathcal{L}^{n+1}$ with respective to $\boldsymbol{q}_{i}^{n+1}$. The action of this Hessian on an arbitrary increment $\delta\boldsymbol{q}^{s}$ is given as follows:

where $\mathbb{A}_{p}^{s}$ is given by:

We linearize the implicit system with one step of Newton’s method, which provides the following symmetric system for $\boldsymbol{\hat{v}}_{i}^{n+1}$:

where $\boldsymbol{I}$ is the identity matrix and $\boldsymbol{\hat{v}}_{i}^{n+1}$ is given in equation (23). We update rasterized positions at grid nodes as shown below:

We project $\boldsymbol{q}_{i}^{n+1}$ to the global grid and process grid-based collisions (Stomakhin et al., 2013) to compute $\boldsymbol{v}_{i}^{n+1}$, which is transferred back to particles using the APIC method (Jiang et al., 2015). The specific updates for $\boldsymbol{q}_{p}^{n+1}$ and $\boldsymbol{v}_{p}^{n+1}$ are given below:

We first update particle deformation gradients with respect to the latest configuration at time $t_{s}$ and use the following MLS-based gradient operator (see equation (17)):

Substituting $\boldsymbol{q}_{i}^{n+1}=\boldsymbol{q}_{i}^{n}+\Delta t\boldsymbol{v}_{i}^{n+1}$ and $\boldsymbol{q}_{p}^{n+1}=\boldsymbol{q}_{i}^{n}+\Delta t\boldsymbol{v}_{p}^{n+1}$ to equation (29) gives:

where $\nabla_{s}\boldsymbol{v}_{p}^{n+1}$ is given by:

Using the chain rule, the update for $\boldsymbol{F}_{p}^{0(n+1)}$ is given by:

We designed the configuration update criterion based on an intuitive assumption that more severe deformation is accompanied with large change of $J$ in the next time step $t^{n+1}$ with respect to the latest configuration $t^{s}$. Using this observation, our criterion is defined as a measure for the amount of particle deformation as follows:

where $J_{p}^{s{n+1}}=\text{det}(\boldsymbol{F}_{p}^{s(n+1)})$ and $J_{p}^{ss}=1$. We mark particles with $\delta J_{p}\geq\epsilon$ as indicators where large deformation is taking place and count the total amount of marked particles ($n_{mp}$). If $n_{mp}/n_{p}\geq\eta$, where $n_{p}$ is the total number of particles, we update $W_{ip}^{s}$, $\mathbb{K}_{p}^{s}$, and $\boldsymbol{F}_{p}^{0s}$. Otherwise, these variables remain the same until a new update occurs. $\epsilon$ and $\eta$ are user-defined parameters to adjust the update frequency.

While our $A$-ULMPM framework is new to computer graphics, setting $s=n$ for the configuration map shares some similarities with MLS-EMPM (Hu et al., 2018) and APIC (Jiang et al., 2015). We prove that the MLS-gradient of velocity is equivalent to the APIC-gradient of velocity by leveraging the properties of interpolation weights $\sum_{i}\boldsymbol{r}_{pi}W_{pi}=\boldsymbol{0}$ (Jiang et al., 2017b) as follows:

Our internal force in equation (23) has the same pattern as that in original MLS-EMPM (Hu et al., 2018), which is derived from the weak-form element-free Galerkin (EFG) framework. The $\mathbb{K}_{p}$ matrices are simplified by $\frac{4}{\Delta x^{2}}\boldsymbol{I}$ for quadratic $W_{ip}$ and $\frac{3}{\Delta x^{2}}\boldsymbol{I}$ for cubic $W_{ip}$. This simplification comes from the properties of splines, and is not generally true for other interpolation functions.

We imposed torsion and stretch boundary conditions at the two ends of a hyperelastic beam to showcase that $A$-ULMPM allows large deformations in solids without non-physical fractures. As shown in Figure 4, traditional EMPM fails to preserve the shape of the beam during the twisting and pulling, while $A$-ULMPM can readily handle the challenging invertible elasticity when one end of the beam is released after twisting. Figure 7 shows that $A$-ULMPM is capable of robustly recovering the beam shape after extreme elastic deformations, while particles in EMPM cluster into one irreversible (or plastic) thin string blocking the “recovery” of elasticity.

Next, we tossed a stretchy Stanford bunny yo-yo with gravity forces. Our $A$-ULMPM hyperelastic bunny demonstrates rich elastic responses and realistic bouncing dynamics, while EMPM breaks due to severe numerical fracture (see Figure 10). Our $A$-ULMPM framework can perfectly handle hyperelastic deformation under severe bending (see Figure 13 (left)) and stretching (see Figure 13 (right)).

We simulated the fluid-like behavior of different materials using our $A$-ULMPM framework, as described in Section 5, such as elastoplastic snow (Stomakhin et al., 2013) and weakly compressible water. We dropped two copies of the snow Stanford bunny with different orientations to a solid wedge, as shown in Figure 8. $A$-ULMPM captures the vivid snow smashing and scattering after the bunnies fall on the wedge, similar to its Eulerian counterparts proposed in prior works (Stomakhin et al., 2013) (see side-by-side comparisons in our video). We simulated a water Stanford bunny falling inside a spherical container (see Figure 9), showing the extreme large deformations and energetic splashes. $A$-ULMPM does not update the configuration maps at each time step, in contrast to EMPM, so it is naturally more efficient in fluid-like simulations (see Table 3).

Our $A$-ULMPM framework can also simulate realistic fluid-structure interaction problems, as shown in Figure 11.

We dropped a water Stanford bunny to a hyperelastic bowl. $A$-ULMPM captures rich fluid-structure interactions, displaying advantages of both TLMPM and EMPM (see Figure 14 (right)), and highlighting that $A$-ULMPM can adaptively handle both fluid and solid simulations without having to switch between TLMPM and EMPM. In contrast, EMPM suffers from numerical fractures. As shown in Figure 11, the hyperelastic bowl fractures, causing the water to spill below.

Our model has generated a large number of compelling examples, but there remains much work to be done. Parameters to adjust the configuration update frequency were tuned by hand, and it would be interesting to calibrate them to measured models. Since each object has it own configuration map, we only briefly investigated contact by projecting particle positions to a global grid and processing grid-based collisions (Stomakhin et al., 2013). By doing this, we observed slight self-penetration in the 3D twisting beam example (see Figure 4). This could be addressed by further introducing contact algorithms, such as (Jiang et al., 2017a; Han et al., 2019). While we did not experience a need for excessively small time steps given the low stiffness of the materials we considered, deforming materials with high wave speed, such as steel, could benefit from a fully implicit discretization. Though side-to-side comparisons with Eulerian MLS-MPM shows that our $A$-ULMPM framework produces similar fluid-like dynamics as EMPM and captures rich interactions, a slight energy loss in $A$-ULMPM framework was observed. This could be addressed by introducing more accurate configuration update criteria. Finally, while our focus was on the material responses of water, snow, and hyperelastic solids, it would be interesting to investigate other materials and multi-physics coupling problems with $A$-ULMPM.

We proposed an arbitrary updated Lagrangian Material Point Method ($A$-ULMPM), which combines advantages of Total Lagrangian frameworks (de Vaucorbeil and Nguyen, 2021; de Vaucorbeil et al., 2020) and Eulerian frameworks (Stomakhin et al., 2013; Hu et al., 2018) in an adaptive fashion. $A$-ULMPM avoids the cell-crossing instability and numerical fracture in solid simulations, while still allowing for large deformations that arise in fluid-like simulations. It can be easily integrated with any existing MPM framework and builds a foundation for devising various MPM schemes, such as PolyPIC (Fu et al., 2017), for enhanced accuracy and visual vividness.