【GAMES201学习笔记】线性有限元

1. 形变参数

1.1 形变映射

形变映射 $\varphi$ ：静止材料位置 $\mapsto$ 形变材料位置

$${\mathbf{x}}_{\mathbf{deformed}}=\varphi ({\mathbf{x}}_{\mathbf{rest}})$$

1.2 形变梯度

形变梯度 $\mathbf{F}$ ：

$$\mathbf{F}:=\frac{\partial {\mathbf{x}}_{\mathbf{deformed}}}{\partial {\mathbf{x}}_{\mathbf{rest}}}$$

Note ：

• 平移过程中形变梯度保持不变：${\varphi }_1=\varphi ({\mathbf{x}}_{\mathbf{rest}})$ 和 ${\varphi }_1=\varphi ({\mathbf{x}}_{\mathbf{rest}})+c$ 有相同的形变梯度

1.3 变形/静止体积比

变形/静止体积比 $J$ ：

$$J=\det (\mathbf{F})$$

Note ：

• 三维空间中，矩阵行列式的性质，即体积变化倍数

2. 弹性体

2.1 超弹性体（Hyperelasticity）

超弹性材料：应力 – 应变关系 由 应变能密度函数 定义 ：

$$\Psi =\Psi (\mathbf{F})$$

直观理解：$\Psi$ 是惩罚形变的势函数。

应力 ：材料的内部弹性力。

应变 ：现在用 形变梯度 $\mathbf{F}$ 代替。

Note ：

• $\Psi$ 是应变能量密度函数
• $\varphi$ 是形变映射

2.2 应力张量

应力：代表材料微元施加在其周围为材料微元的内力。

2.2.1 三种常用的应力张量

• Piola - Kirchhoff 应力张量 （PK1）：
  $$\mathbf{P}(\mathbf{F})=\frac{\partial \Psi (\mathbf{F})}{\partial \mathbf{F}}$$
  （计算简单，是在静止空间计算，需要变换到形变空间）
   
• 基尔霍夫应力（Kirchhoff stress） ：$\tau$
   
• 柯西应力张量（Cauchy stress tensor） ：$\sigma$
  （对称，因为角动量守恒）

2.2.2 关系式

• $\tau =J\sigma ={\mathbf{PF}}^T$
• $\mathbf{P}=J\sigma {\mathbf{F}}^{-T}$

Note：

• ${\mathbf{F}}^{-T}$ ：补偿材料变形
• 用 ${\mathbf{F}}^{-T}$ 替代 ${\mathbf{F}}^{-1}$ ：因为变换的是法向量 $\mathbf{n}$ 而不是 $\mathbf{x}$

2.2.3 牵引力

• $\mathbf{t}={\sigma }^T\mathbf{n}$

直观来说，将法向量乘以应力张量即可获得材料向周围微元施加的力。

2.3 弹性模量（各向同性材料）

• 杨氏模量 ：应力张量与应变张量的比值
  $$E=\frac{\sigma }{ϵ}$$

• 体积模量 ：压强关于体积的变化，常用于液体
  $$K=-V\frac{dP}{dV}$$

• 泊松比 ：
  $$\nu \in [0.0,0.5)$$

Note ：

• 辅助项具有负泊松比；
• 泊松比为 0 时，拉长物体时，截面积不会发生变化；
• 泊松比较大时，物体会尽量保持体积不变，例如在拉长物体时，物体会变细。

拉梅常数（Lamé parameters） ：

• Lamé 第一参数： $\lambda$
  表示材料的压缩性，等价与体弹性模量或者杨氏模量
   
• Lamé 第二参数： $\mu$
  表示材料的剪切模量，用 G 表示

换算公式 ：

• $K=\frac{E}{3(1-2\nu )}$ （常用于可压缩液体）
   
• $\lambda =\frac{E\nu }{(1+\nu )(1-2\nu )}$
   
• $\mu =\frac{E}{2(1+\nu )}$

广义胡克定律 ：均匀和各向同性的材料

$$\sigma =2\mu ϵ+\lambda tr(ϵ)$$

3. 常见的超弹性模型

• 经典的 MPM 方法将每个粒子看做材料的一个微元，材料的本构模型会有一个能量密度函数 $\Psi$ ；
• 对能量密度函数 $\Psi$ 关于整个模型求积分，得到整个材料的势能；
• 势能对材料点的形变梯度进行求导：$P(\mathbf{F})=\frac{\partial \Psi }{\partial \mathbf{F}}$ ；
• 物理意义上来说，势能对位置求导的结果就是力，$P(\mathbf{F})$ 可以看做材料点的受力。

3.1 Neo-Hookean 模型

适用于各向同性材料

• $\Psi (\mathbf{F})=\frac{\mu }{2}{\sum }_i[({\mathbf{F}}^T\mathbf{F}{)}_{ii}-1]-\mu \log (J)+\frac{\lambda }{2}{\log }^2(J)$
   
• $P(\mathbf{F})=\frac{\partial \Psi }{\partial \mathbf{F}}=\mu (\mathbf{F}-{\mathbf{F}}^T)+\lambda \log (J){\mathbf{F}}^{-T}$

因为 Neo-Hookean 模型容易造成能量流失，这时可以考虑 Corotated 模型。

FEM 中应用 Neo-Hookean 模型的示例代码 ：

dim = 2

E, nu = 1000, 0.3
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

element_V = 0.01

pars = ti.Vector(dim, dt=real, shape=n_nodes, needs_grad=True)
vels = ti.Vector(dim, dt=real, shape=n_nodes)

total_energy = ti.var(dt=real, shape=(), needs_grad=True)

# 计算势能
@ti.kernel
def compute_total_energy():
    for i in range(n_elements):
        # get F：见 4.2
        ......
        # NeoHookean
        I1 = (F @ F.transpose()).trace()
        # 防止 J 的值过小引起错误
        J = max(0.2, F.determinant())
        element_energy_density = 0.5 * mu * (I1 - dim) - mu * ti.log(J) + 0.5 * la * ti.log(J)**2
        total_energy[None] += element_energy_density * element_V

# 渲染循环
while