有限元碰撞，实现前所未有的精确、稳定、无穿透仿真效果

有限元碰撞，实现前所未有的精确、稳定、无穿透仿真效果 | SIGGRAPH 2020

2020 年 7 月 22 日 AI科技评论

作者 | 李旻辰

编辑 | 丛末

IPC，全称 Incremental Potential Contact。

其中 Incremental Potential是基于数学优化的有限元物理仿真方法在每一个时间歩要最小化的目标函数；contact包括碰撞与摩擦，在仿真中通常需要被描述成额外的约束条件。IPC将物体的contact巧妙地融合到目标函数Incremental Potential中，能在不直接引入约束条件的情况下达到前所未有的精确、稳定、无穿透仿真效果：

这一切得益于IPC对所有碰撞情景的严谨数学定义、对障碍函数法的降维加速、对静动摩擦转换的可控建模、以及对摩擦算子的巧妙离散。

论文地址：https://ipc-sim.github.io/file/IPC-paper-fullRes.pdf

开源地址：http://ipc-sim.github.io/

碰撞的严谨定义

首先，在不考虑碰撞时，弹性物体的所有采样点的位置

在每个时间歩可通过最小化Incremental Potential

来更新：（使用隐式欧拉时间积分）

这里

，其中

为采样点速度，

为外力（比如重力），

为采样点的质量矩阵，

为时间歩长。

是弹性势能。最小化Incremental Potential是一个已经解决的很好的问题，用牛顿迭代辅助line search就能保证稳定、精确地求解。但如果考虑碰撞呢？其实也很简单：

加一个“无穿透”的约束条件就好，但如何为离散曲面定义无穿透的数学表达式，以及如何求解就显得尤为尴尬了。

之前大家试过很多方法，比如定义物体表面距离很近的“点-三角”对和“边-边”对（接触元对）组成的四面体的有符号体积要恒大于0；或者在每个时间歩开始时，找到距离近的每个接触元对上距离最近的点对，限定这些点对的距离大于0等。这些定义的问题在于，如果物体在某个时间歩移动比较多，他们可能无法正确描述无穿透这件事。另外这类方法大多基于sequential quadratic programming优化方法，没法像无约束条件时的牛顿迭代那样方便地通过line search保证稳定收敛。

另一种思路是把物体外的空间也分割成四面体微元，作为另一些有限元“物体”来模拟，通过保证这些“物体”的每个微元不反转（体积恒正）来保证无穿透。这样做的问题是需要不断更新外界四面体网格的拓扑来防止被卡住，使得模拟非常昂贵，而且这种“变硬的空气”还可能阻碍弹性体本来应有的形变。

而我们观察到，如果准确计算物体表面每一个“点-三角”对和“边-边”对的距离，即在不同相对位置关系的情况下都准确找出距离最近的点对并求出距离（注意这个距离是非负的），那么如果初始状态下无重叠，只要保证所有这些距离在模拟的每一个时间歩 t 中的优化的每一次迭代 i 的更新轨迹上每一处都不为0，物体就不可能穿透！

写出我们所需的距离的表达式就是：

第 k 个“点 -三角 ”对：

第 k 个“边 -边 ”对：

这样“无穿透”三个字就可以在任意情况下准确严谨地表达为：

局部障碍函数

虽然我们的

和

其实都可以根据相对位置关系写成连续分段解析函数，我们对“无穿透”三个字的数学描述还是很昂贵的：a 可以取0到1之间无穷个实数，k 的个数则是物体表面点、边、三角形数量的二次方级别的！而IPC的解决方案，就是巧用连续碰撞检测（CCD）来对付a，并构造光滑的局部障碍函数来对 k 降维！

对于每次迭代更新，物体位置其实是从

线性渐变成

，那么只要我们在这两个状态间做CCD，找出第一次恰好产生穿透的那个 a，然后保证这一步迭代的更新步长不超过这个找出的 a 就好了，而无需对迭代更新过程进行更多的采样。（注意这里说的是迭代步长而不是时间步长。）这其实是一种line search filtering，常在内点法中使用。

配合上述CCD line search，我们把防止穿透的力定义为延距离函数梯度方向，随距离减小而增加，并在距离趋近于0时趋近于无穷，类似于用障碍函数法处理不等式约束条件。也就是说现在我们的问题转化成了

没有显式的约束条件，配合CCD就能精确近似原问题！这里 b 就是障碍函数，常用的比如log函数，它能为物体提供防穿透所需的任意大小的接触力；

是接触硬度，在IPC中通过数值分析自动调节。

传统的障碍函数当距离很远时力很小，对仿真结果的影响可忽略不计。但是接触元对的数量是随采样点二次方增加的，都处理的话仍然会非常昂贵。很自然地，我们会想去在计算中忽略那些距离远的接触元对，但是直接忽略会造成目标函数的不连续性，从而影响一切梯度下降方法的收敛性。IPC的解决方案则是，设计一个随距离增加光滑地下降为0并保持的函数：

就是IPC接触力产生的临界距离，一般设成毫米至微米级别。注意我们的函数 b 在

处也是C2连续的，画出来就是下图中的C2曲线：

IPC用来近似真实碰撞模型(discontinuous)的障碍函数

比起使用C0或者C1的截断函数，IPC的C2障碍函数更利于梯度下降方法的收敛，

设置的越小，IPC对真实碰撞模型的近似也就会越精确，需要在优化中计算的接触力和力贾科比矩阵的接触元对也越少，配合空间哈希快速找出需要考虑的接触元对，就能实现高效的仿真了！

摩擦力

现实中，物体间的接触难免会产生摩擦。传统的摩擦力通常根据库仑模型由NCP问题定义，也是个很难解的带约束条件的问题。为什么摩擦力要由约束条件来描述呢？简单来说，这是因为静摩擦力在两个有着完全相同物体位置与速度的场景中，还是可能由于外力的不同而不同，这就使得摩擦力系统不得不引入更多变量，并由约束条件来描述变量间的关系，就像原本的碰撞问题一样。

IPC则通过可控的平滑近似来描述静、动摩擦力与速度的关系，使得摩擦力的大小方向可以唯一由物体接触元对的切向相对速度决定，从而避免了更多变量以及约束条件的引入。区分静、动摩擦力的速度值越小，近似就越精确：（这个临界速度通常设为毫米每秒至微米每秒之间）

即使这样，很遗憾，摩擦系统还是无法写出Incremental Potential，因为它并不是保守力。IPC在对摩擦算子的时间离散中，巧妙的对法向碰撞力和摩擦切平面进行显式离散，将两者的误差保证在

内，从而构造出了近似的摩擦力势能，并融合到Incremental Potential中，实现了完整的Incremental Potential Contact：

这里

是所有上一个时间歩中的法向接触力不为0的接触元对，

和

是它们在上一个时间歩中的法向接触力和切平面算子，

是近似摩擦势能（

为近似摩擦力）。更多细节见原文。

精度验证

可能大家固有思维是我们计算机图形学里的有限元都是花架子，上面说了那么多“可控近似”、“精确近似”，IPC这样做的精度到底怎么样？有没有convergence under refinement？那么我们一起来看一看文章里的一个实验对比：

这段高尔夫球高速撞墙视频是我们在Youtube上找到的，我们在网上查阅了高尔夫球的材料参数以及这段视频中球的初速度并设置到场景中，使用能量守恒效果更好的Newmark时间积分，加上了一点lagged Rayleigh阻尼，就成功得到了这段仿真，视觉上完美重现高尔夫球撞墙后的形变，以及弹性波的变化。一定程度上验证了IPC碰撞处理的准确性。

计算机图形学的确不关注convergence under refinement，对于IPC来说，接触元对的距离、阻碍函数、以及接触硬度

都是resolution independent，所以refine后随着接触元对数量增加，结果肯定是不收敛的，这点我们在对IPC的后续研究中也发现了，并会有所改进。至于其它力，我们为了效率直接使用的是线性微元和lumped mass matrix，收敛性和能量守恒当然也是一般了，而且在大形变中肯定也有locking，这都是我们没有关注的点，而且后者和IPC关注的有限元碰撞也算是相对独立。

所以计算机图形学和机械、材料科学中有限元的区别一目了然，计算机图形学更关注稳定性和效率，机械、材料科学更关注数值精度，其实是各有千秋。我读过一些IJNME和CMAME的文章，它们为了精度，时间歩长即使是隐式积分也只取1e-4~1e-5秒左右，所以不会遇到大时间歩长牛顿迭代的不稳定性（由泰勒展开精度下降导致），所以也不会去研究或尝试在优化中加凸函数近似和line search。我们文章的三作，泰西欧.施耐德，就是机械工程背景，他第一次看到我仿真的neo-Hookean弹性材料能在这么大的时间歩长、这么大的形变下也保持稳定，直接惊呆了，你们可以发邮件和他求证。当然了，计算机图形学由于只关注视觉效果，所以基本也不会有动力去使用更高阶的有限微元，更不太会去关注isogeometric analysis这种高端的东西，因为贵啊！还【看】不出区别。但是呢，计算机图形学和机械、材料科学的有限元也是相互借鉴、互相促进的。它们就像科学、宗教与艺术，都有自己对这个世界的独特解构。

开源项目

IPC已在Github上全面开源，详细使用手册如下：

https://github.com/ipc-sim/IPC/wiki

只需简单几行的配置文件，就能与世界说你好：

  
  
    
   
   
     shapes input 2
input/tetMeshes/cube.msh 0 3 0  0 0 0  1 1 1
input/tetMeshes/cube.msh 0 1 0  0 0 0  1 1 1

selfFric 0.1

ground 0.1 0

还能方便地设置材料属性、边界条件、以及任意维度的预设物体：

  
  
    
   
   
     shapes input 2
input/tetMeshes/cube.msh 0 3 0  0 0 0  1 1 1 material 3000 1e8 0.4
input/tetMeshes/cube.msh 0 1 0  0 0 0  1 1 1

selfFric 0.1

ground 0.1 0

  
  
    
   
   
     shapes input 2
input/tetMeshes/cube.msh 0 3 0  0 0 0  1 1 1
input/tetMeshes/cube.msh 0 1 0  0 0 0  1 1 1  DBC -0.1 -0.1 -0.1  0.1 1.1 0.1  -0.2 0 -0.2  0 0 0  DBC 0.9 -0.1 0.9  1.1 1.1 1.1  0.2 0 0.2  0 0 0

selfFric 0.1

ground 0.1 0

  
  
    
   
   
     shapes input 2
input/tetMeshes/cube.msh 0 3 0  0 0 0  1 1 1
input/triMeshes/triangle.seg 0 1 0  0 0 0  2 2 2 angularVelocity 10 90 0

selfFric 0.1

ground 0.1 0

zoom 0.3

目前我们正在为IPC开发Houdini插件，以及进一步提升IPC代码的可读性与可扩展性，预计将会与一个新的开源物理仿真项目共同发布。同时我们也在研究如何让IPC更高效地适用于更广泛的应用场景。我们希望为广大研究人员、工程师、设计师、同学们、亲人们、朋友们提供可控、易用、可微分的物理仿真方法，让大家远离调参，直接使用真实物理参数得到预想结果，为这个世界创造更多美好！

如果一定要给最后一个视频加一个期限，我希望它是100s：