【GAMES201学习笔记】01- 拉格朗日视角

1. 拉格朗日视角 vs 欧拉视角

拉格朗日视角：其中的元素跟随一起移动。

欧拉视角 ：其中的元素固定在对应的地方，网格的点对应的速度。

2. 弹簧质点系统

胡克定律

$$\begin{array}{rl}{\mathbf{f}}_{ij}& =-k(\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }_2{l}_{ij})\widehat{{\mathbf{x}}_{ij}}\\ {\mathbf{f}}_i& =\sum _j^{j\ne i}{\mathbf{f}}_{ij}\end{array}$$

牛顿第二定律

$$\begin{array}{rl}\frac{\partial {\mathbf{v}}_i}{\partial t}& =\frac{1}{m_i}{\mathbf{f}}_i\\ \frac{\partial {\mathbf{x}}_i}{\partial t}& ={\mathbf{v}}_i\end{array}$$

2.1 时间积分

前欧拉（显式）积分器 ：根据现有状态推测以后的状态

$$\begin{array}{rl}{\mathbf{v}}_{t+1}& ={\mathbf{v}}_t+\Delta t\frac{{\mathbf{f}}_t}{m}\\ {\mathbf{x}}_{t+1}& ={\mathbf{x}}_t+\Delta t{\mathbf{v}}_t\end{array}$$

半隐式欧拉积分器 ：根据现有状态推测以后的状态，速度根据推测出的速度计算

$$\begin{array}{rl}{\mathbf{v}}_{t+1}& ={\mathbf{v}}_t+\Delta t\frac{{\mathbf{f}}_t}{m}\\ {\mathbf{x}}_{t+1}& ={\mathbf{x}}_t+\Delta t{\mathbf{v}}_{t+1}\end{array}$$

2.2 代码

1. 计算速度：${\mathbf{v}}_{t+1}={\mathbf{v}}_t+\Delta t\frac{{\mathbf{f}}_t}{m}$
2. 与地面碰撞（防止计算位移之后陷入地下）
3. 计算位移：${\mathbf{x}}_{t+1}={\mathbf{x}}_t+\Delta t{\mathbf{v}}_{t+1}$

@ti.kernel
def substep():
    # Compute force and new velocity
    n = num_particles[None]
    for i in range(n):

        v[i] *= ti.exp(-dt * damping[None]) # 阻尼
        total_force = ti.Vector(gravity) * particle_mass
        
        for j in range(n):
            # 两个粒子有弹簧时
            if rest_length[i, j] != 0:
                x_ij = x[i] - x[j]
                total_force += -spring_stiffness[None] * (x_ij.norm() - rest_length[i, j]) * x_ij.normalized()
        v[i] += dt * total_force / particle_mass
        
    # Collide with ground
    for i in range(n):
        if x[i].y < bottom_y:
            x[i].y = bottom_y
            v[i].y = 0

    # Compute new position
    for i in range(num_particles[None]):
        x[i] += v[i] * dt

NOTE ：

• 其中 $\Delta t$ 需要满足一下情况：
  $$\Delta t=c\sqrt{\frac{m}{k}}(c\sim 1)$$

2.3 质点系统

隐式时间积分 ：

$$\begin{array}{rl}{\mathbf{x}}_{t+1}& ={\mathbf{x}}_t+\Delta t{\mathbf{v}}_{t+1}\\ {\mathbf{v}}_{t+1}& ={\mathbf{v}}_t+\Delta t{\mathbf{M}}^{-1}\mathbf{f}({\mathbf{x}}_{t+1})\end{array}$$

NOTE ：其中 $ \bold{M} $ 是质量矩阵。

消除 ${\mathbf{x}}_{t+1}$ ：

$${\mathbf{v}}_{t+1}={\mathbf{v}}_t+\Delta t{\mathbf{M}}^{-1}\mathbf{f}({\mathbf{x}}_t+\Delta t{\mathbf{v}}_{t+1})$$

线性化（牛顿法的一步）：

$${\mathbf{v}}_{t+1}={\mathbf{v}}_t+\Delta t{\mathbf{M}}^{-1}\left[\mathbf{f}({\mathbf{x}}_t)+\frac{\partial \mathbf{f}}{\partial \mathbf{x}}({\mathbf{x}}_t)\Delta t{\mathbf{v}}_{t+1}\right]$$

3. 形变

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

$${\mathbf{x}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}=\varphi ({\mathbf{x}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}})$$

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

$$\mathbf{F}:=\frac{\partial {\mathbf{x}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}}{\partial {\mathbf{x}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}}$$

Note ：

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

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

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

Note ：

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

4. 弹性体

4.1 超弹性体（Hyperelasticity）

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

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

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

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

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

Note ：

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

4.2 应力张量

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

4.2.1 三种常用的应力张量

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

4.2.2 关系式

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

Note：

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

4.2.3 牵引力

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

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

4.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(ϵ)$$

5. 常见的超弹性模型

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

6.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

x = ti.Vector(dim, dt=real, shape=n_nodes, needs_grad=True)
v = 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
        # 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 True:
    for s in range(30):
        # 调用 taichi 的自动微分器
        # 定义的损失函数：compute_total_energy()
        # 函数值：total_energy
        # 变量值：x
        # 微分结果(i 点上的力) ：x.grad[i]
        # f = - \partial (total_energy) / \partial x
        with ti.Tape(total_energy):
            compute_total_energy()
        ......

Note ：

• 在 FEM 中的势能（后续有详细说明）：

$$U(e)={\int }_e\mathbf{\Psi }(\mathbf{F}(\mathbf{x}))\mathrm{d}\mathbf{x}=V_e\mathbf{\Psi }({\mathbf{F}}_{\mathbf{e}})$$

6.2 (Fixed)Corotated 模型

• $\Psi (\mathbf{F})=\mu {\sum }_i({\sigma }_i-1{)}^2+\frac{\lambda }{2}(J-1{)}^2$
   
• $P(\mathbf{F})=\frac{\partial \Psi }{\partial \mathbf{F}}=2\mu (\mathbf{F}-\mathbf{R})+\lambda (J-1)J{\mathbf{F}}^{-T}$

Note ：

• ${\sigma }_i$ 是 $\mathbf{F}$ 的奇异值。

6.3 MPM 教程

《The Material Point Method for Simulating Continuum Materials》

6. 有限元基础

有限元法 ：建立离散模型的 Galerkin 离散格式。

6.1 线性四面体（三角形）有限元法

线性四面体有限元（用于弹性）假设 形变映射 $phi$ 是一个仿射变换，因此 形变梯度 $\mathbf{F}$ 在单个四面体单元内是恒定的，对单个元素来说：

$${\mathbf{x}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}=\mathbf{F}{\mathbf{x}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}+\mathbf{p}$$

对于每个元素 $e$ ，对能量密度函数求体积的积分，可以计算其弹性势能：

$$U(e)={\int }_e\mathbf{\Psi }(\mathbf{F}(\mathbf{x}))\mathrm{d}\mathbf{x}=V_e\mathbf{\Psi }({\mathbf{F}}_{\mathbf{e}})$$

• 其中 形变梯度 $\mathbf{F}$ 在元素 $e$ 上是个常数 ${\mathbf{F}}_{\mathbf{e}}$ ，即 $\mathbf{\Psi }({\mathbf{F}}_{\mathbf{e}})$ 也为常数，因此可以直接得到积分结果。

6.2 计算形变梯度

$${\mathbf{x}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}=\mathbf{F}{\mathbf{x}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}+\mathbf{p}$$

在 2D 三角形元素（三维空间中是四面体元素）中，设：

• 静止时的顶点位置： ${\mathbf{a}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}$ ，${\mathbf{b}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}$ ，${\mathbf{c}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}$ ；
• 变形后的顶点位置：${\mathbf{a}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}$ ，${\mathbf{b}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}$ ，${\mathbf{c}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}$ 。

因为在线性三角形元素中 $\mathbf{F}$ 是常数，则有：

$$\begin{array}{rl}{\mathbf{a}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}& =\mathbf{F}{\mathbf{a}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}+\mathbf{p}\\ {\mathbf{b}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}& =\mathbf{F}{\mathbf{b}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}+\mathbf{p}\\ {\mathbf{c}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}& =\mathbf{F}{\mathbf{c}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}+\mathbf{p}\end{array}$$

消除 $\mathbf{p}$ ：

$$\begin{array}{rl}({\mathbf{a}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}-{\mathbf{c}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}})& =\mathbf{F}({\mathbf{a}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}-{\mathbf{c}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}})\\ ({\mathbf{b}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}-{\mathbf{c}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}})& =\mathbf{F}({\mathbf{b}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}-{\mathbf{c}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}})\end{array}$$

现在 ${\mathbf{F}}_{2\times 2}$ 有四个线性约束：

$$\begin{array}{rl}\mathbf{B}& =[{\mathbf{a}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}-{\mathbf{c}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}|{\mathbf{b}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}-{\mathbf{c}}_{\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t}}{]}^{-1}\\ \mathbf{D}& =[{\mathbf{a}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}-{\mathbf{c}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}|{\mathbf{b}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}-{\mathbf{c}}_{\mathbf{d}\mathbf{e}\mathbf{f}\mathbf{o}\mathbf{r}\mathbf{m}\mathbf{e}\mathbf{d}}]\\ \mathbf{F}& =\mathbf{D}\mathbf{B}\end{array}$$

其中 $\mathbf{B}$ 在整个物理过程中是常数。因此可进行预计算。

B = ti.Matrix(dim, dim, dt=real, shape=n_elements)
vertices = ti.var(dt=ti.i32, shape=(n_elements, 3))

@ti.func
def compute_D(i):
    a = vertices[i, 0]
    b = vertices[i, 1]
    c = vertices[i, 2]
    return ti.Matrix.cols([x[b] - x[a], x[c] - x[a]])

@ti.kernel
def compute_B():
    for i in range(n_elements):
        B[i] = compute_D(i).inverse()

@ti.kernel
def compute_total_energy():
    for i in range(n_elements):
        D = compute_D(i)
        F = D @ B[i]
        ......

6.3 显式时间积分

$$\begin{array}{rl}{\mathbf{v}}_{t+1,i}& ={\mathbf{v}}_{t,i}+\Delta t\frac{{\mathbf{f}}_{t,i}}{m_i}\\ {\mathbf{x}}_{t+1,i}& ={\mathbf{x}}_{t,i}+\Delta t{\mathbf{v}}_{t+1,i}\end{array}$$

• ${\mathbf{v}}_{t,i}$ 和 ${\mathbf{x}}_{t,i}$ 存储在顶点中。

弹性势能对位置求导，结果的相反值即为顶点的受力：

$$\begin{array}{rl}{\mathbf{f}}_{t,i}& =-\frac{\partial U}{\partial {\mathbf{x}}_i}\\ & =-\sum _e\frac{\partial U(e)}{\partial {\mathbf{x}}_i}\\ & =-\sum _e{V}_e\frac{\partial \Psi ({\mathbf{F}}_e)}{\partial {\mathbf{F}}_e}\frac{\partial {\mathbf{F}}_e}{\partial {\mathbf{x}}_i}\\ & =-\sum _e{V}_e\mathbf{P}({\mathbf{F}}_e)\frac{\partial {\mathbf{F}}_e}{\partial {\mathbf{x}}_i}\end{array}$$

@ti.kernel
def integrate():
    for p in x:
        ...... # 碰撞检测
        # 显式时间积分
        v[p] = (v[p] + ((-x.grad[p] / node_mass) + ti.Vector([0, -10])) * dt) * math.exp(dt * -6)
        x[p] += dt * v[p]

• 其中：

$${\mathbf{v}}_{t+1,i}=\left({\mathbf{v}}_{t,i}+\frac{{\mathbf{f}}_{t,i}+m_{i}g}{m_i}\Delta t\right)e^{-6\Delta t}$$

6.4 隐式时间积分

$$\left[\mathbf{I}-\Delta t^2{\mathbf{M}}^{-1}\frac{\partial \mathbf{f}}{\partial \mathbf{x}}({\mathbf{x}}_t)\right]{\mathbf{v}}_{t+1}={\mathbf{v}}_t+\Delta t{\mathbf{M}}^{-1}\mathbf{f}({\mathbf{x}}_t)$$

力的微分计算：

$$\frac{\partial \mathbf{f}}{\partial \mathbf{x}}=-\frac{{\partial }^2\mathbf{\Psi }}{\partial {\mathbf{x}}^2}$$