同一组粒子速度经过网格往返，PIC 可能把旋转磨平，FLIP 可以保住粒子原有速度；那么 APIC 多保存一个矩阵，到底解决了什么问题？关键不只是“回传时少损失一些速度”，而是网格能否看见粒子所代表的局部速度结构。

速度回传篇已经解释 PIC 的覆盖、FLIP 的增量和投影耗散。本文在此基础上比较三种方法的双向传递，以 Jiang 等人 2015 年的 APIC 原论文第 5 节为主要依据，并给出可复现的二维纯传递算例。读者需要了解形函数插值、矩阵乘法和外积；速度梯度与变形梯度的物理含义可参照运动学篇。

这里采用速度各分量存放在同一节点的规则网格，不讨论 MAC 交错网格。文中的算例只运行粒子与网格之间的映射，不推进位置、不施加内力，也不求解压力。实验室的 PIC／FLIP／APIC 对比面板提供相同算例的交互版本。它们是传递算子的验证，不是完整 MPM 动力学模拟；弹性杆实验仍只提供 PIC／FLIP。

一、先把三种方法放在同一张图里

方法粒子携带的速度状态P2G：粒子到网格G2P：网格到粒子
PIC平移速度加权汇集粒子动量用网格插值速度覆盖
FLIP平移速度与这里的 PIC 相同保留旧粒子速度，加上网格速度增量
APIC平移速度与局部仿射状态汇集平移与仿射速度贡献从网格重建平移与仿射状态

PIC 让速度信息完整地经过网格筛选。FLIP 给旧粒子速度保留了一条不经过网格重建的路径，因此网格没有分辨出来的粒子模式也能留下。APIC 则丰富粒子的表示：网格回传的信息不必全部压缩成一个速度向量，还可以留下局部线性变化。

这也解释了一个容易混淆的地方：APIC 的平移速度回传与 PIC 相同，不代表两者是同一种方法。APIC 接收到的网格场可能已经不同，而且它还回传一个矩阵，供下一次 P2G 使用。

二、统一符号与传递阶段

用表示粒子、表示节点。粒子质量为，位置为，节点位置为，权重和相对位置分别为

无量纲，位置单位为米。本文的推导要求使用完整的插值支撑域，并满足常量与线性再现：

讨论加权平均的耗散时，还要求权重非负。所有除以节点质量的操作只在的节点进行。

为避免把本步映射值与上一时刻残留网格值混淆，用上标表示本步 P2G 刚得到的节点速度，用表示本步节点动力学及约束处理后的速度。一次往返中，G2P 仍使用本次 P2G 的位置与权重。下文省略时间步编号，但不意味着可以混用两个时刻的权重。

PIC 与 FLIP 的共同入口

两者在这里都使用集中质量映射：

节点经过力、压力或边界处理后得到，两种回传分别为

按照这里的阶段定义，增量包括映射之后施加的约束修正。实际程序若先约束再保存旧网格速度，或只回传加速度，必须重新说明增量口径；不能把不同实现的边界处理悄悄视为同一公式。

纯传递算例中没有任何节点更新，因此，FLIP 的粒子速度严格不变。这是一条代数结论，不需要靠长时间模拟证明，也不能据此宣布 FLIP 已经准确表示了网格速度场。

三、APIC 多保存的矩阵是什么

APIC 用每个粒子附近的仿射速度表示代替局部常量表示：

此处左边表示粒子所代表的局部速度函数，右边的是该函数在粒子中心的值。为矩阵，是空间维数，单位为；它可以表示局部旋转、剪切和伸缩，而不只是角速度。

但是，首先是传递状态。通常用于本构和变形更新的物理速度梯度是

其中列向量外积使分量满足。则是无量纲的变形梯度，是材料变形历史的状态量。不能因为和都像“速度导数”，就在任意插值、任意更新顺序下互换，更不能用其中一个替代。

对于下面这组公式和张量积线性插值，在非退化位置可利用权重恒等式使重建的等于相应的插值速度梯度；这是一种特定关系，不是对所有 APIC 实现的通用赋值规则。

1. P2G：把局部仿射速度送到节点

节点质量不变，动量映射改为

这不是直接给网格加旋转力，而是改变粒子向各节点贡献的速度。由式，同一粒子的仿射贡献在所有节点求和为零，因此不会凭空增加总线动量。

2. G2P：同时重建速度与矩阵

引入位置二阶矩和速度—位置矩：

的单位为，的单位为。在可逆且位置固定的一次往返内，

原论文以作为携带状态，并在 P2G 中使用。本文的固定位置算例存与存等价；若位置随时间变化，也可能改变，必须明确矩阵对应哪个几何状态，不能机械照搬固定位置循环。

与加权最小二乘拟合的矩阵相同：重建的是网格速度相对粒子中心的线性部分。因为相对位置加权和为零，中也可以用代替。

3. 不要对所有形函数硬编码同一个逆矩阵

在完整支撑的均匀笛卡尔网格上，令网格间距为，标准张量积二次 B 样条满足，三次 B 样条满足。这时求逆可以化为常数缩放，是单位矩阵。

线性插值则不同。一维粒子在单元内的位置写成，，其二阶矩为

它随位置变化，在节点处为零；多维张量积线性插值的相应方向也可能退化。原论文第 5.3 节给出利用权重梯度恒等式避免显式求逆的处理思路。工程实现应选定一致的表示及退化处理，而不是遇到奇异矩阵就随手加一个小常数，之后仍声称保留了原来的再现和守恒性质。

四、为什么 APIC 能再现仿射速度场

设真实速度是全局仿射场

其中是常量速度，是常量矩阵。初始化粒子时，令、。那么每个粒子对节点的局部预测都相同：

因此有质量的节点在 P2G 后得到精确仿射速度。在没有节点更新的 G2P 中，线性再现给出精确粒子中心速度，且

这就证明了式与式的仿射再现性。它要求粒子初始矩阵与目标场一致、权重满足前述矩条件、二阶矩可逆或采用适当的等价处理。

如果初始，第一次 P2G 就仍与 PIC 一样，不能期待 APIC 从少量中心速度中自动恢复已经丢失的信息。后续虽能重建矩阵，也不等于第一次损失从未发生。

五、三个可以复现的纯传递算例

取边长的正方形单元，四个节点位于四角，四个等质量粒子位于

每个粒子质量为。采用双线性权重，四个节点质量均为，每个粒子的二阶矩都是。四角节点不施加固定边界；这里是完整权重支撑，不是一个固定四边的力学问题。

用单元中心定义三个速度场：

算例速度场APIC 初始矩阵
均匀平移
简单剪切，其余为零
刚体旋转，其余为零

取剪切率、逆时针角速度。三种方法具有相同的粒子中心速度；APIC 的矩阵按解析场初始化，它多获得的仿射信息是本次再现测试的已知输入，不是从数据拟合出来的额外成绩。

同一旋转场下，目标节点速度与 APIC 映射一致，PIC 和 FLIP 映射到节点的速度只有目标的四分之一

图中方形为同一个单元，圆点为四个粒子，箭头位于四个网格节点；三幅图使用相同箭头比例，方向均为逆时针。粒子中心速度相同，不代表网格收到的场相同。

将所有采样点速度分量拼接为向量，用欧氏范数定义误差：

这里的表示全部节点组成的数据集合，而不是某个单独节点。这三个场的分母均不为零。直接运行下节代码得到：

算例方法P2G 节点场误差往返后粒子速度误差
平移PIC／FLIP／APIC0%0%
剪切PIC75%75%
剪切FLIP75%0%
剪切APIC0%0%
旋转PIC75%75%
旋转FLIP75%0%
旋转APIC0%0%

APIC 三个算例回传的矩阵也都与初始矩阵一致。剪切和旋转中，PIC／FLIP 映射得到的节点速度只有目标值的四分之一；PIC 回传后，粒子中心速度也缩为原来的四分之一。FLIP 因节点增量为零而保留原值，但没有让本次网格场变准确。

这张表比较的是仿射再现，不是完整动力学精度排名。如果只看，FLIP 和 APIC 会被判为相同；加上，才能看到 APIC 丰富局部表示的意义。即使两者在这里都不改变粒子速度，下一步计算内力或压力时面对的网格状态也可能不同。

在交互实验中，选择“刚体旋转”、初始 C 为“解析矩阵”、粒子距边界为 0.25 m，并回到第一轮，即可复现表中旋转结果。先看“网格速度”，再切换“回传粒子速度”，观察 FLIP 两种误差的差别；点击“再往返一轮”，观察 PIC 在同一固定几何上继续损失速度。将初始 C 改成零，可以验证 APIC 第一次映射与 PIC 相同、后续重建矩阵不能自动补回原始仿射场。往返轮数不是物理时间，图中的粒子不会沿箭头运动。

六、旋转守恒时，究竟在统计什么

普通点粒子关于单元中心的轨道角动量为

APIC 还携带局部仿射贡献。对于本文二维固定位置表示，总角动量应写成

原因可由式展开得到：局部贡献是的出平面分量。它不是额外施加的物理转矩，而是粒子仿射表示所携带的角动量。

旋转算例的结果如下，单位统一为：

方法及统计口径初始粒子态P2G 后网格态G2P 后粒子态
PIC：粒子只计轨道项0.50.50.125
FLIP：粒子只计轨道项0.50.50.5
APIC：粒子计轨道与仿射项222

APIC 初始轨道项同样为，但仿射项为，所以初始总量为。不能把表中的与当成相同初始总角动量下的性能比较。这张表首先说明各自表示内的传递守恒；APIC 的额外状态改变了总量的统计定义。

在前述矩条件和相容的传递公式下，APIC 的角动量守恒不依赖这四个粒子的对称巧合。反过来，集中质量 FLIP 在这里保持角动量，只因为节点根本没有更新；它不能证明任意力更新和边界条件下都严格守恒。外力矩、接触、约束和时间推进的误差也必须另算。

七、完整复现代码

以下代码需要 Python 3 和 NumPy。它使用米、秒、千克为单位，直接构造四粒子算例，没有调用博客实验室，也没有预设结果曲线。代码中的质量都为一，因此部分质量因子可以省略。输出百分误差、APIC 矩阵误差，以及旋转时的三个角动量值；浮点输出的末位差异不影响表中结果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import numpy as np

p = np.array([[.25, .25], [.75, .25], [.25, .75], [.75, .75]])
g = np.array([[0., 0.], [1., 0.], [0., 1.], [1., 1.]])
r = g[None, :, :] - p[:, None, :]
w = np.prod(np.where(g[None, :, :] == 1,
                     p[:, None, :], 1 - p[:, None, :]), axis=2)
D = np.einsum('pi,pij,pik->pjk', w, r, r)
invD = np.linalg.inv(D)  # 本例非退化且几何固定，只需计算一次


def transfer(v, C, mode):
    local = np.broadcast_to(v[:, None, :], r.shape)
    if mode == 'APIC':
        local = local + np.einsum('pab,pib->pia', C, r)
    vg = np.einsum('pi,pia->ia', w, local) / w.sum(axis=0)[:, None]
    # 无节点更新，所以 FLIP 的网格速度增量为零。
    vp = v.copy() if mode == 'FLIP' else w @ vg
    B = np.einsum('pi,ia,pib->pab', w, vg, r)
    Cp = B @ invD if mode == 'APIC' else np.zeros_like(C)
    return vg, vp, Cp


def orbital(x, v):
    return np.sum((x[:, 0] - .5) * v[:, 1]
                  - (x[:, 1] - .5) * v[:, 0])


def affine(C):
    B = C @ D
    return np.sum(B[:, 1, 0] - B[:, 0, 1])


fields = {
    'translation': (np.zeros((2, 2)), np.array([1., -.5])),
    'shear': (np.array([[0., 1.], [0., 0.]]), np.zeros(2)),
    'rotation': (np.array([[0., -1.], [1., 0.]]), np.zeros(2)),
}
for name, (A, b) in fields.items():
    v = (p - .5) @ A.T + b
    exact = (g - .5) @ A.T + b
    for mode in ['PIC', 'FLIP', 'APIC']:
        C = np.broadcast_to(A, (4, 2, 2)).copy()
        vg, vp, Cp = transfer(v, C, mode)
        eg = 100 * np.linalg.norm(vg - exact) / np.linalg.norm(exact)
        ep = 100 * np.linalg.norm(vp - v) / np.linalg.norm(v)
        print(name, mode, 'grid %:', eg, 'particle %:', ep)
        if mode == 'APIC':
            print('matrix error:', np.linalg.norm(Cp - C))
        if name == 'rotation':
            initial = orbital(p, v) + (affine(C) if mode == 'APIC' else 0)
            final = orbital(p, vp) + (affine(Cp) if mode == 'APIC' else 0)
            print('angular momentum:', initial, orbital(g, vg), final)

八、守恒、耗散和实现代价要分开判断

仿射再现不等于任意场无耗散

APIC 能精确传递上述仿射场，但一般速度场还包含更细的变化，局部仿射状态不能表示所有模式。原论文第 8 节的能量实验与第 9 节的局限讨论也表明，APIC 仍然会耗散，FLIP 在一些情况下更能保住总能量；部分被保住的能量可能位于网格不能分辨的模式中。

因此不能把“APIC 保角动量”推导为“APIC 严格保动能”。也不要直接拿只统计的中心速度动能，去证明扩充了仿射状态的两种表示之间完全没有能量损失。需要先说明能量度量，再分开测量 P2G、网格更新和 G2P 的变化。

FLIP 保留细节，也会保留网格看不见的模式

如果某组粒子速度贡献映射到节点后相互抵消，网格就无法基于这部分速度作出对应的动力学响应。FLIP 不会仅因网格看不见它而将其删除；这既可能保留有用细节，也可能让不受网格控制的模式持续存在。APIC 丰富的是局部仿射表示，并没有消除所有未解析模式或所有数值不稳定性。

PIC／FLIP 混合不是 APIC 的另一种写法

定义 FLIP 比例，常见混合为

它调节覆盖与保留旧速度的比例，没有新增仿射矩阵，也没有自动改变式所描述的 P2G 表示。原论文也讨论 APIC／PIC 混合，但那涉及缩放仿射贡献，不能与这个公式混用。

多一个矩阵，成本不是固定倍数

相对于共同的粒子速度状态，完整 APIC 仿射矩阵在二维增加个标量、三维增加个标量，并增加局部矩阵—向量乘法和矩积汇集。存或是实现选择，通常不需要把两者都永久存一份。

实际运行时间还受插值支撑节点数、缓存布局、粒子数量、网格求解器和时间步数影响。FLIP 通常需要保留旧网格速度或等价的增量状态；不同代码的数据布局也不同。本文没有完整三方法求解器的计时实验，因此不报告性能倍数。

九、在 MPM 中怎样选择下一步检查

把问题按所在环节拆开，比直接给三种方法排序更有用：

  • 往返后平移或仿射场不对：先查质量归一化、矩条件、外积顺序、矩阵初始化及新旧网格速度。
  • 旋转明显衰减：同时查 G2P 和角动量统计口径，不要漏掉 APIC 仿射项。
  • 能量看似保持，但粒子出现细尺度乱动：检查网格未解析模式，不能只看总能量。
  • 粒子跨单元时应力异常：检查形函数梯度、支撑域和粒子积分；改用 APIC 并不自动消除网格穿越误差。
  • 传递检查通过，但动力学仍不准：继续用弹性杆解析验证一类问题检查相位、时间步、应力更新和边界。

PIC／FLIP／APIC 主要回答“速度信息怎样在粒子和网格之间传递”；uGIMP／CPDI 主要涉及粒子域与插值处理，B 样条则是基函数选择。它们处在不同的设计维度，可以在一致的数学条件下组合，但不能把这里的二阶矩常数和守恒证明不加检查地搬到另一套权重上。

本文的结论不是“从 PIC 升级到 FLIP，再升级到 APIC”。更准确的理解是：PIC 通过覆盖进行筛选，FLIP 通过增量保留旧粒子信息，APIC 通过扩充局部表示改善网格往返的信息损失。知道每种方法保留了什么、丢弃了什么，以及网格究竟看到了什么，才能为具体物理问题作出选择。

参考资料

  1. Jiang, C., Schroeder, C., Selle, A., Teran, J., & Stomakhin, A. (2015). The Affine Particle-In-Cell Method. ACM Transactions on Graphics, 34(4), Article 51. DOI: 10.1145/2766996。本文传递公式对应原文第 5.3 节式（8）—（11），能量与局限讨论参照第 8、9 节。
  2. Disney Animation：作者发布页及原论文。本文的四粒子数值表和示意图为独立构造的验证，不是论文图表的复刻，也不是对原论文全部动力学实验的复现。
本文含有 AI 生成内容