GAMES103笔记 Lecture 7、8 有限元方法（Finite Element Method）

力学概念汇总

正式开始前先把模拟中涉及到的物理概念汇总放在这，没学过力学的同学直接看公式肯定一脸懵逼，可以结合 PPT 来看，看的过程中对概念有疑惑可回来做个速查。

研究连续介质力学需要使用张量语言，本文中尽量使用矩阵而不是张量（张量在某一组基向量下的投影），对物理概念也不做详细的讲解，如果想进一步学习公式对应的物理本质需学习张量计算和连续介质力学的相关资料。

形变 Deformation

变形梯度 Deformation Gradient 形变后位置对形变前位置的梯度矩阵 $\mathbf{F}$ ，包含形变和一次旋转

应变 Strain Tensor 形变的定量描述，本文中使用格林应变 $\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_{uv} & \varepsilon_{vv} \end{bmatrix}$

格林应变是对称的，$\varepsilon_{uu}$ 和 $\varepsilon_{vv}$ 表示沿坐标轴方向的正应变（Normal Strain）， $\varepsilon_{uv}$ 表示非坐标轴方向的切应变（Shear Strain）

应变能量密度函数 Strain Energy Density Function $W(\mathbf{G})$ ，形变时单位空间的能量，是应变的函数，直接决定材料的性质，不同材料的 $W(\mathbf{G})$ 不同，根据定义有 $E=\int W(\mathbf{G})dA$

应力 Stress Tensor 单位空间所承受的作用力，单位空间的能量对应变求梯度得到应力（类似于力乘位移等于能量，这里是单位空间的力乘应变等于单位空间的能量）

第二皮奥拉-基尔霍夫应力（Second Piola–Kirchhoff Stress Tensor） $\mathbf{S}=\frac{\partial W(\mathbf{G})}{\partial \mathbf{G}}=2\mu\mathbf{G}+\lambda \mathop{trace}(\mathbf{G})\mathbf{I}$

**Traction Vector，**Traction vector 表示力的密度，单位和压强一样，在二维空间中为单位长度的力，在三维空间中为单位面积的力。应力乘法向量等于 traction vector，$\mathbf{T}=\mathbf{SN}$， $\mathbf{T}$ 和 $\mathbf{N}$ 分别是形变前的 traction vector 和法向量。

如果更换为形变后的法向量和形变后的梯度，相应的应力 $\mathbf{\sigma}$ 称为柯西应力（Cauchy Stress Tensor），又称真实应力（True Stress Tensor） $\mathbf{t} =\mathbf{\sigma} \mathbf{n}$

此外还有第一皮奥拉-基尔霍夫应力（First Piola–Kirchhoff Stress Tensor） $\mathbf{P}$ ，对应形变前的法向量和形变后的 traction vector $\mathbf{t} =\mathbf{P} \mathbf{N}$

FEM

有限元方法（Finite Element Method）把物体分割一个个有体积的单元来模拟。

特别地，线性有限元方法在二维空间中把物体分割成三角形，在三维空间中把物体分割成四面体。

接下来讨论二维空间上的有限元方法。

称三角形静止时的状态为 Reference 状态。

假设三角形内部的形变是均匀的，三角形内任意点 $\mathbf{X}$ 形变以后的位置 $\mathbf{x} = \mathbf{FX+c}$

注意大写表示形变之前的状态，小写表示形变之后的状态。

则任意向量 $\mathbf{X}_{ba}$ 形变后变成 $\mathbf{x}_{ba}=\mathbf{FX}_{ba}$

矩阵 $\mathbf{F}=\frac{\partial\mathbf{x}}{\partial\mathbf{X}}$ 称为变形梯度（Deformation Gradient）

用形变前后三角形的两个边向量就能求出 $\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}$ 蕴含了对形变的表示，但也包含了额外的刚体运动。

做奇异值分解 $\mathbf{F}=\mathbf{U}\mathbf{D}\mathbf{V}^T$ ，之前提过奇异值分解得到的三个矩阵可以从右到左看做是依次生效的的旋转、缩放、旋转， $\mathbf{V}^T$ 指示缩放的角度， $\mathbf{D}$ 指示每个维度的缩放量， $\mathbf{D}\mathbf{V}^T$ 合起来描述了物体的形变，最后的 $\mathbf{U}$ 只是再转一次物体，去掉 $\mathbf{U}$ 不影响物体的形状。

如果每次都做奇异值分解来去掉 $\mathbf{U}$ 还是比较麻烦，考虑用 $\mathbf{F}^T$ 乘 $\mathbf{F}$ 消去 $\mathbf{U}$

Green-Lagrange Strain（格林应变）：$\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})$

应变（Strain）定量地描述了形变

没有形变时， $\mathbf{G}=\mathbf{0}$ ，形变增加时， $\|\mathbf{G}\|$ 也会变大

根据定义容易得知 $\mathbf{G}$ 是一个对称矩阵，记为 $\begin{bmatrix} \varepsilon_{uu} & \varepsilon_{uv} \\ \varepsilon_{uv} & \varepsilon_{vv} \end{bmatrix}$

引入应变能量密度函数 Strain Energy Density Function $W(\mathbf{G})$ ，表示形变时单位空间的能量，是应变的函数，根据定义有 $E=\int W(\mathbf{G})dA$

$W(\mathbf{G})$ 直接决定材料的性质

用 Saint Venant-Kirchhoff (StVK) 模型继续推导公式

StVK 模型的 $W(\mathbf{G})=W(\varepsilon_{uu}, \varepsilon_{vv}, \varepsilon_{uv})=\frac{\lambda}{2}(\varepsilon_{uu}+\varepsilon_{vv})^2+\frac{\mu}{2}(\varepsilon_{uu}^2+\varepsilon_{vv}^2+\varepsilon_{uv}^2)$

类似于能量对位移求导得到力，单位空间的能量对应变求导得到单位空间所承受的作用力，称为应力（Stress Tensor）

注意这里都是相对于静止状态的，求导得到第二皮奥拉-基尔霍夫应力（Second Piola–Kirchhoff Stress Tensor）： $\mathbf{S}=\frac{\partial W(\mathbf{G})}{\partial \mathbf{G}}$

如果不是相对于静止状态的会得到其他形式的应力。

StVK 模型的第二皮奥拉-基尔霍夫应力 $\mathbf{S}=\frac{\partial W(\mathbf{G})}{\partial \mathbf{G}}=2\mu\mathbf{G}+\lambda \mathop{trace}(\mathbf{G})\mathbf{I}$ ，其中 $\mathop{trace}(\mathbf{G})$ 是矩阵的迹（三维空间中也是这个形式）

能量对位置求导得到力，经过一番推导最终得到这个公式可以一次算出顶点 $\mathbf{x}_1 $ 和 $\mathbf{x}_2$ 受到的力。

根据合力为0，顶点 $\mathbf{x}_0$ 收到的力 $\mathbf{f}_0=-\mathbf{f}_1-\mathbf{f}_2$

在三维情况下前面介绍的公式依然能用，只需要给各向量和矩阵各增加一个维度。

有了力之后用显式积分法或隐式积分法更新顶点的速度和位置完成模拟即可。

FVM

引入一个新的物理量 **Traction Vector。**Traction vector 在二维空间中为单位长度的力 $\mathbf{T}=\frac{d \mathbf{F}}{dl}$ ，在三维空间中为单位面积的力， $\mathbf{T}=\frac{d \mathbf{F}}{dA}$ 。Traction vector 表示力的密度，单位和压强一样。

FEM（Finite Volume Method） 用能量对顶点位置求导得到力，FVM 方法对 traction vector 积分得到力。

根据定义可以知道应力乘法向量等于 traction vector，$\mathbf{T}=\mathbf{SN}$， $\mathbf{T}$ 和 $\mathbf{N}$ 分别是形变前的 traction vector 和形变前的法向量

如果更换参考系，使用形变后的法向量和形变后的梯度，相应的应力 $\mathbf{\sigma}$ 称为柯西应力（Cauchy Stress）， $\mathbf{t} =\boldsymbol{\sigma} \mathbf{n}$

还有第一皮奥拉-基尔霍夫应力（First Piola–Kirchhoff stress） $\mathbf{P}$ ，对应形变前的法向量和形变后的 traction vector， $\mathbf{t} =\mathbf{P} \mathbf{N}$

这几个表达式的含义完全相同，只是更换了参考的状态

不同参考状态的应力之间相互转换的方法如上表所示

traction vector 是一个矢量，可以直接用变形梯度来转换，所以 $\mathbf{t}=\mathbf{FT}$

又因为 $\mathbf{t} =\mathbf{P} \mathbf{N}$ 和 $\mathbf{T}=\mathbf{SN}$

易得 $\mathbf{P}=\mathbf{FS}$

一个区域受到的合力等于在表面上 traction 的积分

有限体积方法假设应力在一个三角形内是常数

取所有三角形边的中点，用 $\mathbf{x}_0$ 代表周围一圈三角形的边中点连线形成的小块

则绿色三角形对 $\mathbf{x}_0$ 的力 $\mathbf{f}_0=\oint_L\mathbf{t}dl=\oint_L \boldsymbol{\sigma}\mathbf{n}dl$

又因为 $\oint_L \boldsymbol{\sigma}\mathbf{n}dl+ \oint_{L_{20}} \boldsymbol{\sigma}\mathbf{n}dl+ \oint_{L_{10}} \boldsymbol{\sigma}\mathbf{n}dl= \boldsymbol{\sigma} (\oint_{L} \mathbf{n}dl+ \oint_{L_{20}} \mathbf{n}dl+ \oint_{L_{10}} \mathbf{n}dl)=0$ （封闭曲线上对法向量积分等于0，这里 $L_{20}$ 表示 $\mathbf{x}_0$ 和 $\mathbf{x}_2$ 的中点到 $\mathbf{x}_0$ ， $L_{10}$ 表示 $\mathbf{x}_0$ 和 $\mathbf{x}_1$ 的中点到 $\mathbf{x}_0$）

所以中点间的连线形状其实无所谓，最后都有 $\mathbf{f}_0=-\oint_{L_{20}} \boldsymbol{\sigma}\mathbf{n}dl- \oint_{L_{10}} \boldsymbol{\sigma}\mathbf{n}dl= -\boldsymbol{\sigma}(\frac{\|\mathbf{x}_{20}\|}{2}\mathbf{n}_{20}+\frac{\|\mathbf{x}_{10}\|}{2}\mathbf{n}_{10})$

类似地，在三维空间中取顶点相邻四面体上三角形面积的 $\frac{1}{3}$ 连起来形成一个小块

一个四面体对顶点的力 $\mathbf{f}_0=-\oint_{\Omega}\boldsymbol{\sigma}\mathbf{n}dA= -\boldsymbol{\sigma}(\frac{A_{012}}{3}\mathbf{n}_{012}+\frac{A_{023}}{3}\mathbf{n}_{023}+\frac{A_{031}}{3}\mathbf{n}_{031})$

把三角形面积和法向量改成用边的叉乘来表示，化简后

$\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})$

再把柯西应力换成第一皮奥拉-基尔霍夫应力，法向量也换成静止状态的法向量，traction vector 依旧是静止状态下的，积分得到的力不变

再转换成第二皮奥拉-基尔霍夫应力（前面说了 $\mathbf{P}=\mathbf{FS}$ 如何推导）

得到 $\mathbf{f}_1=-\frac{\mathbf{FS}}{6}\mathbf{b}_1$

其中 $\mathbf{b}_1=(\mathbf{X}_{10} \times \mathbf{X}_{20} + \mathbf{X}_{20} \times \mathbf{X}_{30} + \mathbf{X}_{30} \times \mathbf{X}_{10})$

记住 $\mathbf{f}_1$ 是与顶点 $\mathbf{x}_1$ 相邻的一个四面体对 $\mathbf{x}_1$ 的力，$\mathbf{S}$ 是应变能量密度函数的导数

这两张 PPT 依据叉乘和矩阵运算化简运算，把一个四面体上三个顶点的 $\mathbf{b_1}, \mathbf{b_2}, \mathbf{b_3}$ 一同求出来。

Vol 指四面体的体积。

三个顶点受到的力也可一同算出，剩下一个顶点只要根据合力为 $\mathbf{0}$ 对其他三个顶点受力相加取反即可。

至此由形变前后四面体顶点位置求该四面体对其四个顶点的力的整个流程已经打通，接下来用显式积分法更新顶点的速度和位置完成模拟即可。

发现 FVM 方法最后得到的公式和 FEM 完全相同。当划分的单元为三角形或四面体时，FEM 和 FVM 等效。

超弹性模型

前面提到应变能量密度函数 Strain Energy Density Function$W(\mathbf{G})$ 决定了材料的性质，最后再介绍除 StVK 外几种不同的模型。

注意 FVM 方法中 $\mathbf{f}_0$ 是对 traction vector 积分得到的，最终得到的公式右边也是在拿第一皮奥拉-基尔霍夫应力 $\mathbf{P}$ 乘法向量得到 traction vector。

而 $\mathbf{P}$ 由应力能量密度函数 $W (\mathbf{G})$ 决定， $\mathbf{G}$ 又由 $\mathbf{F}$ 得到

所以 $\mathbf{P}$ 也可以看作 $\mathbf{F}$ 的函数 $\mathbf{P} (\mathbf{F})$

再对 $\mathbf{F}$ 做奇异值分解

有 $\mathbf{t}=\mathbf{P}(\mathbf{U}\mathbf{D}\mathbf{V}^T)\mathbf{N}$

上图左边为应用 $\mathbf{D}\mathbf{V}^T$ 形变后的 traction vector

traction vector 是一个向量，所以可以直接用 $\mathbf{U}$ 再旋转一次得到上图右边最终的 $\mathbf{t}$

$\mathbf{U}\mathbf{P}(\mathbf{D}\mathbf{V}^T)\mathbf{N}=\mathbf{t}=\mathbf{P}(\mathbf{U}\mathbf{D}\mathbf{V}^T)\mathbf{N}$

所以 $\mathbf{P}(\mathbf{U}\mathbf{D}\mathbf{V}^T)=\mathbf{U}\mathbf{P}(\mathbf{D}\mathbf{V}^T)$ ，把 $\mathbf{U}$ 提到括号外边不影响应力 $\mathbf{P}$ 的值

如果材料的形变在各个方向上都是均匀的，则第一个旋转 $\mathbf{V}^T$ 也能提到外边

称这种材料为各向同性材料

对于各向同性材料，把 $\mathbf{P}(\mathbf{D})$ 的自变量写成对角矩阵 $\mathbf{D}$ 的特征值

对角矩阵 $\mathbf{D}$ 的特征值是 $\mathbf{F}$ 的奇异值，又称主拉伸（Principal Stretches）

记 $\mathbf{C}=\mathbf{D}^T\mathbf{D}$ ，称为柯西-格林变形张量

$I_\mathbf{C}=trace(\mathbf{C})=\lambda_0^2+\lambda_1^2+\lambda_2^2$

$II_\mathbf{C}=trace(\mathbf{C^2})=\lambda_0^4+\lambda_1^4+\lambda_4^2$

$III_\mathbf{C}=\frac{1}{2}(trace^2(\mathbf{C})-trace(\mathbf{C^2}))=\lambda_0^2\lambda_1^2+\lambda_0^2\lambda_2^2+\lambda_1^2\lambda_2^2$

则四个弹性体的模型可以记为：

绿框部分代表物体抵抗被拉伸的趋势，蓝色部分代表物体抵抗体积改变的趋势

Fung 模型用来模拟人体组织

可用主拉伸来表示第一皮奥拉-基尔霍夫应力 $\mathbf{P} $ ：

这个方法和之前 $\mathbf{P}=\mathbf{FS}=\mathbf{P}\frac{\partial W}{\partial\mathbf{G}}$ 的方法是等效的

StVK 和 Neo Hookean 模型的应力随形变变化趋势：

StVK 的缺点：

• 压缩的时候，压缩幅度变大，抵抗力反而变小
• 四面体被反转（手性改变）之后，力也被反转了

Neo Hookean 没有这些问题

用四个模型做模拟的表现：

上传视频封面

好的标题可以获得更多的推荐及关注者

拉伸时中间会凹下去的现象称为 Poison Effect。

总结

本篇介绍了有限元方法和有限体积方法。有限元方法由能量对位置求导得到力，有限体积方法对 traction vector 积分得到力。当划分的单元为三角形或四面体时，有限元方法和有限体积方法等效。

后面又介绍了四种各向同性材料的应变能量密度函数。对于各向同性材料，其应变能量密度函数可以通过主拉伸（变形梯度的奇异值）来表示，可以通过对主拉伸求导的方式算出应力。

行文过程中花了大量的篇幅做公式推导，其实不关注具体的推导过程也没问题。建议学习时注意把握整体上的思路。当心过于陷入局部细节。经过复杂的推导得到的最终公式其实形式上非常简单，只要理解了关键性的概念，代码写起来不难（但需要调参 ╮(╯﹏╰）╭）。

得到力之后还需要用显式积分或隐式积分的方法来更新速度和位置，显式积分法容易发生数值不稳定模拟的物体直接炸掉的情况，可以通过平滑方法让模拟更稳定一些。隐式积分法的公式本篇没有推导，GAMES 103 有相应的 PPT，有兴趣的同学可以自行查看。

相关论文

Volino et al. 2009. A simple approach to nonlinear tensile stiffness for accurate cloth simulation. TOG

Teran et al. 2003. Finite Volume Methods for the Simulation of Skeleton Muscles. SCA.

Xu et al. 2015. Nonlinear Material Design Using Principal Stretches. TOG (SIGGRAPH).

Reference

1. Continuum Mechanics. https://www.continuummechanics.org

GAMES103课程主页：GAMES103：基于物理的计算机动画入门 - 计算机图形学与混合现实研讨会

上一篇：GAMES103笔记 Lecture 6 基于约束的布料模拟方法（Constrained Approaches）

下一篇：GAMES103笔记 Lecture 9 碰撞处理（Collision Handling）

GAMES103 系列笔记目录