第三讲：刚体动力学模拟

作者: 离原春草 | 来源:发表于2022-04-09 12:04 被阅读0次

上一讲中介绍了单个block或者两个block的弹簧系统在受力情况下的运动过程，其中介绍了三种数值积分的方法：

显式欧拉
隐式欧拉
辛欧拉

不过前面的内容中物理系统的自由度较少，在实际情况中我们遇到的情况比这个复杂的多，因此今天这里来介绍刚体在受力情况下的运动过程，刚体指的是在任何受力情况下都不会发生形变的物体，这是一种理想情况，实际情况中并不存在刚体这种物体。

不过这里值得一提的是，将物体看成刚体的假设实际上也会引入新的问题，我们后面会看到，刚体的计算复杂度比软体更高，这也是为什么现有的刚体模拟中，没有哪种方法是真正完全准确的原因。另外，作为一般性，我们会发现，有时候我们在解决问题中的假设，试图将问题进行简化，但实际上并没有得到很好的效果，这里的刚体假设是一个案例，另一个案例就是后面介绍的流体动力学，在流体中，我们给出的一个假设是流体不可压缩，而这个假设同时也加剧了流体计算的复杂度。

用数学公式来表示刚体，可以给出如下：
$|x_i - x_j|_2 = FixedDistance$
即对物体上的任意两点，其欧式距离都是一个固定的常量。

在上面的情况下，我们可以用一个很长的向量来表示刚体：
$[x_1, x_2, ... , x_n]$

要计算这个很长的向量就会变得复杂，因此一个很直观的方法就是采用reduced coordinate（广义坐标）来表达，简单描述就是，由于任意两点之间的固定位置关系，因此可以用一些较少的自由度来完成对上述长向量的表达，后面我们会介绍，对于一个刚体来说，我们通常只需要六个自由度就能够完成表达。

而广义坐标方法的使用策略给出如下：

定义广义坐标
”immersion“映射，通过表达式将定义的广义坐标映射到full coordinate
根据表达式，就可以定义刚体的动能
根据上面的内容，就可以写出外力做功的表达式
最终输出一个与牛顿第二定律 $f = m \cdot a$ 等价的广义坐标方程

根据上面的广义坐标方程，就可以利用上一讲的内容利用数值积分完成求解，上面策略中前两步称为kinematics，后面三步称为dynamics。

从上面的策略，我们看到，这里是从能量角度而非力的角度（因为在实际的情况中，力的情况已经十分复杂且不清楚了）出发推导出一个广义坐标方程，而这就是经典物理中的一个重要里程碑拉格朗日力学的最重要的思想——第一次从能量的角度来看待物理世界。这样的角度的好处是，能量是跟坐标系无关的，所以推导跟建模的过程就无需考虑广义坐标的具体形式是什么样子。

下面会介绍两块内容：

kinematics在广义坐标系下是什么样的
时间不够，这里会跳过拉格朗日力学的具体推导过程，即一个经典的变分思想，直接给出牛顿第二定律在广义坐标下的形式

这一讲的重点在于外力，而游戏中的外力99%的情况下来自于外界的碰撞，因此这里会主要介绍如何解决碰撞。

1. 广义坐标使用策略细节

首先来定义一下计算的坐标系：

上图给出了一个简单的坐标系示意图，对于刚体上的任意一点，我们可以用一个3D向量 $\vec{x_i}$ 来表示，而在这个表示方法下，我们可以定义一个刚体的质心(center of mass)为：
$\vec{x_{cm}} \equiv \frac{\sum_i \vec{x_i} \cdot m_i}{\sum_i m_i}$

相当于对每一点的坐标进行了加权平均，而权重是当点质量的百分比。

这里来定义广义坐标系，而广义坐标系的原点，我们就选取在物体初始状态时的质心位置，即 $\vec{x_{cm}} = 0$ ，这是大家常用的做法，在后面的计算中按照这种做法会使得计算得到一定的简化。

刚体在初始状态受力之后，位置会发生变化，我们定义初始状态的刚体上某点位置叫 $\vec{r_i}$ ，而变化后的位置可以用如下公式计算得到：
$\vec{x_i} = R(\theta) \cdot \vec{r_i} + \vec{x_{cm}}$

其中 $R(\theta)$ 指的是围绕质心旋转某个角度 $\theta$ 的旋转矩阵，也就是说，整个刚体的运动可以拆解成平移与旋转，平移通过质心的坐标变化来表示，而旋转都是绕着质心完成的。

刚体上每个质点的动能为 $1/2m_i|\vec{v_i}|^2_2$ ，其中 $\vec{v_i}$ 为质点 $\vec{x_i}$ 的速度，而这个速度可以通过如下公式计算得到：
$\vec{v_i} = \vec{\dot{x_i}} = \dot{R}(\theta)r_i + \dot{\vec{x_{cm}}}$

可以看到，质点的速度可以由质心的速度（上面等式右边的求和中的右边一项）与质点的旋转速度表示。

在2D的情况下，旋转矩阵是一个2x2的矩阵：
$R(\theta) = \left[ \begin{matrix} cos \theta && -sin\theta \\ sin \theta && cos\theta \end{matrix} \right]$
导数就可给出为：
$\dot{R}(\theta) = \left[ \begin{matrix} -sin\theta && -cos \theta && \\ cos\theta && sin \theta \end{matrix} \right] \dot{\theta} = R(\frac{\pi}{2})R(\theta) \dot{\theta}$

3D坐标系下情况要复杂一点， $R$ 是一个三维旋转矩阵，属于Lie Group（李群），Lie Group表示的某一类矩阵，而Lie algebra（李代数）是另外一个代数结构，表示另外一类矩阵，这里用A表示。

我们知道，三维旋转可以由旋转轴与旋转角度表示，旋转轴可以表示为：
$\vec{\phi} = (x, y, z)$ ， $|\vec{\phi}| = 1$ ，旋转角度这里就还是用 $\theta$ 表示。根据旋转轴跟旋转角度，我们可以构建出一个Lie Algebra的矩阵：
$A(\theta) = \left[ \begin{matrix} 0 && -z\theta && y\theta \\ z \theta && 0 && -x\theta \\ -y \theta && x\theta && 0 \end{matrix} \right]$

这是一个斜对称矩阵（Skew-symmetric matrix），即其转置等于矩阵本身取负： $A^T = -A$

注意，这里的A不是旋转矩阵，但是跟我们的旋转矩阵有一个很好的映射关系： $R=e^A$ ，这里的矩阵指数计算跟标量指数的计算是一样的，不同的是矩阵指数是将指数应用到矩阵上的每个元素上。

假设某个质点在当前的旋转状态为 $R_t$ ，其旋转角速度为 $\vec{\omega}$ ，角速度可以用两个量进行表达，即旋转轴与角速度大小，令 $\vec{\omega} = (\omega_x, \omega_y, \omega_z)$ ，这个三维坐标可以看成是旋转坐标乘上旋转的角速度大小，那么经过一个短暂时间 $dt$ 的作用之后，其旋转的角度就变成了 $\omega dt$ ，其对应的Lie algebra就等于：
$R(dt) = \left[ \begin{matrix} 0 && -\omega_zdt && \omega_ydt \\ \omega_z dt && 0 && -\omega_xdt \\ -\omega_y dt && \omega_xdt && 0 \end{matrix} \right]$

那么在 $R_t$ 状态下，经过了时间dt之后的旋转状态为 $R_t \cdot R(dt)$ ，那么R_t的导数就可以表示为：
$\dot R(t)=\lim_{dt \rightarrow 0} \frac{R(t+dt) - R_t}{dt} =\lim_{dt \rightarrow 0} \frac{R_tR(dt) - R_t}{dt} =R_t\lim_{dt \rightarrow 0} \frac{R(dt) - I}{dt} =R_t\lim_{dt \rightarrow 0} \frac{e^{A(dt)} - I}{dt}$
对 $e^{A(dt)}$ 进行泰勒展开，即 $e^x \approx 1+x+O(x^2)$ ，省去后面的子项，上面公式可以变成：
$\dot R(t) \approx \frac{A(dt)}{dt} R_t$
将A(dt)代入，就得到了：
$\dot R(t) \approx \left[ \begin{matrix} 0 && -\omega_z && \omega_y \\ \omega_z && 0 && -\omega_x \\ -\omega_y && \omega_x && 0 \end{matrix} \right] R_t = \Omega^\star R_t$
而：
$\vec{\dot x_i }= \dot R(t) r_i + \vec{\dot{x}_{cm}} = \Omega^\star R_t r_i + \vec{\dot{x}_{cm}}$

另外，需要说一下的是， $\Omega \vec{v} = \omega \times \vec{v}$ 。

有了速度，我们就可以计算出动能：
$E_k = \frac{1}{2}\sum_i m_i|v_i|_2^2 = \frac{1}{2}\sum_i m_i v_i^T v_i = \frac{1}{2} \sum_i m_i(\dot{x}_{cm}^T\dot{x}_{cm} + r_i^T R^T_t (\Omega^\star)^T \Omega^\star R_t r_i + 2\dot x_{cm}^T\Omega^{\star}R_tr_i)$
而
$\sum_i m_i \dot x_{cm}^T\Omega^{\star}R_tr_i = \dot x_{cm}^T\Omega^{\star}R_t\sum_i m_i r_i$
又因为 $\sum_i m_i r_i$ 为初始质心坐标，即坐标原点，因此这个求和项结果为0，所以：
$E_k = \frac{1}{2} \sum_i m_i(\dot{x}_{cm}^T\dot{x}_{cm} + r_i^T R^T_t (\Omega^\star)^T \Omega^\star R_t r_i) = \frac{1}{2} \sum_i m_i \dot{x}_{cm}^T\dot{x}_{cm} + \frac{1}{2} \sum_i m_i r_i^T R^T_t (\Omega^\star)^T \Omega^\star R_t r_i = E_{k_t}+E_{k_r}$
$E_{k_t}$ 跟 $E_{k_r}$ 分别称为Translational Energy与Rotational Energy。用矩阵形式来表达：
$E_{k_t} = \dot x_{cm}^T \cdot M_{k_t} \cdot \dot x_{cm}$
这里 $M_{k_t}$ 是一个对角矩阵，矩阵中的每个元素都等于 $\sum_i m_i$

$E_{k_r} = \vec{\omega}^T \cdot M_{k_r} \cdot \vec{\omega}$

这里的 $M_{k_t}$ 被称为3D moment of inertia（通常使用inertia的首字母 $I$ 表示，inertia是惯性，也可以看成是另一种形式的质量，这里的inertia是旋转时的质量），这里会有一个简单的推导过程，时间关系，这里就直接省掉了。

回到微分方程：
$\frac{d}{dt} \left[ \begin{matrix} q \\ \dot q \end{matrix} \right] = \left[ \begin{matrix} \dot q \\\frac{f}{m} \end{matrix} \right]$

令动量 $P = m \dot q$ ，两边乘以m，我们有：
$\frac{d}{dt} \left[ \begin{matrix} q \\ P \end{matrix} \right] = \left[ \begin{matrix} \frac{P}{m} \\ f \end{matrix} \right]$

这是质点移动的微分方程，总共包含了6个自由度（ $q$ 跟 $\dot q$ 各三个），如果我们要为刚体建立一个类似的表达式，我们就需要将前面表示动能的所有自变量总共12个自由度都考虑进去：

质点位置 $\vec{x_{cm}}$ （后面考虑省略上面的箭头以简化书写）三个自由度
旋转轴 $\phi$ 也是三个自由度（虽然多了个 $\theta$ ，但是由于旋转轴方向长度为1，所以又减去了一个自由度，刚好还是3个）,
刚体线速度 $\vec{v}$ 三个自由度
刚体角速度 $\vec{\omega}$ 三个自由度）

同样，这里用动量来取代速度： $P = m v$ 为线动量， $L$ 为角动量，因此我们就可以得到：
$\frac{d}{dt} \left[ \begin{matrix} x_{cm} \\ \phi \\ P \\ L \end{matrix} \right] = \left[ \begin{matrix} \frac{P}{m} \\ \Omega^{\star} R_t \\ F \\ \tau \end{matrix} \right]$

这里来逐项解释一下：

第一项表示质心的位移的导数，可以理解为刚体的线速度，因此可以转化为动量除以质量，这里的质量是整个刚体的质量
第二项是角度的导数，对应的就是角速度，而角速度 $\dot R(t)$ 前面已经求解过，可以直接写入
第三项动量的导数，也就是 $mv$ 的导数，就是 $ma$ 也就是质心的受力
第四项角动量的导数，那就是质心的力矩了，力矩可以在网上搜到求解方法，这里就不做赘述。

这就是刚体下的牛顿第二定律。

2. 碰撞

游戏中最常见的物理模拟就是刚体的碰撞，碰撞主要有两个计算阶段：

碰撞检测，判断是否发生碰撞，以及碰撞发生的位置
受力计算，碰撞发生后会对产生碰撞的两个物体施加怎样的力和力矩

碰撞检测的方法在网上可以搜到大量资料，这里因为时间关系，这里就不做展开了，接下来的内容主要介绍如何计算受力，而这个问题到现在为止还没有一个完美的解决方案，现有的物理引擎提供了一系列的solver，允许用户进行参数调节，但是这些solver都有各自的缺陷。

两种求解方法：

Impulse-based方法：大概思路是通过冲量来完成对碰撞影响的模拟，简单来说，设碰撞前某点的速度为 $v_i^-$ ，碰撞后的速度为 $v_i^+$ ，在一个time step（通常是一个物理帧间隔） $t$ 之内，碰撞对二者产生了一个冲量 $f \cdot t$ ，经过这个冲量导致两者的速度发生了变化
Penalty-based方法，这是游戏引擎中常用的方法，其大致思路是在碰撞发生的时候为双方施加一个penalty，可以理解为在碰撞的时候，两者发生挤压，从而产生弹力，而弹力就会对应于弹性势能 $Q(h)$ 这里的h理解为两者碰撞的时候挤压的幅度（想象弹簧的形变幅度），而通过对势能求导我们可以得到对应的受力： $F =\nabla Q(h)$

先来看下Penalty-based方法，在Penalty-based方法中：

假设我们就用弹性势能作为碰撞时候产生的势能，那么我们就有： $Q(h) = \frac{1}{2}kh^2$ ， $F=kh$ ，通过这个受力来计算碰撞发生后两者的速度与位置，而这种方法的一个问题在于，如果发生碰撞时的速度过快，在给定的k值下，可能没有办法将两者弹开，而是可能会导致一个物体从另一个物体中穿过去的异常现象；

如上图所示，前面的弹性势能对应的曲线用红线绘制，可以看到势能增长的幅度与挤压的幅度的二次方程正比，而如果我们能够给出一种penalty function（称之为barrier function）使得势能在接近某个阈值的时候迅速逼近无穷大就能够有效抑制前面的问题。不过这种算法理想是美好的，现实却很残酷，比如在数值计算中没有办法处理这种无穷大的逻辑

在Imulse-based方法中，有 $v^+ = v^- + \alpha cdot n$ 其中 $n$ 是碰撞点的法线（方向朝外）， $\alpha$ 是当前的冲量。为了避免出现前面Penalty-based方法中物体穿过碰撞体的情况，这里通常会给碰撞施加一个约束，比如发生碰撞的时候的一个质点的位置为q，那么我们就有：
$g(q) \geq 0$
这里g(q)可以表示为当前点高于碰撞面往外方向的长度，举个最简单的例子，某个刚体与平面发生碰撞，我们可以有：
$(\vec q - \vec x) \cdot \vec n \geq 0$
上面公式中 $\vec x$ 表示的是碰撞面上的点， $\vec n$ 则是碰撞面上的法线。在Impulse-based方法中可以通过冲量很容易就能保证这种约束条件的达成。

那么如何求取这里的冲量呢（或者换句话说，如何计算输出的速度）？这里有很多模型来解决这个问题，最经典的做法是根据能量守恒定律，物体碰撞前跟碰撞后的能量保持不变：那么不考虑碰撞时的能量损耗，可以直接将输出速度的幅度设置为与输入速度一样，方向则可以根据反射的方式求得（当然，这是在碰撞没有考虑力矩作用的前提下的结论）；实际情况下，我们为了保证模拟结果的精度，通常会考虑碰撞过程中的能量损耗（如摩擦导致的热能损耗等），这种时候我们会添加一个叫做restitution coefficient的系数来控制能量保留比例，如果为0，则表示能量完全损耗。

上面介绍的是单点碰撞的情况，但我们有时候需要处理多点碰撞（multiple contact）的问题，这种问题就要复杂得多，这里就先不讲了，看看后面课程会不会有相关的内容。

第三讲：刚体动力学模拟

1. 广义坐标使用策略细节

2. 碰撞

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读