Stability of Extrinsic Cohesive-Zone Model with Penalty-Based Contact in Explicit Dynamic Fragmentation Simulations

This study identifies that combining extrinsic cohesive-zone models with penalty-based contact in explicit dynamic fragmentation simulations causes severe unphysical energy growth and artificial fragmentation due to stiffness discontinuities and switching errors, ultimately concluding that this approach is unsuitable for long-term, energy-consistent simulations despite the partial mitigation offered by adaptive penalty strategies.

The Gist

The Big Picture: Simulating a Shattering Glass

Imagine you are trying to predict exactly how a glass window will shatter when hit by a rock. You want to know not just that it breaks, but how many pieces it makes, how big they are, and how fast they fly. To do this, scientists use computer simulations.

This paper investigates a specific type of computer simulation used for high-speed explosions or impacts. The researchers found that their simulations were "lying" to them. Instead of showing a stable break, the computer was creating fake, endless shattering and inventing energy out of thin air.

They set out to find out: Why is the computer glitching, and how do we fix it?

The Setup: The "Glue" and the "Spring"

To simulate breaking, the researchers used two main tools in their computer model:

1. The "Glue" (Cohesive Zone Model): Imagine the material is made of tiny Lego bricks. Between the bricks, there is invisible, stretchy glue. When you pull the bricks apart, the glue stretches and eventually snaps. This models how a crack starts and grows.
2. The "Spring" (Penalty Contact): Once the glue snaps and the bricks separate, they might bounce back and hit each other. To stop them from passing through one another (which is physically impossible), the computer uses a "spring" rule. If two bricks try to overlap, the spring pushes them apart. The stiffer the spring, the harder it is to overlap.

The Problem: The "Bouncy Castle" Effect

When they ran the simulation, the computer started behaving like a bouncy castle that never stops bouncing.

• The Symptom: The total energy in the simulation kept growing higher and higher, even though no new energy was being added.
• The Result: The computer thought the material was breaking into millions of tiny, impossible pieces. The "fragment count" (number of pieces) kept rising forever, which is physically impossible.

The researchers asked: Is the glue too weak? Is the spring too stiff? Or is the math itself broken?

The Investigation: Three Suspects

The team tested three possible reasons for the glitch, like a detective ruling out suspects.

Suspect 1: The "Brand New Glue" (Diverging Initial Stiffness)

The Theory: When a piece of glue is first created (before it stretches at all), it is incredibly stiff. Theoretically, it's infinitely stiff.
The Test: They checked if this super-stiffness was making the computer calculations unstable.
The Verdict: Not the main culprit. While it can cause problems, in their specific test, the glue didn't get stiff enough to break the simulation. It was a red herring.

Suspect 2: The "Softening" (Gradual Weakening)

The Theory: As the glue stretches and breaks, it gets weaker (softens). Maybe this change in strength confused the computer.
The Test: They analyzed the math of the glue getting weaker.
The Verdict: Innocent. The math showed that when the glue weakens, the energy lost is perfectly balanced by the energy used to create the new crack surface. This part of the simulation was actually working correctly.

Suspect 3: The "Switch" (Cohesive-Contact Transition) — THE GUILTY PARTY

The Theory: This is the real problem. Imagine a piece of the material is vibrating. It stretches (glue mode), then snaps back and touches another piece (contact mode).

• In Glue Mode, the material acts like a specific type of spring.
• In Contact Mode, the material acts like a different type of spring (the penalty spring).

The problem is that the computer has to instantly switch from one spring rule to the other the moment the pieces touch. It's like driving a car that suddenly switches from "Gas" to "Brake" every time you hit a bump.
The Result: Every time the material switches between "glue" and "contact," the computer makes a tiny math error. It accidentally adds a tiny bit of energy.

• The Analogy: Imagine a child on a swing. Every time they reach the top, you accidentally give them a tiny, invisible push. You don't notice it at first, but after 1,000 swings, the child is flying so high they hit the ceiling.
• The Reality: In the simulation, these tiny energy errors piled up over millions of steps, causing the "fake energy" explosion and the endless shattering.

The Proposed "Fix" and Why It's Not a Real Fix

The researchers tried a clever trick to stop the glitch. They made the "Contact Spring" change its stiffness to match the "Glue Spring" exactly.

• The Result: The sudden "switch" disappeared. The energy stopped growing. The simulation became stable.
• The Catch: To make the springs match, the "Contact Spring" had to become very weak when the glue was damaged. This meant the pieces of the material were allowed to overlap (interpenetrate) significantly.
• The Conclusion: While this fixed the math glitch, it broke the physics. You can't have a simulation where solid pieces pass through each other just to make the numbers work. So, this "fix" is useful for diagnosing the problem, but it's not a solution for real-world engineering.

The Final Takeaway

The paper concludes that using "penalty springs" to handle contact in high-speed breaking simulations is fundamentally flawed for long-term accuracy.

• The Trade-off: You can't have it all. If you make the contact spring very stiff to stop pieces from overlapping, you force the computer to take tiny, slow steps. If you make it softer to speed things up, you get energy errors and fake shattering.
• The Future: The authors suggest that instead of using "soft springs" (penalty methods), we need to use "hard rules" (nonsmooth mechanics) that treat contact like a strict law rather than a spring. This would stop the energy leaks and allow for accurate, long-term simulations of how things shatter.

In short: The computer was hallucinating a never-ending explosion because it got confused every time a broken piece bounced off another piece. The "spring" method used to stop them from overlapping was the cause of the confusion, and the only way to truly fix it is to change the rules of how the computer handles collisions entirely.

Technical Summary

Technical Summary: Stability of Extrinsic Cohesive-Zone Model with Penalty-Based Contact in Explicit Dynamic Fragmentation Simulations

Problem Statement
Dynamic fragmentation simulations are critical for predicting material response under high strain rates, such as hypervelocity collisions in orbital debris scenarios. While explicit dynamic simulations combining an extrinsic cohesive-zone model (CZM) with penalty-based contact are popular for their scalability, they frequently exhibit severe numerical instabilities. In benchmark tests, these instabilities manifest as exponential energy growth and artificial, unphysical fragmentation, invalidating statistical outputs like fragment mass and velocity distributions. Despite the prevalence of this issue, the specific root causes of the instability in the coupling of extrinsic CZM and penalty contact within explicit time integration have not been systematically isolated or quantified.

Methodology
The authors employ a numerical framework implemented in the open-source code Akantu, utilizing an extrinsic CZM with a linear traction-separation law (TSL) coupled with a penalty-based contact formulation. The simulations use an explicit central-difference (CD) time-integration scheme.

To diagnose the instability, the authors:

1. Established a Benchmark: Simulated the dynamic fragmentation of a 2D square plate (analogous to an exploding hollow sphere) using AD-995 alumina properties, with cohesive strengths sampled from a Weibull distribution.
2. Isolated Instability Mechanisms: Systematically analyzed three potential sources of instability using simplified models:
   • Diverging Initial Stiffness: Theoretically analyzing the effect of newly inserted cohesive elements having near-zero damage (infinite stiffness) on the stable time step.
   • Cohesive–Contact Switching: Modeling the abrupt transition between cohesive stiffness ($k^+$) and contact penalty stiffness ($k^-$) using a single-degree-of-freedom (SDOF) system to quantify energy errors per step.
   • Softening: Analyzing the energy balance during the progressive reduction of cohesive stiffness (damage evolution).
3. Quantified Energy Drift: Used analytical error estimates, phase-space diagnostics, and energy growth metrics to determine how small per-step errors accumulate into long-term energy drift.
4. Tested Mitigation: Assessed an adaptive penalty strategy where the contact stiffness ($k^-$) is dynamically tied to the evolving cohesive stiffness ($k^+$) to eliminate the stiffness discontinuity.

Key Contributions and Results

• Identification of Three Mechanisms: The study isolates and quantifies three distinct sources of instability:

  1. Diverging Initial Cohesive Stiffness: As damage approaches zero, cohesive stiffness diverges, constraining the stable time step. However, in the specific high-strain-rate benchmark used, this mechanism was not the primary cause of the observed global instability, as cohesive elements rapidly overshot the unstable region.
  2. Discontinuous Stiffness Jumps (Primary Cause): The abrupt switch between cohesive stiffness ($k^+$) and a fixed contact penalty ($k^-$) at the interface creates a discontinuity in the traction-separation response. The authors demonstrate that repeated switching between these states generates alternating energy injections and dissipation. Under specific parameter ranges (stiffness ratios and time steps), these small per-step errors accumulate, leading to unbounded energy growth and artificial fragmentation, even when the standard Courant–Friedrichs–Lewy (CFL) stability condition is satisfied.
  3. Discontinuity from Softening: While softening introduces a non-smooth path, the study finds that the resulting algorithmic energy loss is exactly balanced by the physical fracture dissipation. Therefore, softening events alone are numerically conservative and do not trigger instability.

• Quantification of the "Cohesive–Contact" Instability: Using an SDOF model, the authors show that the classical stability bound (derived for linear time-invariant systems) is insufficient for piecewise-linear systems with switching stiffness. They map regions of stability and instability in the parameter space of the stiffness ratio ($\alpha_{ref} = k^+/k^-$) and normalized time step. The results indicate that for frequent switching with significant stiffness contrast, the time step must be reduced significantly below the standard limit (e.g., $\Delta t \lesssim 0.2 \Delta t_c$) to maintain stability.

• Evaluation of Adaptive Penalty: The authors tested a diagnostic modification where the contact penalty stiffness is tied to the cohesive stiffness ($k^- = k^+(d)$). This modification successfully eliminates the stiffness discontinuity, restores energy conservation, and stabilizes the fragment count. However, it introduces significant non-physical interpenetration as damage increases, rendering it unsuitable as a definitive physical remedy but effective as a diagnostic tool to confirm the root cause.

Significance and Claims
The paper concludes that the coupling of extrinsic CZM with standard penalty-based contact in explicit dynamics is not a viable approach for long-term, energy-consistent fragmentation simulations requiring physically meaningful fragment statistics. The primary limitation is the inherent trade-off in penalty methods: high stiffness prevents interpenetration but restricts the time step and introduces energy artifacts via discontinuities, while lower stiffness allows larger time steps but permits interpenetration.

The authors assert that the root cause of the observed unphysical energy drift is the repeated cohesive–contact switching, which accumulates numerical errors. They argue that heuristic fixes (e.g., local time-step reduction, mass redistribution) only mitigate symptoms. Instead, the paper advocates for the adoption of nonsmooth time-stepping schemes that enforce contact at the velocity level (using impulses and complementarity constraints). Such methods would eliminate the artificial stiffness jump, shift stability limits back to bulk and cohesive dynamics, and allow for robust, long-duration simulations without the energy drift associated with penalty-based contact.