Improved Implementation of Approximate Full Mass Matrix Inverse Methods into Material Point Method Simulations

This paper presents a revised and simplified implementation of the FMPM(k) approximate full mass matrix inverse method for Material Point Method simulations, which resolves compatibility conflicts with lumped-mass features and addresses stability and efficiency concerns at higher orders.

The Gist

Imagine you are trying to simulate how a pile of sand, a piece of rubber, or even a crashing car behaves on a computer. This is what the Material Point Method (MPM) does. It's like a digital playground where tiny particles (the material) move around, but the computer calculates the forces and speeds on a fixed grid (like a chessboard) underneath them.

The problem is that the standard way of doing these calculations (called FLIP) is a bit like a noisy, jittery video game. It's fast, but the results can look "fuzzy" or unstable, especially when things crash or move very fast.

To fix this, scientists developed a better way called FMPM(k). Think of this as a "high-definition" mode. Instead of just taking a quick guess at how the grid moves, it does a few extra, more precise calculations (the "k" stands for how many extra steps you take). The more steps you take, the clearer and more stable the simulation becomes.

However, the original "High-Definition" mode had three big headaches:

1. It clashed with other features: It didn't play nice with things like walls, cracks, or when two different materials touched.
2. It was too sensitive: If you tried to use too many steps (high "k"), the simulation would crash or become unstable.
3. It was slow: Doing all those extra steps took a lot of computer power.

This paper by John Nairn is like a mechanic's manual for upgrading that engine. Here is the simple breakdown of the fixes:

1. The "Incremental" Fix (The Staircase Analogy)

The Old Way: Imagine you want to climb a 10-story building. The old method tried to calculate the exact position of the 10th floor all at once, then the 9th, then the 8th, using a complex formula that got messier the higher you went. If you wanted to add a safety railing (a boundary condition) or check for a hole in the floor (contact), it was a nightmare to fit it into that complex formula.

The New Way: The paper suggests climbing one step at a time.

• Step 1: Take a quick guess (the standard method).
• Step 2: Look at where you are, and calculate only the difference needed to get to the next level of accuracy.
• Step 3: Calculate the difference for the next level, and so on.

Why it's better: Because you are only calculating the "difference" (the tiny step up), you can easily stop at any point to check if you hit a wall, if two materials are touching, or if you need to enforce a rule. It's like having a modular staircase where you can easily add a handrail or a safety net at any step without rebuilding the whole building. This solves the "clashing with other features" problem.

2. The "Stability" Fix (The Tightrope Walker)

The Problem: The more steps (higher "k") you take to get a perfect image, the more the simulation wants to wobble and fall off the tightrope. It required such tiny time steps that the simulation would take forever to run.

The Solution: The author found two tricks to keep the tightrope walker steady:

• Blending: Imagine mixing a super-precise but wobbly camera with a slightly blurry but very steady one. By mixing them (blending the math), you get a result that is still very clear but much less likely to crash.
• Periodic Checks: Instead of doing the super-precise calculation every single frame of the movie, you do it every few frames. This keeps the movie smooth and stable without overworking the computer.

3. The "Smart" Fix (The Dynamic Order)

The Problem: You don't always need 10 steps. Sometimes 2 steps are enough. But the computer was blindly doing 10 steps every time, wasting energy.

The Solution: The paper introduces a "Dynamic" mode. It's like a smart thermostat.

• When the simulation is calm (like a block of sand sitting still), the computer says, "I only need 2 steps," and stops early.
• When things get chaotic (like a crash or a shockwave), it says, "Okay, I need 10 steps!" and keeps going until the math settles down.

The Catch: The author admits that finding the perfect "thermostat setting" (knowing exactly when to stop) is tricky. Sometimes the computer gets confused and keeps calculating even when it doesn't need to. So, for now, the best advice is to just pick a "Goldilocks" number (like 4 or 5 steps) that works well for most situations without needing the complex smart thermostat.

The Bottom Line

This paper takes a powerful but finicky tool (FMPM) and makes it robust, compatible, and easier to use.

• Before: You could use the high-precision mode, but it would break if you touched a wall or two materials, and it was hard to tune.
• After: You can use the high-precision mode on complex problems (cracks, collisions, moving walls) without it breaking. It's more stable, and you have options to make it faster or more accurate depending on what you need.

The author has already built this into his software (NairnMPM), so engineers and scientists can now simulate complex crashes and material failures with much higher clarity and fewer headaches.

Technical Summary

1. Problem Statement

The Material Point Method (MPM) is a powerful computational technique for simulating large deformations, but standard implementations often suffer from null-space noise (spurious oscillations) and energy dissipation. To address this, Approximate Full Mass Matrix (FMPM(k)) methods were developed, where $k$ represents the order of approximation. These methods improve the calculation of grid velocities from grid momenta by approximating the inverse of the full mass matrix, reducing noise and enhancing stability.

However, prior implementations of FMPM(k) faced three critical issues that limited their practical application:

1. Incompatibility with Standard MPM Features: Many essential MPM features (grid-based velocity boundary conditions, multimaterial contact, crack contact, and imperfect interfaces) rely on lumped mass matrix calculations. Mixing these with FMPM(k) caused conflicts, leading to degraded results or instability as the order $k$ increased.
2. Temporal Stability Concerns: Replacing the lumped mass matrix with an approximate full mass matrix alters the system dynamics. There was a lack of clarity regarding the stability limits (maximum allowable time step) as $k$ increases, raising fears of ill-conditioning or the need for impractically small time steps.
3. Computational Cost: The cost of FMPM(k) increases linearly with $k$. Without a clear convergence metric, it was difficult to determine the optimal order $k$ for a specific simulation, potentially making high-order calculations inefficient.

2. Methodology

The paper proposes a revised implementation of FMPM(k) based on an incremental approach rather than the previous recursive expansion.

A. The Revised "FMPM Loop"

Instead of summing terms of a Taylor series expansion directly, the new method calculates the velocity increment $\Delta v^{(\ell)}$ between FMPM($\ell$) and FMPM($\ell-1$).

• Mathematical Formulation: The grid velocity $v^{+(k)}$ is expressed as a sum of increments:
  $$v^{+(k)} = \sum_{\ell=1}^{k} \Delta v^{(\ell)}$$
  where $\Delta v^{(\ell)} = A \Delta v^{(\ell-1)}$ and $A = I - S^+S$.
• Advantage: This formulation is independent of the total order $k$ during the recursion. It allows for the insertion of corrective steps (for boundary conditions and contact) at every increment $\ell$, rather than just at the final step.

B. Resolving Feature Conflicts

The incremental nature of the loop allows the integration of standard MPM features into every step of the FMPM calculation:

• Velocity Boundary Conditions (BCs): Instead of applying BCs only to the final velocity (which fails for $k>2$), the new method applies BC constraints to each velocity increment $\Delta v^{(\ell)}$. This ensures that the cumulative velocity satisfies the BCs while maintaining the full mass matrix coupling.
• Multimaterial Contact: Two new incremental contact methods were derived:
  • "Evolving" Method: Restricts current evolving velocities to tangential motion.
  • "Net" Method: Restricts net uncorrected velocities to tangential motion.
  • Result: The "Net" method was found to be superior, correctly reverting to single-material behavior when contact forces are zero and handling non-linear friction laws (like Coulomb friction) without generating artifacts.

C. Stability Analysis and Improvements

The author analyzed the temporal stability using a freely vibrating 1D bar problem:

• Stability Limits: As $k$ increases, the maximum stable Courant number ($C$) decreases but plateaus (e.g., stabilizing around $C \approx 0.24$ for high $k$), proving that FMPM(k) does not become ill-conditioned.
• Blending Options: To improve stability without excessive energy dissipation, the paper proposes blending the full mass matrix with lower-order matrices:
  • Blending with FMPM(2) (rather than the lumped mass FMPM(1)) significantly improves stability while maintaining low energy dissipation.
• Periodic FMPM(k): Running the FMPM loop periodically (every $N$ time steps) rather than every step can extend stability limits and reduce dissipation issues associated with small time steps.

D. Dynamic FMPM(k)

The paper explores Dynamic FMPM, where the order $k$ is determined dynamically during the simulation based on convergence criteria (e.g., when $\Delta v^{(\ell)}$ becomes sufficiently small). While promising for steady-state problems, the study found that for dynamic impact problems with persistent stress waves, the method often required high orders anyway, limiting efficiency gains.

3. Key Contributions

1. Revised Algorithm: A mathematically identical but computationally superior implementation of FMPM(k) using an incremental loop that simplifies coding and enables dynamic adaptation.
2. Compatibility Framework: A robust solution for integrating FMPM(k) with velocity boundary conditions and multimaterial contact, resolving previous conflicts that degraded results for $k > 2$.
3. Stability Characterization: A definitive analysis showing that FMPM(k) stability plateaus at high orders and does not suffer from ill-conditioning, providing safe guidelines for time-step selection.
4. Optimization Strategies: Introduction of blending with FMPM(2) and periodic execution as practical methods to balance stability, accuracy, and computational cost.
5. Open Source Implementation: The algorithms are implemented and validated in the OSParticles and NairnMPM (v19+) codes.

4. Results

• Boundary Conditions: In a moving wall test, the new method maintained low velocity errors (400x lower than FLIP) for all $k$, whereas the previous "combined" method degraded significantly for $k > 2$.
• Contact Mechanics: In shock wave simulations passing through material interfaces, the new "Net" contact method produced results identical to single-material simulations, eliminating the spurious wave reflections seen in prior methods.
• Energy Conservation: In disk impact simulations, FMPM(4) showed significantly lower energy dissipation than FLIP or FMPM(1). The "Net" contact method preserved energy conservation even with frictional interfaces.
• Stability: Stability tests confirmed that FMPM(k) requires smaller time steps than FLIP (limiting $C$ to $\approx 0.24$ for high $k$), but this is a fixed, manageable limit. Blending with FMPM(2) raised this limit to $C \approx 0.40$ without sacrificing energy conservation.
• Dynamic Efficiency: Dynamic FMPM showed mixed results; while it could reduce costs for steady-state problems, it struggled to converge efficiently in high-frequency dynamic events, suggesting that a fixed, modest order (e.g., $k=4$ or $5$) is often the most efficient practical choice.

5. Significance

This paper resolves the primary barriers preventing the widespread adoption of high-order FMPM(k) methods in complex MPM simulations. By decoupling the full mass matrix benefits from the conflicts with standard MPM features (contact and BCs), it enables:

• Higher Accuracy: Reduction of null-space noise and improved shock wave resolution.
• Robustness: Reliable simulation of multimaterial interactions and complex boundary conditions without numerical artifacts.
• Practicality: Clear guidelines on stability limits and the realization that moderate orders ($k=4$ or $5$) often provide the "sweet spot" between accuracy and computational cost.

The work effectively transitions FMPM(k) from a theoretical improvement to a robust, production-ready tool for computational mechanics, particularly for problems involving contact, interfaces, and high-frequency dynamics.