BFEMP: Interpenetration-Free MPM-FEM Coupling with Barrier Contact

Given trial function $\mathbf{Q}(\cdot,t):\Omega^{0}\rightarrow\mathbb{R}^{3}$, the corresponding Lagrangian weak form of Equation (4) is

where $T_{\alpha}=P_{\alpha\beta}N_{\beta}$ (with $\mathbf{N}(\mathbf{X})$ being the material space normal) is the traction field at the domain boundary $\partial\Omega^{0}$, on which one could presribe traction boundary conditions as needed.

While Total Langrangian FEM typically discretizes Equation (8), MPM usually adopts the Updated Lagrangian view and consequently discretizes an Eulerian weak form instead [8]. Correspondingly the stress derivatives are discretized at the current configuration $\Omega^{t}$. We can either push forward Equation (8) or directly integrate Equation (2) to reach

where $\mathbf{q}(\mathbf{x},t)=\mathbf{Q}(\Phi^{-1}(\mathbf{x},t),t)$ is the push forward of $\mathbf{Q}$, i.e., an Eulerian trial function, $\mathbf{a}(\mathbf{x},t)=\mathbf{A}(\Phi^{-1}(\mathbf{x},t),t)=\frac{D\mathbf{v}}{Dt}(\mathbf{x},t)=\frac{\partial\mathbf{v}}{\partial t}+\mathbf{v}\cdot\nabla\mathbf{v}$, and $\mathbf{t}=\bm{\sigma}\mathbf{n}$ is the traction at $\partial\Omega^{t}$ with $\mathbf{n}(\mathbf{x})$ being the outward pointing normal.

To enable the development of efficient optimization-based time integrators, we follow the variational treatment of the time-discretized incremental problem [53, 67]. Concretely, for a broad family of time discretization schemes (we will focus on backward Euler and the Newmark-$\beta$ family in this paper), the solution of Equation $(\ref{eqn:lagrangian-weak-form})$ during an incremental time step, i.e., the advancement from $\phi^{n}=\phi(\mathbf{X},t^{n})$ to $\phi^{n+1}=\phi(\mathbf{X},t^{n+1})$), for an hyperelastic material is given by minimizing the following functional:

where

encodes the inertia term and is a constant field in the minimization problem. The velocity update rule is

Here we have slightly modified the energy proposed by Radovitzky and Ortiz [67] to let it be compatible with backward Euler since their version was explicitly built for Newmark’s algorithm.

Note that in the case of backward Euler ($\alpha=1$, $\beta=\frac{1}{2}$, $\gamma=1$), it can be easily shown that the Euler-Lagrangian equation of functional (10) gives back the time-discretized momentum balance

Similarly, for Middle-point Newmark ($\alpha=\frac{1}{2}$, $\beta=\frac{1}{4}$, $\gamma=\frac{1}{2}$), we would get

Both of them match the results from temporally discretizing the Lagrangian form of Equation (2).

In the proposed framework, the Lagrangian form is discretized on the FEM domain $\Omega_{F}$ (Section 3.1), while the Eulerian form is discretized on the MPM domain $\Omega_{M}$ (Section 3.2). We assume

and their corresponding deformation map $\mathbf{x}_{M}=\phi_{M}(\mathbf{X}_{M},t)$ and $\mathbf{x}_{F}=\phi_{F}(\mathbf{X}_{F},t)$ are fully governed by their own variational forms in the absence of any coupling mechanism. In other words, without coupling, evolving the two domains independently (with two solvers) is equivalent to minimizing

by directly letting $\Omega=\Omega_{F}\cup\Omega_{M}$.

The coupling between FEM and MPM domains are then modeled via imposing the non-interpenetration constraints:

between the two domains. Note that the time is discretized with $t=t^{0},t^{1},t^{2},\dots,t^{n}$. We use $\Omega^{n}$ to denote $\Omega^{t^{n}}$ for notational simplicity. See Figure 1 for an illustration. Thus the final variational MPM-FEM coupling problem can be described as minimizing Equation (16) under the equality constraint described by Equation (17).

The feasible region described by Equation (17) can be equivalently expressed as

Moreover, if we define

to describe the Euclidean proximity between $\Omega_{M}^{t}$ and $\Omega_{F}^{t}$, Equation (18) can be further converted to a strict inequality constraint

Note that in Equation (15), we have assumed that the undeformed domains are non-overlapping. Thus the minimization problem starts with a strictly feasible solution at $t=0$.

In Section 4.1 we describe a barrier method that results in a contact pressure for enforcing Equation (15). See Section 4.3 for extra components on incorporating tangential frictional effects.

For the FEM domain, nodal positions and velocities are stored on mesh vertices and updated directly. The nodal masses are kept constant, i.e., $m^{n}_{i}=m_{i}$.

Inside any simplex element, the material and world space coordinate of an arbitrary location $\mathbf{X}$ are expressed using linear interpolation kernels $N_{i}(\mathbf{X})$ of node $\mathbf{X}_{i}$ as

Then according to the definition of $\mathbf{F}$ (Equation (3)), with $N_{i}$ being the linear hat function, the deformation gradient is piecewise constant. Inside a linear simplex element $e$ it is directly evaluated as a function of $x$:

where $\mathbf{T}_{e}(x)$ is the current triangle basis of element $e$ and $\mathbf{B}_{e}$ is the triangle basis of element $e$ in material space.

For the MPM domain, the nodal positions $x^{n}_{\text{M}}$ are the uniform Cartesian grid coordinates at each time step. The grid velocity $\mathbf{v}_{i}^{n}$ and mass $m_{i}^{n}$ are transferred from particles. The nodal movements are conceptual. $\tilde{x}^{n+1}_{\text{M}}$ and $\tilde{v}^{n+1}_{\text{M}}$ will be transferred back to particles for advection.

Similar to FEM, each MPM grid node $i$ is associated with a kernel function $N_{i}(\mathbf{x})$ for the grid to represent the continuous field. Note that the kernel is defined in terms of $\mathbf{x}$ rather than $\mathbf{X}$ because the grid is essentially a discretization of $\Omega_{M}^{n}$ – a direct consequence of adopting Updated Lagrangian kinematics. When $N_{i}$ is evaluated at a particle location $\mathbf{x}_{q}^{n}$, a shorter notation $N_{i}(\mathbf{x}_{q}^{n})=w_{iq}^{n}$ is from now on used instead. Here $N_{i}$ directly takes the current particle location $\mathbf{x}_{q}^{n}$ as input as opposed to FEM because there is no globally defined reference configuration in MPM and the deformation is evolved over time steps rather than recomputed using a rest shape. More specifically, the deformation gradient of a particle $q$ is defined as

In this paper, without loss of generality, we adopt the quadratic B-spline kernel for $N_{i}(\mathbf{x})$ to avoid MPM’s cell-crossing instability [69]. Other kernels based on the NURBS [70], Generalized Interpolation Material Point Method (GIMP) [9], Convective Particle Domain Interpolation (CPDI) [31, 71, 72], or the Dual Domain Material Point (DDMP) [73, 74] can also be directly used in our framework.

To transfer information between the particles and the grid, we implemented options including the Affine Particle-In-Cell (APIC) method [75, 34] (Table 1), Particle-In-Cell (PIC) method [3] (Table 2), and the Fluid-Implicit Particle (FLIP) method [4] (Table 3). Note that other transfer schemes such as XPIC [76] can also be applied in our framework in a straightforward manner.

Recently for Lagrangian FEM, Li et al.[57, 59] proposed a consistent variational contact model that smoothly approximates the nonsmooth contact phenomena with bounded error, and demonstrated its convergence under refinement for piecewise linear boundary discretization. Here we customize the FEM contact potential [59] to the FEM-MPM coupling setting by defining the inter-surface contact potential between surfaces $\partial\Omega_{\text{M}}$ and $\partial\Omega_{\text{F}}$ to be

where $d^{PP}(\mathbf{x}_{m},\mathbf{x}_{f})=\|\mathbf{x}_{m}-\mathbf{x}_{f}\|$ is the point-point distance function, and

is a smoothly clamped barrier function that serves as the contact energy density with $d$ the input distance, $\hat{d}$ a small distance threshold below which contact activates, and $\kappa$ in $Pa$ the barrier stiffness [57, 59], which scales the magnitude of contact forces at a certain distance.

Intuitively, $\min_{\mathbf{x}_{f}\in\partial\Omega_{\text{F}}^{t}}d^{PP}(\mathbf{x}_{m},\mathbf{x}_{f})$ is the distance between a material point $\mathbf{x}_{m}\in\partial\Omega_{\text{M}}^{t}$ and surface $\partial\Omega_{\text{F}}^{t}$, and Equation (26) can be viewed as an integration of the point($\mathbf{x}_{m}$)-surface($\partial\Omega_{\text{F}}^{t}$) contact energy density over surface $\partial\Omega_{\text{M}}^{0}$. The barrier function $b$ smoothly increases from $0$ to infinity as the input distance decreases from $\hat{d}$ to $0$, providing arbitrarily large repulsion to ensure no interpenetration and at the same time bound the contact gap error within $\hat{d}$ (Figure 2). As $\hat{d}\rightarrow 0$, the approximation error between the contact energy density function $b$ and the real contact phenomenon described in Equation (20) decreases, which also makes $\partial\Omega_{\text{M}}$ and $\partial\Omega_{\text{F}}$ interchangeable in the limit. Note that the integration is performed in the material space while the distance is evaluated in the world space.

After applying FEM and MPM discretization schemes, assuming 2D, the contact potential Equation (26) is discretized to be

where $\mathcal{Q}$ is the set of all MPM particles, $\omega_{q}$ is the integration weight (equivalently, the boundary area) of MPM particle $q$, $\mathcal{B}$ is the set of FEM boundary edges, and $d^{PE}(\mathbf{x}_{q},e)$ is the point-edge distance between particle $\mathbf{x}_{q}$ and edge $e$. Note that the MPM grid nodal positions $\tilde{x}_{\text{M}}$ to be solved and the particle positions $\mathbf{x}_{q}$ after advection using $\tilde{x}_{\text{M}}$ are linearly related through the particle-grid transfer kernel (Figure 3 right), and so the contact force on the MPM nodal degrees of freedom can be calculated by applying the chain rule:

where

with $\alpha,\beta=1,2,...,d$ and $d=2$ or $3$ the spatial dimension since $\mathbf{x}_{q}=\sum_{i}w_{iq}^{n}\mathbf{x}_{i}$.

Ideally, $\omega_{q}$ should be zero for interior particles. It should reveal the proportional surface area of boundary particles. Since FEM boundary elements are always outside the MPM domain and the barrier function $b$ is only activated at a small distance, the activation of $b$ can be applied to conveniently decide whether an MPM particle is at the boundary or the interior without explicitly identifying the MPM domain boundary in each time step (Figure 3 left). Then assuming a close to uniform particle distribution to be maintained throughout the simulation, $\omega_{q}$ can be set to $2\sqrt{V_{q}^{0}/\pi}$ in 2D and $\pi\left(3V_{q}^{0}/\left(4\pi\right)\right)^{\frac{2}{3}}$ in 3D for all particles, which is the area of the largest cross section of a spherical particle $q$.

The minimization operator in the potential (Equation (28)) helps to compute the point-polyline distance from point-edge distances. As the barrier function $b$ is monotonically decreasing, the potential can be rewritten as

Due to the existence of a $\max$ operator, it is only $C^{0}$ continuous, making the incremental potential challenging to be efficiently minimized by gradient-based optimization methods like Newton’s method. Since the barrier function $b$ is with local support around each boundary element, and that it maps a majority of the activate distances to tiny potential values (Figure 4), the maximization of the potential field can be well approximated by summation, with the duplicate potential around FEM boundary nodes compensated by subtraction as proposed by Li et al. [59]:

where $\hat{\mathcal{B}}$ is the set of all FEM boundary nodes and $\eta_{k}$ is the number of FEM boundary edges incident to node $k$. For closed manifold domains in 2D, $\eta_{k}=2$ for all $k$.

Similarly, in 3D, the discretized contact potential becomes

where $\mathcal{B}$ is now the set of all FEM boundary triangles, $\hat{\mathcal{B}}$ is the set of all edges on the FEM boundary with $\eta_{e}$ the number of FEM boundary triangles incident to edge $e$, $\tilde{\mathcal{B}}$ the set of all nodes that are in the interior of the FEM boundary surface mesh, and $d^{PT}(\mathbf{x}_{q},t)$ the point-triangle distance between particle $\mathbf{x}_{q}$ and triangle $t$. For closed manifold domains in 3D, $\eta_{e}=2$ for all $e$.

Adding the contact potential into the incremental potential, the minimization problem for time integration is now fully unconstrained:

Since the distance values measured for contacting MPM particle - FEM simplex pairs are all unsigned, Problem (33) may contain a local optimum at configurations with intersections. To be consistent with the continuous constraint Equation (20), it is also constrained that the iterates always stay in the feasible region on one side of the barrier without crossing. This is achieved by applying the interior-point filter line-search algorithm [77] with continuous collision detection (CCD) [78, 60].

To model frictional contact, local frictional forces $F_{k}$ can be added for every active contact pair $k$. For each such pair $k$, at the current state $x$, a consistently oriented sliding basis $T_{k}(x)\in\mathbb{R}^{dm\times(d-1)}$ can be constructed, where $m$ is the total number of colliding nodes and $d$ is the dimension of space, such that $\mathbf{u}_{k}=T_{k}(x)^{T}(\Delta tv^{n}+\Delta t^{2}((1-\gamma)a^{n}+\gamma a^{n+1}))\in\mathbb{R}^{d-1}$ provides the local relative sliding displacement in the frame orthogonal to the distance gradient. The corresponding sliding velocity is then $\mathbf{v}_{k}=\mathbf{u}_{k}/\Delta t\in\mathbb{R}^{d-1}$.

Maximizing dissipation rate subject to the Coulomb constraint defines friction forces variationally [79, 80]

where $\lambda_{k}$ is the contact force magnitude and $\mu$ is the local friction coefficient. This is equivalent to

with $\mathbf{s}(\mathbf{u}_{k})=\frac{\mathbf{u}_{k}}{\|\mathbf{u}_{k}\|}$ when $\|\mathbf{u}_{k}\|>0$, while $\mathbf{s}(\mathbf{u}_{k})$ takes any unit vector when $\|\mathbf{u}_{k}\|=0$. The friction magnitude function, $f$, is nonsmooth with respect to $\mathbf{u}_{k}$ since $f(\|\mathbf{u}_{k}\|)=1$ when $\|\mathbf{u}_{k}\|>0$, and $f(\|\mathbf{u}_{k}\|)\in[0,1]$ when $\|\mathbf{u}_{k}\|=0$. This nonsmoothness would severely slow and even break convergence of gradient-based optimization.

To enable efficient and stable optimization, the friction-velocity relation in the transition to static friction can be mollified by replacing $f$ with a smoothly approximated function. Following Li et al. [57], we use

where $f_{1}^{\prime}(\Delta t\epsilon_{v})=0$ and a velocity magnitude bound $\epsilon_{v}$ (in units of $m/s$) below which sliding velocities $\mathbf{v}_{k}$ are treated as static is defined for bounded approximation error (Figure 5).

Note that the velocity used in our friction model on the MPM side is the interpolated grid velocity at particle quadrature locations, rather than the particle velocity after grid-to-particle transfer. This makes the velocity seen by frictional forces independent from the choice of the particle-grid transfer scheme. This is important because for example in FLIP, the particle velocity does not reflect its displacement ($\mathbf{v}_{q}^{n+1}\neq(\mathbf{x}_{p}^{n+1}-\mathbf{x}_{p}^{n})/\Delta t$) and thus should not be used to define friction in an implicit solve.

However, challenges remain on incorporating friction into the optimization time integration. A major problem is that friction is not a conservative force and there is no well-defined potential such that taking the opposite of its gradient produces the frictional force. Therefore, following Li et al. [57], we fix the friction constraint set $\mathcal{F}$ along with the normal force magnitude $\lambda$ and the tangent operator $T$ during the nonlinear optimization to the last updated value $\mathcal{F}^{j}=\mathcal{F}(x^{j})$, $\lambda^{j}=\lambda(x^{j})$, and $T^{j}=T(x^{j})$, which then makes the lagged friction force integrable with the pseudo-potential

where $\mathcal{F}^{j}$ is the set of all contact pairs with nonzero $\lambda_{k}^{j}$, $f^{\prime}_{0}(y)=f_{1}(y)$, $\bar{\mathbf{u}}_{k}=(T_{k}^{j})^{T}(\Delta tv^{n}+\Delta t^{2}((1-\gamma)a^{n}+\gamma a^{n+1}))$ and so we have $-\nabla D(x)=-\sum_{k\in\mathcal{F}^{j}}\mu\lambda^{j}_{k}T^{j}_{k}f_{1}(\|\bar{\mathbf{u}}_{k}\|)\>\mathbf{s}(\bar{\mathbf{u}}_{k})$, which is a semi-implicit discretization of the frictional force with lagged variables $\lambda^{j}_{k}$ and $T^{j}_{k}$. Then we can iteratively alternate between the nonlinear optimization with fixed $\mathcal{F}$, $\lambda$, and $T$ given as

and friction update until convergence (Algorithm 1). Although the friction convergence is not guaranteed for arbitrarily large time step sizes due to the nonlinearity and asymmetry of the problem, we have confirmed that all our experiments converge with the practical time step sizes applied (Section 6).

In the BFEMP framework, an experimental setup with a subset or all of the FEM nodes prescribed with Dirichlet boundary conditions on their displacements can be applied to model irregular boundaries for MPM. This can not only resolve detailed boundary geometries even when the MPM grid is relatively coarse (Section 6.2), but also provide accurate and controllable friction on the boundary (Section 6.4).

The collision between two elastic rings is simulated to verify the momentum and energy behavior of BFEMP and to demonstrate the robustness of our framework in handling large deformation. This example is modified from the MPM-MPM contact version in [85].

The experiment setup is shown in the left-middle subfigure of Figure 6. The two rings are identical except that the left ring is discretized with MPM, and the right one is discretized with FEM. The inner radius of the ring shape is $3m$, and the outer radius of the ring shape is $4m$. The Young’s modulus is $E=10^{8}Pa$, the Poisson’s ratio is $\nu=0.2$, and the density is $\rho=1000kg/m^{2}$. The MPM ring is discretized by 20098 particles, where the grid spacing is $0.1m$. The FEM ring is discretized by 1830 vertices and 3310 triangles. The gravitational force and the frictional forces are not included. The two rings are placed $2m$ apart and then move towards each other with an initial speed of $40m/s$. The contact active distance and the contact stiffness are set to $\hat{d}=10^{-2}m$ and $\kappa=10^{7}Pa$ respectively. To minimize numerical dissipation, we use the Newmark time integrator with time step size $\Delta t=2\times 10^{-4}s$. The Newton tolerance is $\epsilon_{d}=10^{-6}m/s$. We also compare the APIC and FLIP transfer schemes. Their differences in displacement and stress are small, so only the stress wave propagation with APIC is shown in Figure 6. The energies and momenta over time are plotted for both APIC and FLIP in Figure 7.