【质点弹簧】如何做一个不会崩溃的质点弹簧

【质点弹簧】如何做一个不会崩溃的质点弹簧

演示视频见：https://www.bilibili.com/video/BV15vPie8EEd/

在绳索、布料、软体等软性物质的模拟上，质点弹簧绝对是最流行的一种物理模型，相关资料在网上非常多。但无一例外的都绕不过一个痛点：动不动就崩溃给你看。那有没有一种能实现绝对不会崩溃的质点弹簧模型，或者说我们能始终确切的知道它崩溃的边缘在哪里，而不是和传统质点弹簧模型一样，总是在调参。

传统质点弹簧的三个缺陷

传统质点弹簧为什么那么容易崩溃？其实总结起来就那几个原因：

• 迭代时间没设置好（时间太短，导致累计误差大于弹簧收敛范围）
• 弹力系数没设置好（弹性过大，导致误差增幅快速超过了弹簧收敛范围）
• 弹簧阻力没设置好（阻力与弹力不平衡，无法限制过大的弹力，甚至火上浇油）

这三点其实就是对应了质点弹簧的三个参数。这三个参数不仅难调，而且其数值含义也比较模糊，产生的效果也只能大概描述，不精准预判它们的效果（例如你只知道弹力系数能让弹簧更硬，但不能精确说明其对质点间位置关系的影响，以及阈值在哪里）。

说到底，这三个参数从设计上就会引发崩溃，所以想要解决问题，就必须对这三个家伙优化。

去除迭代时间

迭代时间是真的可以从根本上解决，因为我们实际上有不依赖时间的位置积分方式。

确定质点位置的积分方式

总所周知质点弹簧模型里一般也就两种实现质点位置积分方式：

• 欧拉积分：$p(t+\mathrm{\Delta }t)=p(t)+v(t)\mathrm{\Delta }t+a(t)\mathrm{\Delta }{t}^2$（这里用常见的半隐式欧拉）
• Verlet 积分：$p(t+\mathrm{\Delta }t)=p(t)+(p(t)-p(t-\mathrm{\Delta }t))+a(t)\mathrm{\Delta }{t}^2$

可以观察到与位移有关项基本都依赖 $\mathrm{\Delta }t$，但很巧的是 Verlet 积分中的第二项 $p(t)-p(t-\mathrm{\Delta }t)$ 却是与时间无关。而由于加速度实际上可以转换为位移，所以我们可以使 $a(t)$永远等于 0，因此第三项我们是完全可以去除的。

划下来我们便能得到一种特别的位置积分方式，一种与时间完全无关的方法：

$$p(t+\mathrm{\Delta }t)=p(t)+(p(t)-p(t-\mathrm{\Delta }t))$$

这里的 $p(t)-p(t-\mathrm{\Delta }t)$ 其实对应的就是欧拉积分中的 $v(t)\mathrm{\Delta }t$，虽然具体数值和含义不同，但效果是一样的，因此后续我们将这一项也简称为“速度”。

确定力对质点的作用入口

该公式中提供了我们两个影响质点的入口：

• $p(t)$（当前位置）
• $p(t-\mathrm{\Delta }t)$（上次位置）

修改上次位置显然不太适合，因为这会使其含义与实际值不匹配，例如如果后续要做 ccd（连续碰撞检测），这个本应该可以用上的参数就完全废掉了。

那便只能修改当前位置了。从公式的效果上来看，这样的结果就相当于力使质点立即发生的位移，并将速度累计了下来（当前位置和上次位置差值变化了）。

（此外这还隐藏了一些额外的好处，我们后续再说。）

确定质点位置积分的时序

传统的积分方式是在质点当前的所有力都施加完毕后，即每帧的末尾再对质点的位置进行积分。但该积分流程如果用在我们的积分方案上就会存在问题，因为我们每次施加力是直接位移，然后累加下速度，接着积分通过速度再次位移，算下来一帧就位移了两次，即多了一次。

所以我们需要调整积分时序，在每帧的开头时进行积分，释放上一帧累加的速度，而后续的力计算，因为是直接作用于位移，因此也不依赖位置积分。

（此外这种迭代方式还隐藏了一些其他好处，我们后续再说。）

总结

新的质点弹簧模型，我们将使用如下方式迭代质点：

1. 每帧开始时对质点积分，积分方式采用：$p(t+\mathrm{\Delta }t)=p(t)+(p(t)-p(t-\mathrm{\Delta }t))$。
2. 对质点积分后再进行弹簧等力计算，结算力的方式为直接修改质点当前位置。
3. 完成帧迭代，后续直接用质点当前位置进行渲染等操作。

优化弹力系数

胡可定律表明：弹力=弹力系数*距离*方向。我相信该公式在现实世界的正确性，但在一个存在时间误差的模拟系统中，弹力系数并不是一个可以随便设置的值，而且由于上述“去除迭代时间”的操作，我们也无法直接使用力、速度（这里指欧拉积分里的常规速度概念）等与时间相关的参数。

弹力系数到底是什么？说到底，实现弹力现象的关键有两点：

1. 是一种力，所以会改变物体的速度。
2. 这种力是距离约束，却不会立即将物体拉到规定的距离上。

第一点我们已经实现，因为当前使用的积分方法中，修改位置就会累计速度。而第二点中，在“非规定距离”的选择上隐含了一些很重要的限制条件：

若当前质点距离弹簧约束的最佳距离位置为 $a$，那么质点下一次可移动到的距离只能在 $[0,2a]$ 间选择，因为超过这些距离后，质点每次迭代只会离目标越来越远，最终崩溃。考虑弹簧还会累计速度，若不存在阻力的话，那距离选择更是被限制在了 $[0,a]$。这一很关键的限制在传统的弹力系数上却没有体现，这就是为什么调大了就会崩溃。

恰好我们这里使用位移施加力，因此我们不再使用传统的弹力系数，而是改成一个 $[0,2]$ 的比例值，对应的就是上述的 $[0,2a]$。这样我们才能更清晰更稳定的控制弹簧弹力的大小。

处理弹簧阻力

很多传统质点弹簧除了弹簧本身的力之外，往往还会多一个步骤处理弹簧阻力，为什么？

观察实验现象很容易发现一点：

• 仅有一个弹簧时，不处理弹簧阻力也不容易崩溃。
• 但一旦多个弹簧共同作用，不施加阻力，就很容易崩溃了。

说到底原因就是多弹簧的弹力累加后使质点一次移动超过了 $2a$，最终弹簧发散崩溃，所以必须要额外的步骤实现阻力，来解决累计弹力过大的问题。

一种常见的阻力处理方式是检测在弹簧方向上的弹力的总和，然后用额外的阻力系数计算阻力。从原理上看这很类似与传统的弹力实现，只是这是基于力约束而不是距离约束。但也因此存在和传统弹力系数一样，阻力系数不够明确的问题。

• 最佳的阻力计算应该是什么样的？

  实际就是就是要考虑弹簧间的相互作用，确保在上一个弹簧处理后，再加上当前的弹簧力，始终不会使质点位移距离不超过 $2a$。

• 这种阻力有这种实现方式吗？

  有，另一种类似半隐式欧拉思想的阻力实现上，会始终利用质点在结算当前持有力后的位置来计算弹簧约束，从而确保了在与多个弹簧的相互作用后，质点相对当前弹簧的位移距离依然在约束范围内。但代价就是使原本每帧一次质点积分，变成了每次计算弹簧都要进行质点积分。

那有没有一种既能实现这种半隐式欧拉的弹簧阻力，同时又不用花费精力去不断积分的方法？很幸运的是，由于我们在积分方式上埋下的伏笔，我们弹簧的每一次力计算都是基于这种半隐式欧拉的：

1. 我们是先进行质点积分，再计算力，所以质点的当前位置始终是结算了当前速度的。
2. 我们的每一次力都是直接作用在位置上，所以质点的当前位置始终是考虑了力间相互作用的。

所以总结下来，若已按前两步“去除迭代时间”和“优化弹力系数”，那么这一步我们什么都不用做，因为我们的质点弹簧默认就已经考虑了阻力问题。

但实际上弹簧真的需要阻力吗？放在现实里，若一个物体同时被两个弹簧拉拽，显然它应该受到两个弹簧的合力，而不是一个弹簧力会变小。因此上述这种考虑力相互作用的弹簧是非真实的，结果就是这样的弹簧在摆球实验中，不能满足能量守恒定律。不过因此换来的稳定性确实是巨大的，相比之下，显然还是这种假弹簧更适用。

具体代码

下面是在 Unity 中基于上述质点弹簧模型实现的一个 Demo 的简化代码片段，你可以直观看到这种质点弹簧的实现方式，完整的代码见：https://www.cnblogs.com/BDFFZI/p/18732684

void Update()
{
    foreach (Transform point in allPoints)
    {
        Vector3 position = point.position;
        point.position += (position - lastPositions[point]);
        lastPositions[point] = position;
    }
    foreach (Spring spring in allSprings)
    {
        Transform pointA = spring.pointA;
        Transform pointB = spring.pointB;
        Vector3 positionA = pointA.position;
        Vector3 positionB = pointB.position;

        //胡克定律：弹力=弹力系数*距离*方向
        Vector3 vector = positionA - positionB;
        float distance = vector.magnitude - spring.length; //距离
        Vector3 direction = vector.normalized; //方向
        Vector3 move = elasticity * distance * direction;//弹力系数为[0,2]，无阻力时应小于1

        pointB.position += 0.5f * move;
        pointA.position += 0.5f * -move;
    }
}