第五讲：弹性体模拟-数值方法
www.zhiqu.org 时间： 2024-06-16
这一讲是紧接着上一讲的内容，开始讨论如何通过数值方法来完成对弹性体形变的模拟，另外，顺带提一句，这部分内容跟游戏开发中的物理引擎具有较高关联。
上一讲降到，我们的弹性势能可以表示为：
值得澄清的是，这里是省略了高阶项之后的弹性势能，即这个公式是在线性形变的情况下成立的，如果要考虑更为复杂的形变，后面还要加上更多的高阶项。
下面来介绍，我们在计算机中要如何计算出弹性体形变之后产生的internal force，大致思路是先计算出形变梯度矩阵（Deformation Gradient），之后根据这个矩阵算出Green Strain Tensor，之后再算出势能表达式，最后对势能求导，就得到了对应的conservative force。
首先来看下，我们要如何计算形变梯度矩阵，这里先考虑3D软体的情况，2D软体的情况如布料等在后面再考虑。
首先，要想在计算机中进行模拟，就需要将3D软体离散化，2D表面离散化是通过将之拆解成三角形来完成，而3D Volume的离散则是将之拆分成一个个的四面体。
形变发生后，单个四面体就会从material space映射到deformed space，即四面体本身发生了形变，如下图所示：
紧接着给出两个假设：
由于四面体足够小，因此 满足微分中的一个极小的局部变化的要求，可以采用泰勒展开：
在线性形变的情况下，我们可以省略高阶项，我们可以得到：
上述公式中的 我们是知道的（最开始的时候， ，所以自然是知道的，之后，经过一个timestep迭代，最新的位置经过后续的计算，自然也是知道的），而形变前的数据 我们也知道，因此就可以计算出F，另外，我们前面也说过，F是一个与位置有关的变量，要精确考虑的话，四面体上的每一点对应的F都应该是不同的，不过这里我们可以假设：
在上述假设下，我们可以有更多的等式：
将前面的三组公式放在一起，我们就得到了下面的矩阵乘法形式，从而推导出F：
简化一下：
跟 都是一个3x3的矩阵，根据上面公式可以推导出：
这里是其中一个四面体的F，其他四面体也会有对应的F，根据F我们就可以计算Green Strain Tensor，这也是每个四面体一份的：
以及势能 ，对势能进行求导（具体来说，consititutive model有公式的话，可以手动推导其弹性力的表达式），我们就能得到internal force（弹性力）：
这里，势能是物体共有的一个参数，是一个标量，而 则是物体上所有四面体的顶点组成的一个向量 ，那么经过求导后得到的弹性力也自然是一个n维的向量，每一维代表的是每个顶点的受力。
上面是内部的力，外面来自于重力或者碰撞的力，我们用外力 表示。
其中的质量m是一个很大的vector，在上面的公式下，我们就可以用之前说过的implicit euler或者explicit euler完成时间步长积分，从而算出更新的位置向量。
上面这种用离散的四面体来模拟软体的弹性形变的方法叫做Finite Element Method（FEM），这种方法可以通过不断提高分辨率（缩小四面体的尺寸）来保证准确性，不过同时也会导致计算复杂度增加，所以游戏或者很多电影特效模拟不是用的这种方法，而是基于优化的方法，Optimization-Based Methods，这种方法可以实现快速模拟。
先来回顾一下隐式欧拉方法，速度跟位置向量的演变可以用下面公式来表达：
这里的两个向量都是包含了多个四面体的相关数据的，是一个3xn的矩阵。
前面说过隐式欧拉最终会变成一个高维方程，很难求解，这里考虑将之转换成一个优化问题。
首先，internal force是势能的导数给出如下：
根据上面的公式，我们有：
这边希望构造出一个函数 ，并且使得这个函数的极值：
达成的条件跟前面的公式相一致，也就是说，希望这个函数的极值对应的解 就正好是上面一个等式的解，而这个函数可以给出为（比较容易推导？）：
其中
从而我们可以知道：
前面说到，这个函数提出的目的是为了求得前面implicit timestep等式的等价解，这里的一个疑惑就是，同样是解方程，前面implicit方程不好解，这里的方程同样也不好解，那么这种变化的意义在哪里呢？
这里的作用在于公式中的后面部分 ，这个部分表示的是系统的势能，前面说过弹性形变的势能公式 是通过一系列的物理推导得出的，是精确解，这种求解起来会比较麻烦，但是在游戏中，我们可以使用非精确解，只要保证模拟的结果看起来真实就好，这样一来，我们就能够选择一些能够让这个求解过程变得容易的公式代入这里的势能项中，从而达到简化问题的目的。
那么这个势能公式要怎么来设计呢？我们知道，弹性形变导致的势能总体来说，可以分成两种，分别是体积变化导致的形变（如从体积为1变成体积为2），或者是扭曲导致的形变，如一个cube变成平行六面体等。
前者势能可以表示为：
det是矩阵的行列式，这个公式的意义在于，当没有形变的时候，矩阵就是单位矩阵，这个势能的结果就是0
后者势能可以表示为：
tr（trace）是矩阵的迹，一个n×n矩阵A的主对角线（从左上方至右下方的对角线）上各个元素的总和被称矩阵A的迹（也就等于所有特征值的和）。同样，当没有形变时，这个F就是单位矩阵，所以迹就是3，那么势能依然为0。
总结可以看出，我们需要的势能公式其实就是一个约束条件，这个约束条件或者说约束函数在无形变的时候输出为0，在有形变的时候输出的是一个正值（或者说在无形变的时候取得极小值0），比如其形式为：
在这个条件下，我们就可以定义势能为：
上式中的A是一个对角阵，C(q)是一个向量，对于其中的每个元素：
这里希望有：
这里的i用来指示四面体的序号，这里翻译一下，就是说我们希望每个四面体的形变梯度矩阵的行列式为1，也就是说每个四面体都是体积不变的。
下面来介绍一下Position-based Dynamics方法，回到前面的公式：
这里，我们将势能用脉冲函数（实际上是脉冲函数取反，不过为了描述方便，这里用脉冲函数来表示）来表示：
其中脉冲函数：
之所以这样设计，就是为了满足前面说过的，我们希望函数J取得最小值，那么就需要保证 ，这里的i依然是四面体的序号。
继续下去，Position-based Dynamics方法会采用Gauss-Siedel Iteration（假设1）进行求解，即对于后面的众多四面体而言，我们会先只考虑一个四面体作用下的求解，并将结果代入原公式之后考虑第二个四面体作用下的求解，循环往复。如果只考虑单个四面体作用，原公式就变成了：
再定义：
公式就变成：
将公式转换一下，我们想要求解J的最小值，就等价于在条件：
下求解
的最小值。
之后，再假设 足够小（假设2），那么前面的条件就可以做泰勒展开为：
而这个就可以用上一讲的拉格朗日乘子来求解，也就是说，在上面条件下的最小值求解就变成了：
而Gauss-siedel Iteration的问题是，需要进行多次迭代才能收敛，从而导致计算比较慢，性能较差，这也是Position-based Dynamics（PBD）的缺点，如果Iteration数目不够就会导致一些模拟效果质量较低，比如布料会看起来像弹簧一样生硬。
我们来回顾下之前的内容：
那么很自然的问题就是，是否有一种方法是居于两者之间的，一方面我们不需要求解一个十分复杂的方程，而是只需要求解一个简单的方程，且这个方程可以通过少数的迭代就能输出结果的。
这个方法就是物理界常用的Projective Dynamics，实际上这种方法跟PBD是处于同一框架下得，不同的是PBD这里的internel energy函数是一个脉冲函数之和，而Projective Dynamics方法的internal energy函数则是：
这个式子中一共引入了两个新的符号，一个是变量p，这个变量代表的是某个粒子（如前面的四面体）的一个目标状态（比如，我们可以选择四面体不发生形变时的状态为一个目标状态，当然我们也可以选择其他状态作为目标状态，具体如何选，后面再说），而d则是表达p跟q差异性的函数，可以是一个非常简单的函数，比如就取两者的距离作为函数。
整个函数是什么意思呢，我们这里来举个例子说明一下，假如我们选择某个函数（simple distance measure）作为上面的函数d：
这里的F跟前面一样，还是deformation gradient， 是一个没有形变，但是经过一个旋转的状态，两者都是3x3的矩阵，其中 可以通过一个叫做Shape Matching的算法求得。
另外一项 是一个脉冲函数，这个脉冲函数的取值受到约束条件c的影响，当 的时候，取值为0，否则取值为无穷大，而由于这个函数的存在，要想使得整体的internal energy最小，就必须得满足这里的约束条件，从而使得脉冲函数取值为0。
在这个例子中，我们选取了p为经过某个旋转 后的状态，而在这个状态下，可以满足上面的脉冲函数结果为0，那么最终的internal energy的计算就变成了一个非常简单的矩阵的线性计算。
这里，我们还可以选择其他的p，比如：
这里的 也是一个deformation gradient，表示的是最接近于当前形变状态的deformation gradient，但是在这个gradient的作用下，形变后的四面体的体积跟原始四面体保持一致（这种对于动物形体的模拟比较有用，比如肚子在受击之后体积保持不变）， 的计算同样可以通过脉冲函数求得，比如只需要找到体积不变的形变的约束条件，就能得到对应的deformation gradient。
在这种internal energy的作用下，最优解的求解就变成了如下的函数形式：
可以看到，p只出现在后两项，q则只出现在前两项，从而我们可以通过有限的迭代来求解这个最优化问题。
而由于脉冲函数的存在，p的求解只需要保证脉冲函数为0即可，因此整个问题可以快速变成q变量下的最优解，这个过程叫做local projection，即找到限制条件c（constraint）来保证脉冲函数为0，从而将问题从两个变量降维成一个变量。
在求得p之后，就可以来求解q以获得整个函数的最优解：
而由于上面的函数实际上都是q的一个二次函数，而其导数则是q的一次函数，从而只需要将导数等于0，就能得到最终的解，这一步叫做global solve。
前面说过，projective Dynamics是一个迭代的过程，上面展示的local projection跟global solve就是这个过程中的仅有两个步骤，这两个步骤需要循环往复的去执行（因为p是local projection的，也就是一个局部解，这个解并不一定就对应着最优解，而是需要不断调整q来逼近最优解）。
这个方法的迭代数目比PBD会少很多，目前UE、Houndini等软件使用的是X-PBD方案，是在PBD方案上的一种更新方案，但是其迭代次数依然比Projective Dynamics方法要多，之所以不用Projective Dynamics方法，是因为线性系统的求解在维数很大（比如q是一个很大的向量）的时候会存在比较大的困难。
~
#酆茜丽# MATLAB怎么进行数学建模? -
(15782337865)： 一、数学建模的一般步骤 数学建模并不是新东西,粗略地说, 数学建模是一个多次迭代的过程,每一次 迭代大体上包括:实际问题的抽象、简化, 做出假设,明确变量和参数;形成明确的 数学问题;以解析形式或者数值形式求解 该数学模型...
#酆茜丽# 什么是排除假设法 (建模方法) -
(15782337865)： 解: 这个方法就是排除假设法 (建模方法) 先把不可能的情况排除在外 再把可能的情况进行假设,再排除 最后得到正确的答案.
#酆茜丽# 有数字0,1,2,3,4,5,(1)可组成多少个没有重复数字的六?
(15782337865)： (1) 六个数字全排列:P(6,6)= 6! = 720个但第一位数字为0的不行,这种数字有P(5,5)= 120个所以能排出720-120 = 600个(2)全排列(包括第一个数字为0)的720个数,...
上一讲降到，我们的弹性势能可以表示为：
值得澄清的是，这里是省略了高阶项之后的弹性势能，即这个公式是在线性形变的情况下成立的，如果要考虑更为复杂的形变，后面还要加上更多的高阶项。
下面来介绍，我们在计算机中要如何计算出弹性体形变之后产生的internal force，大致思路是先计算出形变梯度矩阵（Deformation Gradient），之后根据这个矩阵算出Green Strain Tensor，之后再算出势能表达式，最后对势能求导，就得到了对应的conservative force。
首先来看下，我们要如何计算形变梯度矩阵，这里先考虑3D软体的情况，2D软体的情况如布料等在后面再考虑。
首先，要想在计算机中进行模拟，就需要将3D软体离散化，2D表面离散化是通过将之拆解成三角形来完成，而3D Volume的离散则是将之拆分成一个个的四面体。
形变发生后，单个四面体就会从material space映射到deformed space，即四面体本身发生了形变，如下图所示：
紧接着给出两个假设：
由于四面体足够小，因此 满足微分中的一个极小的局部变化的要求，可以采用泰勒展开：
在线性形变的情况下，我们可以省略高阶项，我们可以得到：
上述公式中的 我们是知道的（最开始的时候， ，所以自然是知道的，之后，经过一个timestep迭代，最新的位置经过后续的计算，自然也是知道的），而形变前的数据 我们也知道，因此就可以计算出F，另外，我们前面也说过，F是一个与位置有关的变量，要精确考虑的话，四面体上的每一点对应的F都应该是不同的，不过这里我们可以假设：
在上述假设下，我们可以有更多的等式：
将前面的三组公式放在一起，我们就得到了下面的矩阵乘法形式，从而推导出F：
简化一下：
跟 都是一个3x3的矩阵，根据上面公式可以推导出：
这里是其中一个四面体的F，其他四面体也会有对应的F，根据F我们就可以计算Green Strain Tensor，这也是每个四面体一份的：
以及势能 ，对势能进行求导（具体来说，consititutive model有公式的话，可以手动推导其弹性力的表达式），我们就能得到internal force（弹性力）：
这里，势能是物体共有的一个参数，是一个标量，而 则是物体上所有四面体的顶点组成的一个向量 ，那么经过求导后得到的弹性力也自然是一个n维的向量，每一维代表的是每个顶点的受力。
上面是内部的力，外面来自于重力或者碰撞的力，我们用外力 表示。
其中的质量m是一个很大的vector，在上面的公式下，我们就可以用之前说过的implicit euler或者explicit euler完成时间步长积分，从而算出更新的位置向量。
上面这种用离散的四面体来模拟软体的弹性形变的方法叫做Finite Element Method（FEM），这种方法可以通过不断提高分辨率（缩小四面体的尺寸）来保证准确性，不过同时也会导致计算复杂度增加，所以游戏或者很多电影特效模拟不是用的这种方法，而是基于优化的方法，Optimization-Based Methods，这种方法可以实现快速模拟。
先来回顾一下隐式欧拉方法，速度跟位置向量的演变可以用下面公式来表达：
这里的两个向量都是包含了多个四面体的相关数据的，是一个3xn的矩阵。
前面说过隐式欧拉最终会变成一个高维方程，很难求解，这里考虑将之转换成一个优化问题。
首先，internal force是势能的导数给出如下：
根据上面的公式，我们有：
这边希望构造出一个函数 ，并且使得这个函数的极值：
达成的条件跟前面的公式相一致，也就是说，希望这个函数的极值对应的解 就正好是上面一个等式的解，而这个函数可以给出为（比较容易推导？）：
其中
从而我们可以知道：
前面说到，这个函数提出的目的是为了求得前面implicit timestep等式的等价解，这里的一个疑惑就是，同样是解方程，前面implicit方程不好解，这里的方程同样也不好解，那么这种变化的意义在哪里呢？
这里的作用在于公式中的后面部分 ，这个部分表示的是系统的势能，前面说过弹性形变的势能公式 是通过一系列的物理推导得出的，是精确解，这种求解起来会比较麻烦，但是在游戏中，我们可以使用非精确解，只要保证模拟的结果看起来真实就好，这样一来，我们就能够选择一些能够让这个求解过程变得容易的公式代入这里的势能项中，从而达到简化问题的目的。
那么这个势能公式要怎么来设计呢？我们知道，弹性形变导致的势能总体来说，可以分成两种，分别是体积变化导致的形变（如从体积为1变成体积为2），或者是扭曲导致的形变，如一个cube变成平行六面体等。
前者势能可以表示为：
det是矩阵的行列式，这个公式的意义在于，当没有形变的时候，矩阵就是单位矩阵，这个势能的结果就是0
后者势能可以表示为：
tr（trace）是矩阵的迹，一个n×n矩阵A的主对角线（从左上方至右下方的对角线）上各个元素的总和被称矩阵A的迹（也就等于所有特征值的和）。同样，当没有形变时，这个F就是单位矩阵，所以迹就是3，那么势能依然为0。
总结可以看出，我们需要的势能公式其实就是一个约束条件，这个约束条件或者说约束函数在无形变的时候输出为0，在有形变的时候输出的是一个正值（或者说在无形变的时候取得极小值0），比如其形式为：
在这个条件下，我们就可以定义势能为：
上式中的A是一个对角阵，C(q)是一个向量，对于其中的每个元素：
这里希望有：
这里的i用来指示四面体的序号，这里翻译一下，就是说我们希望每个四面体的形变梯度矩阵的行列式为1，也就是说每个四面体都是体积不变的。
下面来介绍一下Position-based Dynamics方法，回到前面的公式：
这里，我们将势能用脉冲函数（实际上是脉冲函数取反，不过为了描述方便，这里用脉冲函数来表示）来表示：
其中脉冲函数：
之所以这样设计，就是为了满足前面说过的，我们希望函数J取得最小值，那么就需要保证 ，这里的i依然是四面体的序号。
继续下去，Position-based Dynamics方法会采用Gauss-Siedel Iteration（假设1）进行求解，即对于后面的众多四面体而言，我们会先只考虑一个四面体作用下的求解，并将结果代入原公式之后考虑第二个四面体作用下的求解，循环往复。如果只考虑单个四面体作用，原公式就变成了：
再定义：
公式就变成：
将公式转换一下，我们想要求解J的最小值，就等价于在条件：
下求解
的最小值。
之后，再假设 足够小（假设2），那么前面的条件就可以做泰勒展开为：
而这个就可以用上一讲的拉格朗日乘子来求解，也就是说，在上面条件下的最小值求解就变成了：
而Gauss-siedel Iteration的问题是，需要进行多次迭代才能收敛，从而导致计算比较慢，性能较差，这也是Position-based Dynamics（PBD）的缺点，如果Iteration数目不够就会导致一些模拟效果质量较低，比如布料会看起来像弹簧一样生硬。
我们来回顾下之前的内容：
那么很自然的问题就是，是否有一种方法是居于两者之间的，一方面我们不需要求解一个十分复杂的方程，而是只需要求解一个简单的方程，且这个方程可以通过少数的迭代就能输出结果的。
这个方法就是物理界常用的Projective Dynamics，实际上这种方法跟PBD是处于同一框架下得，不同的是PBD这里的internel energy函数是一个脉冲函数之和，而Projective Dynamics方法的internal energy函数则是：
这个式子中一共引入了两个新的符号，一个是变量p，这个变量代表的是某个粒子（如前面的四面体）的一个目标状态（比如，我们可以选择四面体不发生形变时的状态为一个目标状态，当然我们也可以选择其他状态作为目标状态，具体如何选，后面再说），而d则是表达p跟q差异性的函数，可以是一个非常简单的函数，比如就取两者的距离作为函数。
整个函数是什么意思呢，我们这里来举个例子说明一下，假如我们选择某个函数（simple distance measure）作为上面的函数d：
这里的F跟前面一样，还是deformation gradient， 是一个没有形变，但是经过一个旋转的状态，两者都是3x3的矩阵，其中 可以通过一个叫做Shape Matching的算法求得。
另外一项 是一个脉冲函数，这个脉冲函数的取值受到约束条件c的影响，当 的时候，取值为0，否则取值为无穷大，而由于这个函数的存在，要想使得整体的internal energy最小，就必须得满足这里的约束条件，从而使得脉冲函数取值为0。
在这个例子中，我们选取了p为经过某个旋转 后的状态，而在这个状态下，可以满足上面的脉冲函数结果为0，那么最终的internal energy的计算就变成了一个非常简单的矩阵的线性计算。
这里，我们还可以选择其他的p，比如：
这里的 也是一个deformation gradient，表示的是最接近于当前形变状态的deformation gradient，但是在这个gradient的作用下，形变后的四面体的体积跟原始四面体保持一致（这种对于动物形体的模拟比较有用，比如肚子在受击之后体积保持不变）， 的计算同样可以通过脉冲函数求得，比如只需要找到体积不变的形变的约束条件，就能得到对应的deformation gradient。
在这种internal energy的作用下，最优解的求解就变成了如下的函数形式：
可以看到，p只出现在后两项，q则只出现在前两项，从而我们可以通过有限的迭代来求解这个最优化问题。
而由于脉冲函数的存在，p的求解只需要保证脉冲函数为0即可，因此整个问题可以快速变成q变量下的最优解，这个过程叫做local projection，即找到限制条件c（constraint）来保证脉冲函数为0，从而将问题从两个变量降维成一个变量。
在求得p之后，就可以来求解q以获得整个函数的最优解：
而由于上面的函数实际上都是q的一个二次函数，而其导数则是q的一次函数，从而只需要将导数等于0，就能得到最终的解，这一步叫做global solve。
前面说过，projective Dynamics是一个迭代的过程，上面展示的local projection跟global solve就是这个过程中的仅有两个步骤，这两个步骤需要循环往复的去执行（因为p是local projection的，也就是一个局部解，这个解并不一定就对应着最优解，而是需要不断调整q来逼近最优解）。
这个方法的迭代数目比PBD会少很多，目前UE、Houndini等软件使用的是X-PBD方案，是在PBD方案上的一种更新方案，但是其迭代次数依然比Projective Dynamics方法要多，之所以不用Projective Dynamics方法，是因为线性系统的求解在维数很大（比如q是一个很大的向量）的时候会存在比较大的困难。
~
#酆茜丽# MATLAB怎么进行数学建模? -
(15782337865)： 一、数学建模的一般步骤 数学建模并不是新东西,粗略地说, 数学建模是一个多次迭代的过程,每一次 迭代大体上包括:实际问题的抽象、简化, 做出假设,明确变量和参数;形成明确的 数学问题;以解析形式或者数值形式求解 该数学模型...
#酆茜丽# 什么是排除假设法 (建模方法) -
(15782337865)： 解: 这个方法就是排除假设法 (建模方法) 先把不可能的情况排除在外 再把可能的情况进行假设,再排除 最后得到正确的答案.
#酆茜丽# 有数字0,1,2,3,4,5,(1)可组成多少个没有重复数字的六?
(15782337865)： (1) 六个数字全排列:P(6,6)= 6! = 720个但第一位数字为0的不行,这种数字有P(5,5)= 120个所以能排出720-120 = 600个(2)全排列(包括第一个数字为0)的720个数,...
