GAMES103课程笔记07-Finite Element Method I
这个系列是GAMES103：基于物理的计算机动画入门(GAMES 103: Intro to Physics Based Animation)的同步课程笔记。本课程是计算机动画的入门课程，着重介绍各种基于物理的动画仿真模拟技术。本节课主要介绍有限元在物理仿真中的应用。
Linear Finite Element Method
Deformation Gradient
线性有限元(linear FEM)是有限单元法(FEM)中最基本的形式。在linear FEM中我们使用三角形网格来表示物体，并且假设三角形内部任意点的运动可以使用一个线性变换来表示：
\[\mathbf{x} = \mathbf{F} \mathbf{X} + \mathbf{c}\]式中\(\mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}}\)称为变形梯度(deformation gradient)。同时不难发现对于三角形内部的任意两点\(\mathbf{X}_a\)、\(\mathbf{X}_b\)，它们构成的矢量在变形前后满足：
\[\mathbf{x}_{ba} = \mathbf{x}_b - \mathbf{x}_a = (\mathbf{F} \mathbf{X}_b + \mathbf{c}) - (\mathbf{F} \mathbf{X}_a + \mathbf{c}) = \mathbf{F} \mathbf{X}_{ba}\]
因此我们可以利用三角形的两条边来计算变形梯度：
\[\left\{ \begin{aligned} & \mathbf{F} \mathbf{X}_{10} = \mathbf{x}_{10} \\ & \mathbf{F} \mathbf{X}_{20} = \mathbf{x}_{20} \end{aligned} \right. \Rightarrow \mathbf{F} \begin{bmatrix} \mathbf{X}_{10} & \mathbf{X}_{20} \end{bmatrix} = \begin{bmatrix} \mathbf{x}_{10} & \mathbf{x}_{20} \end{bmatrix}\] \[\mathbf{F} = \begin{bmatrix} \mathbf{x}_{10} & \mathbf{x}_{20} \end{bmatrix} \begin{bmatrix} \mathbf{X}_{10} & \mathbf{X}_{20} \end{bmatrix}^{-1}\]这样我们就可以利用三角形在变形前后的坐标来计算变形梯度。但需要注意的一点是这样计算出的变形梯度实际上同时包含了变形和旋转：即使三角形只进行了旋转而没有发生变形，我们仍然可以计算出一个变形梯度，此时\(\mathbf{F}\)等价于旋转矩阵。
Green Strain
显然我们希望求解出的变形能够与刚体变换无关。因此可以对\(\mathbf{F}\)进行SVD来分离出旋转：
\[\mathbf{F} = \mathbf{U} \mathbf{D} \mathbf{V}^T\]其中\(\mathbf{V}^T\)和\(\mathbf{D}\)控制了变形的方向和大小，而\(\mathbf{U}\)为变形后的旋转与变形本身无关。
因此我们希望可以从\(\mathbf{F}\)中分离出\(\mathbf{U}\)。直接的方法是对\(\mathbf{F}\)进行SVD，但更常用的做法是利用\(\mathbf{F}\)构造出一个与\(\mathbf{U}\)无关的矩阵：
\[\mathbf{G} = \frac{1}{2} (\mathbf{F}^T \mathbf{F} - \mathbf{I}) = \frac{1}{2} (\mathbf{V} \mathbf{D}^2 \mathbf{V}^T - \mathbf{I}) = \begin{bmatrix} \varepsilon_{uu} & \varepsilon_{uv} \\ \varepsilon_{vu} & \varepsilon_{vv} \end{bmatrix}\]上式中\(\mathbf{G}\)称为Green应变张量(Green strain)。显然在没有变形的情况下，Green应变张量为0；变形越大，Green应变张量越大。因此Green应变张量描述了物体的变形而且与旋转无关。同时注意到Green应变张量是一个对称阵，因此它只具有3个自由度。
Strain Energy Density Function
有了变形后就可以定义系统的能量。我们定义面积微元上的变形能量为能量密度(energy density)\(W(\mathbf{G})\)，则三角形上的能量等于能量密度的积分：
\[E = \int W(\mathbf{G}) d A\]对于linear FEM整个三角形上的变形是相同的，因此三角形的能量等于能量密度乘以三角形的面积：
\[E = \int W(\mathbf{G}) d A = A^{\text{ref}} W(\varepsilon_{uu}, \varepsilon_{uv}, \varepsilon_{vv})\]在图形学中一般使用Saint Venant-Kirchhoff Model(StVK)来计算能量密度：
\[W(\varepsilon_{uu}, \varepsilon_{uv}, \varepsilon_{vv}) = \frac{\lambda}{2} (\varepsilon_{uu} + \varepsilon_{vv})^2 + \mu (\varepsilon_{uu}^2 + \varepsilon_{vv}^2 + 2 \varepsilon_{uv}^2)\]其中\(\lambda\)和\(\mu\)称为拉梅参数(Lamé parameters)。
对能量密度求导可以发现：
\[\begin{aligned} \frac{\partial W}{\partial \mathbf{G}} &= \begin{bmatrix} \frac{\partial W}{\partial \varepsilon_{uu}} & \frac{\partial W}{\partial \varepsilon_{uv}} \\ \frac{\partial W}{\partial \varepsilon_{vu}} & \frac{\partial W}{\partial \varepsilon_{vv}} \\ \end{bmatrix} \\ &= \begin{bmatrix} 2 \mu \varepsilon_{uu} + \lambda (\varepsilon_{uu} + \varepsilon_{vv}) & 2 \mu \varepsilon_{uv} \\ 2 \mu \varepsilon_{vu} & 2 \mu \varepsilon_{vv} + \lambda (\varepsilon_{uu} + \varepsilon_{vv}) \\ \end{bmatrix} \\ &= 2 \mu \mathbf{G} + \lambda \text{trace}(\mathbf{G}) \mathbf{I} \\ &= \mathbf{S} \end{aligned}\]其中\(\mathbf{S}\)称为第二类Piola-Kirchhoff应力张量(second Piola-Kirchhoff stress tensor)。
Forces
通过对能量求导就可以计算出作用在每个顶点上的力：
\[\mathbf{f}_i = -\frac{\partial E}{\partial \mathbf{x}_i} = -A^{\text{ref}} \ \frac{\partial W}{\partial \mathbf{x}_i} = -A^{\text{ref}} \ \bigg(\frac{\partial W}{\partial \varepsilon_{uu}} \frac{\partial \varepsilon_{uu}}{\partial \mathbf{x}_i} + \frac{\partial W}{\partial \varepsilon_{vv}} \frac{\partial \varepsilon_{vv}}{\partial \mathbf{x}_i} + 2 \frac{\partial W}{\partial \varepsilon_{uv}} \frac{\partial \varepsilon_{uv}}{\partial \mathbf{x}_i} \bigg)\]其中\(\frac{\partial W}{\partial \varepsilon}\)项已经在\(\mathbf{S}\)中，因此我们只需要再计算应变对位置的导数。通过繁琐的数学推导可以得到：
最后整理得到顶点上的力：
\[\begin{bmatrix} \mathbf{f}_1 & \mathbf{f}_2 \end{bmatrix} = -A^{\text{ref}} \ \mathbf{F} \ \mathbf{S} \begin{bmatrix} \mathbf{r}_{10} & \mathbf{r}_{20} \end{bmatrix}^{-1}\]利用单元上合力为0可以得到\(\mathbf{x}_0\)顶点上的力：
\[\mathbf{f}_0 = - (\mathbf{f}_1 + \mathbf{f}_2)\]Finite Volume Method
除了FEM外在图形学中还经常使用有限体积法(finite volume method, FVM)来求解弹性体的变形问题，而且对于线性单元FVM和FEM实际上是等价的。
Traction and Stress
首先我们定义内力在单位长度(面积)上的分布为traction，记为\(\mathbf{t}\)。这样作用在物体上的力可以表示成traction的积分：
\[\mathbf{f} = \int_L \mathbf{t} \ d l\]同时，traction可以表示为应力张量\(\boldsymbol{\sigma}\)与法向的内积：
\[\mathbf{t} = \boldsymbol{\sigma} \mathbf{n}\]因此我们可以得到应力张量表达的内力：
\[\mathbf{f} = \int_L \boldsymbol{\sigma} \mathbf{n} \ d l\]
不难发现，在FVM中我们更多地是从积分的角度来看待物体的力与变形。对于三角网格，我们定义每个面片对于顶点力的贡献为：
\[\mathbf{f}_0 = \int_L \boldsymbol{\sigma} \mathbf{n} \ d l\]其中\(L\)为顶点在面片上对应区域的边界。为了保证三角形对三个顶点的贡献是相同的，我们规定每个顶点对应的区域都要通过它相邻边的中点。对于线性单元，其应力张量在三角形内都是相同的，因此根据散度定理有：
\[\int_L \boldsymbol{\sigma} \mathbf{n} \ d l + \int_{L_{20}} \boldsymbol{\sigma} \mathbf{n}_{20} \ d l + \int_{L_{10}} \boldsymbol{\sigma} \mathbf{n}_{10} \ d l = \mathbf{0}\]这样就得到了边界\(L\)对\(\mathbf{x}_0\)节点力的贡献：
\[\begin{aligned} \mathbf{f}_0 &= - \bigg( \int_{L_{20}} \boldsymbol{\sigma} \mathbf{n}_{20} \ d l + \int_{L_{10}} \boldsymbol{\sigma} \mathbf{n}_{10} \ d l \bigg) \\ &= -\boldsymbol{\sigma} \bigg( \frac{\Vert \mathbf{x}_{20} \Vert}{2} \mathbf{n}_{20} + \frac{\Vert \mathbf{x}_{10} \Vert}{2} \mathbf{n}_{10} \bigg) \end{aligned}\]
对于三维的情况只需要把曲线边界推广成面积即可：
\[\mathbf{f}_0 = -\frac{\boldsymbol{\sigma}}{6} ( \mathbf{x}_{10} \times \mathbf{x}_{20} + \mathbf{x}_{20} \times \mathbf{x}_{30} + \mathbf{x}_{30} \times \mathbf{x}_{10} )\]
Different Stresses
从上面的推导可以发现计算顶点力的核心是应力张量\(\boldsymbol{\sigma}\)，但需要说明的是这里的应力张量和FEM中的应力张量是完全不同的概念。在FEM中单元的几何信息都是定义在参考系(变形前)上的，而在FVM中则是定义在当前系(变形后)的。
实际上根据参考系的不同应力张量有不同的定义，在FEM中使用的是第二类Piola-Kirchhoff应力张量，而在FVM中则使用的是柯西应力张量(Cauchy stress)。还有第一类Piola-Kirchhoff应力张量(first Piola-Kirchhoff stress)，它表示法向在变形前但traction在变形后的应力张量：
\[\mathbf{P} = \mathbf{F} \mathbf{S}\]
第一类Piola-Kirchhoff应力张量和第二类Piola-Kirchhoff应力张量之间只相差一个变形梯度矩阵，而第一类Piola-Kirchhoff应力张量和柯西应力张量之间的变换关系可以推导如下：
总结一下，不同应力张量之间的转换关系可参见下表：
这样我们就可以通过前面计算出的\(\mathbf{S}\)来计算柯西应力张量\(\boldsymbol{\sigma}\)，进而计算出节点力。根据节点力公式，顶点\(\mathbf{x}_1\)处的力可以表示为：
\[\begin{aligned} \mathbf{f}_1 &= -\frac{\boldsymbol{\sigma}}{6} ( \mathbf{x}_{01} \times \mathbf{x}_{21} + \mathbf{x}_{21} \times \mathbf{x}_{31} + \mathbf{x}_{31} \times \mathbf{x}_{01} ) \\ &= -\frac{\mathbf{P}}{6} ( \mathbf{X}_{01} \times \mathbf{X}_{21} + \mathbf{X}_{21} \times \mathbf{X}_{31} + \mathbf{X}_{31} \times \mathbf{X}_{01} ) \\ &= -\frac{\mathbf{F} \mathbf{S}}{6} \mathbf{b}_1 \end{aligned}\]其中\(\mathbf{b} = \mathbf{X}_{01} \times \mathbf{X}_{21} + \mathbf{X}_{21} \times \mathbf{X}_{31} + \mathbf{X}_{31} \times \mathbf{X}_{01}\)为参考系下的向量叉乘，因此是固定值可以事先计算出来。
接下来考虑\(\mathbf{b}_1\)叉乘的性质：
整理得到：
最终FVM的求解过程如下：
Hyperelastic Models
First Piola–Kirchhoff Stress
我们可以把第一类Piola-Kirchhoff应力张量\(\mathbf{P}\)看做是关于变形梯度\(\mathbf{F}\)的函数，它把参考系下的法向\(\mathbf{N}\)映射成当前系下的traction：
\[\mathbf{t} = \mathbf{P}(\mathbf{U} \mathbf{D} \mathbf{V}^T) \mathbf{N}\]
同时应力张量\(\mathbf{P}\)具有旋转不变性，显然变形后的旋转不会改变traction：
\[\mathbf{P}(\mathbf{U} \mathbf{D} \mathbf{V}^T) = \mathbf{U}(\mathbf{P} \mathbf{D} \mathbf{V}^T)\]
Isotropic Materials
假设弹性体是各向同性的(isotropic)，我们还可以进一步把变形前的旋转分离出去：
\[\mathbf{P}(\mathbf{D} \mathbf{V}^T) = \mathbf{P} (\mathbf{D}) \mathbf{V}^T\]
因此\(\mathbf{P}\)对\(\mathbf{F}\)的作用可以表示为：
\[\mathbf{P} (\mathbf{F}) = \mathbf{P} (\mathbf{U} \mathbf{D} \mathbf{V}^T) = \mathbf{U} \mathbf{P}(\mathbf{D}) \mathbf{V}^T = \mathbf{U} \mathbf{P}(\lambda_0, \lambda_1, \lambda_2) \mathbf{V}^T\]其中\(\lambda_0\)、\(\lambda_1\)、\(\lambda_2\)为变形梯度的三个奇异值，表示在三个坐标轴上的伸展。在很多文献中还使用其它的方式来表示\(\mathbf{P}(I_\mathbf{C}, {II}_\mathbf{C}, {III}_\mathbf{C})\)，其中：
\[I_\mathbf{C} = \text{trace} (\mathbf{C}) = \lambda_0^2 + \lambda_1^2 + \lambda_2^2\] \[{II}_\mathbf{C} = \text{trace} (\mathbf{C^2}) = \lambda_0^4 + \lambda_1^4 + \lambda_2^4\] \[{III}_\mathbf{C} = \frac{1}{2} (\text{trace}^2 (\mathbf{C}) - \text{trace} (\mathbf{C^2})) = \lambda_0^2 \lambda_1^2 + \lambda_0^2 \lambda_2^2 + \lambda_1^2 \lambda_2^2\]其中\(\mathbf{C} = \mathbf{F}^T \mathbf{F}\)称为右Cauchy-Green变形张量(right Cauchy-Green deformation tensor)。
我们可以使用\(I_\mathbf{C}, {II}_\mathbf{C}, {III}_\mathbf{C}\)来表示弹性体的能量密度，如StVK和neo-Hookean模型的能量密度公式如下：
其中第一项表示材料抵抗剪切的能力，而第二项则表示材料保持原有体积(面积)的能力。其它常用模型的能量密度计算公式如下：
当我们获得了能量密度后就可以通过对能量密度求导来计算应力张量：
\[\mathbf{P}(\lambda_0, \lambda_1, \lambda_2) = \begin{bmatrix} \frac{\partial W}{\partial \lambda_0} & & \\ & \frac{\partial W}{\partial \lambda_1} & \\ & & \frac{\partial W}{\partial \lambda_2} \end{bmatrix}\] \[\mathbf{P} = \mathbf{U} \mathbf{P}(\lambda_0, \lambda_1, \lambda_2) \mathbf{V}^T\]这样通用模型的弹性体节点力计算过程如下：
