PBD（Position Based Dynamics）学习笔记

符号说明

仿真物体包括 $N$ 个节点和 $M$ 个约束。

每个节点 $i\in [1,...,N]$ 包含的参数有质量 $m_i$ 、位置 ${\mathbf{x}}_{\mathbf{i}}$ 和速度 ${\mathbf{v}}_{\mathbf{i}}$ ；

每个约束 $j\in [1,...,M]$ 包含的参数有：

• 基数 $n_j$ ；
• 约束函数 $C_j:{\mathbb{R}}^{3n_j}\to \mathbb{R}$ ；
• 节点编号集合 {i1,...,inj}, ik∈[1,...,N]\{i_1,...,i_{n_j}\}, \ i_k\in[1,...,N]{i1​,...,inj​​}, ik​∈[1,...,N] ；
• 刚度参数 $k_j\in [0,1]$ ，描述约束的强度；
• 约束类型：等式 或 不等式；
  • 等式表示 $C_j({\mathbf{x}}_{{\mathbf{i}}_1},...,{\mathbf{x}}_{{\mathbf{i}}_{{\mathbf{n}}_{\mathbf{j}}}})=0$ ；
  • 不等式表示 $C_j({\mathbf{x}}_{{\mathbf{i}}_1},...,{\mathbf{x}}_{{\mathbf{i}}_{{\mathbf{n}}_{\mathbf{j}}}})\ge 0$ ；

算法流程

首先初始化 ${\mathbf{x}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}^0$，${\mathbf{v}}_{\mathbf{i}}={\mathbf{v}}_{\mathbf{i}}^0$，${\omega }_i=1/m_i$ ，然后在每个 $\Delta t$ 时间内重复一次以下过程：

1. 处理外力：${\mathbf{v}}_{\mathbf{i}}\leftarrow {\mathbf{v}}_{\mathbf{i}}+\Delta t{w}_i{f}_{ext}({\mathbf{x}}_{\mathbf{i}})$ ；
2. 速度阻尼：$dampVelocities({\mathbf{v}}_1,...,{\mathbf{v}}_{\mathbf{N}})$ ；
3. 下一个位置预测：${\mathbf{p}}_{\mathbf{i}}\leftarrow {\mathbf{x}}_{\mathbf{i}}+\Delta t{\mathbf{v}}_{\mathbf{i}}$ ；
4. 碰撞检测，生成 $M_{coll}$ 个碰撞约束：$generateCollisionConstraints({\mathbf{x}}_{\mathbf{i}}\to {\mathbf{p}}_{\mathbf{i}})$ ；
5. 若干次迭代：每次迭代处理一遍所有约束 $projectConstraints(C_1,...,C_{M+M_{coll}},{\mathbf{p}}_1,...,{\mathbf{p}}_{\mathbf{N}})$ ；
6. 速度、位置更新：${\mathbf{v}}_{\mathbf{i}}\leftarrow ({\mathbf{p}}_{\mathbf{i}}-{\mathbf{x}}_{\mathbf{i}})/\Delta t$ ，${\mathbf{x}}_{\mathbf{i}}\leftarrow {\mathbf{p}}_{\mathbf{i}}$ ；
7. 碰撞节点的速度更新：$velocityUpdate({\mathbf{v}}_{\mathbf{i}},...,{\mathbf{v}}_{\mathbf{N}})$ ；

流程详细说明

2

由于物体内部的摩擦力、粘性等特性，节点速度变化时会有阻尼。PBD中速度阻尼的计算方法如下：

首先计算系统的质心、角速度：

• 质心 ${\mathbf{x}}_{\mathbf{c}}=({\Sigma }_i{\mathbf{x}}_{\mathbf{i}}{m}_i)/({\Sigma }_i{m}_i)$ ；
• 质心速度 ${\mathbf{v}}_{\mathbf{c}}=({\Sigma }_i{\mathbf{v}}_{\mathbf{i}}{m}_i)/({\Sigma }_i{m}_i)$ ；
• 系统角动量 $\mathbf{L}={\Sigma }_i{\mathbf{r}}_{\mathbf{i}}\times (m_i{\mathbf{v}}_{\mathbf{i}})$ ，其中 ${\mathbf{r}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}-{\mathbf{x}}_{\mathbf{c}}$ ；
• 系统转动惯量 $\mathbf{I}={\Sigma }_i\widetilde{{\mathbf{r}}_{\mathbf{i}}}{\widetilde{{\mathbf{r}}_{\mathbf{i}}}}^T{m}_i$ ；
• 系统角速度 $\mathbf{\omega }={\mathbf{I}}^{-1}\mathbf{L}$ ；

然后设定一个全局的参数 $k_{damping}\in [0,1]$ ，对于每个节点的速度 ${\mathbf{v}}_{\mathbf{i}}$，做如下变换：

1. $\Delta {\mathbf{v}}_{\mathbf{i}}={\mathbf{v}}_{\mathbf{c}}+\mathbf{\omega }\times {\mathbf{r}}_{\mathbf{i}}-{\mathbf{v}}_{\mathbf{i}}$；
2. ${\mathbf{v}}_{\mathbf{i}}\leftarrow {\mathbf{v}}_{\mathbf{i}}+k_{damping}\Delta {\mathbf{v}}_{\mathbf{i}}$ ；

从中可看出 $k_{damping}$ 描述了阻尼强度，如果 $k_{damping}=1$，那么系统为刚体；如果 $k_{damping}=0$，那么每个节点自由运动，不存在阻尼。

4

当某个物体的一个节点从 $\mathbf{x}$ 运动到 $\mathbf{p}$ 时，可能与其它物体发生碰撞，这里碰撞有两种类型（原文描述的两种类型为 continuous collision 和 static collision，但我不明白这么命名的原因）：

• $\mathbf{x}$ 在另一个物体外部而 $\mathbf{p}$ 在另一个物体内部；
• $\mathbf{x}$ 和 $\mathbf{p}$ 都在另一个物体内部；

对于这两种类型，我们都希望 $p$ 能够被约束到另一个物体外部，因此构造如下约束：

• 对于第一种类型，构造不等式约束 $C(\mathbf{p})=(\mathbf{p}-{\mathbf{q}}_{\mathbf{c}})\cdot {\mathbf{n}}_{\mathbf{c}}$ ，刚度参数 $k=1$；
• 对于第二种类型，构造不等式约束 $C(\mathbf{p})=(\mathbf{p}-{\mathbf{q}}_{\mathbf{s}})\cdot {\mathbf{n}}_{\mathbf{s}}$ ，刚度参数 $k=1$；

还有其它很多碰撞检测的约束方法，比如运动节点和运动三角形的。

5

这一步要迭代若干次，因为其本质是用了类似**牛顿迭代法（Newton’s Method / Newton-Raphson Method）**的方法。

牛顿迭代法，用于求解方程 $f(x)=0$，迭代函数为 $x_{k+1}=x_k-\frac{f(x_k)}{f^{{}^{\prime }}(x_k)}$，迭代足够多次时 $x$ 接近零点。

而每一次迭代，采取了 Gauss-Seidel method 的思想，即独立地一个个地处理每个约束（Gauss-Seidel method 在求解 $Ax=b$ 时，一个个地求取 $x_i$，并且当前求取的 $x_i$ 会用到后面分量的求取中）。

因此，第5步的过程就是：循环若干次，每次循环时一个个处理所有约束。

约束有两类，一类是固有的 $M$ 个约束，另一类是每次碰撞检测生成的 $M_{coll}$ 个约束，前者需要满足动量和角动量守恒，而后者不用。由于 $C(\mathbf{p})$ 的值与刚体运动（平移、旋转）无关，也就是说刚体运动对应 $\mathbf{p}$ 的变化位于 $C(\mathbf{p})$ 的等值线上，因此梯度 ${\nabla }_{\mathbf{p}}C(\mathbf{p})$ 与刚体运动正交，如果 $\mathbf{p}$ 的变化 $\Delta \mathbf{p}$ 位于 ${\nabla }_{\mathbf{p}}C(\mathbf{p})$ 的方向上，那么自然能保持动量守恒（这里存疑……）。

假设所有节点质量相等，下面说明如何处理一个等式约束 $C(\mathbf{p})=C({\mathbf{p}}_1,...,{\mathbf{p}}_{\mathbf{n}})=0$：

对于当前的 $\mathbf{p}$ ，我们需要计算一个 $\mathbf{\Delta }\mathbf{p}$，使得 $C(\mathbf{p}+\mathbf{\Delta }\mathbf{p})\approx C(\mathbf{p})+{\nabla }_{\mathbf{p}}C(\mathbf{p})\cdot \mathbf{\Delta }\mathbf{p}=0$ ，由于 $\Delta \mathbf{p}$ 在梯度的方向上，因此可设 $\Delta \mathbf{p}=\lambda {\nabla }_{\mathbf{p}}C(\mathbf{p})$ ，联立求解可得
$$\Delta \mathbf{p}=-\frac{C(\mathbf{p})}{|{\nabla }_{\mathbf{p}}C(\mathbf{p}){|}^2}{\nabla }_{\mathbf{p}}C(\mathbf{p})\text{\hspace{1em}(1)}$$
上式就是牛顿迭代法中的一步变化 $\mathbf{p}\leftarrow \mathbf{p}+\Delta \mathbf{p}$，而一次迭代需要对所有约束都作这样的计算。

我们将 $\Delta \mathbf{p}$ 进行拆分，可以得到与 (1) 等价的式子
$$\Delta {\mathbf{p}}_i=-s{\nabla }_{{\mathbf{p}}_i}C(\mathbf{p})\text{\hspace{1em}(2)}$$
其中
$$s=\frac{C(\mathbf{p})}{|{\nabla }_{\mathbf{p}}C(\mathbf{p}){|}^2}=\frac{C(\mathbf{p})}{{\Sigma }_j|{\nabla }_{{\mathbf{p}}_j}C(\mathbf{p}){|}^2}\text{\hspace{1em}(3)}$$
这么做是为了处理节点质量不相等的情况，当质量不相等时，${\mathbf{p}}_i$ 的变化应当与质量成反比。因此考虑质量的倒数 ${\omega }_i=1/m_i$ ，可得
$$\Delta {\mathbf{p}}_i=-s\frac{n{\omega }_i}{{\Sigma }_j{\omega }_j}{\nabla }_{{\mathbf{p}}_i}C(\mathbf{p})\text{\hspace{1em}(4)}$$
对于不等式约束，我们只在 $C(\mathbf{p})<0$ 时进行迭代变换。

如果再考虑每个约束的刚度参数 $k$，假设需要 $n_s$ 轮迭代（这里或许是：假设还需要 $n_s$ 轮迭代？），我们将位置修正 $\Delta \mathbf{p}$ 乘上系数 $1-(1-k{)}^{1/n_s}$ ，存在线性误差 $\Delta \mathbf{p}(1-k)$ 。

7

这一步是用来更新发生了碰撞的节点的速度的。

如果一个节点发生了碰撞，假设碰撞位置的法向量为 $\mathbf{n}$ ，那么该节点速度垂直于 $\mathbf{n}$ 的分量将因为摩擦力而减少，并且由于反弹会产生与 $\mathbf{n}$ 同方向的速度分量。