弹性体模拟

9 弹性体模拟

9.1 连续介质的守恒定律

• 连续方程: \[ \pdv{\rho}{t} = -\nabla\cdot(\rho v). \]

• 线动量守恒 (Cauchy 第一运动方程): \[ \rho \dv{v}{t} = \nabla\cdot\sigma + f_{\sf body}. \]

  • \(\sigma\) 为 Cauchy 应力张量. 假如把物体沿着一个小曲面 (单位法向量为 \(n\)) 分为两部分, 则这两部分就会通过这张曲面相互施加作用力 \(t\) (应力, 即单位面积上的内力). 应力张量描述了 \(t\) 与 \(n\) 的线性关系: \(t=\sigma n\).
  • \(f_{\sf body}\) 为外力, 比如重力 \(\rho g\).

• 角动量守恒 (Cauchy 第二运动方程): \[ \sigma = \sigma\T. \]

  • \(\sigma\) 是一个对称的 \((0,2)\) 型张量.

Newton 流体的 Cauchy 应力张量为 \[ \sigma = -pI + \mu(\nabla v + (\nabla v)\T). \]

9.2 模拟框架

• Shape Changes:
  • Deformation Measure
  • Strain
• Elastic Energy
• Forces: \(\rho\dd{v}/\dd{t}=\nabla\cdot\sigma+f_{\sf body}\).
• Dynamics:
  • Explicit Time Integration.
  • Implicit Time Integration.

10 Linear FEM

参考资料: FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction (SIGGRAPH 2012 Course).

后文会使用一些张量指标记号, 上标为反变指标, 下标为协变指标, 并遵从 Einstein 求和约定. 在欧氏空间 \(\R^3\) 中, 可以使用欧氏度量 \(g\) 升降指标 (music isomorphism); \((0,2)\) 型或 \((2,0)\) 型张量的迹 \(\tr\) 定义为其与 \(g\) (或 \(g^{-1}\)) 缩并.

Linear FEM (线性有限元) 将形变近似为分片线性映射, 其计算流程如下: \[ \textsf{形变} \xrightarrow{\textsf{超弹性模型}} \textsf{应力张量} \longrightarrow \textsf{内力} \longrightarrow \textsf{时间积分}. \]

10.1 形变场

设 \(\R^3\) 中的物体 \(\Omega\) 经过了连续可微的形变 \(\varphi:\Omega\to\Omega'\), \(X\mapsto\varphi(X)=x\). 变形前的 \(\Omega\) 称为参考状态 (reference state), 变形后的 \(\Omega'\) 称为当前状态 (current state).

形变梯度 (deformation gradient) 定义为切映射 \(T_X\varphi:\R^3\to\R^3\), 它是一个 \((1,1)\) 型张量 \[ F^i_J := \pdv{x^i}{X^J}, \] 称为形变张量 (deformation tensor). 注意参考状态和当前状态的指标用大小写区分开.

形变梯度给出了形变的局部线性近似. 在计算机模拟中, 我们通常不能模拟一个连续的场, 只能离散地近似它们. 在线性有限元方法 (FEM) 中, 我们用一个分片线性映射 \(\tilde\varphi\) 去近似形变场 \(\varphi\). 将区域 \(\Omega\subset\R^3\) 剖分为若干四面体, 设 \(\tilde\varphi\) 在每个四面体 \({\cal T}\) 上都是线性映射, 即 \({\cal T}\) 按照 \(\tilde\varphi\) 的形变张量 \(F\) (在 \({\cal T}\) 上常值) 进行线性变换.

如下图, 设 \({\cal T}\) 的四个顶点为 \(X_0,X_1,X_2,X_3\), \({\cal T}'\) 的对应顶点为 \(x_0,x_1,x_2,x_3\), 则有 \[ F\underbrace{\pmqty{E_1 & E_2 & E_3}}_E = \underbrace{\pmqty{e_1 & e_2 & e_3}}_e,\quad\textsf{即}\; F = eE^{-1}, \] 其中 \(E_I=X_I-X_0,e_i=x_i-x_0\) 为边向量, \(E,e\) 称为局部标架矩阵.

10.2 应力张量

下面介绍几个应力张量.

• 回顾前文, Cauchy 应力张量 \(\sigma\) 给出了当前状态下的应力 \(t\) 与法向 \(n\) 的关系: \[ t^i = \sigma^{ij} n_j. \]

• 第二 Piola-Kirchhoff 应力张量 \(S\) 给出了参考状态 (\({\cal T}\)) 中应力 \(T\) 与法向 \(N\) 的关系: \[ T^I = S^{IJ} N_J. \]

• 第一 PK 应力张量 \(P\) 给出了当前状态 (\({\cal T}'\)) 下的应力 \(t\) 与参考状态下法向 \(N\) 的关系: \[ t^i = P^{iJ} N_J. \] 结合 \(t^i=F^i_JT^J\), 可知第一和第二 PK 应力张量满足 \(P^{iJ}=F^i_I S^{IJ}\).

这三个应力张量都是对称的 \((2,0)\) 型张量.

下面看看 Cauchy 张量和 PK 张量的关系. 取 \({\cal T}\) 中的面元 \(N\dd{S}\), 它经过 \(\tilde\varphi\) (\(F\)) 变成 \({\cal T}'\) 中的面元 \(n\dd{S'}\), 则由 Nanson 公式可得 \[ \dd{S'} = (\det{F}) \|NF^{-1}\| \dd{S}. \] 内力平衡要求 \(t^i\dd{S'}=F^i_JT^J\dd{S}\), 于是 \[ (\det{F}) \|NF^{-1}\| t^i = F^i_J T^J. \tag{a} \] 另一方面, 结合法向量的变换关系, 我们有 \[ t^i = \sigma^{ij} n_k = \sigma^{ij} \frac{N_I (F^{-1})^I_j}{\|NF^{-1}\|}. \tag{b} \] 将 \(\textrm{(b)}\) 代入 \(\textrm{(a)}\), 整理可得 \[ T_I = \underbrace{ (\det{F}) \cdot \sigma^{ij} (F^{-1})^I_i (F^{-1})^J_j }_{S^{IJ}} \cdot N_J. \]

10.3 超弹性模型

超弹性模型 (hyperelastic model) 认为物体的形变总是弹性形变, 不论形变有多大. 超弹性模型假设应力张量 (\(S,P\) 等) 可以由应变能量密度 \(\Psi\) 的导数给出.

应变能量密度是关于形变 \(F\) 的非负实值函数 \(\Psi=\Psi(F)\), 而总能量就是能量密度的积分: \[ \mathcal{E}(F) = \int_{\Omega} \Psi(F)\dd{S}. \] 第一 PK 应力张量是 \(\Psi\) 关于 \(F\) 的导数: \[ P^{iJ} = g^{ij} \pdv{\Psi}{F^j_J}. \] 在线性 FEM 中, \(F\) 在四面体上常值, 故 \(\Psi(F)\) 在四面体上常值. 在四面体 \({\cal T}\) 上有 \(\mathcal{E}(F)={\rm Vol}({\cal T})\Psi(F)\). 能量的负梯度为内力 \(f\): \[ \Align{ f^i_j = -\pdv{\mathcal{E}}{e_i^j} &= -{\rm Vol}({\cal T}) \pdv{\Psi}{e_i^j} \\ &= -{\rm Vol}({\cal T}) \pdv{\Psi}{F^k_K} \pdv{F^k_K}{e_i^j} \\ &= -{\rm Vol}({\cal T}) (g_{km}P^{mK}) \cdot \pdv{e_i^j} \pqty{ e_l^k (E^{-1})^l_K } \\ &= -{\rm Vol}({\cal T}) g_{ij} (E^{-1})^i_I P^{jI}, } \] 其中 \(f^i_j\) 表示第 \(i\) 个力的第 \(j\) 个分量, \(e_i^j\) 表示第 \(i\) 个基向量的第 \(j\) 个分量.

10.3.1 St. Venant-Kirchhoff 模型

Saint Venant-Kirchhoff 模型 (StVK) 的能量密度 \(\Psi\) 与 Green 应变张量有关.

Green 应变张量 (Green strain tensor) 是一个对称的 \((0,2)\) 型张量, 定义为 \[ G(X,Y) = \frac12(g(FX,FY) - g(X,Y)). \]

容易看出如果 \(F\) 复合一个正交变换 \(Z\) (变成 \(ZF\)) 的话, \(G\) 不会改变, 所以 Green 应变张量是旋转不变的 (与旋转无关). 应变张量 \(G\) 的分量为 \((G_{IJ})=\frac12(g_{kl} F^k_I F^l_J - g_{IJ})\).

Green 应变描述了某个向量在 \(\tilde\varphi\) 下的长度变化, 可以反映形变的弹性势能. 设切向量 \(X\in\R^3\), 它在切映射 \(F\) 下的像为 \(x=FX\), 则 \[ \Align{ G(X,X) &= \frac12 ( g(FX,FX) - g(X,X) ) \\ &= \frac12 (\|x\|^2 - \|X\|^2). } \]

StVK 的能量密度表达式为 \[ \Psi(F) = \frac\lambda2\tr^2(G) + \mu \|G\|_F^2. \] 其中 \(\lambda,\mu\) 为 Lamé 参数. Green 应变的旋转不变性导致 StVK 的旋转不变性, 即 \(\Psi(ZF)=\Psi(F)\) 对任意正交变换 \(Z\). 在 StVK 下, 第一 PK 张量为 \[ \Align{ P^{iJ} = g^{ij} \pdv{\Psi}{F^j_J} &= \pdv{\Psi}{G_{KL}} \pdv{G_{KL}}{F^j_J} \\ &= F^i_I \cdot \underbrace{(\lambda \tr(G) g^{IJ} + 2\mu G^{IJ})}_{S^{IJ}}, } \] 比较上式和 \(P=FS\), 可以发现括号里的项 (即 \(\partial\Psi/\partial G\)) 就是 \(S\).

Lamé 参数与物理量的联系: 设 Poisson 比为 \(\nu\) (a measure of incompressibility), 杨氏模量为 \(E\) (a measure of stretch resistance), 则 \[ \Align{ \mu &= \frac{E}{2(1+\nu)}, & \lambda &= \frac{E\nu}{(1+\nu)(1-2\nu)}. } \]

10.3.2 Neohookean 模型

StVK 在大形变以及翻转的情形下效果不够好. 考虑一维的情形 (形变 \(\tilde\varphi:\R\to\R\), 形变梯度 \(F\in\R\)), 则 StVK 的第一 PK 应力为 \[ P(F) = \frac12F[\lambda(F^2-1) + 2\mu(F^2-1)] = \frac{2\mu+\lambda}2 F(F+1)(F-1), \] 这是一个三次函数. 当形变梯度很大时, 应力也很大, 带来数值不稳定性. 当形变梯度接近零时 (极度压缩), 应力接近零甚至为正, 使得物体应变朝着另一个稳定点 \(-1\) 演化, 造成物体翻转.

沿袭 StVK 的思路, 我们希望能量密度 \(\Psi\) 是旋转不变的. 考虑正交不变量:

• (右 Cauchy-Green 不变量) \(C(X,Y):=g(FX,FY)\),

它衍生出三种不变量:

• \({\rm I}_C:=\tr(C)=\lambda_1^2+\lambda_2^2+\lambda_3^2\). (反映长度的形变)
• \({\rm II}_C:=\dfrac12\bqty{(\tr{C})^2-\tr(C^2)}=(\lambda_1\lambda_2)^2+(\lambda_2\lambda_3)^2+(\lambda_3\lambda_1)^2\). (反映面积的形变)
• \({\rm III}_C:=\det(C)=(\lambda_1\lambda_2\lambda_3)^2\). (反映体积的形变)

其中 \(\lambda_1,\lambda_2,\lambda_3\) 为 \(F\) 的奇异值 (\(\lambda_1^2,\lambda_2^2,\lambda_3^2\) 为 \(C\) 的特征值). Neohookean 的能量密度是这三种不变量的函数: (假定 \(\det{F}>0\))

\[ \Psi(F) = \frac\lambda8 [\ln({\rm III}_C)]^2 + \frac\mu2 ({\rm I}_C-3) - \frac\mu2 \ln({\rm III}_C), \] 第二 PK 应力为 \[ P^{iJ} = g^{ij} \pdv{\Psi}{F^j_J} = \mu\!\bqty{ F^i_I g^{IJ} - F^i_I (C^{-1})^{IJ} } + \frac\lambda2 (\ln\det{C}) F^i_J (C^{-1})^{IJ}. \]

把 \(C\) 展开, 上式可化简为 \[ P^{iJ} = \mu[ F^i_I g^{IJ} - (F^{-1})^J_j g^{ij} ] + \lambda (\ln\det{F}) (F^{-1})^J_j g^{ij}. \] 如果要模拟嵌入在三维空间中的二维物体 (比如布料), 则张量 \(F:\R^2\to\R^3\) 不可逆, 则应按照原始形式 (包含张量 \(C\)) 计算.

Neohookean 可以比较好的保持体积; 一个缺点是 \(\det{F}\leq0\) (定向翻转) 时未定义.

一维情形的可视化: 我们把 Neohookean 和 StVK 的第二 PK 应力 \(P\) 随着 \(F\in\R\) 的图像画出来:

它们在 \(F=1\) 处是相切的 (对高维也成立), 故 StVK 和 NH 在小形变下几乎等价.

一个弹性模型称为各向同性的 (isotropic), 如果它关于右乘正交变换不变: \(\Psi(FZ)=\Psi(F)\). 如果 \(\Psi\) 同时为旋转不变的和各向同性的, 将 \(F\) 奇异值分解可以得到 \[ \Psi(F) = \Psi(U\Sigma V\T) = \Psi(\Sigma), \] 即 \(\Psi\) 只与 \(F\) 的奇异值有关, \(\Psi(F)=\Psi(\lambda_1,\lambda_2,\lambda_3)\). 右 Cauchy-Green 不变量不是各向同性的, 但 \({\rm I}_C,{\rm II}_C,{\rm III}_C\) 是各向同性的. 由此衍生出的 Neohookean 等超弹性模型也是各向同性的.

10.3.3 More...

一些常用的超弹性模型的能量密度, in terms of invariants:

• StVK: \[ \Psi_{\sf StVK}(F) := \frac\lambda2({\rm I}_C-3)^2 + \frac\mu4 ({\rm I}_C^2 - 2{\rm I}_C-2{\rm II}_C+3). \]

• Neohookean: \[ \Psi_{\sf NH}(F) := \frac\lambda2 \ln^2({\rm III}_C^{-1/2}) + \frac\mu2 ({\rm I}_C-3) - \mu \ln({\rm III}_C^{-1/2}). \]

• Mooney-Rivlin: \[ \Psi_{\sf MR}(F) := \sum_{r,s\geq0} \mu_{r,s}({\rm I}_C-3)^r({\rm II}_C-3)^s. \]

• Arruda-Boyce: \[ \Psi_{\sf AB}(F) := \mu\sum_{i=0}^n \frac{C_i}{N^{i-1}} ({\rm I}_C^i - 3^i). \]

• Fung: \[ \Psi_{\sf Fung}(F) := \frac12\bqty{a({\rm I}_C-3) + b (\exp({\rm I}_C-3)-1)}. \]

• Yeoh: \[ \Psi_{\sf Yeoh}(F) := \sum_{i=0}^n C_i({\rm I}_C - 3)^i. \]

小形变常用 StVK (快), 橡胶常用 Neohookean, 人体组织常用 MR 或 Fung...

10.4 Linear FEM 框架

把上面的推导总结一下 (将分量式转化为矩阵乘法), 可以得到如下的模拟框架.

• (初始化) 预计算每个四面体局部标架的逆 \(E^{-1}=\pmqty{X_1-X_0&X_2-X_0&X_3-X_0}^{-1}\) 和体积 \(\mathrm{Vol}(\mathcal{T})=\det(E)/6\).

• (时间步) 对每个四面体 \({\cal T}\):

  1. 计算当前局部标架 \(e=\pmqty{x_1-x_0&x_2-x_0&x_3-x_0}\).

  2. 计算形变梯度 \(F=eE^{-1}\).

  3. 计算第一 PK 应力 \(P\).

     • StVK 模型:
       1. Green 应变 \(G=(F\T F-I)/2\).
       2. 第二 PK 应力 \(S=2\mu G+\lambda\tr(G)I\).
       3. 第一 PK 应力 \(P=FS\).
     • Neohookean 模型:

       1. 右 Cauchy-Green 张量 \(C=F\T F\).

       2. 第一 PK 应力 \(P=\mu(F-FC^{-1})+\lambda(\ln\det{C})FC^{-1}/2\).

          或者, 当 \(F\) 为方阵时, \(P=\mu(F-F^\mathsf{-T})+\lambda(\ln\det{F})F^\mathsf{-T}\).

  4. 计算内力 \(\pmqty{f_1&f_2&f_3}=-{\rm Vol}({\cal T})PE^{\sf-T}\) 以及 \(f_0=-f_1-f_2-f_3\).

  5. 时间积分.