The semi-explicit nonsmooth Newmark time integrator for robust unilateral contact in dynamic fragmentation simulations

本文介绍并验证了一种半显式非光滑 Newmark (NSN) 时间积分方案，该方案通过强力执行约束，能够稳健地处理动态破碎模拟中的单边接触，从而在稳定性与准确性方面优于基于惩罚的方法，并揭示了接触耗散可以通过改善损伤局部化，从而反常地增加碎片数量。

要点

想象一下，你正在尝试模拟一个固体物体（比如陶瓷盘或岩石）在受到剧烈撞击时破碎成数千个微小碎片的场景。这不仅仅是简单的破裂，而是一场混乱的爆炸：碎片四处飞溅、互相碰撞、在墙壁上反弹并相互磨损。

这篇论文介绍了一种新的计算机程序（一种“时间积分器”），旨在模拟这种混乱现象，而不会导致计算机崩溃或产生荒谬的结果。以下是其工作原理的分解，通过简单的概念进行说明：

1. 问题所在：“弹簧”陷阱

为了模拟破碎，科学家通常使用一种方法，即假装材料是由微小的弹簧组成的。当材料断裂时，弹簧会断开。当碎片互相碰撞时，他们使用“惩罚弹簧”（penalty springs）将它们推开，以防止它们互相穿透。

类比： 想象试图用一根橡皮筋来阻止一个保龄球。

• 如果橡皮筋太松（刚度低），球会直接穿过它（不符合物理规律）。
• 如果为了完美阻止球而使用极紧的（高刚度）橡皮筋，它就像一面坚硬的墙。但如果你把它弄得太紧，计算机就必须采取极其微小的步长来计算反弹过程，导致模拟运行得异常缓慢。
• 论文的观点： 旧的方法（使用这些紧绷的弹簧）是不稳定的。它会导致计算机出现漂移、能量损失或崩溃，尤其是在同时发生数百万次碰撞时。

2. 解决方案：“交通警察”（非光滑新马克法）

作者创建了一种名为**非光滑新马克法（Nonsmooth Newmark, NSN）**的新方法。该方法不再使用橡皮筋来推开碎片，而是像一个在繁忙路口工作的严格交通警察。

类比：

• 主体（汽车）： 物体的主体部分自由且平滑地移动。计算机预测如果没有障碍物，汽车将会走向何方。这部分计算速度非常快（显式计算）。
• 接触（路口）： 如果汽车撞到了墙或其他汽车，交通警察就会介入。他不是用弹簧把车推回去，而是瞬间宣布：“停！你不能去那里。”他执行一条硬性规则：禁止穿透。
• 神奇之处： 这种方法将“禁止穿透”规则视为物理学的硬性定律，而非软性的弹簧。这使得计算机可以采取更大的时间步长，因为它不必担心橡皮筋变得过紧。

3. “双重人格”方法

论文将这种方法描述为“半显式”的。可以将其看作是一个两步走的舞蹈：

1. 步骤 A（预测）： 计算机猜测在下一个时刻所有物体的位置，忽略碰撞。
2. 步骤 B（修正）： 如果预测显示有两个碎片发生了重叠，计算机会立即修正它们的运动速度和位置以修复重叠，就像台球撞击另一颗球并瞬间改变方向一样。

这使得模拟既具有快速的特性（像预测阶段），又具有准确性和稳定性（像修正阶段）。

4. 他们的发现（实验）

作者使用三种场景，将这种新的“交通警察”方法与旧有的“橡皮筋”方法进行了对比测试：

• 弹跳球： 一个在地面上弹跳的简单球体。新方法与现有的最佳方法一样准确，但处理弹跳时不会出现能量损失或抖动。
• 撞击棒： 一根金属棒撞击墙壁。旧方法在处理冲击速度方面表现挣扎，但新方法完美处理了这种“挤压”过程，保持了能量计算的正确性。
• 破碎棒： 一根已经产生裂纹的棒子撞击墙壁。旧方法为了保持稳定需要如此微小的步长，以至于运行速度极其缓慢。新方法可以采取巨大的步长，运行速度快了 27 倍，且更加准确。

5. 令人惊讶的发现：受限破碎

论文中最有趣的部分涉及一个“受限”实验。想象一根棒子是在一个小盒子里而不是在开放空间中破碎。

• 旧有的直觉： 你可能会认为，如果碎片在撞击墙壁时会损失能量（耗散），那么剩下的能量就会更少，从而导致产生的碎片更少、更大。
• 论文的发现： 情况恰恰相反。当碎片在墙壁上反弹并损失一些能量（接触耗散）时，材料实际上破碎成了更多、更小的碎片。
• 原因： 作者解释说，“反弹”起到了过滤器的作用。在一个完全弹性（极具弹性）的世界里，应力波会疯狂地来回反弹，导致材料变得“混乱”，产生许多细小但无法完全分离的微弱裂纹。当墙壁吸收了部分能量时，这些波动会变得平缓。这使得应力能够集中在特定点上，驱动裂纹贯穿始终，从而产生干净、独立的碎片。

总结

这篇论文提出了一种新的数学工具，通过将碰撞视为硬性的、瞬时的规则，而非软性的弹簧，来模拟物体的破碎。这使得计算机模拟：

1. 更稳定： 不会崩溃或产生漂移。
2. 更快： 它可以采取更大的时间步长。
3. 更准确： 它能正确预测物体会破碎成多少块。

作者得出结论，该工具已准备好用于复杂的 3D 模拟，例如理解太空碎片如何破碎，或岩石如何在雪崩中破碎，因为它为处理数百万个碰撞碎片的混乱舞蹈提供了一种稳健的方法。

技术摘要

技术摘要：用于动态破碎模拟中鲁棒单边接触的半显式非光滑 Newmark 时间积分器

问题陈述
动态破碎的数值模拟面临双重挑战：既要捕捉固体连续体破碎成离散碎片的快速拓扑演化，又要解决碎片之间以及裂纹面之间发生的密集、同步接触问题。虽然采用外在内聚带模型（CZM）的有限元法（FEM）能有效处理断裂，但传统的接触处理方法面临显著困难。常用于显式时间积分的惩罚法存在根本性的权衡：低刚度会导致非物理性的相互穿透，而高刚度则会将稳定时间步限制在计算上难以承受的极低水平。此外，断裂与正则化接触之间的相互作用引入了非光滑性，导致持续的数值能量漂移和不稳定，从而掩盖了底层的物理现象。

方法论
作者提出并验证了一种新型半显式时间积分方案：非光滑 Newmark（NSN）法。该方法通过将整体动力学与接触相互作用解耦，将非光滑接触动力学（NSCD）框架应用于动态破碎：

1. 半显式公式化： 该方案将动力学分为“光滑”部分（整体变形、外部载荷和内聚断裂）和“非光滑”部分（接触反力）。
   • 光滑动力学使用二阶 Newmark-$\beta$ 方案（$\beta=0, \gamma=1/2$）进行显式积分，在自由飞行阶段保持高精度。
   • 非光滑动力学（单边接触）在速度层面上通过一阶隐式欧拉法进行隐式处理。这强力执行了 Signorini 间隙约束，且不引入人工柔度。
2. 算法实现： 接触问题被表述为一个线性互补问题（LCP），并通过将其限制在活跃接触集上的凸二次规划（QP）问题来求解。这实现了高效、可并行的求解，其规模随活跃接触数量而非系统总规模进行扩展。
3. 修正牵引-分离法则（TSL）： 为了解决显式方案中损伤趋于零时内聚刚度趋于无穷大的数值挑战，作者引入了一个刚度上限（$\tilde{k}$）。这种正则化防止了时间步趋于零，同时保持了接触问题的凸性。
4. 验证与应用： 该框架在开源 Akantu FEM 库中实现。通过碰撞球和冲击棒等基准测试，该方法通过了针对解析解及成熟非光滑方法（Moreau-Jean 和 CD-Lagrange）的验证。随后，该方法被应用于自由膨胀和受限膨胀下的 1D 动态破碎研究。

核心贡献与结果

• 精度与稳定性： 基准测试表明，NSN 方案达到了与成熟非光滑方法（Moremore-Jean 和 CD-Lagrange）相当的精度，并显著优于基于惩罚的方法。在冲击棒测试中，NSN 在位移方面保持了二阶精度，并在速度精度方面优于 CD-Lagrange，尤其是在弹性阶段。
• 计算效率： 尽管 NSN 方法每步的计算成本较高（由于 QP 求解，约为纯显式惩罚法的 5 倍），但其鲁棒的稳定性允许使用显著更大的时间步。对于 1D 基准测试，其总求解时间与纯显式方法相当或更好，且能量守恒误差降低了数个数量级（接近机器精度）。
• 受限破碎中的物理洞察：
  • 能量预算转移： 在受限膨胀场景下（碎片与刚性壁碰撞），断裂的能量预算从单个碎片的局部动能转向整个系统的全局动能。与自由膨胀相比，这导致平均碎片尺寸大幅减小，并增加了总断裂能量耗散。
  • 接触耗散的反直觉作用： 研究发现，引入接触耗散（非弹性碰撞，$e < 1$）会减少断裂的总能量耗散，但会增加最终的碎片数量。这是因为接触耗散起到了低通滤波器的作用，抑制了高频应力波振荡。这种滤波促进了更清晰的应力波传播和更好的损伤定位，驱动裂纹走向完全分离，而不是让损伤停留在部分损伤区域（$d < 1$）。

意义
本文确立了 NSN 方案作为在高密度接触场景下生成高保真破碎统计数据的鲁棒工具。通过在保留显式效率处理整体动力学的同时，严格地将接触视为非光滑过程，该方法克服了惩罚法固有的数值不稳定。作者声称，该框架对于模拟复杂的多体相互作用（如轨道碎片场景）至关重要，其中频繁的碎片间相互作用驱动了二次破碎。这项工作为将高保真破碎分析扩展到更高维度迈出了关键一步，利用了该方法在不依赖全局隐式求解的高昂代价下保持物理保真度的能力。