物质点法（Material Point Method, MPM）是一种混合欧拉-拉格朗日数值方法，在计算机图形学中被广泛用于模拟各种材料的大变形行为，包括雪、沙、泥浆、弹性体等。MPM结合了粒子方法的灵活性和网格方法的计算效率，特别适合处理拓扑变化、碎裂、混合等复杂现象。本章将深入探讨MPM的理论基础、移动最小二乘MPM（MLS-MPM）的改进，以及高级技术与应用。
在介绍MPM之前，我们需要深入理解连续介质力学的基本概念。考虑一个连续体 $\Omega$，其运动可以用映射 $\phi: \Omega_0 \times [0,T] \rightarrow \mathbb{R}^d$ 描述，其中 $\Omega_0$ 是初始（参考）构型，$d$ 是空间维度（通常为2或3）。
运动学描述：
变形梯度定义为： \(\mathbf{F} = \frac{\partial \phi}{\partial \mathbf{X}} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}}\)
变形梯度 $\mathbf{F}$ 是一个二阶张量，包含了所有的局部变形信息。从几何角度看，$\mathbf{F}$ 将参考构型中的无穷小线元 $d\mathbf{X}$ 映射到当前构型中的 $d\mathbf{x} = \mathbf{F} d\mathbf{X}$。
变形梯度的性质：
极分解定理： 任何可逆的变形梯度可以唯一分解为旋转和拉伸的组合： \(\mathbf{F} = \mathbf{R}\mathbf{U} = \mathbf{V}\mathbf{R}\) 其中：
计算极分解的算法通常使用SVD：$\mathbf{F} = \mathbf{U}{\text{svd}} \boldsymbol{\Sigma} \mathbf{V}{\text{svd}}^T$，则 $\mathbf{R} = \mathbf{U}{\text{svd}} \mathbf{V}{\text{svd}}^T$。
应变度量： 应变量化了变形中的拉伸部分，不同的应变张量适用于不同的变形程度和分析需求：
| 假设：$ | \mathbf{H} | \ll 1$（小变形、小转动） |
应变率和自旋张量： 速度梯度张量 $\mathbf{L} = \nabla\mathbf{v} = \dot{\mathbf{F}}\mathbf{F}^{-1}$ 可分解为： \(\mathbf{L} = \mathbf{D} + \mathbf{W}\) 其中：
应力度量： 应力描述内力的分布，不同的应力张量对应不同的参考构型：
守恒方程： 连续介质必须满足基本的物理守恒定律：
对于不可压缩材料：$\nabla \cdot \mathbf{v} = 0$
其中物质导数：$\frac{D(\cdot)}{Dt} = \frac{\partial(\cdot)}{\partial t} + \mathbf{v} \cdot \nabla(\cdot)$
| 总能量形式：$\rho\frac{D}{Dt}(e + \frac{1}{2} | \mathbf{v} | ^2) = \nabla \cdot (\boldsymbol{\sigma} \cdot \mathbf{v}) - \nabla \cdot \mathbf{q} + \rho(\mathbf{b} \cdot \mathbf{v} + r)$ |
边界条件和初始条件： 完整的边值问题需要指定：
虚功原理和弱形式： 强形式的平衡方程等价于虚功原理，对任意运动学可容许的虚位移 $\delta\mathbf{u}$（满足 $\delta\mathbf{u} = \mathbf{0}$ 在 $\partial\Omega_u$ 上）： \(\delta W = \delta W^{\text{int}} - \delta W^{\text{ext}} = 0\)
展开为： \(\int_{\Omega} \rho \ddot{\mathbf{u}} \cdot \delta\mathbf{u} \, dV + \int_{\Omega} \boldsymbol{\sigma} : \delta\boldsymbol{\epsilon} \, dV = \int_{\Omega} \rho\mathbf{b} \cdot \delta\mathbf{u} \, dV + \int_{\partial\Omega_t} \mathbf{t} \cdot \delta\mathbf{u} \, dA\)
这是有限元方法和MPM的理论基础，通过分部积分将二阶导数降为一阶，降低了对解的光滑性要求。
MPM巧妙地结合了拉格朗日粒子法和欧拉网格法的优势。物质点（粒子）携带所有的历史信息和材料属性，而背景网格仅用于求解动量方程，避免了网格畸变问题。
MPM的核心思想：
理论基础： MPM可以看作是有限元方法的粒子化版本。考虑弱形式： \(\int_\Omega \rho \frac{D\mathbf{v}}{Dt} \cdot \delta\mathbf{v} \, dV + \int_\Omega \boldsymbol{\sigma} : \nabla\delta\mathbf{v} \, dV = \int_\Omega \rho\mathbf{b} \cdot \delta\mathbf{v} \, dV\)
MPM使用粒子进行积分近似： \(\int_\Omega (\cdot) \, dV \approx \sum_p V_p (\cdot)_p\)
离散化策略： 将连续体 $\Omega$ 离散为 $N_p$ 个粒子，每个粒子 $p$ 代表一小块材料：
粒子初始化： 粒子的初始分布影响数值精度：
基本变量：
MPM的计算循环（单个时间步 $[t^n, t^{n+1}]$）：
粒子到网格的传输（P2G）：
质量传输： \(m_i = \sum_p w_{ip} m_p\)
这确保了总质量守恒：$\sum_i m_i = \sum_p m_p$
动量传输： \(\mathbf{p}_i = \sum_p w_{ip} m_p \mathbf{v}_p\)
其中插值权重：$w_{ip} = N_i(\mathbf{x}_p)$，$N_i$ 是与节点 $i$ 关联的形函数。
速度计算： \(\mathbf{v}_i = \begin{cases} \frac{\mathbf{p}_i}{m_i} & \text{if } m_i > \epsilon_m \\ \mathbf{0} & \text{otherwise} \end{cases}\)
其中 $\epsilon_m$ 是防止除零的小量（如 $10^{-10}$）。
力的计算：
内力（来自应力散度）： \(\mathbf{f}_i^{\text{int}} = -\sum_p V_p \boldsymbol{\sigma}_p \nabla w_{ip}\)
这是虚功原理的离散形式。对于任意虚速度场 $\delta\mathbf{v} = \sum_i \delta\mathbf{v}_i N_i$： \(\delta W^{\text{int}} = -\int_\Omega \boldsymbol{\sigma} : \nabla\delta\mathbf{v} \, dV \approx -\sum_p V_p \boldsymbol{\sigma}_p : \sum_i \delta\mathbf{v}_i \otimes \nabla w_{ip}\)
整理得到节点力： \(\mathbf{f}_i^{\text{int}} = -\frac{\partial W^{\text{int}}}{\partial \mathbf{v}_i}\)
外力：
阻尼力（可选）： \(\mathbf{f}_i^{\text{damp}} = -c_d m_i \mathbf{v}_i\) 其中 $c_d$ 是阻尼系数。
网格上的动量更新：
显式欧拉（一阶精度）： \(\mathbf{v}_i^{n+1} = \mathbf{v}_i^n + \frac{\Delta t}{m_i} (\mathbf{f}_i^{\text{int}} + \mathbf{f}_i^{\text{ext}})\)
辛欧拉（保持相空间体积）： \(\mathbf{p}_i^{n+1} = \mathbf{p}_i^n + \Delta t \mathbf{f}_i\) \(\mathbf{v}_i^{n+1} = \mathbf{p}_i^{n+1} / m_i\)
速度Verlet（二阶精度）： \(\mathbf{v}_i^{n+1/2} = \mathbf{v}_i^n + \frac{\Delta t}{2m_i} \mathbf{f}_i^n\) \(\mathbf{v}_i^{n+1} = \mathbf{v}_i^{n+1/2} + \frac{\Delta t}{2m_i} \mathbf{f}_i^{n+1}\)
边界条件处理：
Dirichlet边界（位移/速度约束）：
Neumann边界（力/应力约束）： 已包含在 $\mathbf{f}_i^{\text{traction}}$ 中
接触边界（单侧约束）：
网格到粒子的传输（G2P）：
PIC更新（Particle-In-Cell，数值耗散大）： \(\mathbf{v}_p^{n+1} = \sum_i w_{ip} \mathbf{v}_i^{n+1}\)
FLIP更新（Fluid-Implicit-Particle，保持细节但有噪声）： \(\mathbf{v}_p^{n+1} = \mathbf{v}_p^n + \sum_i w_{ip} (\mathbf{v}_i^{n+1} - \mathbf{v}_i^n)\)
PIC/FLIP混合（平衡耗散和噪声）： \(\mathbf{v}_p^{n+1} = (1-\alpha)\mathbf{v}_p^{\text{PIC}} + \alpha\mathbf{v}_p^{\text{FLIP}}\) 典型取 $\alpha \in [0.95, 0.99]$。
位置更新： \(\mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \Delta t \mathbf{v}_p^{n+1}\)
或使用中点法提高精度： \(\mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \Delta t \sum_i w_{ip} \frac{\mathbf{v}_i^{n+1} + \mathbf{v}_i^n}{2}\)
变形梯度更新：
速度梯度计算： \(\mathbf{L}_p = \sum_i \mathbf{v}_i^{n+1} \otimes \nabla w_{ip}\)
一阶更新（简单但可能引入误差）： \(\mathbf{F}_p^{n+1} = (\mathbf{I} + \Delta t \mathbf{L}_p) \mathbf{F}_p^n\)
指数映射（保持正定性和体积）： \(\mathbf{F}_p^{n+1} = \exp(\Delta t \mathbf{L}_p) \mathbf{F}_p^n\)
对于小时间步，使用Padé近似： \(\exp(\mathbf{A}) \approx (\mathbf{I} - \frac{1}{2}\mathbf{A})^{-1}(\mathbf{I} + \frac{1}{2}\mathbf{A})\)
体积更新： \(J_p^{n+1} = \det(\mathbf{F}_p^{n+1})\) \(V_p^{n+1} = J_p^{n+1} V_p^0\)
应力更新：
根据本构模型计算新的应力： \(\boldsymbol{\sigma}_p^{n+1} = \mathcal{C}(\mathbf{F}_p^{n+1}, \{\text{history}_p\})\)
其中 $\mathcal{C}$ 是本构关系，history包含：
返回映射算法（塑性）：
算法特点与优势：
数值考虑：
粒子与网格之间的信息传递是MPM的核心，形函数的选择直接影响数值精度、稳定性和计算效率。本节深入探讨形函数理论、实现细节和性能优化。
形函数的数学基础： 形函数 $N_i(\mathbf{x})$ 是定义在网格上的基函数，必须满足以下性质：
| 紧支性：$N_i(\mathbf{x}) = 0$ 当 $ | \mathbf{x} - \mathbf{x}_i | > r_{\text{support}}$ |
理论基础： 从重构理论角度，任意函数 $f(\mathbf{x})$ 可以近似为： \(f(\mathbf{x}) \approx \sum_i f_i N_i(\mathbf{x})\)
重构误差与形函数的逼近阶数相关： \(||f - f_h||_{L^2} = O(h^{k+1})\) 其中 $k$ 是多项式精度阶数。
常用形函数族：
线性形函数（帐篷函数/三线性插值）：
一维基函数： \(N^1(\xi) = \begin{cases} 1 - |\xi| & |\xi| < 1 \\ 0 & \text{otherwise} \end{cases}\)
导数： \(\frac{dN^1}{d\xi} = \begin{cases} -\text{sign}(\xi) & |\xi| < 1 \\ 0 & \text{otherwise} \end{cases}\)
多维张量积（分离变量）： \(N(\boldsymbol{\xi}) = \prod_{d=1}^D N^1(\xi_d)\)
梯度（链式法则）： \(\nabla N = \frac{1}{h}\begin{bmatrix} \frac{dN^1}{d\xi_1}(\xi_1) \prod_{d=2}^D N^1(\xi_d) \\ N^1(\xi_1) \frac{dN^1}{d\xi_2}(\xi_2) \prod_{d=3}^D N^1(\xi_d) \\ \vdots \end{bmatrix}\)
特点：
二次B样条：
一维基函数（中心化）： \(N^2(\xi) = \begin{cases} \frac{3}{4} - \xi^2 & |\xi| \leq \frac{1}{2} \\ \frac{1}{2}(\frac{3}{2} - |\xi|)^2 & \frac{1}{2} < |\xi| \leq \frac{3}{2} \\ 0 & |\xi| > \frac{3}{2} \end{cases}\)
导数： \(\frac{dN^2}{d\xi} = \begin{cases} -2\xi & |\xi| \leq \frac{1}{2} \\ -\text{sign}(\xi)(\frac{3}{2} - |\xi|) & \frac{1}{2} < |\xi| \leq \frac{3}{2} \\ 0 & |\xi| > \frac{3}{2} \end{cases}\)
特点：
三次B样条：
一维基函数： \(N^3(\xi) = \frac{1}{6}\begin{cases} (2-|\xi|)^3 & 1 < |\xi| \leq 2 \\ 4 - 6\xi^2 + 3|\xi|^3 & |\xi| \leq 1 \\ 0 & |\xi| > 2 \end{cases}\)
导数： \(\frac{dN^3}{d\xi} = \frac{1}{6}\begin{cases} -3\text{sign}(\xi)(2-|\xi|)^2 & 1 < |\xi| \leq 2 \\ -12\xi + 9\xi|\xi| & |\xi| \leq 1 \\ 0 & |\xi| > 2 \end{cases}\)
特点：
GIMP（Generalized Interpolation Material Point）形函数：
核心思想：考虑粒子的有限尺寸，通过卷积获得形函数： \(S_{ip} = \int_{\Omega_p} N_i(\mathbf{x}) \chi_p(\mathbf{x}) d\mathbf{x}\)
其中 $\chi_p$ 是粒子的特征函数（形状函数）。
对于轴对齐矩形粒子，一维GIMP权重： \(S^1_{ip}(\xi, l) = \begin{cases} 1 - \frac{(|\xi|-l)^2}{4l} & l < |\xi| < 1+l \\ 1 - |\xi| & |\xi| < l \\ 0 & |\xi| > 1+l \end{cases}\)
其中 $l = l_p/h$ 是归一化粒子半宽度。
梯度（使用Leibniz积分法则）： \(\nabla S_{ip} = \int_{\Omega_p} \nabla N_i(\mathbf{x}) \chi_p(\mathbf{x}) d\mathbf{x}\)
优势：
变种：
权重和梯度的高效计算：
给定粒子位置 $\mathbf{x}_p$ 和网格节点 $i$ 的位置 $\mathbf{x}_i$：
基节点索引: base = floor(x_p / h)
局部坐标: ξ = (x_p - x_i) / h
相对索引: offset = i - base
梯度计算技巧： \(\nabla w_{ip} = \frac{1}{h} \begin{bmatrix} \frac{\partial N}{\partial \xi_1} \frac{N(\xi_2)N(\xi_3)}{N(\boldsymbol{\xi})} \\ \frac{\partial N}{\partial \xi_2} \frac{N(\xi_1)N(\xi_3)}{N(\boldsymbol{\xi})} \\ \frac{\partial N}{\partial \xi_3} \frac{N(\xi_1)N(\xi_2)}{N(\boldsymbol{\xi})} \end{bmatrix} w_{ip}\)
避免重复计算基函数值。
误差分析与收敛性：
插值误差估计： 对于 $k$ 阶精度的形函数： \(||f - \Pi_h f||_{L^2(\Omega)} \leq C h^{k+1} ||f||_{H^{k+1}(\Omega)}\)
其中 $\Pi_h$ 是插值算子。
积分精度： 数值积分误差： \(\left|\int_\Omega f dV - \sum_p V_p f_p\right| \leq C h^{k+1} ||f||_{W^{k+1,1}(\Omega)}\)
条件数分析： 质量矩阵条件数：$\kappa(M) = O(1)$（对角占优） 刚度矩阵条件数：$\kappa(K) = O(h^{-2})$（与FEM相同）
数值稳定性考虑：
梯度爆炸预防： 限制速度梯度： \(||\mathbf{L}_p|| \leq L_{\max} = \frac{c_{\max}}{h}\)
其中 $c_{\max}$ 是材料波速上界。
高效实现策略：
struct GridNode {
float mass;
vec3 momentum;
vec3 force;
// 填充到缓存行边界
};
struct Particle {
vec3 position;
vec3 velocity;
mat3 F; // 变形梯度
mat3 stress;
float volume;
// 按访问模式组织
};
// 避免重复计算
float wx = N1(xi.x);
float wy = N1(xi.y);
float wz = N1(xi.z);
float wip = wx * wy * wz;
// 利用对称性
grad_wip.x = dN1(xi.x) * wy * wz / h;
自适应形函数选择：
根据局部物理状态动态选择形函数：
应变率准则： \(\text{order} = \begin{cases} 1 & ||\mathbf{D}|| > D_{\text{thresh}} \\ 2 & \text{otherwise} \end{cases}\)
形函数选择决策树：
材料类型？
├─ 弹性固体
│ └─ 二次B样条（C¹连续性）
├─ 流体
│ ├─ 低雷诺数 → 三次B样条
│ └─ 高雷诺数 → 二次B样条
├─ 颗粒材料
│ └─ 线性 + GIMP（防穿透）
└─ 大变形塑性
└─ 自适应（应变率相关）
MPM中常用的时间积分方案包括：
CFL条件： \(\Delta t \leq C \frac{h}{c}\) 其中 $c = \sqrt{E/\rho}$ 是材料中的声速，$C$ 是CFL数（通常取0.1-0.5）。
MPM可以处理各种材料模型：
线弹性： \(\boldsymbol{\sigma} = \lambda \text{tr}(\boldsymbol{\epsilon})\mathbf{I} + 2\mu\boldsymbol{\epsilon}\)
超弹性（Neo-Hookean）： \(\boldsymbol{\sigma} = \mu(\mathbf{B} - \mathbf{I}) + \lambda\ln(J)\mathbf{I}\) 其中 $\mathbf{B} = \mathbf{F}\mathbf{F}^T$，$J = \det(\mathbf{F})$。
塑性模型：
应力计算通常涉及：
传统MPM虽然在处理大变形问题上有优势，但存在一些固有问题：
MLS-MPM的核心思想是使用移动最小二乘（Moving Least Squares）方法在粒子周围重构速度场，从而获得更精确的变形梯度更新。
速度场重构： 假设在粒子 $p$ 附近的速度场可以表示为： \(\mathbf{v}(\mathbf{x}) = \mathbf{c}_0 + \mathbf{C}_1 (\mathbf{x} - \mathbf{x}_p)\)
其中 $\mathbf{c}_0$ 是常数项，$\mathbf{C}_1$ 是一阶导数矩阵。
最小二乘拟合： 通过最小化加权误差： \(E = \sum_i w_{ip} m_i ||\mathbf{v}_i - \mathbf{c}_0 - \mathbf{C}_1(\mathbf{x}_i - \mathbf{x}_p)||^2\)
求解得到： \(\mathbf{c}_0 = \sum_i w_{ip} \mathbf{v}_i\) \(\mathbf{C}_1 = \left(\sum_i w_{ip} \mathbf{v}_i \otimes (\mathbf{x}_i - \mathbf{x}_p)\right) \mathbf{D}_p^{-1}\)
其中 $\mathbf{D}p = \sum_i w{ip} (\mathbf{x}_i - \mathbf{x}_p) \otimes (\mathbf{x}_i - \mathbf{x}_p)$ 是惯性张量。
APIC矩阵： MLS-MPM引入仿射粒子-网格（APIC）矩阵 $\mathbf{C}p$ 来存储速度梯度信息： \(\mathbf{C}_p = \mathbf{B}_p \mathbf{D}_p^{-1}\) 其中 $\mathbf{B}_p = \sum_i w{ip} \mathbf{v}_i \otimes (\mathbf{x}_i - \mathbf{x}_p)$。
P2G传输（改进版）： \(m_i = \sum_p w_{ip} m_p\) \(m_i \mathbf{v}_i = \sum_p w_{ip} m_p [\mathbf{v}_p + \mathbf{C}_p(\mathbf{x}_i - \mathbf{x}_p)]\)
网格动量更新（与传统MPM相同）： \(\mathbf{v}_i^{n+1} = \mathbf{v}_i^n + \frac{\Delta t}{m_i} \mathbf{f}_i\)
G2P传输（改进版）： \(\mathbf{v}_p^{n+1} = \sum_i w_{ip} \mathbf{v}_i^{n+1}\) \(\mathbf{C}_p^{n+1} = \sum_i w_{ip} \mathbf{v}_i^{n+1} \otimes (\mathbf{x}_i - \mathbf{x}_p) \cdot \mathbf{D}_p^{-1}\)
惯性张量的高效计算： 对于规则网格和对称核函数： \(\mathbf{D}_p = \frac{4h^2}{3}\mathbf{I}\)（二维情况） \(\mathbf{D}_p = \frac{h^2}{2}\mathbf{I}\)（三维情况）
体积守恒： 通过修正变形梯度保持体积： \(\mathbf{F}_p^{\text{corrected}} = J^{1/d} \mathbf{F}_p / \det(\mathbf{F}_p)^{1/d}\)
数值稳定性增强：
MLS-MPM可以看作是APIC（Affine Particle-In-Cell）方法的推广：
优势对比：
实际应用中经常需要模拟多种材料的相互作用，如流体-固体耦合、多相流等。
多材料表示方法：
界面处理：
接触算法： \(\mathbf{v}_i^{k,\text{final}} = \mathbf{v}_i^k + \frac{\Delta t}{m_i^k} \mathbf{f}_i^{\text{contact}}\)
接触力通过罚函数或拉格朗日乘子计算。
混合理论： 应力通过体积加权平均： \(\boldsymbol{\sigma} = \sum_k \phi_k \boldsymbol{\sigma}_k\)
表面张力： 对于流体界面，需要考虑表面张力： \(\mathbf{f}_{\text{surface}} = \gamma \kappa \mathbf{n}\) 其中 $\gamma$ 是表面张力系数，$\kappa$ 是曲率，$\mathbf{n}$ 是法向量。
为了提高计算效率和精度，可以使用自适应技术：
分裂准则： \(\text{split if } \det(\mathbf{F}_p) > \theta_{\text{split}}\)
合并准则： \(\text{merge if } ||\mathbf{x}_p - \mathbf{x}_q|| < r_{\text{merge}} \text{ and } \det(\mathbf{F}_p) < \theta_{\text{merge}}\)
时间自适应： 根据局部CFL条件调整时间步长： \(\Delta t_p = \min\left(C_{\text{CFL}} \frac{h}{||\mathbf{v}_p|| + c_p}, \Delta t_{\text{max}}\right)\)
MPM天然适合GPU并行化，关键优化策略包括：
典型加速比：
本章深入介绍了物质点法（MPM）及其现代变种MLS-MPM，这是计算机图形学中模拟大变形材料的强大工具。
关键概念：
重要公式：
算法选择指南：
形函数性质 证明二次B样条形函数满足分割单位性：$\sum_i N_i^2(\mathbf{x}) = 1$。
提示：考虑一维情况，利用B样条的递归定义。
变形梯度的物理意义 给定变形梯度 $\mathbf{F} = \begin{bmatrix} 2 & 0.5 \ 0 & 1 \end{bmatrix}$，计算： a) 体积比 $J$ b) 右Cauchy-Green张量 $\mathbf{C}$ c) Green-Lagrange应变 $\mathbf{E}$
提示：使用 $J = \det(\mathbf{F})$，$\mathbf{C} = \mathbf{F}^T\mathbf{F}$。
CFL稳定性分析 对于杨氏模量 $E = 10^6$ Pa、密度 $\rho = 1000$ kg/m³ 的材料，网格间距 $h = 0.01$ m，计算最大允许时间步长（取CFL数 $C = 0.3$）。
提示：声速 $c = \sqrt{E/\rho}$。
能量守恒分析 证明在无外力和无耗散的情况下，显式MPM的总能量（动能+势能）在一个时间步内的变化为 $O(\Delta t^2)$。
提示：考虑辛结构和中点法则。
MLS-MPM的理论分析 推导MLS-MPM中APIC矩阵 $\mathbf{C}_p$ 的更新公式，并说明其如何保持角动量守恒。
提示：从最小二乘问题出发，考虑惯性张量的作用。
多材料耦合 设计一个算法处理流固耦合问题，其中固体使用Neo-Hookean模型，流体使用不可压缩条件。讨论界面条件的处理。
提示：考虑分离式求解和耦合约束。