第四讲：弹性体模拟-弹性力学基础

Posted 2023-04-21

参考技术A 这一讲的主题是弹性体模拟中的弹性力学基础，首先来看一个问题，就是当一个球跟多个球同时发生碰撞，如下图所示（称之为Bernoulli Problem）：

或者当第一个球跟第二个球碰撞的瞬间，第二个球同时与第三个球碰撞，如下图所示（称之为Newton's Cradle）：

这些问题的求解都可以归结为线性互补问题Linear Complementary Problem（简称LCP）的求解，线性互补问题的定义可以参考 物理引擎之约束求解(一)——线性互补问题 。假设有两个刚体（分别标注为1号跟2号），由于每个刚体都有六个自由度（translation+rotation），那么这两个刚体的位置就可以用下述公式来表示：

其中

同样 也有对应的变量表达，而由于刚体是不可穿插形变的，因此q还会有一些额外的约束：

这里的i指的是第i个碰撞点， 可以表示快要碰撞时候两个物体在第i个碰撞点处的signed distance。

如果用上一讲中的impulse-based方法来解方程，那么就意味着我们需要在碰撞的瞬间为两者施加一个冲量，而要想产生冲量，也就意味着两者产生了穿插，前面的条件将不再满足，即：

而如果我们希望在下一个timestep回复到正确的状态，那么上面这个函数的导数就需要是为正的：

这个表达式的物理意义是，等式后面表达式中的前面部分 表达的是碰撞点处的表面的法线（signed distance的梯度就是signed distance下降最快的方向，即表面法线方向），后面部分是刚体的速度，两者相乘大于等于0表示的是这两个向量的夹角小于九十度，也就是说在这个速度下，会使得两者的signed distance不断变大。

在impulse的作用下，速度会从 变成 ，即速度从碰撞变成背离。在Houdini中有个叫做restitution coefficient（恢复系数）的参数 ，这个参数是做什么的呢？用一个例子来说明，一个球在重力的作用下下坠，碰到地表，如果这个参数为0，那么这个小球在碰撞后，就会贴在地表上，如果是1的话，碰撞之后就会回弹到起始的高度，也就是说，这个参数表达的是经过碰撞之后能量的残留比例。

将恢复系数跟表面法线结合起来看，我们可以得到如下公式（这个条件我们称之为一号条件）：

如下图所示：

入射方向为 ，出射方向为 ， 可以简单看成是是出射方向的速度的长度与入射方向的速度长度的比例， 是法线方向，那么上面这个公式自然就是成立的，毕竟能量是守恒的，出射速度的长度不可能超过入射速度的长度。

上面这个公式是针对单个碰撞点（contact point）而言的，当我们有多个碰撞点的时候， 就变成了一个矩阵，这个矩阵我们称之为 （这个矩阵的物理含义为，将刚体的速度映射到碰撞点上的沿着法线方向的速度？不是特别理解，如果单个碰撞点的情况，对应的是碰撞点的法线，那么这里就对应的是多个碰撞点法线组成的矩阵，既然是多个碰撞点，也就没有一个整体的法线概念了），这里的转置后面会解释。

另外，这里还有另外一个公式：

这个公式中 表示的是刚体的质量， 是一个矩阵，这个矩阵就是上面梯度向量组成的矩阵的转置（推导这里就不说了，抛开多个碰撞的情况不说，单个碰撞点的情况下，冲量的作用方向肯定是垂直于碰撞表面的，也就是说经过冲量作用后的加速度应该也是与碰撞表面的法线方向保持一致的，从主观上来说，这个结果是可以理解的）， 则是前面说的碰撞时的冲量impulse，如果有多个碰撞点的话，这个变量也就是一个向量（几个碰撞点就是几维的）。矩阵 跟冲量的乘法得到的是一个力，而作用力除以质量就得到了加速度，通过这个加速度（在时间的作用下）就能完成速度的一个变化。

这个公式中，冲量（每个分量）必须要大于等于0（当刚体碰撞时的速度在表面法线方向上的速度非0，且速度法线方向的分量与法线方向相反，那么就取＞0，否则就取等于0），否则在碰撞后，刚体就直接装进另一个刚体中，从而发生穿插，这个条件我们称之为二号条件：

一号条件跟二号条件组成linear complementary condition（即如果 （对应的是碰撞点碰撞之后没有施加冲量，也就是说碰撞时候的速度与表面法线正好垂直，这种情况下，表面法线与速度正交，一号条件不等式前面的结果就是0）的时候，就不需要考虑一号条件；而如果 就必须要考虑一号条件，我们将这种情况称之为互补），最终我们需要求解的问题就可以表示为，在一二号条件：

约束下，求解：

即LCP，这个问题如果我们能够求解出冲量向量，那么我们就能得到需要的解，但是在多个碰撞点的情况下，这个会比较复杂。

Bullet物理引擎用的是一种叫做Gauss Siedal Method，一个个碰撞点去考虑，比如先考虑一号碰撞点的impulse是多少，之后再看一号impulse对二号碰撞点的速度的改变，之后再看二号碰撞点的impulse为多少并加上一号碰撞点的影响，同时算出对三号碰撞点的影响，以此类推，直到最后一个碰撞点对一号碰撞点的影响，不断迭代，经历多轮迭代最终达到一个平衡，这种方法是串行计算，对GPU不友好，性能较差。

另一种方法则是将LCP（这是一种timestep的方法）转换成一个优化问题，这种思想是物理模拟中的一种十分重要的思想，在很多问题的求解上都有应用。

先来说说拉格朗日乘子（ lagrange multiplier ），这是数学上的一种优化策略中的术语，这个优化策略用以在给定的一些等式约束下，求得某个函数的局部极值（极大极小值），这个算法的基本思想可以叙述为，对于一个需要求解极小值（极大值可以转换为极小值）的函数f(x)，已知需要满足g(x) = 0条件（或者说，已知极值是在g(x) = 0条件下满足），那么拉格朗日函数可以表示为：

我们知道，如果不考虑条件函数g(x) = 0的话，f(x)的极值可以直接通过 来求得，而在加上条件之后，问题就会变得复杂一点，上面的拉格朗日函数是自变量x跟 （这个就叫拉格朗日乘子）的函数，而原函数f(x)的极点则肯定会出现在拉格朗日函数的saddle point（极点）上，所以只需要对拉格朗日函数求偏导，并令各个偏导结果为0进行求解，就能得到取得极值的坐标点。

而如果对上面的情况进行泛化处理，比如将条件从等式变成不等式，那么我们就需要使用 KKT（Karush Kuhn Tucker） 条件方法，同样的，KKT会将原始的函数f(x)，与等式约束 以及不等式约束 写到一个式子里：

f(x)取得极值的条件为：

联合这些等式，就能得到最终的极值点的坐标。这里值得一提的是，上面第三个条件中，由于h(x)是小于等于0的，因此要想条件成立，要么h(x) = 0，要么 ，也就是说，这个条件可以拆成两个条件：

且这两个条件是complementary（互补）的，如果移除g(x)的干扰，这个条件就跟前面LCP问题的求解的形式完全一致了，也就是说，LCP的求解就变成了KKT（优化问题）的求解，而这种问题的求解方法就十分丰富了，这里就不做展开，否则可以直接讲一个学期。

但是我们发现，用LCP来求解Bernoulli Problem可以得到正确的结果，但是用来求解Newton Cradle问题的时候结果却是错误的，这是因为理论上LCP问题的求解是不唯一的，也就是说我们得到的解有多个，其中选取的那个可能并不是正确的解。

这里需要注意的是，上面的求解是没有考虑碰撞时的摩擦力的，实际上如果添加了摩擦力的话，问题会变得更为复杂。但是有意思的是，即使加上摩擦力，依然可以表达为LCP问题（不过更为复杂），转化成的优化问题，其限制条件将会变成二次函数，因此其求解也会变得更为复杂。

最后做个总结，在目前情况下的物理引擎，对刚体碰撞的模拟依然是存在较大的精度问题，比如说大多数都没有办法解决上面的Bernoulli跟Newton Cradle问题，即使是不考虑摩擦力的情况，目前都没有物理引擎能给出完善的解决方案，更何况是更为精确的考虑摩擦力的情况。

Continuum Mechanics中描述了形变的相关理论，包括弹性形变与塑性形变，我们这里着重考察弹性形变。

最简单的弹性形变就是前面说的弹簧+质点的粒子系统，当质点不断移动从而压缩弹簧的话，弹簧就会不断的积累势能，当手松开之后，势能可能就会转换为动能，此时弹簧产生的力叫做conservative force。

但是弹性形变不只是局限于弹簧系统，普通的物件也会发生此类形变，但是这里有亮点需要明确：

可以知道，弹性体的形变是一个连续的变化，即形变前相邻的点在形变后依然是相邻的。我们可以利用微分思想来对这个连续的规律进行考虑，假设形变前物体上某一点X（处在material space）经过形变后变成了点x（处在deformed space），那么形变就可以表示成如下的公式：

对这个映射关系采用泰勒展开：

由于 是一个三维的变量（三维空间中的点），因此上面公式中的偏微分项 就是一个3x3的矩阵，这个矩阵我们叫做deformation gradient（形变梯度矩阵），通常用F来表示。

从这个公式我们可知道，这个矩阵的第一列进行乘法运算时对应的是 的第一个元素，换句话说，第一列就是在material space中沿着x方向移动一个距离时，在deformed space中就需要沿着某个方向移动一个相同大小的距离，这个方向就用第一列向量来表示，同理，其他两列就分别对应y/z方向上移动时对应的deformed space中的移动方向。

接下来我们看看要如何用数学公式表示物体的弹性形变长度，假设在material space中某个点移动的位移为 ，那么移动的长度的平方就可以表示为 ，而在deformed space中对应点的移动位移就可以用 来表示，其长度则为 ，那么具体的形变量就可以用后者开平方后减去前者的开平方，不过这里为了避免因为开平方导致的高额计算消耗，可以先直接用平方之差来表示形变的测度：

上面等式中最后的 称之为Green（格林） Strain（形变） Tensor（本质上是个矩阵），这个矩阵可以用来描述微分意义下物体形变的规律，可对于推导或者计算出这个 的点X而言，在其周边较小的一块区域（微分算子）里的所有点都可以通过这个矩阵来计算其形变量，但是，如果超出这个区域范围， 可能会不同，也就是说， 准确来说是物件上某点的一个函数，可以写成

接下来，就是如何将前面的 变成势能，势能可以表示为如下形式的函数：

我们知道，这个函数有如下的性质：

根据上面两个性质，我们可以推断，这个函数可以大致用下图来表示：

至少是一个两阶以上的模型，应用泰勒展开，可以得到：

而由于 跟 都是0，那么上面公式就变成了：

由于 是个矩阵，那么这个公式就可以改写成如下的形式，这个公式我们称之为constitution model：

到目前为止，所有的公式都是通过数学推导给出的，没有掺杂物理相关的概念与理论，接下来，为了能够得到上面的公式的具体形式，就需要给出系数 ，而物理学中针对这个系数提出了众多的模型，甚至还有人考虑通过神经网络训练得到这个系数，这就造成了不同的势能物理模型。

要知道这个公式是如何在物理中进行应用的，我们就放到下一讲中进行介绍了。

AEJoy —— 表达式之弹性（韧性）模拟详解JS

效果图

弹性（elastic）

表达式代码与注释

var p = 0.6; ///< 弹性期限
var a = 140; ///< 弹性振幅

/// @note 弹性函数
function outElastic(t, b, c