这个系列是GAMES103：基于物理的计算机动画入门(GAMES 103: Intro to Physics Based Animation)的同步课程笔记。本课程是计算机动画的入门课程，着重介绍各种基于物理的动画仿真模拟技术。本节课主要介绍刚体动力学的相关内容。
Rigid Body
刚体(rigid body)是指不会发生形变的物体。尽管现实世界中并不存在理想刚体，但只要可以忽略物体的变形我们都可以将其视为是理想刚体。
同时在游戏以及动画中刚体模型也都拥有非常多的应用：
刚体的运动包括平移(translation)和旋转(rotation)两部分，在进行仿真时需要计算每个时刻刚体的位姿。
由于刚体不存在变形，只要我们知道刚体的位姿就可以计算出刚体上任意点的位置。
Translational Motion
Integration Methods
我们首先来考虑刚体的平移，根据牛顿定律刚体的位置可以通过在时间上积分来计算：
\[\mathbf{v} (t^{[1]}) = \mathbf{v} (t^{[0]}) + \mathbf{M}^{-1} \int_{t^{[0]}}^{t^{[1]}} \mathbf{f} (\mathbf{x}(t), \mathbf{v}(t), t) dt\] \[\mathbf{x} (t^{[1]}) = \mathbf{x} (t^{[1]}) + \int_{t^{[0]}}^{t^{[1]}} \mathbf{v} (t) dt\]
积分项通常是通过数值积分来计算，计算方法包括显式欧拉(explicit Euler)、隐式欧拉(implicit Euler)、中点法(mid-point)等：
- Explicit Euler: \(\int_{t^{[0]}}^{t^{[1]}} \mathbf{v} (t) dt \approx \Delta t \mathbf{v} (t^{[0]})\)
- Implicit Euler: \(\int_{t^{[0]}}^{t^{[1]}} \mathbf{v} (t) dt \approx \Delta t \mathbf{v} (t^{[1]})\)
- Mid-point: \(\int_{t^{[0]}}^{t^{[1]}} \mathbf{v} (t) dt \approx \Delta t \mathbf{v} (t^{[0.5]})\)
直观来说，显式欧拉是使用当前时刻的速度进行积分，隐式欧拉是使用下一时刻的速度进行积分，而中点法则是使用中间时刻的速度进行积分：
此外我们还可以将显式和隐式混合起来，在计算速度时使用显式欧拉法而在计算位置时使用隐式欧拉法：
\[\mathbf{v}^{[1]} = \mathbf{v}^{[0]} + \Delta t \mathbf{M}^{-1} \mathbf{f}^{[0]}\] \[\mathbf{x}^{[1]} = \mathbf{v}^{[0]} + \Delta t \mathbf{v}^{[1]}\]这种方法一般称为半隐式(semi-implicit)方法，也称为leapfrog method。
Types of Forces
在进行仿真时我们还需要考虑物体的受力，一般情况下我们需要考虑重力：
\[\mathbf{f}_\text{gravity}^{[0]} = \mathbf{M} \mathbf{g}\]同时当物体进行自由落体时还需要考虑空气阻力：
\[\mathbf{f}_\text{drag}^{[0]} = -\sigma \mathbf{v}^{[0]}\]式中\(\sigma\)为阻力系数。除了直接计算空气阻力外，实际中有时会使用速度衰减来近似阻力：
\[\mathbf{v}^{[1]} = \alpha \mathbf{v}^{[0]}\]式中\(\alpha\)为衰减系数。我们把力的计算与刚体位置的更新结合到一起就得到了刚体平移的仿真流程：
Rotational Motion
Representation
接下来讨论刚体的旋转。我们知道旋转可以用一个正交矩阵来表示，因此最基本的旋转表示方法是使用旋转矩阵：
\[\mathbf{R} = \begin{bmatrix} r_{00} & r_{01} & r_{02} \\ r_{10} & r_{11} & r_{12} \\ r_{20} & r_{21} & r_{22} \\ \end{bmatrix}\]但使用旋转矩阵来进行表示有一些缺点：比如说旋转矩阵过于冗余了，它需要9个参数才能表示旋转但实际上我们只需要3个参数就可以了；同时旋转矩阵非常不直观，对它进行求导也很不方便。
另一种表示旋转的常用方法是使用欧拉角(Euler angle)，它将旋转分解为绕三个轴的转动因此只需要3个参数就能表示旋转。欧拉角的优点在于它非常直观，很容易让人理解旋转的过程，因此在Unity中就使用了欧拉角来表示刚体的旋转。但需要注意的是欧拉角存在万向锁(gimbal lock)的问题，在某些情况下按欧拉角进行旋转会导致自由度的损失。同时对欧拉角进行求导也是比较困难的，因此在刚体仿真中也不太常用欧拉角来表示旋转。
在刚体仿真中较为常用的旋转表示方法是四元数(quaternion)。四元数可以看成是复数的推广，它具有一个实部和\(\mathbf{i}\)、\(\mathbf{j}\)、\(\mathbf{k}\)三个虚部。虚部之间的乘法法则如下：
我们可以使用一个向量来表示四元数\(\mathbf{q} = [s, \mathbf{v}]\)，四元数之间的运算法则包括：
- 数乘：\(a \mathbf{q} = [as, a\mathbf{v}]\)
- 加减法：\(\mathbf{q}_1 \pm \mathbf{q}_2 = [s_1 \pm s_2, \mathbf{v}_1 \pm \mathbf{v}_2]\)
- 乘法：\(\mathbf{q}_1 \times \mathbf{q}_2 = [s_1 s_2 - \mathbf{v}_1 \cdot \mathbf{v}_2, s_1 \mathbf{v}_2 + s_2 \mathbf{v}_1 + \mathbf{v}_1 \times \mathbf{v}_2]\)
- 模：\(\Vert \mathbf{q} \Vert = \sqrt{s^2 + \mathbf{v} \cdot \mathbf{v}}\)
假设物体绕\(\mathbf{v}\)轴旋转了\(\theta\)度，则旋转的过程可以使用一个单位四元数(unit quaternion)来表示：
使用单位四元数对矢量进行旋转可以表示为：
\[[0, \mathbf{r}'] = \mathbf{q} \ [0, \mathbf{r}] \ \mathbf{q}^{-1}\]其中\(\mathbf{r}\)为位置矢量；\(\mathbf{q}^{-1} = [s, -\mathbf{v}]\)为单位四元数\(\mathbf{q}\)的逆。我们还可以把上式的四元数乘法改成矩阵乘法，这样得可以到四元数与旋转矩阵的转换公式：
\[\mathbf{R} = \begin{bmatrix} s^2 + x^2 - y^2 - z^2 & 2(xy - sz) & 2(xz + sy) \\ 2(xy + sz) & s^2 - x^2 + y^2 - z^2 & 2(yz + sx) \\ 2(xz - sy) & 2(yz + sx) & s^2 - x^2 - y^2 + z^2 \\ \end{bmatrix}\]Torque and Inertia
类似于质量与加速度，在计算刚体旋转时我们有力矩(torque)和转动惯量(inertia)的概念。力矩是刚体角速度发生变化的原因，它定义为位置矢量与力的叉乘：
\[\boldsymbol{\tau}_i = (\mathbf{R} \mathbf{r}_i) \times \mathbf{f}_i\]
转动惯量则类似于平移运动中的质量，它表示刚体抵抗力矩的能力。但需要注意的是转动惯量不仅与刚体的质量分布有关，也与刚体旋转的轴有关。因此转动惯量不是一个常数，而是一个张量(矩阵)。通常情况下我们把刚体初始状态的转动惯量记为\(\mathbf{I}_\text{ref}\)，计算公式为：
\[\begin{aligned} \mathbf{I}_\text{ref} &= \begin{bmatrix} \int_V (\mathbf{x}_y^2 + \mathbf{x}_z^2) \rho(\mathbf{x}) dV & -\int_V (\mathbf{x}_x \mathbf{x}_y) \rho(\mathbf{x}) dV & -\int_V (\mathbf{x}_x \mathbf{x}_z) \rho(\mathbf{x}) dV \\ -\int_V (\mathbf{x}_x \mathbf{x}_y) \rho(\mathbf{x}) dV & \int_V (\mathbf{x}_x^2 + \mathbf{x}_z^2) \rho(\mathbf{x}) dV & -\int_V (\mathbf{x}_y \mathbf{x}_z) \rho(\mathbf{x}) dV \\ -\int_V (\mathbf{x}_x \mathbf{x}_z) \rho(\mathbf{x}) dV & -\int_V (\mathbf{x}_y \mathbf{x}_z) \rho(\mathbf{x}) dV & \int_V (\mathbf{x}_x^2 + \mathbf{x}_y^2) \rho(\mathbf{x}) dV \\ \end{bmatrix} \\ &= \begin{bmatrix} \sum_i (y_i^2 + z_i^2) m_i & -\sum_i (x_i y_i) m_i & -\sum_i (x_i z_i) m_i \\ -\sum_i (x_i y_i) m_i & \sum_i (x_i^2 + z_i^2) m_i & -\sum_i (y_i z_i) m_i \\ -\sum_i (x_i z_i) m_i & -\sum_i (y_i z_i) m_i & \sum_i (x_i^2 + y_i^2) m_i \\ \end{bmatrix} \\ &= \sum_i m_i (\mathbf{r}_i^T \mathbf{r}_i - \mathbf{r}_i \mathbf{r}_i^T) \end{aligned}\]
同时需要注意的是在每一时刻我们都需要对刚体当前的转动惯量进行更新：
\[\begin{aligned} \mathbf{I} &= \sum_i m_i (\mathbf{r}_i^T \mathbf{R}^T \mathbf{R} \mathbf{r}_i - \mathbf{R} \mathbf{r}_i \mathbf{r}_i^T \mathbf{R}^T) \\ &= \sum_i m_i (\mathbf{r}_i^T \mathbf{r}_i - \mathbf{R} \mathbf{r}_i \mathbf{r}_i^T \mathbf{R}^T) \\ &= \sum_i m_i (\mathbf{R} \mathbf{r}_i^T \mathbf{r}_i \mathbf{R}^T - \mathbf{R} \mathbf{r}_i \mathbf{r}_i^T \mathbf{R}^T) \\ &= \mathbf{R} \mathbf{I}_\text{ref} \mathbf{R}^T \end{aligned}\]
这样就可以得到计算刚体的旋转时在时间上进行积分的公式：
\[\boldsymbol{\omega}^{[1]} = \boldsymbol{\omega}^{[0]} + \Delta t (\mathbf{I}^{[0]})^{-1} \boldsymbol{\tau}^{[0]}\] \[\mathbf{q}^{[1]} = \mathbf{q}^{[0]} + [0, \ \frac{\Delta t}{2} \boldsymbol{\omega}^{[1]}] \times \mathbf{q}^{[0]}\]
Rigid Body Simulation
最后我们把平移和旋转的更新结合到一起就得到了刚体运动的仿真流程：
Reference
- Lecture 3: Rigid Body Dynamics
- Foundations of Physically Based Modeling and Animation, Chapter 9: Rigid Body Dynamics
- Foundations of Physically Based Modeling and Animation, Appendix E: Quaternions