Implicit Bonded Discrete Element Method with Manifold Optimization

TL;DR

This work addresses efficient, accurate fracture simulation using Bonded Discrete Element Methods (BDEM) by introducing an optimization-based implicit integrator and a manifold optimization framework for quaternion constraints. A constant-mass variational formulation with $M_s$ enables a stable, second-order symplectic update, while a nullspace-based, second-order manifold solver dramatically speeds up quaternion handling. The approach yields substantial speedups (up to 12×) over state-of-the-art methods and demonstrates strong scale-consistency compared with FEM/MPM in fragmentation tasks, supported by tests on ceramic plates and fabric tearing. The method’s practical impact lies in enabling fast, realistic fracture simulations with coarse-to-fine resolution and robust convergence, while remaining extensible to collision handling and surface reconstruction.

Abstract

This paper proposes a novel approach that combines variational integration with the bonded discrete element method (BDEM) to achieve faster and more accurate fracture simulations. The approach leverages the efficiency of implicit integration and the accuracy of BDEM in modeling fracture phenomena. We introduce a variational integrator and a manifold optimization approach utilizing a nullspace operator to speed up the solving of quaternion-constrained systems. Additionally, the paper presents an element packing and surface reconstruction method specifically designed for bonded discrete element methods. Results from the experiments prove that the proposed method offers 2.8 to 12 times faster state-of-the-art methods.

Figures (11)

• Figure 1: showcases the experimental results of our method on a chocolate drop impact simulation. The figure illustrates the moment when a chocolate drop hits the ground and cracks into several pieces and small fragments. The simulation shows the effectiveness of our approach in accurately modeling high-speed object collisions and fracture phenomena.
• Figure 2: A ceramic plate dropped on a hard floor and broken into pieces. Our simulation captures highly detailed fractures using only 219K particles, with a 2.4x speedup over the state-of-the-art method.
• Figure 3: This scene illustrates the simulation of a woven fabric undergoing progressive stretching and eventual rupture using our method. With only 165K particles, our approach effectively captures the realistic behavior of yarn-level cloth and the phenomenon of fabric tearing, demonstrating a convincing rupture process and rich details at the tear propagation. Compared to current explicit methods, our technique achieves a threefold acceleration.
• Figure 4: The figure, from left to right, illustrates the performance of the FEM, MPM, and our approach across varying scales. The upper, middle, and lower rows represent the fracture behavior of the three methods under identical loads at different scales. It is observable that FEM exhibits completely distinct fracture occurrences across the three scales. MPM shows similar fracture timings at medium and high scales but deviates at lower scales. In contrast, our method consistently demonstrates fracture at the same load across all scales, indicating superior scale consistency.
• Figure 5: We conducted a three-point beam bending experiment with identical settings for our method and rotation-vector formulation in geilinger2020add. The horizontal axis represents the time step size used in the simulation, with curves of different colors indicating various material Young's moduli. The vertical axis denotes the average number of iterations required during the simulation process. The dashed lines in the figure represent the rotation vector formulation, while the solid lines represent our approach. It is observable that, across different time steps and stiffness coefficients, our method needs fewer iterations to achieve the same convergence criteria.