A Total Lagrangian Finite Element Framework for Multibody Dynamics: Part II -- GPU Implementation and Numerical Experiments

1. Problem Statement

The paper addresses the computational challenges of simulating finite-deformation flexible multibody dynamics (FMBD) involving large rotations and deformations. While the theoretical framework (Total Lagrangian formulation, element kinematics, and material models) was established in a companion paper [1], this work focuses on the numerical implementation and GPU acceleration required to achieve real-time or faster-than-real-time simulation rates.

Key challenges addressed include:

• Nonlinear Solvers: Efficiently solving the coupled, nonlinear residual equations arising from implicit time integration and bilateral constraints.
• Scalability: Handling large-scale meshes (up to ~1 million degrees of freedom) with complex contact interactions.
• Constraint Enforcement: Managing bilateral kinematic constraints (e.g., joints) without the sensitivity issues associated with penalty methods or the complexity of index reduction.
• Collision Detection: Performing robust, high-throughput collision detection for deformable bodies on GPUs, avoiding the bottlenecks of traditional tree-based structures.

2. Methodology

A. Mathematical Formulation

• Time Integration: Uses a velocity-based implicit Backward-Euler scheme. The equations of motion are reduced to a nonlinear residual in velocity ($v$).
• Constraint Enforcement: Employs an Augmented Lagrangian Method (ALM). This creates a two-level structure:
  • Outer Loop: Updates Lagrange multipliers ($\lambda$) to enforce constraints ($c(q)=0$).
  • Inner Loop: Solves for velocity $v$ minimizing the augmented cost function $\Phi_\rho$.
• Material Models: Supports hyperelasticity (St. Venant–Kirchhoff and Mooney–Rivlin) with optional Kelvin–Voigt viscous damping.

B. GPU Implementation Strategy

The authors introduce a highly optimized, GPU-native framework with several key architectural decisions:

1. Fixed-Sparsity Matrix Strategy:

   • The sparsity patterns of the mass matrix, constraint Jacobian, and Newton Hessian are determined once at initialization based on the mesh topology.
   • This eliminates repeated symbolic analysis during the simulation.
   • The cuDSS library (NVIDIA's GPU-resident direct sparse solver) is used to perform numerical refactorization on the fixed pattern at every Newton iteration, significantly reducing overhead compared to full factorization.

2. Two-Stage Parallelization for Internal Forces & Stiffness:

   • Stage 1 (Stress Evaluation): Assigns one thread per (element, quadrature point) pair to compute the First Piola–Kirchhoff stress tensor ($P$). This is fully data-parallel with no write conflicts.
   • Stage 2 (Assembly): Assigns one thread per (element, node) pair. Threads accumulate local contributions into the global force vector using atomic additions.
   • The tangent stiffness assembly reuses the cached stress tensors from Stage 1 to avoid redundant constitutive evaluations.

3. Dual Inner Solvers:

   • First-Order (AdamW): A fully parallel optimizer requiring only the gradient. Suitable for coarse meshes or when Hessian assembly is too costly.
   • Second-Order (Newton): Assembles the sparse global Hessian ($H = \frac{1}{h}M + hK_t + h^2\rho J^T J$) and solves the linear system using cuDSS. This offers superlinear convergence.

4. BVH-Free Collision Detection:

   • Broad-Phase: Instead of Bounding Volume Hierarchies (BVH), which require expensive tree updates for deformable bodies, the algorithm uses a bin-based spatial partitioning on "triangle soups." This maps naturally to GPU thread blocks and avoids per-step tree restructuring.
   • Asynchronous Pipelining: A two-thread architecture is used:
     • Kinematics Thread: Updates contact candidates (Active Contact Set) using the bin-based method.
     • Dynamics Thread: Performs physics integration and narrow-phase force calculation.
     • The threads run asynchronously, overlapping collision detection costs with physics integration.
   • Narrow-Phase: Uses projection-based triangle-triangle overlap tests and a "patch-level" reduction to smooth contact normals for stable friction modeling.

3. Key Contributions

1. GPU-Resident Sparse Direct Solving: Demonstrates the viability of using cuDSS for implicit nonlinear FEA on GPUs by leveraging a fixed-sparsity strategy to enable efficient refactorization.
2. Scalable Parallelization: A novel two-stage parallelization strategy for internal force and tangent stiffness evaluation that scales to arbitrary mesh sizes without restructuring.
3. BVH-Free Collision Detection: A robust, GPU-native collision algorithm that replaces hierarchical tree traversal with bin-based partitioning and asynchronous execution, eliminating the bottleneck of dynamic tree updates.
4. Dual-Solver Framework: Provides both a first-order (AdamW) and second-order (Newton) solver within the same ALM framework, allowing users to trade off convergence speed against assembly cost.
5. Comprehensive Validation: Rigorous unit tests validating frictional contact (stick-slip and oblique impact) and bilateral constraints (revolute and spherical joints) against analytical rigid-body solutions.

4. Results and Performance

The framework was benchmarked on an NVIDIA GeForce RTX 5090 against CPU baselines (FEniCS for T10 elements, Project Chrono for ANCF elements).

• Scaling Benchmarks:

  • Tested on three element types: T10 (tetrahedral), ANCF3243 (beam), and ANCF3443 (shell).
  • Tested across six mesh resolutions (up to ~500k DOFs for T10, ~6M DOFs for shells).
  • Performance Gain: The Newton solver achieved approximately one order of magnitude (10x) reduction in Real-Time Factor (RTF) compared to CPU baselines at the largest resolutions.
    • Example (T10, RES32): GPU Newton RTF $\approx$ 1,250 vs. CPU FEniCS RTF $\approx$ 27,597 (meaning the GPU is ~22x faster).
  • The AdamW solver performed well at coarse resolutions but struggled with convergence at very fine resolutions due to first-order convergence rates.

• Large-Scale Demonstrations:

  • Geometrically Complex Meshes: Successfully simulated a Utah Teapot, Stanford Bunny, and Deformable Tire with complex topologies.
  • Mixed-Item Dropping: A stress test involving 9 deformable tires, a rigid container, and a flexible ANCF shell beam (~934k free DOFs). The system handled ~67% active contact time with stable convergence, validating the pipeline's robustness in chaotic, many-body scenarios.

• Validation:

  • Frictional Contact: Quasi-static brick sliding and dynamic oblique impact tests matched closed-form rigid-body analytical predictions (Coulomb friction, coefficient of restitution) with high fidelity, especially when material damping was sufficient to suppress elastic oscillations.
  • Constraints: Double-pendulum simulations with revolute and spherical joints showed constraint residuals below $10^{-8}$ and joint reaction forces matching analytical equilibrium predictions.

5. Significance

This paper represents a significant advancement in high-performance multibody dynamics simulation. By moving the entire pipeline—including sparse matrix assembly, factorization, and collision detection—onto the GPU, the authors have demonstrated that implicit, second-order accurate simulations of large-deformation flexible bodies can be performed at speeds approaching or exceeding real-time.

The work bridges the gap between theoretical finite element formulations and practical, large-scale engineering applications (e.g., soft robotics, deployable space structures, and surgical simulation). The elimination of BVH updates and the use of fixed-sparsity refactorization are particularly notable innovations that overcome historical bottlenecks in GPU-based FEA. The framework is implemented within the open-source Project Chrono ecosystem, making these high-performance capabilities accessible to the broader research community.