基于物理的动画:刚体动力学

本文是博主在实习公司制作的课程讲义,内容是基于物理的动画:刚体动力学。

原课程对应的幻灯片可以点击这里下载

动画与物理仿真

概述

动画流水线 - 基于数据的动画

image
image

动画流水线 - 实时物理仿真

image
image

每帧的动画计算

image
image
  • 动画是每帧对场景物体状态,包括视觉状态和非视觉的状态的更新
  • 常见的视觉状态包括

    • 位置
    • 朝向(旋转)
    • 外观
  • 常见的非视觉的物理状态包括

    • 外力
    • 速度
    • 密度
    • ...
  • 离线动画应用仅需存储视觉状态,物理仿真算法则按需存储其他物理状态

物理仿真流水线

\(t^{(k)}\) 为第 \(k\) 帧的时刻,\(\mathbf s^{(k)}\) 为仿真对象在第 \(k\) 帧时的状态。

仿真器在每一帧执行物理演算,并将结果提交给渲染器。

image
image

刚体动力学

预备知识

刚体的性质

image
image
  • 刚体是不会产生形状、大小改变的理想物体。只有整体的平动转动
  • 一般地,物理引擎存储的刚体仿真状态\(\mathbf s = \{\mathbf v, \mathbf x, \mathbf \omega, \mathbf q\}\)是一个四元组,其中 \(\mathbf v\) 为刚体的线速度,\(\mathbf x\)为位移,\(\mathbf \omega\)为角速度,\(\mathbf q\)为表示朝向的四元数 。
  • \(n\)为刚体的顶点数。如果刚体的各顶点质量不一,记顶点\(i\)的质量为\(m_i\),则\(M = (m_1, \dots, m_n)^T\),记\(M^{-1} = (1 / m_1, \dots, 1 / m_n) \in \mathbb R^{1 \times n}\)。此时我们需要让\(\mathbf f = (\mathbf f_1, \dots, \mathbf f_n)^T\)分别储存各顶点的受力情况。

刚体的运动

刚体整体不缩放,只发生平动\(\mathbf x\)(translation)和转动\(\mathbf R\)(rotation)变换。

image
image

平动动力学

平动动力学

image
image

根据牛顿运动定律:

\[ \begin{cases} \mathbf v(t^{(1)}) = \mathbf v(t^{(0)}) + M^{-1} \int_{t^{(0)}}^{t^{(1)}} \mathbf f(\mathbf x, \mathbf v, t) \, \mathrm d t \\ \mathbf x(t^{(1)}) = \mathbf x(t^{(0)}) + \int_{t^{(0)}}^{t^{(1)}} \mathbf v(t) \, \mathrm d t \end{cases} \]

进行离散半隐式(semi-implicit)积分得到

\[ \begin{cases} \mathbf v^{(1)} = \mathbf v^{(0)} + \Delta t M^{-1} \mathbf f^{(0)}, \\ \mathbf x^{(1)} = \mathbf x^{(0)} + \Delta t \mathbf v^{(1)} \end{cases} \]

其中\(\mathbf v\)用上一帧的\(\mathbf f\)数据更新,是显式积分,\(\mathbf x\)用当前帧的\(\mathbf v\)数据更新,是隐式积分。可以证明,如此更新的\(\mathbf x\)相对完全显式法和完全隐式法是二阶精确的。

刚体仿真-仅平动

仿真流程:

  1. 输入状态\(\mathbf s^{(0)} = \{\mathbf v^{(0)}, \mathbf x^{(0)}\}\)
  2. \[ \begin{cases} \mathbf f_{i}^{(0)} = \text{Force}(\mathbf x_i^{(0)}, \mathbf v_i^{(0)}, \dots), \\ \mathbf f^{(0)} = (\mathbf f_1^{(0)}, \dots, \mathbf f_n^{(0)})^T, \\ \mathbf v^{(1)} = \mathbf v^{(0)} + \Delta t M^{-1} \mathbf f^{(0)}, \\ \mathbf x^{(1)} = \mathbf x^{(0)} + \Delta t \mathbf v^{(1)} \end{cases} \]
  3. 输出状态\(\mathbf s^{(1)} = \{\mathbf v^{(1)}, \mathbf x^{(1)}\}\)

注意:

  • 各顶点质量\(M = (m_1, \dots, m_n)^T\)为常向量,表记\(M^{-1} = (1/m_1, \dots, 1/m_n) \in \mathbb R^{1 \times n}\)
  • 时间间隔\(\Delta t\)可自定义,不一定是常量,也不一定与实际物理帧计算的时间间隔相同

转动动力学

旋转的表达 - 旋转矩阵与欧拉角

在讨论三维旋转时,我们只讨论绕某个过原点的旋转轴旋转。若旋转轴不过原点,需要先进行平移变换。

旋转矩阵:

  • 三维空间中的旋转可用一个自由度为3的\(3 \times 3\)矩阵\(\mathbf R\)表示
  • 对空间向量\(\mathbf v \in \mathbb R^{3}\),矩阵-向量乘积\(\mathbf R \mathbf v\)可表示对该向量施加旋转变换
  • 缺点:不够直观;自由度低,浪费空间;难以计算旋转相关的速度

欧拉角:

  • 三维空间中的旋转亦可用直观的三元组------欧拉角(euler angle)表示
  • 给定空间中三个固定正交轴(局部的或全局的),欧拉角是按固定顺序分别绕这些轴的旋转角度的记录。不同的欧拉角三元组,最终可能产生相同的旋转结果
  • 缺点:可能丢失自由度(万向节死锁);难以计算旋转相关的速度

旋转的表达 - 四元数

四元数\(\mathbf q = a + b \textrm i + c \textrm j + d \textrm k \in \mathbb H \leftrightarrow \mathbb R^4\)是对复数域\(\mathbb C \leftrightarrow \mathbb R^2\)的扩展,遵从运算法则

\[ \textrm i^2 = \textrm j^2 = \textrm k^2 = \textrm i \textrm j \textrm k = -1, \]

我们可以表记四元数\(\mathbf q = [s, \mathbf v]\),其中\(s\)\(\mathbf q\)的标量部分,\(\mathbf v\)为向量部分。可以推导四元数的标量积、加法和乘法运算:

\[ \begin{cases} a \mathbf q = [as, a \mathbf v], \\ \mathbf q_1 + \mathbf q_2 = [s_1 + s_2, \mathbf v_1 + \mathbf v_2], \\ \mathbf q_1 \times \mathbf q_2 = \mathbf q_1 \mathbf q_2 = [s_1s_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]. \end{cases} \]

特别地,我们可以让三维向量\(\mathbf v\)也参与上述计算,只需视其为四元数\([0, \mathbf v]\)即可。

旋转的表达 - 四元数运算

正如一个归一化的复数可以表示二维旋转一样,一个归一化的四元数也可以表示一个三维旋转,它满足

\[ ||\mathbf q|| = \sqrt{s^2 + \mathbf v \cdot \mathbf v} = 1. \]

具体地,我们让向量\(\mathbf v\)绕某个过原点的空间轴\(\mathbf u\)旋转\(\theta\)角度,旋转方向使用右手螺旋判断。记

\[ \begin{cases} \mathbf q = \left[\cos \frac{\theta}{2}, \mathbf u \right], \\ ||\mathbf q|| = 1, \end{cases} \]

那么,可以证明,\(\mathbf q \mathbf v \mathbf q^{*}\)就是空间向量\(\mathbf v\)如上旋转后的结果的四元数表示,即\(\mathbf q \mathbf v \mathbf q^{*}\)的标量部分总是为\(0\);其中\(\mathbf q^{*} = \left [\cos \frac{\theta}{2}, -\mathbf u \right]\)\(\mathbf q\)的共轭四元数。

归一化的四元数可以与三维旋转形成一一对应。四元数也可以方便地转换为矩阵和欧拉角,因此成为了物理引擎内部旋转运算的常用数据类型。

转动动力学

image image

\(\mathbf q\)为表示刚体旋转的四元数,\(\mathbf R\)表示该旋转的变换矩阵,\(\mathbf \omega \in \mathbb R^3\)为旋转的角速度,规定:

  • \(\mathbf \omega\)的方向为旋转轴(与旋转方向形成右手螺旋关系),
  • \(\mathbf \omega\)的大小为旋转的角速率。

力矩

转动动力学中的力矩(torque)是平动动力学中的力(force)的等效,定义式为

\[ \mathbf \tau_i = (\mathbf R \mathbf r_i) \times \mathbf f_i. \]

image
image

转动惯量

转动动力学中的转动惯量(inertia)是平动动力学中的质量(mass)的等效,是一个\(3 \times 3\)矩阵。在初始状态下,转动惯量为一常量,定义式为

\[ \mathbf I_{\text{ref}} = \sum m_i (\mathbf r_i^T \mathbf r_i \mathbf 1 - \mathbf r_i \mathbf r_i^T). \]

在刚体经过旋转变换\(\mathbf R\)后,可以证明,此时转动惯量为

\[ \mathbf I = \mathbf R \mathbf I_{\text{ref}} \mathbf R^T. \]

image
image

请注意:一个三维空间中的物体存在无数个可能的转动轴,每个转动轴都对应着一个描述其转动惯性大小的标量\(I\)。但这些量并非完全相互独立,它们可以被整理成一个\(3 \times 3\)的矩阵\(\mathbf I_{\text{ref}}\)

刚体仿真-仅转动

  1. 输入状态\(\mathbf s^{(0)} = \{\mathbf q^{(0)}, \mathbf \omega^{(0)}\}\)
  2. 更新中间物理状态

    \[ \begin{cases} \mathbf R^{(0)} = \text{Matrix.Rotate}(\mathbf q^{(0)}), \\ \mathbf \tau_i^{(0)} = (\mathbf R^{(0)} \mathbf r_i) \times \mathbf f_i^{(0)}, \\ \mathbf \tau^{(0)} = \sum \mathbf \tau_i^{(0)}, \\ \mathbf I^{(0)} = \mathbf R^{(0)} \mathbf I_{\text{ref}} (\mathbf R^{(0)})^T, \end{cases} \]
  3. 更新刚体仿真状态

    \[ \begin{cases} \mathbf \omega^{(1)} = \mathbf \omega^{(0)} + \Delta t (\mathbf I^{(0)})^{-1} \mathbf \tau^{(0)}, \\ \mathbf q^{(1)} = \mathbf q^{(0)} + \left[ 0, \frac{\Delta t}{2} \mathbf \omega^{(1)} \right] \times \mathbf q^{(0)} \end{cases} \]
  4. 输出状态\(\mathbf s^{(1)} = \{\mathbf q^{(1)}, \mathbf \omega^{(1)}\}\)

刚体动画

一般刚体仿真

刚体仿真流水线

  1. 输入状态\(\mathbf s^{(0)} = \{\mathbf v^{(0)}, \mathbf x^{(0)}, \mathbf \omega^{(0)}, \mathbf q^{(0)}\}\)
  2. 更新中间物理状态

    \[ \begin{cases} \mathbf f_{i}^{(0)} = \text{Force}(\mathbf x_i^{(0)}, \mathbf v_i^{(0)}, \dots), \\ \mathbf f^{(0)} = (\mathbf f_1^{(0)}, \dots, \mathbf f_n^{(0)})^T, \\ \mathbf R^{(0)} = \text{Matrix.Rotate}(\mathbf q^{(0)}), \\ \mathbf \tau_i^{(0)} = (\mathbf R^{(0)} \mathbf r_i) \times \mathbf f_i^{(0)}, \\ \mathbf \tau^{(0)} = \sum \mathbf \tau_i^{(0)}, \\ \mathbf I^{(0)} = \mathbf R^{(0)} \mathbf I_{\text{ref}} (\mathbf R^{(0)})^T, \end{cases} \]
  3. 更新刚体仿真状态

    \[ \begin{cases} \mathbf v^{(1)} = \mathbf v^{(0)} + \Delta t M^{-1} \mathbf f^{(0)}, \\ \mathbf x^{(1)} = \mathbf x^{(0)} + \Delta t \mathbf v^{(1)} \\ \mathbf \omega^{(1)} = \mathbf \omega^{(0)} + \Delta t (\mathbf I^{(0)})^{-1} \mathbf \tau^{(0)}, \\ \mathbf q^{(1)} = \mathbf q^{(0)} + \left[ 0, \frac{\Delta t}{2} \mathbf \omega^{(1)} \right] \times \mathbf q^{(0)} \end{cases} \]
  4. 输出状态\(\mathbf s^{(1)} = \{\mathbf v^{(1)}, \mathbf x^{(1)}, \mathbf \omega^{(1)}, \mathbf q^{(1)}\}\)

Shape Matching 刚体仿真

预备知识:奇异值分解与极分解

  • 奇异值分解:任何变换均可被依次分解成三个部分:(可能带翻转的)旋转、缩放和(可能带翻转的)旋转。即

    \[ \mathbf A = \mathbf U \mathbf D \mathbf V^T, \]

    其中\(\mathbf A\)为任意方阵,\(\mathbf U\)\(\mathbf V^T\)为正交矩阵,\(\mathbf D\)为对角矩阵
  • 极分解:任何变换均可被依次分解成两个部分:(可能带翻转的)旋转和其他变形。即

    \[ \mathbf A = \mathbf U \mathbf D \mathbf V^T = (\mathbf U \mathbf V^T) ( \mathbf V \mathbf D \mathbf V^T) = \mathbf R \mathbf S, \]

    其中\(\mathbf R\)为正交矩阵,\(\mathbf S\)为半正定对称矩阵

Shape Matching

  • 基本思想:先对每个顶点做自由模拟,然后再把这些顶点用某种方法变回刚体
  • 变回刚体:让各顶点的位置满足刚体的位置约束,并且最小化各顶点在该过程中的移动距离。旋转变换通过极分解(polar decomposition)得到
  • Shape Matching 思想与可变形体(弹性体)的仿真方案一脉相承,因此适合在弹性体上实现刚体,如衣服中的纽扣等挂件
  • 由于欠缺物理正确性,Shape Matching 不适合处理带摩擦的碰撞情形
image
image

Shape Matching - 公式推导

  • 记仿真得到的刚体的预期质心为\(\mathbf c\),顶点\(i\)的位置为\(\mathbf r_i\),依质心性质有\(\sum_i \mathbf r_i = \mathbf 0\)
  • \(\mathbf y_i\)为自由仿真后的顶点位置,\(\mathbf x_i\)为满足刚体约束的顶点位置,则有\(\mathbf x_i = \mathbf c + \mathbf R \mathbf r_i\),其中\(\mathbf R\)为某旋转矩阵
  • 优化目标为

    \[ \{ \mathbf c, \mathbf R \} = \text{argmin} \sum_{i} ||\mathbf c + \mathbf R \mathbf r_i - \mathbf y_i||^2 \]

image
image

Shape Matching - 公式推导

\[ \{ \mathbf c, \mathbf R \} = \text{argmin} \sum_{i} ||\mathbf c + \mathbf R \mathbf r_i - \mathbf y_i||^2 \]

考虑\(\mathbf A\)为任意矩阵,其中的旋转分量为\(\mathbf R\),那么我们可以优化

\[ \{ \mathbf c, \mathbf A \} = \text{argmin} \sum_{i} ||\mathbf c + \mathbf A \mathbf r_i - \mathbf y_i||^2 \]

求导求极值

\[ \frac{\partial E}{\partial \mathbf c} = 0 \leftrightarrow \sum_{i} \mathbf c + \mathbf A \mathbf r_i - \mathbf y_i = 0 \]

得到

\[ \mathbf c = \frac{1}{n} \sum_{i} \mathbf y_i \]

即变回刚体前后的刚体质心位置不变。

Shape Matching - 公式推导

\[ \{ \mathbf c, \mathbf R \} = \text{argmin} \sum_{i} ||\mathbf c + \mathbf R \mathbf r_i - \mathbf y_i||^2 \]

另一方面

\[ \frac{\partial E}{\partial \mathbf A} = 0 \leftrightarrow \sum_{i} (\mathbf c + \mathbf A \mathbf r_i - \mathbf y_i) \mathbf r_i^T = 0 \]

得到

\[ \mathbf A = \left(\sum_i (\mathbf y_i - \mathbf c) \mathbf r_i^T \right) \left( \sum_i \mathbf r_i \mathbf r_i^T \right )^{-1} \]

最后,我们将\(\mathbf A\)做极分解

\[ \mathbf A = \mathbf R \mathbf S \]

使用其中的旋转部分\(\mathbf R\)即可。

Shape Matching 流水线

在 Shape Matching 刚体仿真流程中,不再以统一的\(\mathbb R^3\)向量更新刚体的位置、速度等仿真状态,而是单独仿真每个顶点;也不再考虑转动动力学。即此时\(\mathbf x = (\mathbf x_1, \dots, \mathbf x_n)^T\)\(\mathbf v = (\mathbf v_1 \dots, \mathbf v_n)^T\)

  1. 输入状态 \(\mathbf s^{(0)} = \{ \mathbf x^{(0)}, \mathbf v^{(0)} \}\)
  2. 对每个顶点单独仿真

    \[ \begin{cases} \mathbf f_i = \text{Force}(\mathbf x_i^{(0)}, \mathbf v_i^{(0)}, \dots), \\ \mathbf v_i = \mathbf v_i^{(0)} + \Delta t m_i^{-1} \mathbf f_i, \\ \mathbf y_i = \mathbf x_i^{(0)} + \Delta t \mathbf v_i, \end{cases} \]
  3. 计算刚体位置约束

    \[ \begin{cases} \mathbf c = \frac{1}{n} \sum_{i} \mathbf y_i, \\ \mathbf A = \left(\sum_i (\mathbf y_i - \mathbf c) \mathbf r_i^T \right) \left( \sum_i \mathbf r_i \mathbf r_i^T \right )^{-1}, \\ \mathbf R = \text{Polar}(\mathbf A), \end{cases} \]
  4. 更新各顶点仿真状态

    \[ \begin{cases} \mathbf x_i^{(1)} = \mathbf c + \mathbf R \mathbf r_i, \\ \mathbf v_i^{(1)} = (\mathbf x_i^{(1)} - \mathbf x_i^{(0)}) / \Delta t, \end{cases} \]
  5. 输出状态 \(\mathbf s^{(1)} = \{ \mathbf x^{(1)}, \mathbf v^{(1)} \}\)

质点碰撞处理

预备知识

有向距离场

有向距离场(signed distance field, SDF)\(\phi(\mathbf x): \mathbb R^3 \to \mathbb R\)针对空间中的位置输出其到某个物体表面的最短距离。若该位置在物体内部,则输出距离为负,否则距离为正。

image
image

\(\phi(\mathbf x)\)的梯度\(\nabla \phi(\mathbf x)\)的方向给出了最快朝外远离物体表面的方向,可以类比法线方向理解。

罚函数法

朴素线性罚函数

  • 质点碰撞检测的一般方法:检测其位置\(\mathbf x\)是否满足\(\phi(\mathbf x) \geq 0\),若否,则意味着该质点位于物体内部,需进行碰撞处理
  • 质点碰撞处理的一般方法:将其沿\(\nabla \phi(\mathbf x)\)的方向,移动到物体外
  • 罚函数法(penalty method)是一类简单的碰撞处理方法,但在宏观物理上不一定正确。它的一般思路是在质点位于物体内部,或接近物体时,给其一个推力,使其远离物体。推力的方向\(\mathbf N\)一般是归一化的\(\nabla \phi(\mathbf x)\)
  • 朴素的罚函数法使用一个正比于\(\phi(\mathbf x)\)的惩罚推力

    \[ \mathbf f = -k \phi(\mathbf x) \mathbf N \]
  • 该方法只有顶点位于碰撞体内部时才能将其推开,存在穿透瑕疵

image
image

有缓冲区的线性罚函数

  • 加入一个宽度为\(\varepsilon\)的缓冲区,可一定程度上缓解穿透问题。此时若检测到\(\phi(\mathbf x) < \varepsilon\),就要启动碰撞处理,施加惩罚推力

    \[ \mathbf f = -k (\phi(\mathbf x) - \varepsilon) \mathbf N \]
  • 此方法可以让质点在穿透之前就有远离物体的倾向,且计算简单。但若质点速率很高,直接穿进了碰撞体,将导致推力过大,产生视觉不正确的"飞天"现象。

image
image

使用log-barrier的罚函数

  • 同样地,考虑在一定距离上就施加惩罚推力,可以设计使用log-barrier的罚函数

    \[ \mathbf f = \rho \frac{1}{\phi(\mathbf x)} \mathbf N \]
  • 问题仍然存在:穿透后,将产生错误的推力方向

image
image

总结:

  • 如今的物理引擎一般都实装一套简单的罚函数法,罚函数一般为有缓冲区的线性罚函数或log-barrier罚函数。它们实现简单,但都可能存在穿透瑕疵和推力可能过大的问题
  • 罚函数法不适合处理摩擦碰撞的情况

冲量法

冲量法

  • 冲量法(impulse method)是一类假设质点的碰撞将导致位置和速度的突变的碰撞处理方法,它基于库仑摩擦定律,在宏观物理上表现得更正确。它的一般思路是在质点位于物体内部时,将质点的位置直接移动到最近的物体表面,并改变其速度方向使其不要再接近物体
  • 对碰撞位置\(\mathbf x\),新的位置\(\mathbf x^{\text{new}}\)的更新方法是

    \[ \mathbf x^{\text{new}} = \mathbf x - \nabla \phi(\mathbf x) \phi(\mathbf x) \]

image
image

冲量法 - 速度更新

  • 仅仅更新位置是不够的,为了避免接下来仍然发生穿透,我们还需要更新质点的速度
  • 对碰撞时的速度\(\mathbf v\)正交分解得到

    \[ \begin{cases} \mathbf v_N = (\mathbf v \cdot \mathbf N) \mathbf N, \\ \mathbf v_T = \mathbf v - \mathbf v_N \end{cases} \]

    image
  • 考虑库仑摩擦,分别修改法向速度和切向速度

    \[ \begin{cases} \mathbf v_N^{\text{new}} = - \mu_N \mathbf v_N, \\ \mathbf v_T^{\text{new}} = a \mathbf v_T \end{cases} \]
  • 最后合成得到

    \[ \mathbf v^{\text{new}} = \mathbf v_N^{\text{new}} + \mathbf v_T^{\text{new}} \]

其中 \(\mu_N, \mu_T\) 是摩擦因数,参数 \(a\) 应当满足质点摩擦的库仑定律:

\[ a = \max \left (1 - \mu_T(1 + \mu_N)\frac{||\mathbf v_N||}{||\mathbf v_T||}, 0 \right ) \]

前项是动摩擦的结果,后项是静摩擦的结果。

刚体碰撞检测与处理

刚体碰撞检测

刚体碰撞检测与惩罚法的碰撞处理

image
image
  • 朴素的刚体碰撞检测方法:对刚体的所有顶点\(i\),逐一检测\(\phi(\mathbf x_i)\)是否为负值
  • 聪明一点的碰撞检测方法:对碰撞不频繁的刚体,预先建立BVH,进行\(O(\log n)\)级别的检测
  • 惩罚法的碰撞处理:对产生碰撞的顶点\(i\),累加相应的惩罚力\(\mathbf f_i\)

冲量法的碰撞处理

冲量法的碰撞处理

image
image
  • 使用冲量法进行刚体碰撞处理时,我们需要修改碰撞顶点的状态包括位置\(\mathbf x_i\)、速度\(\mathbf v_i\)和角速度\(\mathbf \omega_i\)。但刚体为一整体,单独修改顶点的状态是不可行的
  • 为此,我们需要从局部状态的修改,推导出整体状态的修改,思路类似逆向运动学(IK)

冲量法的碰撞处理-从局部到整体

image
image
  • 简单起见,我们先只讨论碰撞点为单个顶点的情况
  • 初始时,\(\mathbf v_i = \mathbf v + \mathbf \omega \times \mathbf R \mathbf r_i\)
  • 顶点\(i\)被施加冲量\(\mathbf j_i = \mathbf f_i \Delta t\)后,有

    \[ \begin{cases} \mathbf v^{\text{new}} = \mathbf v + M^{-1} \mathbf j_i, \\ \mathbf \omega^{\text{new}} = \mathbf \omega + \mathbf I^{-1} (\mathbf R \mathbf r_i \times \mathbf j_i) \end{cases} \]
  • 此时,\(\mathbf v_i^{\text{new}}\)与整体\(\mathbf v^{\text{new}}, \mathbf \omega^{\text{new}}\)的关系为

    \[ \mathbf v_i^{\text{new}} = \mathbf v^{\text{new}} + \mathbf \omega^{\text{new}} \times \mathbf R \mathbf r_i \]

冲量法的碰撞处理 - 推导

代入计算得到

\[ \begin{aligned} \mathbf v_i^{\text{new}} &= \mathbf v^{\text{new}} + \mathbf \omega^{\text{new}} \times \mathbf R \mathbf r_i \\ &= \mathbf v + M^{-1} \mathbf j_i + (\mathbf \omega + \mathbf I^{-1} (\mathbf R \mathbf r_i \times \mathbf j_i)) \times \mathbf R \mathbf r_i \\ &= \mathbf v_i + M^{-1} \mathbf j_i + (\mathbf I^{-1} ( \mathbf R \mathbf r_i \times \mathbf j_i )) \times \mathbf R \mathbf r_i \\ &= \mathbf v_i + M^{-1} \mathbf j_i - \mathbf F_i \mathbf j_i \end{aligned} \]

其中\(\mathbf F_i\)是一个与\(\mathbf R, \mathbf r_i, \mathbf I\)有关的线性算子。于是

\[ \mathbf j_i = \mathbf K_i^{-1} (\mathbf v_i^{\text{new}} - \mathbf v_i) \]

其中

\[ \mathbf K_i = M^{-1} \mathbf 1 - \mathbf F_i. \]

冲量法的碰撞处理 - 流水线

冲量法碰撞处理的一般步骤:

  1. 输入状态\(\mathbf s = \{ \mathbf c, \mathbf v, \mathbf \omega \}\)
  2. 使用SDF检测碰撞,记碰撞点集合为\(U\)。若\(U\)为空集,则无需碰撞处理
  3. 否则,计算每个碰撞点更新后的位置和速度

    \[ \begin{cases} \mathbf v_i = \mathbf v + \mathbf \omega \times \mathbf R \mathbf r_i, \\ \mathbf v_{N, i} = (\mathbf v_i \cdot \mathbf N) \mathbf N, \\ \mathbf v_{T, i} = \mathbf v_i - \mathbf v_{N, i}, \\ a = \max \left (1 - \mu_T(1 + \mu_N)\frac{||\mathbf v_{N, i}||}{||\mathbf v_{T, i}||}, 0 \right ), \\ \mathbf v_{N, i}^{\text{new}} = - \mu_N \mathbf v_{N, i}, \\ \mathbf v_{T, i}^{\text{new}} = a \mathbf v_{T, i}, \\ \mathbf v^{\text{new}}_i = \mathbf v_{N, i}^{\text{new}} + \mathbf v_{T, i}^{\text{new}} \end{cases} \]
  4. 对每个碰撞点,计算使得其速度发生如上变化的冲量\(\mathbf j_i\)

    \[ \begin{cases} \mathbf K_i = M^{-1} \mathbf 1 - \mathbf F_i, \\ \mathbf j_i = \mathbf K_i^{-1} (\mathbf v_i^{\text{new}} - \mathbf v_i), \end{cases} \]
  5. 若有多个碰撞点,求该冲量的平均值;取陷入物体最深的那个碰撞点,作为位移的参考

    \[ \begin{cases} \mathbf j = \frac{\sum \mathbf j_i}{ |U| }, \\ \mathbf x = \min_U(\mathbf x_i) \end{cases} \]
  6. 最后,利用该冲量更新刚体整体的速度和角速度,并更新位置以移除碰撞点

    \[ \begin{cases} \mathbf c^{\text{new}} = \mathbf c + \mathbf R (\mathbf x - \nabla \phi(\mathbf x) \phi(\mathbf x)) , \\ \mathbf v^{\text{new}} = \mathbf v + M^{-1} \mathbf j, \\ \mathbf \omega^{\text{new}} = \mathbf \omega + \mathbf I^{-1} (\mathbf R \mathbf r_i \times \mathbf j) \end{cases} \]
  7. 输出状态\(\mathbf s^{\text{new}} = \{ \mathbf c^{\text{new}}, \mathbf v^{\text{new}}, \mathbf \omega^{\text{new}} \}\)

并行碰撞处理

当多点同时接触时,并行碰撞处理问题,会变成维度更高的线性系统。实际处理时,常进行简单的串行处理。

image
image

Shape Matching 的碰撞处理

Shape Matching 的碰撞处理

  • Shape Matching仿真流程中,只需对每个顶点单独进行碰撞处理即可。这一处理方式与可变形体相同
  • 若使用冲量法进行碰撞处理,则为了保证无穿透,需要在刚体约束之前、约束之后均再执行一遍碰撞检测和处理
image
image

Thanks for your watching!