GAMES103笔记 Lecture5 基于物理的布料模拟（Physics-based Cloth Simulation）

弹簧质点系统

布料模拟中需要使用到弹簧模型，第二次课上老师已经提过弹簧模型，我没放到第二次课的笔记中，与弹簧有关的内容都在本篇中。

先来看一个简单的弹簧，弹簧的一侧固定在原点 $\mathbf{o}$ 上， $L$ 是弹簧原本的长度， $\mathbf{x}$ 是弹簧末端的位置（这是一个三维向量）。

第一个式子计算当前系统的弹性势能 $E(\mathbf{x}) = \frac{1}{2}k(\|\mathbf{x}\|-L)^2$ ，其中 $\|\mathbf{x}\|-L$ 代表弹簧的形变

弹簧的力希望释放掉弹簧系统的弹性势能， $E(\mathbf{x})$ 对 $\mathbf{x}$ 求梯度再取反即可得到弹力。

胡克定律（Hooke’s law）：$\mathbf{f}(\mathbf{x})=-k(\|\mathbf{x}\|-L) \frac{\mathbf{x}}{\|\mathbf{x}\|}$

其中 $\frac{\mathbf{x}}{\|\mathbf{x}\|}$ 是单位向量，指示力的方向。

Tangent stiffness $\mathbf{H}(\mathbf{x})$ 是 $E(\mathbf{x})$ 的 Hessian 矩阵，以后做优化计算时需要用到 Hessian 来判断算法的有效性。

stiffness 是硬度、劲度的意思，对应系数 k。

如果弹簧的两个端点都可以自由移动，则弹簧的形变与两个端点的位置都有关。

这里 $\mathbf{x}_0$ 和 $\mathbf{x}_1$ 是可以自由移动的两个端点。

$E(\mathbf{x})$ 的自变量 $\mathbf{x}=[\mathbf{x}_0, \mathbf{x}_1]^T$ ，是一个6维向量，包括2个顶点的位置。

$\nabla_0E(\mathbf{x})$ 代表 $E(\mathbf{x})$ 对 $\mathbf{x}_0$ 求梯度，$\nabla_1E(\mathbf{x})$ 代表 $E(\mathbf{x})$ 对 $\mathbf{x}_1$ 求梯度。

弹簧两端受到的弹力仍然等于负的能量对端点位置的梯度，弹簧两端受到的力分别等于 $-\nabla_{\mathbf{x_0}E(\mathbf{x})}$ 和 $-\nabla_{\mathbf{x_1}E(\mathbf{x})}$ ，它们大小相同，方向相反。

$\mathbf{H}(\mathbf{x})$ 是 $E(\mathbf{x})$ 的 Hessian 矩阵，维度为 6x6，可以把两个端点的向量分开，写成分块矩阵形式。Lecture2 中证明过，对于这个形状的分块矩阵，如果 $\mathbf{H}_e$ 是正定矩阵，则 $\mathbf{H}(\mathbf{x})$ 一定是半正定矩阵。

如果有多个弹簧，可以独立计算每个弹簧的能量，然后对端点位置求梯度得到每个弹簧的分力。对于多个弹簧的公共端点，直接叠加弹簧的分力即可。

上图中的中心点 $\mathbf{x}_i$ 受到的弹力 $\mathbf{f}_i$ 等于四个弹簧施加的弹力之和。

用弹簧质点系统模拟布料

建模

如果布料（的Mesh）由方方正正的结构化的网格构成，我们需要在上面添加三种弹簧：

1. 在网格原本的边上添加弹簧，抵抗布料产生横向和纵向的的小皱褶
2. 在每个小格子的两个对角线上加两根弹簧，抵抗布料在对角方向上的皱褶
3. 跳过一个顶点连接的弹簧，抵抗布料沿中间的长边产生大的翻折（即Bending，后面可能翻译为弯曲、弯折，都是一个意思）

实际应用中可以简化一些，只在一个小格子的其中一个对角线上添加弹簧。

对于不规则形状的三角形 Mesh 布料，同样在三角网格的每个边上都加一根弹簧，然后跨过所有内部边添加弹簧（称任意两个相邻三角形的公共边为内部边，两个相邻三角形组成一个四边形，新添加的弹簧为这个四边形的对角线）来抵抗弯曲。

拿到一个任意的三角 Mesh 后，需要一个算法来完成添加弹簧的操作。

三角形网格在计算机中存储为两个数组。一个顶点数组依次存储 0 号 - n 号顶点的坐标，每个顶点坐标有3维度。另一个数组存储三角形的顶点在顶点数组中的索引，每 3 个数代表一个三角形，代表顶点数组中的这三个顶点构成一个三角形。

直接在 Mesh 的每个边上都添加弹簧会造成在内部边上重复添加弹簧，因此需要想办法做一些额外处理。

先从原始的 Mesh 数据中提取出所有边的信息构成一个边数组（Triple list），每个边由3个数字组成，前两个数代表它连接的两个顶点，第三个数代表它属于哪一个三角形。允许提取得到的边数组中存在出现两次的边，但它们的第三个一定不相同。

接下来对边数组依据字典序排序，排序之后内部边一定出现在相邻位置，因为它们的第一个数和第二个数都一样，只有第三个数不同。

遍历排序后的边数组时就能判断哪些是内部边，去除内部边可得到不重复的边数组（Edge list）。

对于三角形边上的弹簧，直接遍历 Edge list，依次添加弹簧即可。

对于跨过内部边的弹簧，查询内部边所在三角形的剩余两个顶点，把它们连接起来即可。

显式欧拉法

模拟时把 Mesh 的所有顶点看做一个粒子系统，弹簧被插入到这个粒子系统中。

先遍历所有的弹簧，算出每个弹簧对每个顶点施加的力。

然后遍历每个顶点，把弹簧对这个点的力加起来就是这个点受到的合力。接下来用之前刚体动力学介绍过的方法更新速度和位置即可。

显式积分的问题：overshooting。弹簧受到非常大的力时，会在一个时间步内产生过度形变，不符合能量守恒，导致顶点被弹得越来越远，造成数值不稳定。

原因：真实世界中当弹簧形变量缩小时弹力会逐渐减小，而离散的积分方法假定了一个时间步内受力都相同。当时间步长太大，或者弹性系数 k 太大时就会在一个时间步内产生过度的位移。

解决方法：减小时间步长（会降低计算效率），或者换成隐式积分。

隐式欧拉法

隐式积分法更新速度时用的是新的位置受到的力，也就是说当弹簧缩小时，显式积分用的弹力一定大于这个时间步内的平均弹力，而隐式积分用的力小于一个时间步内的平均弹力，因此隐式积分法在数值上更稳定一些。

隐式积分更新速度时用的是新位置的力 $\mathbf{f}^{[1]}$ ，在弹簧模型中这是一个与新位置 $\mathbf{x}^{[1]}$ 有关的未知量（PPT 中holonomic 指这个力只与位置有关），所以隐式积分法的更新公式无法像显式积分一样简单地算出来。

作业里模拟时还需要考虑重力，重力可看作常量，所以这里可以把力写作关于 $\mathbf{x}^{[1]}$ 的函数 $\mathbf{f}(\mathbf{x}^{[1]})$ 。

把速度更新公式代入位置更新公式消去 $\mathbf{v}^{[1]}$ 得到关于 $\mathbf{x}^{[1]}$ 的非线性方程 ：

$\mathbf{x}^{[1]}=\mathbf{x}^{[0]}+\Delta t \mathbf{v}^{[0]} + \Delta t^2 M^{-1} \mathbf{f}(\mathbf{x}^{[1]})$

解这个方程前，先转换一个视角看待这个方程，把它转换为优化问题。（下面介绍的的牛顿法其实也能直接解，不转化也行，转化成优化问题可以用到更多优化的理论）

老师没讲怎么转换的，我试着推一下

令 $g'(\mathbf{x}) = \mathbf{x} -\mathbf{x}^{[0]} - \Delta t \mathbf{v}^{[0]} - \Delta t^2 M^{-1} \mathbf{f}(\mathbf{x})$

求解原方程等价于求解 $g'(\mathbf{x}) = 0$ ，等价于求 $g'(\mathbf{x})$ 的原函数 $g(\mathbf{x})$ 的极值点（$g(\mathbf{x})$ 可能有多个极值点，后面会分析什么时候 $g(\mathbf{x})$ 只有一个极值，并且就是最小值）。

积分得到

\begin{align} g(\mathbf{x}) &= 2(\mathbf{x} -\mathbf{x}^{[0]}-\Delta t \mathbf{v}^{[0]})^2 + \Delta t^2 M^{-1} \mathbf{E}(\mathbf{x}) \\ &= 2\|\mathbf{x} -\mathbf{x}^{[0]}-\Delta t \mathbf{v}^{[0]}\|^2 + \Delta t^2 M^{-1} \mathbf{E}(\mathbf{x}) \end{align}

两边同乘 $\frac{M}{\Delta t^2}$ ，极值点不变

令 \begin{align} F(\mathbf{x}) = \frac{M}{\Delta t^2} g(\mathbf{x}) &= 2\frac{M}{\Delta t^2}\|\mathbf{x} -\mathbf{x}^{[0]}-\Delta t \mathbf{v}^{[0]}\|^2 + \mathbf{E}(\mathbf{x}) \end{align}

则求解方程 $\mathbf{x}^{[1]}=\mathbf{x}^{[0]}+\Delta t \mathbf{v}^{[0]} + \Delta t^2 M^{-1} \mathbf{f}(\mathbf{x}^{[1]})$ 等价于求 $\mathbf{x^{[1]}} = \arg \min F(\mathbf{x})$

反过来也可以令 $F(\mathbf{x})$ 的梯度等于 0，发现 $\mathbf{x}^{[1]}$ 是 $F(\mathbf{x})$ 的极值点，由此可以证明解上面的非线性方程就等价于最小化 $F(\mathbf{x})$ 。（后面会分析 $F(\mathbf{x})$ 的 Hessian，在某些条件下 $F(\mathbf{x})$ 只存在一个最小值）

牛顿法（Newton-Raphson method）

牛顿法可用来解方程（求函数的零点）。

当函数具有利普希茨连续性，即函数可以被两条直线夹逼时，牛顿法收敛。

例如用牛顿法解方程 $F'(x)=0$ ，初始解为 $x^{(0)}$

每次迭代时，求当前位置切线 $y-F'(x^{(k)})=F''(x^{(k)})(x-x^{(k)})$ 的零点得到下一轮的解 $x^{(k+1)}=x^{(k)}-\frac{F''(x^{(k)})}{F'(x^{(k)})}$

经过一定次数的迭代后，牛顿法收敛到其中一个解。

多维情况下牛顿法的迭代公式也是类似的。

代入弹簧隐式积分的问题，Hessian 里的第二部分 $\mathbf{H}(\mathbf{x}^{(k)})$ 就是之前算过的弹簧系统 $E(\mathbf{x})$ 的 Hessian。

黄色部分等价于解一个解线性方程组，后面会介绍解线性系统的方法。

$\mathbf{H}_e$ 可以分成两个正定矩阵相加，当 $1-\frac{L}{\|\mathbf{x}_{ij}\|}>0$ ，即 $\|\mathbf{x}_{ij}\| > L$ 时，加起来一定还是正定矩阵，所以：

弹簧拉伸时， $\mathbf{H}(\mathbf{x})$ 一定是s.p.d.， $\frac{\partial^2F(\mathbf{x}^{(k)})}{\partial\mathbf{x}^2}$ 一定是s.p.d.， $F(\mathbf{x})$ 只有一个最小值，牛顿法收敛到唯一的解。

弹簧压缩时， $\mathbf{H}(\mathbf{x})$ 可能不是s.p.d.，$\frac{\partial^2F(\mathbf{x}^{(k)})}{\partial\mathbf{x}^2}$ 可能不是s.p.d.。

因此当弹簧拉伸时系统更稳定。

一个直观的解释，在二维和三维空间中，当弹簧压缩时，可能会向多个方向弯曲，因此弹簧压缩时存在不稳定性。而一维的情况下二阶导就等于 $k$ ， $k>0$ ，因此不存在不稳定性。

解线性系统 $\mathbf{Ax}=\mathbf{b}$ 的方法

牛顿法中迭代的那一步需要解线性系统。

解线性系统分为直接法和迭代法。

直接法对 A 的要求小，适合在 CPU 上做，有一定内存开销。

迭代法对 A 有一定的要求，例如要求 A 正定，但可以在 GPU 上实现，有一些加速方法。

Jacobi 方法是一个解线性系统的迭代方法。

用切比雪夫法加速 Jacobi 方法的求解速度。

布料的弯曲（Bending）

用之前的弹簧模型来对布料建模时，如果布料只发生了非常轻微的弯曲，弹簧只发生微小的形变，弹簧产生的抵抗力也非常小，不太符合实际情况，还可以继续改进。

二面角模型

二面角模型用夹角来判断弯曲情况。

记夹角为 $\theta$ ，每个顶点 $\mathbf{x}_i$ 受到的力的大小等于 $f(\theta)$ ，力的方向为 $\mathbf{u}_i$

使 $\mathbf{u}_i$ 满足一些条件：

1. 要让两个三角形恢复共面，因此 $\mathbf{u}_1$ 和 $\mathbf{u}_2$ 分别与两个面的法向量同向。
2. 力不应该拉长或挤压公共边，因此 $\mathbf{u}_4-\mathbf{u}_3$ 应该垂直于公共边，而 $\mathbf{u}_1$ 和 $\mathbf{u}_2$ 也垂直于公共边，所以 $\mathbf{u}_4 - \mathbf{u}_3$ 可以用 $\mathbf{u}_1$ 和 $\mathbf{u}_2$ 的线性组合表示出来。
3. 合力为 0，所以 $\mathbf{u}_1+\mathbf{u}_2+\mathbf{u}_3+\mathbf{u}_4=\mathbf{0}$

由此解出 $\mathbf{u}_i$

力的大小 $f(\theta)$ 见下图，具体的推导过程参考论文 Bridson et al. 2003. Simulation of Clothing with Folds and Wrinkles. SCA.

二面角模型没有定义能量，隐式积分很麻烦，一般用显式积分。

Quadratic Bending Model

二次弯曲模型基于两个假设：（1）静止状态下两个三角形形成一个平面（2）顶点间的拉伸比较小，只考虑弯曲

定义能量 $E(\mathbf{x})$ 如下

其中Q是一个常数矩阵，$E(\mathbf{x})$ 是关于顶点位置的二次函数

最后算出的 $\mathbf{q}^T\mathbf{x}$ 代表了两个三角形的曲率

二次弯曲模型的力和 Hessian 都容易计算，便于隐式积分。

缺点是需要满足两个前提条件，布料本身的拉伸幅度小，静止状态下是平面。

对于非平面的情况，有后续处理方法 Cubic shell model 和 Projective dynamics model，这里没展开讲。

Locking Issues

之前的弹簧模型和弯曲模型把布料的拉伸和弯曲分开考虑，实际也会遇到不能拉伸但可以弯曲的材料比如纸。

但用弹簧模型来模拟不能拉伸的材料时，弹簧的弹性系数非常大，沿着穿着弹簧的边弯曲也会受到很大的阻力，因为弹簧几乎不允许发生形变，也不会跟着弯曲，只要穿过弹簧翻折布料就会受到很大的弹力顶住两端。这被称为 Locking 问题。

Locking 问题在弹簧弹性较弱、Mesh 的分辨率低的时候比较明显。

发生 Locking 问题的本质原因是模拟中自由度不足。

弹簧模型中 Mesh 的顶点可以自由移动，每个顶点有3个自由度，Mesh 之间的边是约束条件，自由度就等于顶点数*3减边数。

根据欧拉定理，Mesh 的顶点数*3减边数等于边界上边的数量加 3。

所以在弹簧模型中布料的自由度相对不足。

解决方法：

1. 减小弹簧压缩时的弹性系数
2. 使弹簧可以在一定范围内自由移动，超出范围才开始计算弹力

相关论文

Baraff and Witkin. 1998. Large Step in Cloth Simulation. SIGGRAPH.

Bridson et al. 2003. Simulation of Clothing with Folds and Wrinkles. SCA.

Choi and Ko. 2002. Stable But Responsive Cloth. TOG (SIGGRAPH)

Bergou et al. 2006. A Quadratic Bending Model for Inextensible Surfaces. SCA.

上一篇：GAMES103笔记 Lecture4 刚体碰撞（Rigid Body Contacts）

下一篇：GAMES103笔记 Lecture6 基于约束的方法（Constrained Approaches）

GAMES103 系列笔记目录