FEM、MPM 与粒子仿真：从连续方程到一个可手算的时间步

1. 先把四层概念分开

因此“MPM 比 FEM 更软”不是正确陈述。两者若使用同一本构、分辨率和积分精度，应逼近同一个连续问题；实际差异来自离散误差、接触、数值耗散和稳定性。

2. 连续介质运动学：$\mathbf F$ 是整篇文章的核心状态

给材料中的每一点一个不随运动改变的参考坐标 $\mathbf X$。运动映射 $\boldsymbol\phi$ 把它送到当前世界坐标：

局部小线段满足 $d\mathbf x=\mathbf Fd\mathbf X$，其中

• $\mathbf F$ 同时包含旋转、拉伸和剪切；它不是应变，但本构通常以它为输入。
• $J$ 是当前体积与参考体积之比：$J=1$ 保体积，$J<1$ 压缩，$J>1$ 膨胀；$J\le 0$ 通常表示单元翻转或数值崩溃。
• 速度梯度 $\mathbf L=\nabla_{\mathbf x}\mathbf v$ 给出 $\dot{\mathbf F}=\mathbf L\mathbf F$，显式一步常写成 $\mathbf F^{n+1}=(\mathbf I+\Delta t\,\mathbf L^n)\mathbf F^n$。

3. 动量方程、边界条件和应力度量

当前构形中的局部动量平衡为

其中 $\boldsymbol\sigma$ 是 Cauchy 应力，$\mathbf b$ 是单位质量体力。位移边界给 $\mathbf u=\bar{\mathbf u}$，力边界给 $\boldsymbol\sigma\mathbf n=\bar{\mathbf t}$。在参考构形中常用第一 Piola–Kirchhoff 应力

这两个应力不能随意混用：$\mathbf P$ 与参考面积、$\nabla_{\mathbf X}$ 配对，$\boldsymbol\sigma$ 与当前面积、$\nabla_{\mathbf x}$ 配对。

4. 本构模型：同一个 $\mathbf F$ 为什么能像橡胶、雪或沙

4.1 从 $E,\nu$ 到 Lamé 参数

$E$ 控制拉伸刚度，$\nu$ 控制横向收缩；三维各向同性材料中 $\nu\to0.5$ 接近不可压缩，此时 $\lambda\to\infty$，显式步长更小，低阶位移 FEM 还可能出现 volumetric locking。

4.2 小形变线弹性

它只对小应变、小旋转可靠。一个刚体旋转也会被线性应变误判为形变，因此软体大转动通常用 corotated 或超弹性模型。

4.3 超弹性：能量先于力

超弹性从单位参考体积的能量 $\psi(\mathbf F)$ 出发，$\mathbf P=\partial\psi/\partial\mathbf F$。可压缩 Neo-Hookean 的常见形式是

fixed-corotated 先做极分解 $\mathbf F=\mathbf R\mathbf S$，用最近旋转 $\mathbf R$ 去掉刚体旋转：

4.4 塑性：把形变拆成可恢复与永久部分

弹性能只读 $\mathbf F_e$。一次典型 return mapping 是：先由速度梯度得到 trial $\mathbf F_e$，再做 SVD，把主伸长或试应力投影回 yield surface，超出的部分写进 $\mathbf F_p$。von Mises 限制偏应力，适合延性金属；Drucker–Prager 的屈服随压力变化，常用于沙土；雪模型会在压实时 hardening。这里的“粒子”价值很明显：$\mathbf F_p$、硬化量等历史变量一直跟着材料点走。

5. FEM：从强形式到可组装的矩阵

直接要求微分方程在每一点成立称为强形式。FEM 乘测试函数 $\mathbf w$ 并分部积分，把应力的空间导数转移到测试函数：

在单元内用形函数插值 $\mathbf u_h=\sum_aN_a\mathbf q_a$，得到半离散系统

线弹性时 $\boldsymbol\sigma=\mathbf D\mathbf B\mathbf q$，于是 $\mathbf K^e=\int\mathbf B^T\mathbf D\mathbf B,dV$。非线性 FEM 不再只有一个恒定 $\mathbf K$，而是在 Newton 迭代中重复计算 residual 和 tangent。

6. FEM 手算：一个二维三角形从位移到节点力

形函数 $N_1=1-x-y,N_2=x,N_3=y$，工程剪应变采用 $\gamma_{xy}$，所以

代入数值：

因为这个单元的 $\mathbf B$ 和应力为常数，积分直接变成乘 $tA=0.05$：

三个节点的 $x$ 力之和、$y$ 力之和都为 0，说明内力满足整体平移平衡。单元刚度为

未施加边界条件前 $\mathbf K^e$ 必须是奇异的，因为刚体平移和旋转不产生应变。若代码组装后它反而满秩，通常是 $\mathbf B$、自由度顺序或边界处理出了问题。

7. FEM 怎样推进时间

弹性波速量级 $c_p\approx\sqrt{(\lambda+2\mu)/\rho}$，显式稳定步长约

网格更细、材料更硬或密度更低都会迫使 $\Delta t$ 变小。substep 不是“白赚精度”，而是在支付稳定性成本。

8. MPM：粒子记历史，网格只活一个时间步

粒子 $p$ 通常保存 $m_p,\mathbf x_p,\mathbf v_p,V_p^0,\mathbf F_p$、仿射速度 $\mathbf C_p$ 和塑性变量。网格节点 $i$ 只保存本步质量、动量和力；下一步先清零，因此不会随材料变形而缠结。

8.1 P2G：把粒子质量和动量散射到网格

令 $w_{ip}=N_i(\mathbf x_p)$，$\nabla w_{ip}=\nabla N_i(\mathbf x_p)$：

$\mathbf C_p=0$ 就是 PIC 的常值粒子速度；APIC 用局部仿射场保存旋转和剪切信息。当前构形下，粒子应力产生的节点内力为

8.2 Grid solve：在规则网格上算力、碰撞和边界

随后把地面、刚体或域边界条件施加到 $\mathbf v_i^*$。这一步是 MPM 易于处理大形变、自接触和拓扑变化的关键：不同粒子可以共享网格动量，不需要维护永久扭曲的材料网格。

8.3 G2P：把更新后的网格速度带回粒子

PIC 覆盖粒子速度，会滤掉网格无法表示的模式，因此稳定但耗散；FLIP 只加网格速度增量，耗散小但会保留 grid 看不见的噪声。工程中常做 $\mathbf v_p=(1-\alpha)\mathbf v_p^{PIC}+\alpha\mathbf v_p^{FLIP}$。APIC 则更新

最后推进位置与形变：

之后做塑性投影并由新 $\mathbf F_e$ 计算下一步应力。不同实现可能交换 stress update 与 transfer 的具体时序，但必须保持使用的是同一时间层。

8.4 权重并非“离粒子最近就给它”

一维均匀网格上，令 base=floor(xp/h-0.5)、$f=x_p/h-\mathrm{base}$，二次 B-spline 的三个权重可写成

多维权重是各维权重的张量积。分片线性 hat function 更容易手算，但梯度在 cell boundary 不连续，粒子穿越网格时容易出现应力噪声；B-spline、GIMP、CPDI 和 MLS-MPM 都在改善这类 grid-crossing artifacts。

9. MPM 手算 A：两粒子看懂 PIC 为什么耗散

P2G 后：

没有力，所以 grid update 前后速度相同。PIC 插值回粒子：

粒子动能从 $1$ J 降到 $0.25$ J；丢失的反向高速模式在网格中无法表示。FLIP 使用网格增量，而本例增量为 0，所以仍得到 $[1,-1]$ m/s，动能不变；这也解释了它为什么更容易把未解析噪声一直留在粒子上。

10. MPM 手算 B：一次应力散射与 $F$ 更新

首先散射质量：$m_0=m_1=0.5$ kg。内力是

左右节点互相拉近，符号与受拉材料恢复原长的直觉一致；总内力为 0。于是 $a=[+40,-40]$ m/s²。取 $\Delta t=0.01$ s，网格速度更新为 $[+0.4,-0.4]$ m/s。

PIC 回传速度为 $v_p=0.5(0.4)+0.5(-0.4)=0$，所以质心不动；但局部速度梯度不为 0：

也就是说粒子位置可以暂时不动，但它代表的材料体积在收缩。这个小例子特别适合查四类 bug：内力负号、shape gradient 单位、lumped mass 除法以及 $\mathbf v\otimes\nabla w$ 的转置。

11. APIC 与 MLS-MPM 到底改了什么

APIC 不把粒子附近速度看成常数，而是 $\mathbf v(\mathbf x)=\mathbf v_p+\mathbf C_p(\mathbf x-\mathbf x_p)$。因此刚体旋转和线性剪切不必在每次 transfer 中消失，粒子—网格往返还能守恒线动量与角动量。MLS-MPM 从 moving least squares / Galerkin 视角重新推导 transfer，使 APIC/PolyPIC 与弱形式离散进入同一框架，并能把应力项直接并入 affine momentum transfer。

这不意味着“MLS-MPM 永远最好”。核函数、边界、接触、plastic projection、粒子密度和时间积分仍决定最终稳定性与精度。

12. FEM 与 MPM 的真正取舍

13. “之类的粒子法”分别在算什么

13.1 SPH：用核函数从邻居重建连续场

它是真正 mesh-free，擅长自由表面流体；代价是邻域搜索、压力噪声、边界缺邻居和不可压缩性。粒子代表的是场的采样，不必是一颗真实沙粒。

13.2 DEM：每颗粒子就是一个物理颗粒

DEM 逐接触解球、块体或颗粒的 Newton–Euler 方程，适合颗粒尺度问题；如果一亿粒沙必须一粒一粒表示，成本会远高于把沙当连续体的 MPM。

13.3 PBD / XPBD：直接修正位置满足约束

对约束 $C(\mathbf x)=0$，XPBD 的一次标量约束更新可写为

它以鲁棒、快速和可控外观见长；compliance 让刚度比传统 PBD 更少依赖步长和迭代数，但这仍不代表任意约束资产都等价于经过验证的连续体材料。

14. 选择方法时不要只看“软体还是颗粒”

15. 数值仿真该报告哪些指标

• 收敛：把 $h$、粒子间距与 $\Delta t$ 减半，位移、应力、接触力是否趋于稳定。
• 守恒：无外力时总线动量/角动量；保守系统中能量漂移与算法耗散。
• 可行性：最小 $J$、翻转单元数、负质量节点、粒子逃逸、接触穿透。
• 离散质量：FEM 看 element aspect ratio/Jacobian；MPM 看 particles per cell、活跃 cell 和每粒子 support。
• 稳定成本：实际 $\Delta t$、substeps、Newton 迭代数、linear solve residual，而不只是 frames/s。
• 物理校准：用拉伸、压缩、剪切、回弹、堆积角和接触实验分别辨识参数，不要靠一个 demo 同时猜 $E,\nu,\mu$ 和 damping。

16. 一份最小 MPM 伪代码

for each step:
    clear(grid.mass, grid.momentum, grid.force)

    for p in particles:                    # P2G
        stress[p] = constitutive(F_e[p], history[p])
        for i in support(x[p]):
            grid.mass[i]     += w(i,p) * mass[p]
            grid.momentum[i] += w(i,p) * mass[p] * (
                velocity[p] + C[p] @ (x_grid[i] - x[p]))
            grid.force[i]    -= volume[p] * stress[p] @ grad_w(i,p)

    for i in active_grid:
        v_old[i] = grid.momentum[i] / grid.mass[i]
        v_new[i] = v_old[i] + dt * (grid.force[i] / grid.mass[i] + gravity)
        v_new[i] = apply_collision_and_boundary(v_new[i])

    for p in particles:                    # G2P
        velocity[p] = interpolate_or_flip(v_new, v_old)
        C[p]        = reconstruct_affine_velocity(v_new)
        grad_v      = sum_i outer(v_new[i], grad_w(i,p))
        x[p]       += dt * velocity[p]
        F_trial     = (I + dt * grad_v) @ F[p]
        F[p], history[p] = plastic_return_map(F_trial, history[p])

真正实现时还要处理 active-grid 稀疏化、atomic scatter、粒子排序、边界核截断、接触多材料速度场、塑性 SVD、autodiff checkpoint 和确定性。公式正确只是第一步。