【GAMES201学习笔记】物质点法入门

1. 流体力学中的两种描述

1.1 拉格朗日描述

1.1.1 什么是拉格朗日描述

拉格朗日描述 的观察对象是组成物体的 质点 ， 粒子（particle） ，或者说 流体微元 ：

• 随着物体变化，对应 “粒子” 的位置也会随之变化，跟踪的观察对象是单个 “粒子” 的属性变化；
• 拉格朗日描述最显著的特征就是观察对象会发生位置或者形状的变换。

1.1.2 常见的使用拉格朗日描述的方法

• 弹簧质点模型（观察质点的变换）
• 有限元法（观察有限微元的变换）

1.1.3 拉格朗日描述擅长的地方

• 拉格朗日描述擅长对独立粒子的操作；
• 物体的变换操作：这时对粒子施加相应的变换即可；
• 粒子的质量和动量守恒是简单的；
• 定义物体的材质：定义单个粒子相关参数。

1.1.4 拉格朗日描述不擅长的地方

由于粒子经过变换之后，对应的位置和形变都是不确定的，因此虽然可以查询之前已经建立了 邻接信息 的粒子之间的关系（当然，在拉格朗日中找邻接粒子本身是一个困难的问题）；但查询未确立 邻接信息 的粒子之间的关系（例如任意两个粒子之间的距离），会异常困难。

这导致当有查询 未邻接的两个粒子 之间关系的需求时，会遇到一定的麻烦：例如在使用有限元进行仿真时，如果没有合理处理自相交问题，会发生模型穿透的情况。

弹性方块的 FEM 模拟 ：

弹性平面的 FEM 模拟 ：

• 由于仅做了邻接粒子之间、粒子与球之间的碰撞检测，所以在平面折叠时发生了穿透。

不过目前在有限元领域，对于碰撞问题，也已经有了对应的解决方案：

• [SIGGRAPH 2020] Incremental Potential Contact (IPC)

1.2 欧拉描述

1.2.1 什么是欧拉描述

欧拉描述 的观察对象是物体所在 空间的场 ， 网格（Grid） ，或者说 流体元胞（Cell） ：

• 网格是在空间中的绝对参照系的微元，随着物体变化，网格的位置和形状不会随之变化，仅物体映射到网格上的信息发生改变；
• 欧拉描述最显著的特征就是观察对象不会发生位置或者形状的变换。

1.2.2 欧拉描述擅长的地方

• 拉格朗日描述擅长处理粒子之间的相对关系；
• 物体性质的变化：质量密度、速度、温度、熵、焓，甚至单位流体中的磁通量；
• 物体内部的压力压强计算，可以很快的获得能量密度函数；
• 边界碰撞检测：欧拉描述可以很快的获得任意粒子之间的相对关系信息，或者说在欧拉视角下物体自然而然地在进行着碰撞检测。

1.2.3 欧拉描述不擅长的地方

• 物体的变换：显然，拉格朗日描述擅长的地方就是欧拉描述不擅长的地方，在物体发生变换的过程中，需要更新网格和相邻网格的信息，这个过程是较为繁琐的。

2. Particle-In-Cell methods(PIC)

Note ：

• 在混合拉格朗日 - 欧拉方法中，“粒子”是一等公民，“网格”是用来计算场信息的中间量。

2.1 粒子信息映射到网格（P2G）

2.1.1 思路

将粒子信息映射到网格上，或者说从一组粒子重建整个场的信息，这里采用了一个 SPH 中处理类似问题的方法：

• 使用一个径向对称的 核函数 （有时称为 基函数 ），来确定一个粒子对周边空间的 “影响（贡献）”，即一种加权计算 ；
• 其中粒子的影响是局部化的，核函数通常选择 有限支撑 ，即存在一个最大支撑距离；
• 核函数是归一化的： ${\oint }_R\hat{N}(\mathbf{r})d\mathbf{r}=1$ ，其中 $R$ 是整个支撑集。

B 样条核函数 ：

在二维空间下，仅考虑粒子周围的九个网格节点。

2.1.2 代码实现

mpm88.py Line 30,42 （除去了 APIC 的部分）：

for p in x:
        Xp = x[p] / dx
        base = int(Xp - 0.5)
        fx = Xp - base
        # Quadratic B-spline
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]

        for i, j in ti.static(ti.ndrange(3, 3)):
            offset = ti.Vector([i, j])
            weight = w[i].x * w[j].y
            grid_v[base + offset] += weight * (p_mass * v[p])
            grid_m[base + offset] += weight * p_mass

Note ：

• 粒子在中间网格里；
• base ：粒子所在周边网格的网格点 ${\mathbf{p}}_{[0,0]}$ ，坐标 (0.5 * dx, 0.5 * dx) ；
• fx ：粒子距离网格点 ${\mathbf{p}}_{[0,0]}$ 的距离；
   
• 距离中间的网格点：$\omega =\frac{3}{4}-|({\mathbf{x}}_p-{\mathbf{x}}_{base})-1{|}^2$
   
• 距离靠左下的网格点：$\omega =\frac{1}{2}(\frac{3}{2}-|{\mathbf{x}}_p-{\mathbf{x}}_{base}|{)}^2$
   
• 距离靠右上的网格点：$\omega =\frac{1}{2}(2-(\frac{3}{2}-|{\mathbf{x}}_p-{\mathbf{x}}_{base}|){)}^2=\frac{1}{2}(|{\mathbf{x}}_p-{\mathbf{x}}_{base}|-\frac{1}{2}{)}^2$

mpm88.py Line 43,46：

# 对速度做一个网格质量的加权平均
for i, j in grid_m:
        if grid_m[i, j] > 0:
            grid_v[i, j] /= grid_m[i, j]
        grid_v[i, j].y -= dt * gravity

2.2 网格信息叠加到粒子（G2P）

2.1.1 思路

在计算完网格上各种场的对应信息之后，将周围网格的信息叠加到原粒子上。

2.2.2 代码实现

mpm88.py Line 55,68 （除去了 APIC 的部分）：

for p in x:
        Xp = x[p] / dx
        base = int(Xp - 0.5)
        fx = Xp - base
        # Quadratic B-spline 
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]

        new_v = ti.Vector.zero(float, 2)

        for i, j in ti.static(ti.ndrange(3, 3)):
            offset = ti.Vector([i, j])
            weight = w[i].x * w[j].y
            new_v += weight * g_v

        v[p] = new_v
        x[p] += dt * v[p]

2.3 PIC 的局限性

主要是在 G2P 的过程中，粒子信息是由网格信息叠加而成的，存在大量的信息缺失，直观来说就是出现了能量损耗，主要在旋转和切变的过程（主要是角动量损耗及其严重）。

• 以速度场为例，粒子拥有 2 个自由度，而一个 3×3 的网格点有 18 个自由度：

3. Affine Particle-In-Cell（APIC）

3.1 APIC 简介

通过 2.3 可知，PIC 在场信息到粒子信息的过程中存在信息缺失的情况，这时不妨在这个过程中存储更多的场的信息：

符 号	含 义
${\mathbf{v}}_0$ 、 ${\mathbf{v}}_1$	均匀常量速度场，描述粒子原本的 2 个自由度的速度信息，不会随着 ${\mathbf{x}}_p$ 位置的变化而变化。
${\mathbf{C}}_{00}$	随着粒子在 x 分量的变化，在 x 分量上线性变化的速度场，场的分量 C[0] 随着 x[p][0] 变化的情况。
${\mathbf{C}}_{10}$	随着粒子在 x 分量的变化，在 y 方向上线性变化的速度场，场的分量 C[1] 随着 x[p][0] 变化的情况。
${\mathbf{C}}_{01}$	随着粒子在 y 分量的变化，在 x 方向上线性变化的速度场，场的分量 C[0] 随着 x[p][1] 变化的情况。
${\mathbf{C}}_{11}$	随着粒子在 y 分量的变化，在 y 方向上线性变化的速度场，场的分量 C[1] 随着 x[p][1] 变化的情况。

3.2 APIC 的角动量是守恒的

数学证明见：

• Tech report: the affine Particle-in-Cell method <5.3.1> <5.3.2>

3.3 代码实现

3.3.1 P2G

mpm88.py Line 30,42（除去压力计算部分）：

x = ti.Vector.field(2, float, n_particles)
v = ti.Vector.field(2, float, n_particles)
C = ti.Matrix.field(2, 2, float, n_particles)

for p in x:
        Xp = x[p] / dx
        base = int(Xp - 0.5)
        fx = Xp - base
        # Quadratic B-spline
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]

        affine = C[p]

        for i, j in ti.static(ti.ndrange(3, 3)):
            offset = ti.Vector([i, j])
            # 网格点位置和粒子位置的偏差值
            dpos = (offset - fx) * dx
            weight = w[i].x * w[j].y

            grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)
            grid_m[base + offset] += weight * p_mass

$$\begin{array}{rl}{\mathbf{g}}_v& =\left[\begin{array}{c}{\mathbf{g}}_v.x\\ {\mathbf{g}}_v.y\end{array}\right]\\ & =\left[\begin{array}{cc}{\mathbf{C}}_{00}& {\mathbf{C}}_{01}\\ {\mathbf{C}}_{10}& {\mathbf{C}}_{11}\end{array}\right]\left[\begin{array}{c}{\mathbf{d}}_{pos}.x\\ {\mathbf{d}}_{pos}.y\end{array}\right]\\ & =\left[\begin{array}{c}{\mathbf{C}}_{00}\cdot {\mathbf{d}}_{pos}.x+{\mathbf{C}}_{01}\cdot {\mathbf{d}}_{pos}.y\\ {\mathbf{C}}_{10}\cdot {\mathbf{d}}_{pos}.x+{\mathbf{C}}_{11}\cdot {\mathbf{d}}_{pos}.y\end{array}\right]\end{array}$$

Note ：

• 其中网格点速度 g_v 、粒子位移的差分值 dpos ，与速度场 C 之间的关系，见 3.1 中的表格。

3.3.2 G2P

mpm88.py Line 30,42（除去压力计算部分）：

for p in x:
        Xp = x[p] / dx
        base = int(Xp - 0.5)
        fx = Xp - base
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]

        new_v = ti.Vector.zero(float, 2)
        new_C = ti.Matrix.zero(float, 2, 2)

        for i, j in ti.static(ti.ndrange(3, 3)):
            offset = ti.Vector([i, j])
            # 网格点位置和粒子位置的偏差值
            dpos = (offset - fx) * dx
            weight = w[i].x * w[j].y
            # 网格点的速度
            g_v = grid_v[base + offset]

            new_v += weight * g_v
            new_C += 4 * weight * g_v.outer_product(dpos) / dx**2

        v[p] = new_v
        x[p] += dt * v[p]
        C[p] = new_C

$$\begin{array}{rl}{\mathbf{C}}_{new}& =\left[\begin{array}{cc}{\mathbf{C}}_{00}& {\mathbf{C}}_{01}\\ {\mathbf{C}}_{10}& {\mathbf{C}}_{11}\end{array}\right]\\ & =\frac{4w}{{\Delta x}^2}{\mathbf{g}}_v\otimes {\mathbf{d}}_{pos}\\ & =\frac{4w}{{\Delta x}^2}\left[\begin{array}{c}{\mathbf{g}}_v.x\\ {\mathbf{g}}_v.y\end{array}\right]\left[\begin{array}{cc}{\mathbf{d}}_{pos}.x& {\mathbf{d}}_{pos}.y\end{array}\right]\\ & =\frac{4w}{{\Delta x}^2}\left[\begin{array}{cc}{\mathbf{g}}_v.x\cdot {\mathbf{d}}_{pos}.x& {\mathbf{g}}_v.x\cdot {\mathbf{d}}_{pos}.y\\ {\mathbf{g}}_v.y\cdot {\mathbf{d}}_{pos}.x& {\mathbf{g}}_v.y\cdot {\mathbf{d}}_{pos}.y\end{array}\right]\end{array}$$

Note ：

• 其中网格点速度 g_v 、粒子位移的差分值 dpos ，与速度场 C 之间的关系，见 3.1 中的表格。

3.4 如果需要更精确的解

APIC 可以类比于泰勒一级展开，如果有追求更高精度的需要（二级展开等），则可以使用 ployPIC ：

• [SIGGRAPH Asia 2017] A polynomial particle-in-cell method

4. Material Point Method（MPM）

4.1 物质点法简介

在解决了欧拉-拉格朗日之间的信息传递问题之后（APIC），就需要真正开始进行模拟，这时需要额外增加一些信息（形变梯度 $\mathbf{F}$ ，体积应变比 $\mathbf{J}$ ，杨氏模量 $E$ 等）

最早的一篇关于 MPM 论文：

• [SIGGRAPH 2013] A material point method for snow simulation

距 2018 年，MPM 的相关论文：

4.2 MLS-MPM

P2G ：

• 使用（仿射 Affine）速度更新粒子；
• 将质量和动量分散到附近的 3×3×3 个节点。

对于每个网格点 ：

• 用动量除以质量得到速度；
• 应用重力和边界条件

G2P ：

• 从 3×3×3 个节点收集 速度/仿射速度（Affine）。

4.3 代码实现

4.3.1 相关参数

mpm88.py Line 2, 23：

import taichi as ti
ti.init(arch=ti.gpu)

n_particles = 8192
n_grid = 128
dx = 1 / n_grid
dt = 2e-4

p_rho = 1
p_vol = (dx * 0.5)**2
p_mass = p_vol * p_rho
gravity = 9.8
bound = 3
E = 400

x = ti.Vector.field(2, float, n_particles)
v = ti.Vector.field(2, float, n_particles)
C = ti.Matrix.field(2, 2, float, n_particles)
# 体积比（当前体积/初始体积）：
# J < 1 ：粒子压缩
# J > 1 ：粒子膨胀
J = ti.field(float, n_particles)

grid_v = ti.Vector.field(2, float, (n_grid, n_grid))
grid_m = ti.field(float, (n_grid, n_grid))

4.3.2 P2G

mpm88.py Line 30, 42：

@ti.kernel
def substep():
    # init grid
    ...... 

    # P2G
    for p in x:
        Xp = x[p] / dx
        base = int(Xp - 0.5)
        fx = Xp - base
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
        stress = -dt * 4 * E * p_vol * (J[p] - 1) / dx**2
        # 动量的仿射信息
        affine = ti.Matrix([[stress, 0], [0, stress]]) + p_mass * C[p]
        for i, j in ti.static(ti.ndrange(3, 3)):
            offset = ti.Vector([i, j])
            dpos = (offset - fx) * dx
            weight = w[i].x * w[j].y
            grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)
            grid_m[base + offset] += weight * p_mass

在 MPM 中，添加了一项 stress ，表示粒子受到压缩之后推开周围粒子的力（实际上是在网格层面操作）：

$$(mv{)}_{stress}=-\frac{4}{\Delta x^2}\Delta tE{\mathbf{p}}_{vol}(J_{\mathbf{p}}-1)$$

4.3.3 边界条件

mpm88.py Line 43, 54：

    for i, j in grid_m:
        if grid_m[i, j] > 0:
            grid_v[i, j] /= grid_m[i, j]
        grid_v[i, j].y -= dt * gravity
        if i < bound and grid_v[i, j].x < 0:
            grid_v[i, j].x = 0
        if i > n_grid - bound and grid_v[i, j].x > 0:
            grid_v[i, j].x = 0
        if j < bound and grid_v[i, j].y < 0:
            grid_v[i, j].y = 0
        if j > n_grid - bound and grid_v[i, j].y > 0:
            grid_v[i, j].y = 0

4.3.4 G2P

mpm88.py Line 55, 72：

    for p in x:
        Xp = x[p] / dx
        base = int(Xp - 0.5)
        fx = Xp - base
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
        
        new_v = ti.Vector.zero(float, 2)
        new_C = ti.Matrix.zero(float, 2, 2)

        for i, j in ti.static(ti.ndrange(3, 3)):
            offset = ti.Vector([i, j])
            dpos = (offset - fx) * dx
            weight = w[i].x * w[j].y
            g_v = grid_v[base + offset]
            new_v += weight * g_v
            new_C += 4 * weight * g_v.outer_product(dpos) / dx**2

        v[p] = new_v
        x[p] += dt * v[p]
        J[p] *= 1 + dt * new_C.trace()
        C[p] = new_C

$$J_{\mathbf{p}-new}=J_{\mathbf{p}}(1+({\mathbf{C}}_{00}+{\mathbf{C}}_{11})\Delta t)$$

Note ：

• ${\mathbf{C}}_{00}$ 和 ${\mathbf{C}}_{11}$ 同样可以用来表示粒子的膨胀（正值）和收缩（负值）。

参考课程

文中相关资源来自课程：

• [Chinagraph 2020] 用太极编写物理引擎 (5) 物质点法