太极图形学——弹性物体仿真 1

1.仿真基础，时间和空间的离散化

渲染对于物体的展现来说是一件很重要的事

但除了渲染之外，还需要物理来控制，形成动画

弹性物体的仿真实际上非常重要，特别是对于头发，衣物的仿真，虚拟手术中软组织的仿真

仿真无非就是求在不同的时间下，物体在什么位置，就像是下面这张图，经典的中学物理题

那么今天就会搭建一个可以形变的物体模型

那么我们需要的就是这个details是怎么实现的，在代码例程中包含的有质点弹簧和线性有限元模型

在时间和空间上做积分，重点将放在如何制作一个仿真器

要用到的一些基础的公式

对于线性材料而言，就是胡克定律，利用此定律可以求出力和加速度，很多时候，现实中的材料都不是线性材料，这种仿真一般都是用于声音的仿真

大部分材料都是非线性材料，弹性仿真无非是用力来更新速度，用速度来更新位置

2.时间积分

从直观的角度上理解，仿真我们所需要的是一个黑盒，这个黑盒的作用是若我们给黑盒输入一个物体当前的位置和速度，我们就能知道下一时刻的该物体的位置和速度，而要实现这一些就是黑盒所做的工作

那这个黑盒要怎么写呢，根据前面所写的公式，我们可以很容易的知道，速度是位置的导数，加速度是速度的导数，又根据所呈现在屏幕上的图像永远是离散化的，那么就有一个时间的步长，我们需要做到的是根据这一时间的位置和速度求出下一时间的位置和速度，结合牛顿第二定律，那么就有位置的变化是当前位置加上位置的变化量，也就是速度对于时间的积分，速度则是当前的速度加上速度的变化量，也就是加速度对于时间的积分，加速度可以由物体的质量和外力和来表示

这个公式的存在，有一个问题，首先就是积分问题本身就不好求解，这是其一，第二是我们事实上并不知道任何在tn，也就是当前时间之后的值，也就是一时刻的力和速度，这也同时就代表着我们无法使用类似蒙特卡洛的方法来进行这些估值计算

这样的话，我们就只能使用数值积分的方法，那么从这里开始就出现了两种经典的方法

方法一：前向欧拉算法（显式积分法）

因为此时我们的问题是不知道在当前时间之后下一时刻的任何值，那么我们就从加速度出发，因为加速度对于时间的积分为速度，且在积分的表现上就是两时刻之间，加速度（力）函数所代表的面积，那么我们就认为此时（当前时间到下一时刻）的力是恒定的，面积的求法直接使用当前时刻的加加速度（力）为高乘以时间的长度，就可以得到速度，那么传统的前向欧拉算法的过程就是使用当前的速度更新我的位移，再使用当前的加速度（力）来更新速度

但这个方法的问题在于，随着时间的增长，这个系统的能量会不断的膨胀，用不了很长的间，速度和位置的更新就会发散

这种前向欧拉算法很容易理解，也很容易实现，但由于能量的问题，很少会用在仿真上，现在我们讲的显示积分法都是先以力更新速度，再用更新过的速度来更新位置，这样做的好处就是整个系统的能量会处于一种振荡波动的状态，也可以认为是一种稳定的状态，保持了动量守恒

在高精度的仿真，或者天体仿真中，这种显示积分法一般是首选

方法二：后向欧拉法（隐式积分法）

所谓的隐式积分法，就是一个解方程的过程，直观的理解就是使用下一时刻的力作为不变的力来更新速度，再用速度来更新位移

隐式积分在使用的过程中，能量会一直减小，虽然会产生数值阻尼，但是，随着能量的减小，也代表着我们可以采用更大的时间步长，可以防止系统的爆炸

这种积分的使用方式一般用在对时间要求比较高的模拟系统中，例如游戏，动画，混合现实这些系统运用中

显式积分做法实际操作

整个对时间积分的显式方法的流水线可以被很简单的概括，首先若是对于弹性物体的模拟，我们需可以先根据能量求出力，再由力来更新速度，再由更新好的速度来更新位置，这个流程从理解上十分的容易

以之前的天体运动来作为例子，首先是根据能量求导得到力

最大的问题是，对于保守力来说，我们一直所需要的就是能量，用这个推导就可以得到之后的一切

但是对于弹性物体呢，弹性物体是一个连续介质物体，能量需要对整个空间的能量密度求积分，那么我们应该怎么对能量密度求解积分呢

空间积分：

需要做到的就是对空间进行积分，那么对于空间进行积分的方法主流的有两种，一种是质点弹簧系统，一种是线性有限元

在进行这两种方法的总结时，最重要的问题就两个，一个是怎样去描述形变的，一个是有了形变之后怎样去描述能量的

方法一：质点弹簧系统（MSM）

这是一种简单但却有效的方法，首先，对于一个连续介质的弹性物体，我们需要将其离散化成很多的顶点，每一个顶点具有质量，根据三角形或者体积的大小将整个物体的质量分配到这些顶点上，然后使用弹簧来连接这些顶点，在物体未形变的时候，弹簧会被回归原长

每一个连接顶点的边都被认为是一个弹簧，那么这种弹簧是怎样去描述形变的呢，很简单，根据弹簧本身的性质，两个质点之间的距离是弹簧的原长，而在形变后质点之间的距离是弹簧的现在的长度，两个长度之间的差值就代表了物体的形变

那么是怎么定义能量和力呢，也很简单，根据胡克定律就可以了，根据弹簧长度的变化就可以求出能量的变化

根据这个能量，就可以求倒数，根据能量对位置的求导就可以求出某一点所受到的力的大小

那么在空间中的积分是怎么做的呢，总能量对其中每一个点的位置进行求导就可以得到每一个点位置的力的大小

代码写起来也很直观，分为空间积分部分和时间积分部分，空间积分部分用于求梯度，也就是力的大小（方向相反），再在时间积分中首先计算加速度，再利用加速度更新速度，再用速度更新位置，使用了显式积分的更新方法

弹簧系统常常用于布料和头发的仿真当中，特别是这种一维的系统

但是在体积的更新方法中，很难使用弹簧系统来进行更新，因为很容易出现翻转

所以对于有面积和体积的物体，大量的研究还是倾向于使用线性有限元的方法

方法二：线性有限元（FEM）

对于形变表达的映射

这个映射我们称为形变图，这个map在很多的方面都很像渲染中的仿射变换，包括位移，旋转以及拉伸

使用微分的方法去思考的话

可以被描述为一个仿射变换，这个F就被称为形变梯度（constant matrix），这一块的内容在game103课程中也有详细的描述，在这里我们是认为在一个很小的范围里，可以认为这个很小的地方是一个仿射变换

举几个例子，首先是平动的例子，矩阵F就是映射对X求导，最后得到的就是一个单位矩阵

如果是旋转，也是一样的，F此时就是R

当然，对拉伸的映射也是

，其实从上图中就可以看到，这个形变梯度只和旋转以及拉伸有关，那么严格的定义中，对于某一个较小的位置，我们都可以认为形变后的位置是一个形变前位置的仿射变换，这个变换中，形变梯度是一个雅克比矩阵

对于能量，我们是怎样去理解的呢，我们希望任何一个非刚体的形变，希望它的形变梯度可以带来一个非0的能量，也就是说物体在做非刚体运动时，能量不置0

那么能量怎么定义呢，通过能量密度来定义，也就是在一个非常小的空间内能量是一个什么样子的

这个能量密度公式中，x可以由X的映射来表示，所以能量可以写成一个X的函数。又因为，事实上能量和平动没有关系，继续化简，x和X只是物体的初始或者静止位置，所以说，能量可以完全定义成只和形变梯度F有关系的东西

这样就完成了对能量密度的定义，一个关于F的函数

这个应该怎么具体来定义呢

需要这个形变的定义是纯粹的，具有旋转不变性和平移不变性，在没有形变的时候，能量需要为0，所以这时候就需要引入应变张量来描述形变，需要满足当为初始状况时，形变应该为0，但有旋转的时候，该张量是不变的，有哪些张量可以用来使用呢，需要找到一些本构模型，其中stvk在图形学中用的是最多的，满足了上述的条件

stvk model使用格林应变张量，这个的好处是使用了一个F的转置，因为如果此时有旋转的话，那么代入旋转矩阵，由于旋转举证的转置等于其逆矩阵，所以具有旋转不变性

此时应变就是用来描述形变的东西，那么有了对于形变的描述，接下来就是如何描述能量，能量的公式一般用第二个，也分为两个部分，第一个部分表示硬度，一个部分用来保体积

那么在能量密度定义后我们需要将其离散化，也就是知道总能量是多少，用能量密度对于整个物体的面积做一次积分

那么此时就轮到FEM登场对物体做离散化了、依然将整个物体分为很多的区块，使用FEM还有一个假设，就是在每一个区块内，所有的点都是平均的，做的是相同的线性映射，也就是F形变梯度是一个常数

这样就将其转换成了一个对空间求和的过程，每一个小的元素的能量，由于在每一块区域的内部，能量密度是一个常数（F不变，应变张量就不变，能量密度就不变），和此时的位置没有关系，就可以提取出来

有了线性元素的存在，可以直观的将F求出来，也就是说在已知形变前和形变后的位置后，能够很容易的求解出形变矩阵的大小

我们要求力的话，首先能量对于位移的导数是力，那么能量密度对应变张量的求导呢

最后得到的是第一皮尔科舍夫应力张量，所谓的应力张量就是在张量空间中，力的扰动对能量的变化，也就是这是张量空间下，力是怎么表示的

第一皮尔科舍夫应力张量

所以，第一皮尔科舍夫应力张量是怎么得来的呢，总的来说，首先需要得到的是位移前的位置X和位移后的位置x,因为已知的位置就知道了边的位置向量，再根据几条边的位置向量我们可以根据公式x = FX + t 得到F（形变梯度），得到形变梯度后，求解得到格林应变张量 Green strain，由格林应变张量我们可以结合空间积分得到的能量公式，结合stvk的本构方程，得到能量密度，再由能量密度对应变求导可以得到应力张量，此时的应力张量其实是第二皮尔科舍夫应力张量，还需要乘以一个F，才是第一