YSYSIMONCODE. ART. SCIENCE. VISION.

应力张量 Stress Tensor

整理 Dirichlet、StVK、ARAP 和 Neo-Hookean 能量对应的应力张量。

Aug 25, 2024约 8 分钟Simulationysysimon
SimulationDeformationStress TensorElasticity

上一章我们已经求出了三阶张量 ∂F∂x\dfrac{\partial\mathbf{F}}{\partial \mathbf{x}} 在非线性有限元分析中，这个三阶张量有时被描述为与几何刚度相关的导数或贡献，因为它捕捉了在几何非线性情况下的变形过程中，当空间位置变化时，变形梯度如何发生变化

我们接着要来看看 ∂Ψ∂F\dfrac{\partial\Psi}{\partial\mathbf{F}} 的情况，在上一章中，我们只是已经分析出它是一个矩阵 (二阶张量)，但是并没有真正计算它，因为对于不同的 弹性势能 定义，会得到不同的结果，所以它也称为 应力张量 Stress Tensor

不同的弹性势能 Ψ\Psi 由变形梯度 F\mathbf{F} 的不同形式作为能量函数 Ψ\Psi 的参数，也就是说，对不同的弹性势能函数求 ∂Ψ∂F\dfrac{\partial\Psi}{\partial\mathbf{F}}，将得到不同的 应力张量 Stress Tensor

狄利克雷能量 Dirichlet Energy - 第一类 Piola-Kirchhoff 应力张量

当我们对 狄利克雷能量 Dirichlet Energy

ΨDirichlet=∥F∥F2\Psi_{\text{Dirichlet}} = \|\mathbf{F} \|_F^2

求 ∂Ψ∂F\frac{\partial\Psi}{\partial\mathbf{F}}，我们将得到第一类 Piola-Kirchhoff 应力张量 可以简写为 PK1，有时它被写为 P(F)\mathbf{P}(\mathbf{F})

∂ΨDirichlet∂F=P(F)=∂∂F∥F∥F2=2F\frac{\partial\Psi_{\text{Dirichlet}}}{\partial\mathbf{F}} = \mathbf{P}(\mathbf{F}) = \frac{\partial}{\partial\mathbf{F}} \|\mathbf{F} \|_F^2 = 2\mathbf{F}

St. Venant Kirchhoff - 第二类 Piola-Kirchhoff 应力张量

我们可以对 StVK 能量的拉伸部分进行分析

ΨStvk,stretch=∥E∥F2\Psi_{\text{Stvk,stretch}} = \|\mathbf{E}\|_{F}^2

其中 E=12(FTF−I)\mathbf{E} = \frac{1}{2}(\mathbf{F}^T\mathbf{F} - \mathbf{I}) 是 Green-Lagrange 应变张量

根据我们学到的经验，我们应该对函数取其直接变量形式的导数，然后把复杂的变量代换放入链式法则中，Stvk 能量的直接变量是 E\mathbf{E}，也就是 Green-Lagrange 应变张量 ，这相当于定义了一个新的中间变量

∂Ψ∂x=∂Ψ∂E∂E∂F∂F∂x\frac{\partial\Psi}{\partial \mathbf{x}} = \frac{\partial\Psi}{\partial \mathbf{E}} \frac{\partial \mathbf{E}}{\partial \mathbf{F}} \frac{\partial \mathbf{F}}{\partial \mathbf{x}}

我们需要计算一个新的基变换张量 ∂E∂x\frac{\partial \mathbf{E}}{\partial \mathbf{x}}，一旦计算出了它，那么计算 第二类 Piola-Kirchhoff 应力张量 将非常简单

∂ΨStvk,stretch∂E=E\frac{\partial\Psi_{\text{Stvk,stretch}}}{\partial\mathbf{E}} = \mathbf{E}

第一类应力张量是 弹性势能 对 变形梯度 求导定义的 ∂Ψ∂F\frac{\partial\Psi}{\partial\mathbf{F}}，而第二类应力张量是 弹性势能 对 应变张量 求导定义的 ∂Ψ∂E\frac{\partial\Psi}{\partial\mathbf{E}}

StVK 是唯一依赖于 应变张量 E\mathbf{E} 的能量，其他所有能都直接由变形梯度 F\mathbf{F} 定义，我们引入这个中间张量可以简化计算

ΨStvk,stretch=14∥FTF−I∥F2=14(∥FTF∥F2+tr⁡I−2tr⁡(FTF))=14(∥FTF∥F2+tr⁡I−2∥F∥F2)∂ΨStvk,stretch∂F=PStvk,stretch(F)=F(FTF−I)=FE\begin{align*} \Psi_{\text{Stvk,stretch}} &= \frac{1}{4}\|\mathbf{F}^T\mathbf{F} - \mathbf{I}\|_F^2 \\ &= \frac{1}{4} \left(\|\mathbf{F}^T\mathbf{F}\|_F^2 + \operatorname{tr}\mathbf{I} - 2\operatorname{tr}(\mathbf{F}^T\mathbf{F}) \right) \\ &= \frac{1}{4} \left(\|\mathbf{F}^T\mathbf{F}\|_F^2 + \operatorname{tr}\mathbf{I} - 2\|\mathbf{F}\|_F^2 \right) \\ \frac{\partial\Psi_{\text{Stvk,stretch}}}{\partial\mathbf{F}} &= \mathbf{P}_{\text{Stvk,stretch}}(\mathbf{F}) \\ &= \mathbf{F}\left(\mathbf{F}^T\mathbf{F} - \mathbf{I}\right) \\ &= \mathbf{F}\mathbf{E} \end{align*}

完整的 StVK 能量

我们一直只研究了 StVK 的拉伸项，现在是时候来看看完整的能量定义了

ΨStVK=μ∥E∥F2+λ2(tr⁡E)2\Psi_{\text{StVK}} = \mu \|\mathbf{E}\|_F^2 + \frac{\lambda}{2}(\operatorname{tr}\mathbf{E})^2

其中 ∥E∥F2\|\mathbf{E}\|_F^2 项是原先的拉伸项，但是新的体积保持项 (tr⁡E)2(\operatorname{tr}\mathbf{E})^2 被添加进来，常量系数，μ\mu 和 λ\lambda 告诉模型需要多少相对拉伸阻力或体积保持

最终 PK1 是

∂ΨStvk∂F=PStvk(F)=μFE+λ(tr⁡E)F\frac{\partial\Psi_{\text{Stvk}}}{\partial\mathbf{F}} = \mathbf{P}_{\text{Stvk}}(\mathbf{F}) = \mu\mathbf{F}\mathbf{E} + \lambda(\operatorname{tr}\mathbf{E})\mathbf{F}

As-Rigid-As-Possible ARAP 能量

ARAP 能量 可以被展开为如下形式

ΨARAP=μ2∥F−R∥F2=μ2(∥F∥F2+∥R∥F2−2tr(FTR))=μ2(∥F∥F2+∥R∥F2−2tr S)\begin{align*} \Psi_{\text{ARAP}} &= \frac{\mu}{2}\|\mathbf{F}-\mathbf{R}\|_F^2 \\ &= \frac{\mu}{2} \left(\|\mathbf{F}\|_F^2 + \|\mathbf{R}\|_F^2 - 2\text{tr} (\mathbf{F}^T\mathbf{R}) \right) \\ &= \frac{\mu}{2} \left(\|\mathbf{F}\|_F^2 + \|\mathbf{R}\|_F^2 - 2\text{tr} \ \mathbf{S} \right) \end{align*}
  • 其中 S=FTR\mathbf{S} = \mathbf{F}^T\mathbf{R} 这是由于 极分解 polar decomposition F=RS\mathbf{F} = \mathbf{R}\mathbf{S} 得到的 S\mathbf{S} 是对称的

我们可以得到 PK1

∂ΨARAP∂F=P(F)ARAP=μ(F−R)\frac{\partial\Psi_{\text{ARAP}}}{\partial\mathbf{F}} = \mathbf{P}(\mathbf{F})_{\text{ARAP}} = \mu (\mathbf{F} - \mathbf{R})

Neo-Hookean 能量

我们知道 Neo-Hookean 能量 是

ΨNeo-Hookean=∥F∥F2−3\Psi_{\text{Neo-Hookean}} = \|\mathbf{F}\|_F^2 - 3

如果我们求它的 PK1

∂ΨNeo-Hookean∂F=PNeo-Hookean(F)=2F\frac{\partial\Psi_{\text{Neo-Hookean}}}{\partial\mathbf{F}} = \mathbf{P}_{\text{Neo-Hookean}}(\mathbf{F}) = 2\mathbf{F}

这看上去似乎和 狄利克雷能量 的应力张量一致，但是在 Mooney 在 1940 年的论文中，变形梯度 F\mathbf{F} 需要满足一个硬性条件: det⁡F=1\det\mathbf{F} = 1 但是这个条件对于 CG 里常见的形变来说过于苛刻

有许多其他形式的 Neo-Hookean 能量，其中比较流行的是 Bonet 和 Wood 在 2008 年提出的

ΨBW08=μ2(∥F∥F2−3)−μlog⁡(J)+λ2(log⁡(J))2\Psi_{\text{BW08}} = \frac{\mu}{2}(\|\mathbf{F}\|_F^2 - 3) - \mu\log(J) + \frac{\lambda}{2}(\log(J))^2
  • F\mathbf{F} 是形变梯度
  • J=det⁡(F)J = \det(\mathbf{F}) 是体积比，表示体积变形的程度
  • μ\mu 和 λ\lambda 是材料的拉梅常数 (材料常数)

在这个能量函数中

  1. μ2(∥F∥F2−3)\frac{\mu}{2}(\|\mathbf{F}\|_F^2 - 3) 这一项类似之前定义的简单 Neo-Hookean 能量 形式，描述了材料的剪切变形
  2. −μlog⁡(J)- \mu\log(J) 这一项涉及体积变形，通过 J=det⁡FJ = \det\mathbf{F} 的对数来描述体积变化，它确保了在无体积变化的情况下，体积变形项不对能量有贡献
  3. λ2(log⁡(J))2\frac{\lambda}{2}(\log(J))^2 这一项进一步控制体积变形，尤其是在大体积变化时，提供了额外的刚度

可以看到，体积保持 volume preserving 是这种能量形式的特色，而它来自于 J=det⁡FJ = \det\mathbf{F} 这告诉我们有多少体积被压缩或没有被压缩，事实上，任何声称“体积保持”但不包含 det⁡F\det\mathbf{F} 的能量都不能真正良好的保持体积

由于 det⁡F\det\mathbf{F} 会大量出现，所以为它引入一个专门的符号 JJ 是值得的，并且这个记号在大量的文献中被广泛使用

Bonet and Wood (2008), Marsden and Hughes (1994), Bower (2009)

ΨBW08\Psi_{\text{BW08}} 的 PK1 是

∂ΨBW08∂F=PBW08(F)=μ(F−1J∂J∂F)+λlog⁡JJ∂J∂F\frac{\partial\Psi_{\text{BW08}}}{\partial\mathbf{F}} = \mathbf{P}_{\text{BW08}}(\mathbf{F}) = \mu \left(\mathbf{F} - \frac{1}{J}\frac{\partial J}{\partial\mathbf{F}} \right) + \lambda \frac{\log J}{J}\frac{\partial J}{\partial\mathbf{F}}

这里新出现的项是 JJ 的梯度 ∂J∂F\frac{\partial J}{\partial\mathbf{F}}，在 3D 情况下，这一项可以用变形梯度 F\mathbf{F} 列的叉积来表示

F=[ ∣∣∣  f0f1f2 ∣∣∣]∂J∂F=[ ∣∣∣  f1×f2f2×f0f0×f1 ∣∣∣]\begin{align*} \mathbf{F} &= \begin{bmatrix} \ | &| &|\\ \ \ \mathbf{f_0} & \mathbf{f_1} & \mathbf{f_2}\\ \ | &| &|\\ \end{bmatrix} \\ \\ \frac{\partial J}{\partial\mathbf{F}} &= \begin{bmatrix} \ | &| &|\\ \ \ \mathbf{f_1} \times \mathbf{f_2} & \mathbf{f_2} \times \mathbf{f_0} & \mathbf{f_0} \times \mathbf{f_1}\\ \ | &| &|\\ \end{bmatrix} \end{align*}