图形学中的有限元方法

Finite Element Method (FEM, 有限元方法) 是一种将物体看作许多有体积的微小单元来进行仿真的方法。 比如将二维图形拆分成若干三角形，将三维体拆分成若干四面体。 这种方法能够很好的模拟弹性体的形变等特点，能够在数学上证明其仿真方法一定收敛到解析解，比MPM等使用粒子模拟的方法更加准确，但同时也更加复杂，运算更慢。

本文将先从弹性体相关的物理量基础定义开始，再讲述各种弹性体的能量函数模型，最后阐述如何用有限元方法离散化进行仿真。

本文主要参考FEM simulation of 3D deformable solids Siggraph 2012 Course，并修正了一些不明确的可能错误的地方

\[\DeclareMathOperator{\tr}{tr}\]

弹性体基础

弹性体指的是能够自己恢复初始状态的材料。

超弹性体（hyperelasticity）是一种特殊的弹性体，指的是弹力与压缩（拉伸）的路径没有关系，只与压缩后的状态有关系。

本文只考虑超弹性体。

形变与能量

定义弹性体的区域为$\Omega$，弹性体上初始时每个点的位置为$X$。 弹性体发生平移、形变等变化后，弹性体上每一个点都将移动到另一个位置$\mathbf{x}$，可以用一个形变映射(deformation map)来表示弹性体的这一个变换$\mathbf{x} = \phi (X)$，如下图：

$\phi(X)$并不能表示形变，因为弹性体整体的移动不产生形变，而相近的点移动的区别产生形变，则可以使用形变映射函数到初始位置的导数作为表示形变的量： \(F = \frac {\partial \phi (X)} {\partial X} = \frac {\partial \mathbf{x}} {\partial X}\) 在三维情况下，这是一个3×3的矩阵，即雅可比矩阵，$F$被称为形变梯度(deformation gradient)。

$ F$形变梯度表示某一点$X$相对其邻域的形变，形变产生弹性势能，所以需要定义针对这一点的能量密度$\Psi(\phi, X)$，弹性体的总能量即为$E(\phi) = \int_{\Omega} \Psi(\phi, X) \mathrm{d} X$. 一般情况下，能量密度只与形变梯度有关，可写为$\Phi(F)$。

注意到$\text{det}(F) = J = \frac {V_{变换后}} {V_{初始}}$，即$F$的行列式可以表示体积的变化。

弹性内力

接下来计算弹性内力，在数学上需要区分内部力和边界力：

对于体积内部，令$f(X)$表示位置$X$处的力密度(单位$N/m^3$)，则一个体积$A \subset \Omega$内的受力为$\mathbf{f}_{agg}(A) = \int_A f(X) \mathrm{d}X$。

对于边界上，令$\tau(X)$表示边界位置$X$处的力密度(单位$N/m^2$)，则一块三维体表面$B \subset \partial \Omega$的受力为$\mathbf{f}_{agg}(B) = \oint_B \tau(X) \mathrm{d}X$。

内力即为弹性势能梯度的相反方向，因为内力总是往使得势能减小最快的方向进行，所以内力可以为$\mathbf{f} = - \frac {\partial E} {\partial \mathbf{x}}$，注意这里是变换后的位置。

PK1数

因为内力是势能的梯度相反方向，所以我们试图用总势能对位置求梯度，看看能得到什么。 \(\begin{align} \delta E & = \delta \left[ \int_\Omega \Psi(F) \mathrm{d}X \right] \\ & = \int_\Omega \delta \left[ \Psi(F) \right] \mathrm{d}X \\ & = \int_\Omega \left[ \frac {\partial \Psi} {\partial F} : \delta F\right] \mathrm{d}X \end{align}\)

公式中的$\left[A : B\right]$表示矩阵$A, B$按位点乘，上述公式是微分链式法则用到了点乘。 之后$[ : ]$均表示按位点乘，有时候还需要按位点乘后再将每个元素加起来，请根据上下文判断。

令$P = \frac {\partial \Psi} {\partial F}$，被称为第一皮奥拉-基尔霍夫应力张量(The first Piola-Kirchhoff stress tensor，简称PK1)，下面将推导$P$的几何意义。

继续推导能量的微分，直接利用斯托克斯公式定理（分部积分）： \(\begin{align} \delta E & = \int_\Omega \left[ P : \delta F\right] \mathrm{d}X \\ & = -\int_\Omega (\nabla \cdot P) \cdot \delta \phi(\mathbf X)\,\mathrm{d}X + \oint_{\partial \Omega} (P \cdot N) \delta \phi(\mathbf X) \, \mathrm{d}X\\ & = -\int_\Omega (\nabla \cdot P) \cdot \delta \mathbf{x} \,\mathrm{d}X + \oint_{\partial \Omega} (P \cdot N) \delta \mathbf{x} \, \mathrm{d}X \\ \end{align}\) 因为出现了对于当前位置$\mathbf{x}$的微分$\delta \mathbf{x} $，根据力的定义，可以把积分前半部分视为力密度。 即 \(f(\mathbf{x}) = \nabla \cdot P, \quad \tau(\mathbf{x}) = - P\cdot N\) 可以理解为PK1数表示的是对该点某一方向微平面的压强。

学习过力学的读者可能记得另一种应力张量：柯西应力张量(Cauchy stress tensor)。 实际上，柯西应力张量是在某一不变的参考坐标系下的应力张量（欧拉视点），而PK1应力则是在形变后材料的坐标系下的张量（拉格朗日视点），两者满足 \(P = \mathbf J F \sigma F^{-T}\) 其中$\mathbf J F$是$F$的雅可比行列式，$F$是形变梯度； 这个变换称为皮奥拉变换。 在无穷小形变的情况下，两个应力张量相等。

弹性能量模型

本节将探讨不同的弹性模型如何从形变梯度定义能量密度，计算出内力。

应变量度

为了计算出能量，我们需要一个衡量形变程度的度量方式，$F$只是一个表示形变的量，并不能衡量形变程度。

注意到形变梯度$F$实际上仅包含变换的旋转和拉伸，矩阵一定可以分解为$F=RS$，其中$R$为旋转矩阵，满足$R^T R=I$；$ S$为拉伸矩阵，一定为对称的。形变的量度可以仅与拉伸有关： \(E=\frac 1 2 (S^2-I) = \frac 1 2 (F^T F - I)\) 注意这里的$E$不是能量，这被称为格林应变张量(Green strain tensor)，是一个3×3的矩阵，用于衡量各个方向的形变程度。 这是一个与旋转、平移无关的量度，但它是非线性的。

尝试简化，我们对$E$在$F=I$进行泰勒展开： \(\begin{align} E(F) & ≈ E(I) + \frac {\partial E} {\partial F}\bigg |_{F=I} : (F-I) \\ & = E(I) +\frac 1 2 \left[ (F-I)^T I + I^T (F-I) \right] \\ & = \frac 1 2 (F+F^T) - I \end{align}\) 其中第一行到第二行的求法是使用微分 \(\frac {\partial E} {\partial F} : \delta F = \delta E = \frac 1 2 (\delta F^T F + F^T \delta F)\) 我们由此定义了另一个量度： \(\epsilon = \frac 1 2 (F+F^T) - I\) 称为微小应变张量(small strain tensor, infinitesimal strain tensor)，是一个线性的应变张量，但其却是与旋转有关的，弹性体在不形变的情况下发生旋转也会改变$\epsilon $。

$\epsilon$一定为对称矩阵，因此其点乘一个反对称矩阵一定为$\mathbf{0}$。

在静力学中，由于$\epsilon$的性质更好，且一般要求微小形变——没人希望自己房子的梁大幅形变——因此一般使用微小形变张量作为研究对象。 图形学中则常常处理大形变的问题，因此常常考虑使用其他形变度量。

线性弹性模型

定义能量密度为： \(\Psi(F) = \mu \tr(\epsilon^2) + \frac \lambda 2 \tr^2(\epsilon)\) 的模型为线性弹性模型（Linear elasticity），其中$\mu = \frac {k} {2(1+\nu)}$, $\lambda = \frac {k \nu} {(1 + \nu) ( 1-2 \nu)}$被称为拉梅系数 (Lamé coefficients) 。 其中的$k$为杨氏模量 (Young’s modulus) ，一般表示拉伸抵抗性。$\nu$为泊松比，一般表示不可压缩性。

通过对$\Psi$微分来求$P$ \(\begin{align} \delta \Psi & = \mu \cdot \tr(\delta \epsilon \cdot \epsilon +\epsilon \delta \epsilon) + \lambda \tr(\epsilon ) \tr(\delta \epsilon) \\ \delta \epsilon & = \frac 1 2 \left( \delta F + \delta F^T \right) \doteq \text{Sym}(\delta F) \ \end{align}\)

1. 任意矩阵都可以拆分为一个对称矩阵和反对称矩阵之和$A = \text{Sym}(A) + \text{Anti-Sym}(A)$ 对称矩阵点乘一个矩阵：$\epsilon : A = \epsilon : \text{Sym}(A)$
2. 矩阵的迹中两矩阵相乘是可以交换的：$\tr(AB) = \tr(BA)$ ，可以展开对角线元素计算证明。还可以是可转置的$\tr (A) = \tr (A^T) $
3. 矩阵的迹的微分可以视为矩阵对角元素的微分，其余为0，则可以写为$\delta \tr(\epsilon) = \tr(\delta \epsilon) = I : \delta \epsilon$

\[\begin{align} \delta \Psi & = 2\mu \cdot \tr(\delta \epsilon \cdot \epsilon ) + \lambda \tr (\epsilon) \left( I : \text{Sym}(\delta F) \right) \\ & = 2\mu \cdot (I : \text{Sym}(\delta F)\epsilon) + \lambda \tr (\epsilon) \left( I : \text{Sym}(\delta F) \right) \\ & = 2\mu \cdot (\epsilon : \text{Sym}(\delta F)^T) + \lambda \tr (\epsilon) \left( I : \text{Sym}(\delta F) \right) \\ & = (2\mu \epsilon + \lambda \tr (\epsilon) I) : \text{Sym}(\delta F) \\ & = (2\mu \epsilon + \lambda \tr (\epsilon) I) : \delta F \\ \Rightarrow P & = 2\mu \epsilon + \lambda \tr (\epsilon) I \\ & = \mu (F + F^T - 2I) + \lambda \tr (F - I) I \end{align}\]

我们得出了线性弹性模型的PK1数： \(P = \mu (F + F^T - 2I) + \lambda \tr (F - I) I\)

这个模型是线性的，容易计算，但是其与旋转相关。

圣维南-基尔霍夫模型

使用格林应变张量($E = \frac 1 2 (F^T F - I)$)定义的非线性模型： \(\Psi (F) = \mu \tr(E^2) + \frac \lambda 2 \tr^2(E)\) 称为圣维南-基尔霍夫模型（St. Venant-Kirchhoff model, StVK）。 注意$E$为对称矩阵。

线性弹性模型实际上就是无穷小形变假设下的StVK模型。

同样，对其微分来求$P$： \(\begin{align} \delta \Psi (F) & = 2\mu \cdot \tr (E\delta E) + \lambda \tr(E) \tr(\delta E) \\ \delta E & = \frac 1 2 (\delta F^TF + F^T \delta F) = \text{Sym}(F^T \delta F) \\ \tr (E\delta E) & = I: E\cdot \text{Sym} (F^T \delta F) = E : \text{Sym} (F^T \delta F)^T = E: F^T \delta F = FE : \delta F \\ \tr (\delta E) & = I: \text{Sym}(F^T \delta F) = I:F^T \delta F = F : \delta F \\ \Rightarrow \delta \Psi(F) & = 2 \mu FE : \delta F + \lambda \tr(E) F:\delta F \\ & = (2\mu FE + \lambda \tr(E) F):\delta F \\ \Rightarrow P &= 2\mu FE + \lambda \tr (E) F = F(2 \mu E + \lambda \tr(E) I) \end{align}\)

1. 上面第三行用了性质$A : B^T C = BA : C$，可以通过展开点乘和矩阵乘法证明；
2. 点乘易错点$A(B:C) \neq (AB) : C$。

于是我们得到了StVK模型的PK1数：$P = F(2 \mu E + \lambda \tr(E) I)$。

这其实是一个关于$F$的三次式，十分复杂，运算复杂度高，优势在于他是一个旋转无关的模型，更符合现实。 其仍有仿真缺陷，当弹性体被压缩到一定程度后，其抵抗压缩的力会突然减小（？），导致弹性体更容易被压缩到一个质点。

共旋转线性模型

线性模型没有旋转无关性，StVK效率很低，于是共旋转线性模型试图折衷考虑两者。 比较复杂，不细讲

能量密度满足： \(\Psi (F) = \mu \| S - I \| ^2_F +\frac \lambda 2 \tr^2(S-I)\) 的模型称为共旋转线性模型(Corotated Linear Elasticity)， 其中$S$为包含在$F$中的拉伸矩阵$F = RS$。

\(P = 2 \mu (F - R) + \lambda \tr (R^T F - I) R\) 相比于线性模型，直观的理解是这个模型手动减去$R$直接去除旋转的影响。

旋转无关、各向同性与新胡克模型

我们接下来关心这些模型的几个性质：

旋转无关（rotational invariance），指弹性体整体旋转不会改变每个位置的能量密度，可以写为$\Psi(RF)=\Psi(F)$，即先变换后，再进行旋转不会影响能量。

各向同性（Isotropy），指弹性体从任意方向进行相同的变换（拉伸旋转平移等）不会影响每个位置的能量密度，可以写为$\Psi(FQ) = \Psi(F)$，即先旋转再进行相同的变换不影响能量。

若模型既是旋转无关也是各向同性的，则能量密度函数满足：$\Psi (RFQ) = \Psi (F)$。

根据这个形式，自然想到将$F$进行奇异值分解(SVD)： $F = U \Sigma V^T$，其中$U$和$V$一定为酉矩阵，即$U^T U =V^T V=I$，与旋转矩阵一致。 而$\Sigma$为对角矩阵，在三维情况下自由度仅为3。

故一个满足旋转无关且各向同性的模型，一定满足$\Psi (F) = \Psi (\Sigma)$，并且一定存在某种定义方式，只需要自由度为3的变量，就能定义能量密度函数$\Psi$。 使用$F$的特征值当然是一种办法，但其非常不方便，所以这里定义了另外3个线性无关的变量，被称作各向同性不变量(isotropic invariants)： \(\begin{align} I_1(F) &= \tr(\Sigma ^2) = \tr(F^T F) \\ I_2(F) &= \tr(\Sigma ^4) = \tr\left[(F^T F)^2\right] \\ I_3(F) &= \text{det}(\Sigma ^2) = (\text{det}F)^2 \\ \end{align}\)

没有下标的$I$均表示单位矩阵，有下标的表示各向同性不变量。

对这3个变量求微分可以得到： \(\begin{align} \delta I_1 & = 2(I:F:\delta F) = 2F : \delta F \\ \delta I_2 & = 4\tr \left( F^TFF^T\delta F \right) = 4 I: F^TFF^T\delta F = FF^TF:\delta F \\ \delta I_3 & = 2\text{det}(F)\delta [\text{det} F] = 2 (\text{det}F)^2(F^T)^{-1} : \delta F \end{align}\) 其中$I_3$用了行列式微分的雅可比公式： \(\delta \det A = \det A \tr(A^{-1} \delta A) = \det A \cdot A^{-T} : \delta A\)

现在可以定义新胡克模型了：

能量密度满足： \(\Psi(I_1, I_3) = \frac \mu 2 \left[I_1-\log(I_3) - 3\right] + \frac \lambda 8 \log^2(I_3)\) 的模型称为新胡克（Neohookean）模型。

此时PK1数很好求： \(\begin{align} P & = \frac {\partial \Psi} {\partial F} \\ & = \frac \mu 2 \left( \frac {\partial I_1} {\partial F} - \frac 1 {I_3} \frac {\partial I_3} {\partial F} \right) + \frac {\lambda \log(I_3)} {4 I_3} \frac {\partial I_3} {\partial F} \\ & = \mu \left( F- F^{-T} \right) + \frac \lambda 2 \log(I_3) F^{-T} \end{align}\) Neohookean模型保证了旋转无关，各向同性，其能处理强压缩，虽然容易数值爆炸，因为$I_3 = J^2$，当压缩力很大时，体积变化分数$J \rightarrow 0$，会导致能量中的$\log$趋于无穷大，能量也非常大。

各种弹性能量模型对比

借用课程上的一张表

名称	$\Psi$	$P$	效率	旋转无关	各向同性	能处理强压缩
线性模型	$\mu \tr(\epsilon^2) + \frac \lambda 2 \tr^2(\epsilon)$	$\mu (F + F^T - 2I) + \lambda \tr (F - I) I$	高	否	否	否
StVK	$\mu \tr(E^2) + \frac \lambda 2 \tr^2(E)$	$F(2 \mu E + \lambda \tr(E) I)$	低	是	是	否
Corotated	$\mu |S - I |^2_F +\frac \lambda 2 \tr^2(S-I)$	$2 \mu (F - R) + \lambda \tr (R^T F - I) R$	中	是	是	否
Neohookean	$\frac \mu 2 \left[I_1-\log(I_3) - 3\right] + \frac \lambda 8 \log^2(I_3)$	$\mu \left( F- F^{-T} \right) + \frac \lambda 2 \log(I_3) F^{-T}$	低	是	是	是

有限元方法离散化

在前面有关弹性体的所有推导中，都是将弹性体视为一个拥有无限点的三维区域，这使得计算机无法进行仿真。 有限元方法的目的就是用有限的元素去拟合一个无穷的区域，在这里，我们将三维弹性体视为若干个四面体元素：

在四面体内部，变换函数将变成线性函数：$\hat{\phi}(X) = AX+b$，因此变换梯度也变得简单：$F = A$。 并且$F$可以使用四面体上的四个点解方程解出来： \(\begin{align} \left[\begin{matrix} \mathbf{x}_1-\mathbf{x}_4 & \mathbf{x}_2-\mathbf{x}_4 & \mathbf{x}_3-\mathbf{x}_4 \end{matrix}\right] & = F\cdot \left[\begin{matrix} X_1-X_4 & X_2-X_4 & X_3-X_4 \end{matrix}\right] \\ D_s & = F \cdot D_m \\ \Rightarrow F & = D_s D_m^{-1} \end{align}\) 令$W$为四面体的体积，可以计算得到$W = \frac 1 6 \left | \text{det} D_m \right |$

则四面体内部$T_i$的弹性势能为： \(E_i = \int_{T_i} \Psi (F_i) \mathrm{d}X = W_i \Psi (F_i)\) 每一个四面体都会对它的4个顶点施加弹力，对每一个顶点需要计算它相关的所有四面体对它贡献的弹力，第$i$个四面体对其第$j$个顶点的贡献为 \(f_i^j = - \frac {\partial E_i} {\partial \mathbf{x}_j}\)

令 \(H_i = \begin{bmatrix} f_i^1 & f_i^2 & f_i^3 \end{bmatrix}\) 则： \(\begin{align} f_i^j & = - \frac {\partial E_i} {\partial \mathbf{x}_j} \\ & = -\frac {\partial E_i} {\partial \Psi} \frac {\partial \Psi} {\partial F_i} \frac {\partial F_i} {\partial \mathbf{x}_j} \\ & = \left[-W_i P(F_i)D_m^{-T}\right]_{ji} \\ \Rightarrow H_i & = -W_i P(F_i)D_m^{-T} \end{align}\) 对于四面体第四个顶点：$f_i^4 = -f_i^1 -f_i^2 -f_i^3$，因为弹性内力要保证整个四面体合力为零。 由此可以算出弹性体分成四面体元素后，内部每一个顶点的受力。

时间积分

前文已经详细讲述了如何计算通过将弹性体视为四面体后计算每个顶点的受力，知道受力后就可以很方便的使用显示时间积分进行仿真，但隐式时间积分仍然比较困难。 本节介绍弹性体仿真时间积分的处理方式，以及一些算法实践上的具体处理。

显示时间积分

以下符号，向量都代表整个弹性体所有点，即$m$个点，如$\mathbf{x}$为$3m$维

对每一个节点进行如下操作： \(\begin{align} \mathbf{x}_{n+1} & = \mathbf{x}_n + \mathbf{v}_{n+1} \Delta t \\ \mathbf{v}_{n+1} & = \mathbf{v}_n + M^{-1} \left[f_{int} (\mathbf{x}_n)+f_{damp}(\mathbf{v}_n) + f_{ext}\right] \Delta t \end{align}\) 其中$f_{int}$表示内力，就是弹性体内部的弹力，需要用前文所述的各种能量计算出来； $f_{damp}$表示阻尼力，与速度有关，用于能量损耗； $f_{ext}$为外力，一般与弹性体形态无关，很容易计算，如重力。

隐式时间积分

需要计算如下方程： \(\begin{align} \mathbf{x}_{n+1} & = \mathbf{x}_n + \mathbf{v}_{n+1} \Delta t \\ \mathbf{v}_{n+1} & = \mathbf{v}_n + M^{-1} \left[f_{int} (\mathbf{x}_{n+1}) + f_{damp}(\mathbf{v}_{n+1}) + f_{ext}\right] \Delta t \end{align}\) 对于阻尼力 $f_{damp}$，如果难以处理，也可以使用$f_{damp} (\mathbf{v}_n )$做半隐式积分。

对于内力$f_{int} (\mathbf{x}_{n+1})$，一般采用泰勒展开近似：(省略了下标) \(f (\mathbf{x}_{n+1}) ≈ f (\mathbf{x}_n) + \frac {\partial f} {\partial \mathbf{x}} \bigg |_{\mathbf{x}_n}(\mathbf{x}_{n+1} - \mathbf{x}_n)\) 为此我们需要求解刚度矩阵：$K(\mathbf{x}_n) = - \frac {\partial f} {\partial \mathbf{x}} \bigg |_{\mathbf{x}_n}$，这个刚度矩阵是一个$3m \times 3m$的矩阵，但我们计算时一般只考虑一个四面体，一个四面体至少需要考虑到3个顶点，那么算法运行时需要用到一个$3 \times 3 \times 3 \times 3$的矩阵，这依然是一个比较大的运算量。 实现时，可以直接将这个巨型矩阵直接算出来（比如VegaFEM）；也可以用隐式计算的方法，直接计算 $K \cdot w$，只需要在求微分$\delta K$之后，把里面的所有$\delta \mathbf{x}$换成$w$即可。

现在来考察一下$K$，通过对力$f$求微分来推导： \(\begin{align} K(\mathbf{x}_n) \delta \mathbf{x} & = \frac {\partial f} {\partial \mathbf{x}} \delta \mathbf{x} = \delta f \\ \delta H & = \begin{bmatrix} \delta f_1 & \delta f_2 & \delta f_3 \end{bmatrix} \\ & = - W \delta P (F) D_m^{-T} \\ \delta P(F) & = \begin{cases} \delta F \left[ 2 \mu E + \lambda \tr(E) I \right] + F \left[ 2 \mu \delta E + \lambda \tr(\delta E) I \right] &\text{if using StVK model} \\ \mu \delta F + \left[ \mu - \lambda \log (J) \right] F^{-T} \delta F^T F^{-T} + \lambda \tr(F^{-1} \delta F) F^{-T} & \text{if using Neohookean model} \\ \end{cases} \\ \delta F & = (\delta D_s) D_m^{-1} \\ \delta D_s & = \begin{bmatrix} \delta \mathbf{x}_1 - \delta \mathbf{x}_4 & \delta \mathbf{x}_2 - \delta \mathbf{x}_4 & \delta \mathbf{x}_3 - \delta \mathbf{x}_4 \end{bmatrix} \end{align}\)

$J = \sqrt {I_3} = \text{det} F$是体积变化分数。

根据以上公式可以对任意弹性模型计算$K(\mathbf{x}_n)\cdot w$，只需要最终在$\delta D_s$中，把$\delta \mathbf{x}$替换为对应位置的$w$即可。 在Siggraph 2012的FEM course中推荐的就是这个计算$K$的方法。

VegaFEM的实现方法

令 \(\Delta \mathbf{x} = \mathbf{x}_{n+1} - \mathbf{x}_n , \Delta \mathbf{v} = \mathbf{v}_{n+1} - \mathbf{v}_n\) 对于阻尼力，令 \(\begin{multline} f_{damp} (\mathbf{v}_{n+1}) \\ = (\alpha_1 \cdot (-K(\mathbf{x}_n)) + \alpha_2 \cdot M) \cdot \mathbf{v}_{n+1} \\ = -D \cdot \mathbf{v}_{n+1} \end{multline}\) 则 \(\begin{align} \Delta \mathbf{v} & = \Delta t M^{-1}\left[ f_{int}(\mathbf{x}_n) - K(\mathbf{x}_n) \cdot \Delta \mathbf{x} - D \cdot \mathbf{v}_{n+1} \right] \\ M \cdot \Delta \mathbf{v} & = \Delta t f_{int} (\mathbf{x}_n) - \Delta t^2 K(\mathbf{x}_n) \cdot \Delta \mathbf{v} - \Delta t D (\mathbf{v}_n + \Delta \mathbf{v}) \\ \left[\Delta t^2 K(\mathbf{x}_n) + \Delta t D + M \right] \cdot \Delta \mathbf{v} & = \Delta t \left[f_{int}(\mathbf{x}_n) - (\Delta t K(\mathbf{x}_n) + D) \mathbf{v}_n\right] \end{align}\) 方程已经变成了$AX=b$的形式，可以使用各种方式迭代解出$\Delta \mathbf{v}$，比如下面Projective Dynamics将会详细讲的一种雅可比迭代法(Jacobi)

Siggraph 2012的FEM course中的方法与此类似，只不过公式的推导过程变成了求解$\Delta \mathbf{x}$而不是$\Delta \mathbf{v}$，与下面的Projective Dynamics也很类似。

Projective Dynamics (PD) 方法

PD算法的主要思路是，对每一个非线性约束，用投影的方式找到满足约束的位置，然后将约束视为到投影位置的某种距离（是个二次式），从而将整个非线性问题转化为一个二次优化问题。

比如弹性体中，四面体网格的弹性内力可以视为约束力，通过投影到最近的零势能位置（rest position），就可以转化为二次优化问题。

但是，仔细考虑就可以发现，这种投影方式本质上就是之前提到过的：将非线性的弹力计算函数用一阶泰勒展开近似。

下面来推导一下： \(\begin{align} \mathbf{x}_{n+1} & = \mathbf{x}_n + \mathbf{v}_{n+1} \Delta t \\ \mathbf{v}_{n+1} & = \mathbf{v}_n + M^{-1} \left[f_{int} (\mathbf{x}_{n+1}) + f_{ext}\right] \Delta t \\ \Rightarrow \mathbf{x}_{n+1} & = \mathbf{x}_n + \Delta t \mathbf{v}_n + \Delta t^2 M^{-1} \left[f_{int}(\mathbf{x}_{n+1}) + f_{ext}\right] \\ \frac 1 {\Delta t^2} &M \left(\mathbf{x}_{n+1} - \mathbf{x}_n - \Delta t \mathbf{v}_n - \Delta t^2 M^{-1} f_{ext} \right) = f_{int}(\mathbf{x}_{n+1})\\ 令 \hat{\mathbf{x}}_n & = \mathbf{x}_n + \Delta t \mathbf{v}_n + \Delta t^2 M^{-1} f_{ext} \\ \frac 1 {\Delta t^2} M \left(\mathbf{x}_{n+1} - \hat{\mathbf{x}}_n \right)& = f_{int}(\mathbf{x}_{n+1}) \approx K(\mathbf{x}^{(k)}) \cdot \left(\mathbf{x}_{n+1} - \mathbf{x}^{(k)} \right) + f_{int} (\hat{\mathbf{x}}_n)\\ \end{align}\)

假设我们有一个迭代求解序列 $\mathbf{x}^{(0)}, …, \mathbf{x}^{(k)}$，使得其越来越逼近于 $\mathbf{x}_{n+1}$，并且初始 $\mathbf{x}^{(0)}=\hat{\mathbf{x}}_n$

我们可以将方程写为：

\[\left(\frac 1 {\Delta t^2} M - K(\mathbf{x}^{(k)})\right) \cdot \mathbf{x}^{(k+1)} = \frac 1 {\Delta t^2} M \hat{\mathbf{x}}_n + f_{int}(\mathbf{x}^{(k)}) - K(\mathbf{x}^{(k)})\mathbf{x}^{(k)}\]

此时求解 $\mathbf{x}^{(k+1)}$已经是$AX=b$的形式了，接下来讲一讲迭代求解。

令 $A=B-C$， 则 \((B-C)X=b \\ BX - CX = b \\ BX = CX + b \\ X = B^{-1} (CX + b)\) 使用一个迭代序列 $X^{(0)}, …, X^{(k)}$，令： \(X^{(k+1)} = B^{-1} (CX^{(k)} + b)\)

一般我们取$B$为$A$的对角元素，因为往往物理方程上面$A$的对角元素都比其他元素大，占主要部分，这样这个序列会更快收敛到方程的解。

我们计算一下$X^{(k+1)}-X^{(k)}$： \(X^{(k+1)} = B^{-1}(B-A)X^{(k)} + B^{-1} b\\ X^{(k+1)}-X^{(k)}=B^{-1}(b-AX)\)

现在我们代入到之前的隐式欧拉积分方程： \(\begin{align} A&=\frac 1 {\Delta t^2} M - K(\mathbf{x}^{(k)})\\ B&=\text{diag}\{A\} \\ X^{(k)}&=\mathbf{x}^{(k)}\\ b&=\frac 1 {\Delta t^2} M \hat{\mathbf{x}}_n + f_{int}(\mathbf{x}^{(k)}) - K(\mathbf{x}^{(k)})\mathbf{x}^{(k)} \end{align}\)

\[\begin{align} \mathbf{x}^{(k+1)} - \mathbf{x}^{(k)} &= B^{-1}\left( \frac 1 {\Delta t^2} M \hat{\mathbf{x}_n} + f_{int} (\mathbf{x}^{(k)})-K(\mathbf{x}^{(k)})\mathbf{x}^{(k)} - \frac 1 {\Delta t^2} M \mathbf{x}^{(k)} + K(\mathbf{x}^{(k)})\mathbf{x}^{(k)}\right) \\ \mathbf{x}^{(k+1)} - \mathbf{x}^{(k)} &= B^{-1}\left( \frac 1 {\Delta t^2} M \left(\hat{\mathbf{x}_n} - \mathbf{x}^{(k)}\right) + f_{int} (\mathbf{x}^{(k)}) \right) \end{align}\]

最终等式右边以全部是已知量，直接迭代求解即可。

还可以参考这篇文章使用切比雪夫多项式加速：A Chebyshev Semi-Iterative Approach for Accelerating Projective and Position-based Dynamic

计算出$\mathbf{x}^{(k+1)} - \mathbf{x}^{(k)}$之后，还可以进一步使用线搜索(line search)算法来计算迭代步长，进一步优化收敛速度

在A Chebyshev Semi-Iterative Approach for Accelerating Projective and Position-based Dynamic的实现中，使用的能量模型并不是前文提到的3中模型之一，而是只用了共旋转线性模型(Corotated Linear Elasticity)的一半，即$\Psi (F) = \mu | S - I | ^2_F = \mu | F - R | ^2_F$. 此时$P(F) = 2 \mu (F-R)$ 并且该文章的实现在计算$B^{-1}$中的$K(\mathbf{x}^{(k)})$的对角元时，忽略了$P(F)$中矩阵分解出来的$R$，故$\frac {\partial P} {\partial F} \approx 2 \mu \cdot I$，可以推出$\delta H \approx - 2\mu W \cdot \delta D_s \cdot D_m^{-1} D_m^{-T}$. 这使得计算出来的$K(\mathbf{x}^{(k)})$与$\mathbf{x}^{(k)}$无关，故可以将其在读取模型的时候进行预处理，从而大幅加速求解速度。 文章的实验显示这种近似是可行的。

Incremental Potential Contact (IPC) 方法

\[\begin{align} \mathbf{x}_{n+1} & = \mathbf{x}_n + \mathbf{v}_{n+1} \Delta t \\ \mathbf{v}_{n+1} & = \mathbf{v}_n + M^{-1} \left[f_{int} (\mathbf{x}_{n+1}, \mathbf{v}_{n+1}) + f_{ext}\right] \Delta t \\ \Rightarrow \mathbf{x}_{n+1} & = \mathbf{x}_n + \Delta t \mathbf{v}_n + \Delta t^2 M^{-1} \left[f_{int}(\mathbf{x}_{n+1}, \mathbf{v}_{n+1}) + f_{ext}\right] \\ \mathbf{x}_{n+1} & - \left( \mathbf{x}_n + \Delta t \mathbf{v}_n + \Delta t^2 M^{-1} f_{ext}\right) - \Delta t^2 M^{-1} f_{int}(\mathbf{x}_{n+1}, \mathbf{v}_{n+1}) = 0 \\ 令 \hat{\mathbf{x}}_n & = \mathbf{x}_n + \Delta t \mathbf{v}_n + \Delta t^2 M^{-1} f_{ext} \\ M(\mathbf{x}_{n+1} & - \hat{\mathbf{x}}_n) - \Delta t^2 f_{int} (\mathbf{x}_{n+1} , \mathbf{v}_{n+1}) = 0 \end{align}\]

我们将其转为一个优化问题，即： \(\min_{\mathbf x} E(\mathbf{x}) = \frac 1 2 \left \| \mathbf{x} - \hat{\mathbf{x}}_n \right \|_M + \Delta t^2 \Psi (\mathbf{x})\) 因为 \(E'(\mathbf{x}) = M(\mathbf{x}_{n+1} - \hat{\mathbf{x}}_n) - \Delta t^2 f_{int} (\mathbf{x}_{n+1} , \mathbf{v}_{n+1})\) 则只需求解 \(\mathbf{x}_{n+1} = \arg \min_{\mathbf{x}} E(\mathbf{x})\) $E(\mathbf{x})$被称为 Increment Potential，可以看出$\frac 1 2 \left | \mathbf{x} - \hat{\mathbf{x}}_n \right |_M$与惯性势能相关，$ \Delta t^2 \Psi(\mathbf{x})$与弹性势能有关。 优化问题的目的就是找到一个势能最小的位置作为下一时间步的状态。 注意到这个$E(\mathbf{x})$还可以添加其他能量包括碰撞，摩擦等。

IPC使用牛顿-拉夫逊方法结合线搜索(line search)迭代优化$E(\mathbf{x})$，即 \(H_E(\mathbf{x}) \cdot v = - D_E(\mathbf{x})\) 其中$H_E(\mathbf{x})$为$E$的二阶导（海塞矩阵），$D_E$为一阶导，这里面涉及到对$\Psi$求二阶导。

求解出$v$后作为搜索的方向，使用线搜索结合连续碰撞检测(CCD)，找到最佳步长进行移动，重复这个步骤直到时间步长加起来达到一个仿真步。

参考资料

FEM simulation of 3D deformable solids Siggraph 2012 Course

A Chebyshev Semi-Iterative Approach for Accelerating Projective and Position-based Dynamic

Projective Dynamics: Fusing Constraint Projections for Fast Simulation

VegaFEM

IPC

北京大学可视计算与交互概论