Frictional-Contact Pipeline in Implicit MPM

• The paper introduces a frictional-contact pipeline that integrates GIMP-based deformable body representation with a monolithic, implicit Newton–Newmark scheme.
• It employs robust contact detection via closest-point search and enforces normal and tangential friction using penalty methods and consistent return mapping.
• Validated on benchmarks like soil-structure interaction, the framework effectively models complex, nonlinear contact evolution between deformable and rigid bodies.

The frictional-contact pipeline for implicit Material Point Method (MPM) provides a robust numerical framework for simulating large-deformation contact between highly deformable materials and rigid bodies with enforceable slip/stick Coulomb friction. This pipeline, as developed in Bird et al. (2024), features generalized interpolation material point (GIMP) domains for representing deformable bodies, explicit triangulation of rigid surfaces, contact localization via closest-point search, consistent stick-slip friction, and tight integration into an implicit Newton–Raphson/Newmark time-integration scheme (Bird et al., 2024). The approach addresses the computational challenges associated with evolving, intersecting geometries under history-dependent, nonlinear material response and complex contact evolution.

1. Kinematic Representation and Contact Detection

The pipeline represents each deformable material point $p$ as a mass $m_p$ and volume $V_p$ augmented by a cuboid domain $\Omega_p$, forming the basis for GIMP interpolation. The union of these domains captures the current spatial extent of the deforming body without explicit surface reconstruction. The rigid body is meshed with planar triangles $\ell$, each supporting a local $(\xi^1, \xi^2)$ affine coordinate system and parameterized by rigid generalized coordinates $\theta$.

Contact detection is performed for each GIMP domain corner, executing a closest-point projection against the rigid mesh. For each candidate triangle, the normal gap $g_N(x_p, \xi) = (x_p - x'(\xi)) \cdot n(\xi)$ is minimized with respect to $\xi$ via Newton iteration. If $g_N < 0$ and the projection lies inside the triangle, the material point is designated in-contact ($p \in P_c$). This procedure is robust to large penetrations and mesh misalignments due to the volumetric character of GIMP domains.

2. Variational Formulation and Discretization

The dynamic evolution of the system is governed by the weak form of force balance over the deformed domain $\phi_t(\Omega)$. After spatial discretization and use of GIMP basis functions $S_{vp}$, the nodal residual at grid node $v$ includes elastic, body-force, inertial, and contact components: $$R_v = \sum_p [ (\nabla S_{vp})^T \sigma_p V_p - (S_{vp})^T b_p V_p + (S_{vp})^T m_p \ddot{x}_p ] - \sum_{p\in P_c}[F_{N,vp}^{(c)} + F_{T,vp}^{(c)}] = 0$$ where $\sigma_p$ is the Cauchy stress, $F_{N,vp}^{(c)}$ and $F_{T,vp}^{(c)}$ denote nodal contact forces in the normal and tangential directions, respectively.

Normal contact is enforced by a penalty method: $p_N = \epsilon_N g_N$, where $\epsilon_N$ is chosen to enforce near-zero gap without inducing ill-conditioning. The associated contact virtual work and residual entries are formulated to admit direct linearization for Newton iteration.

Tangential (frictional) constraints employ a penalty-regularized Coulomb friction law. The relative tangential motion $\Delta g_T$ is split into stick and slip components, with the friction force $p_T$ determined by a trial step, plastic return if yielding is exceeded, and enforcement on the yield surface: $f = \|p_{tr}\| - \mu |p_N| \leq 0$ The method supports both sticking and slip, with consistent return mapping for plastic slip analogous to constitutive frictional models.

3. Assembly, Linearization, and Implicit Solution

At each Newton–Raphson iteration, the global residual $R$ and tangent $K = \partial R/\partial u$ (where $u$ collects all grid nodal displacements) are assembled, including consistent contributions from contact pairs. The stiffness matrix includes elastic, inertial, and detailed contact blocks capturing normal and tangential frictional effects: $$K = K_{\mathrm{int}} + \frac{1}{\beta \Delta t^2} M + K_c$$ with $K_c = K^{(N)} + K^{(T)}$ arising from contact linearization. Linearization of contact contributions uses consistent derivatives of contact force increments with respect to both positions and normals (including orientation changes of the rigid triangles), as detailed in the explicit formulations of $\delta F_{N,vp}$ and $\delta F_{T,vp}$.

4. Implicit Time Integration and Monolithic Update

The time-stepping algorithm leverages an implicit Newmark–Newton integration scheme. Each step initializes the material point and rigid body states, projects velocities and accelerations to the grid, and executes a Newton loop:

• Update grid velocities/accelerations using the Newmark relations.
• Evaluate stress updates, internal force contributions, and consistent tangent operators at the material points.
• Re-detect contact and evaluate the current normal and tangential gaps, forces, and update stick/slip decisions.
• Assemble the full residual and tangent, solve for $\Delta u$, and update unknowns.
• Upon Newton convergence, update particle states, plastic/internal variables, and rigid body positions.
• Remesh or reset the background grid as appropriate, then proceed to the next time step.

Throughout each Newton iteration, contact detection and force computation is repeated, ensuring the compatibility of the residual and linearization with the current iterate. Rigid body dynamics (generalized coordinates $\theta$) are solved monolithically within the global system, incorporating truss inertia and contact forces.

5. Implementation Details and Stabilization

Key parameters include normal and tangential penalty stiffnesses ($\epsilon_N \approx 10-100 E_{\mathrm{local}}$, $\epsilon_T \approx 10-50 E_{\mathrm{local}}$). For small cut elements arising near boundaries, ghost-penalty stabilization terms may be included to improve mass and stiffness matrix conditioning. Large deformations are handled via an updated Lagrangian formulation with multiplicative $F=F_eF_p$ and logarithmic strain updates (no modification to small-strain plasticity solvers is required).

After each time step, GIMP domain positions are "snapped back" by minimizing an energy based on previous and current penalty-weighted gaps, mitigating spurious jumps in contact gap for material points.

6. Contact Model Validation and Application Scope

The pipeline is validated on benchmark problems where analytical solutions exist and in soil-structure interaction studies driven by geotechnical centrifuge experiments (Bird et al., 2024). The method demonstrates accurate normal/tangential traction recovery, robust stick-slip transitions, stable large deformation behavior, and effective simulation of interaction between highly deformable (e.g., soil) and rigid (e.g., structural) bodies. The framework admits extension to multiple rigid and deformable bodies, provided the contact discretization and penalty parameters are tuned to the local problem.

7. Relation to Broader Implicit MPM Contact Research

This pipeline closely aligns with the trend towards monolithic, penalty-regularized, implicit contact solvers in MPM and related multiphysics methods. Notable distinctions include its avoidance of explicit boundary reconstruction for deformable materials, consistent stick-slip frictional treatment with return mapping, and monolithic coupling of rigid dynamics. The methodology is compatible with advances in level-set-based implicit contact (Liu et al., 2020), convex variational formulations (Zong et al., 2024), cone-complementarity solvers (Ménager et al., 2 Feb 2026), and barrier-based particle-mesh coupling (Li et al., 2021), and provides a reference architecture for practical engineering analysis involving complex, highly nonlinear contact scenarios.