第12章：物质点法

物质点法（Material Point Method, MPM）是一种混合欧拉-拉格朗日数值方法，在计算机图形学中被广泛用于模拟各种材料的大变形行为，包括雪、沙、泥浆、弹性体等。MPM结合了粒子方法的灵活性和网格方法的计算效率，特别适合处理拓扑变化、碎裂、混合等复杂现象。本章将深入探讨MPM的理论基础、移动最小二乘MPM（MLS-MPM）的改进，以及高级技术与应用。

章节大纲

12.1 物质点法基础理论

• 连续介质力学基础
• MPM的基本框架
• 粒子-网格交互
• 时间积分方案

12.2 MLS-MPM（移动最小二乘物质点法）

• 传统MPM的局限性
• MLS-MPM的理论推导
• APIC与MLS-MPM的关系
• 实现细节与优化

12.3 高级MPM技术与应用

• 多材料MPM
• 自适应MPM
• GPU加速技术
• 工程应用案例

12.1 物质点法基础理论

12.1.1 连续介质力学基础

在介绍MPM之前，我们需要深入理解连续介质力学的基本概念。考虑一个连续体 $\Omega$，其运动可以用映射 $\phi: \Omega_0 \times [0,T] \rightarrow \mathbb{R}^d$ 描述，其中 $\Omega_0$ 是初始（参考）构型，$d$ 是空间维度（通常为2或3）。

运动学描述：

• 拉格朗日描述：跟踪材料点 $\mathbf{X} \in \Omega_0$ 的运动轨迹，适合追踪材料历史和内部变量演化
• 欧拉描述：关注空间点 $\mathbf{x} \in \Omega_t$ 处的物理量变化，适合流体力学和大变形问题
• 映射关系：$\mathbf{x} = \phi(\mathbf{X}, t)$
• 速度场：$\mathbf{v}(\mathbf{X}, t) = \frac{\partial \phi}{\partial t}(\mathbf{X}, t)$（拉格朗日）或 $\mathbf{v}(\mathbf{x}, t)$（欧拉）

变形梯度定义为： $$\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}$。

变形梯度的性质：

• 可逆性：物理上合理的变形要求 $\det(\mathbf{F}) > 0$（避免材料穿透自身）
• 体积比：$J = \det(\mathbf{F})$ 表示局部体积变化率
• $J > 0$：保持定向（物理上合理）
• $J = 1$：体积守恒（不可压缩）
• $J < 1$：压缩
• $J > 1$：膨胀
• 恒等变形：未变形状态对应 $\mathbf{F} = \mathbf{I}$

极分解定理： 任何可逆的变形梯度可以唯一分解为旋转和拉伸的组合： $$\mathbf{F} = \mathbf{R}\mathbf{U} = \mathbf{V}\mathbf{R}$$ 其中：

• $\mathbf{R}$ 是旋转张量（正交张量）：$\mathbf{R}^T\mathbf{R} = \mathbf{I}$，$\det(\mathbf{R}) = 1$
• $\mathbf{U}$ 是右拉伸张量（对称正定）：描述旋转前的拉伸
• $\mathbf{V}$ 是左拉伸张量（对称正定）：描述旋转后的拉伸

计算极分解的算法通常使用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$。

应变度量： 应变量化了变形中的拉伸部分，不同的应变张量适用于不同的变形程度和分析需求：

1. Green-Lagrange应变张量（拉格朗日描述，大变形）： $$\mathbf{E} = \frac{1}{2}(\mathbf{F}^T\mathbf{F} - \mathbf{I}) = \frac{1}{2}(\mathbf{C} - \mathbf{I})$$ 其中 $\mathbf{C} = \mathbf{F}^T\mathbf{F}$ 是右Cauchy-Green变形张量。

• 物理意义：度量参考构型中线元长度的变化
• 特点：在刚体旋转下不变，适合大变形分析

2. Almansi-Euler应变张量（欧拉描述）： $$\mathbf{e} = \frac{1}{2}(\mathbf{I} - \mathbf{F}^{-T}\mathbf{F}^{-1}) = \frac{1}{2}(\mathbf{I} - \mathbf{B}^{-1})$$ 其中 $\mathbf{B} = \mathbf{F}\mathbf{F}^T$ 是左Cauchy-Green变形张量。

• 物理意义：度量当前构型中线元长度的变化
• 应用：适合更新拉格朗日格式

3. 小应变张量（线性化，工程应变）： 位移场 $\mathbf{u} = \mathbf{x} - \mathbf{X}$，位移梯度 $\mathbf{H} = \nabla\mathbf{u} = \mathbf{F} - \mathbf{I}$ $$\boldsymbol{\epsilon} = \frac{1}{2}(\mathbf{H} + \mathbf{H}^T) = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$$

• 假设：$||\mathbf{H}|| \ll 1$（小变形、小转动）
• 优点：线性理论，计算简单

4. 对数应变（Hencky应变）： $$\mathbf{E}_{\ln} = \ln(\mathbf{V}) = \frac{1}{2}\ln(\mathbf{B})$$

• 特点：对于单轴拉伸给出真实应变
• 计算：需要矩阵对数，通常通过特征值分解

应变率和自旋张量： 速度梯度张量 $\mathbf{L} = \nabla\mathbf{v} = \dot{\mathbf{F}}\mathbf{F}^{-1}$ 可分解为： $$\mathbf{L} = \mathbf{D} + \mathbf{W}$$ 其中：

• 变形率张量（对称）：$\mathbf{D} = \frac{1}{2}(\mathbf{L} + \mathbf{L}^T)$
• 自旋张量（反对称）：$\mathbf{W} = \frac{1}{2}(\mathbf{L} - \mathbf{L}^T)$

应力度量： 应力描述内力的分布，不同的应力张量对应不同的参考构型：

1. Cauchy应力张量 $\boldsymbol{\sigma}$（真实应力）： - 定义：当前构型中单位面积上的力 - 对称性：角动量守恒要求 $\boldsymbol{\sigma} = \boldsymbol{\sigma}^T$ - 牵引力：$\mathbf{t} = \boldsymbol{\sigma} \mathbf{n}$（$\mathbf{n}$ 是当前构型的外法向）

2. 第一Piola-Kirchhoff应力张量 $\mathbf{P}$（名义应力）： $$\mathbf{P} = J\boldsymbol{\sigma}\mathbf{F}^{-T}$$

• 定义：当前力对参考面积的比值
• 非对称：$\mathbf{P} \neq \mathbf{P}^T$（一般情况）
• 功共轭：$\mathbf{P} : \dot{\mathbf{F}} = J\boldsymbol{\sigma} : \mathbf{D}$

3. 第二Piola-Kirchhoff应力张量 $\mathbf{S}$： $$\mathbf{S} = \mathbf{F}^{-1}\mathbf{P} = J\mathbf{F}^{-1}\boldsymbol{\sigma}\mathbf{F}^{-T}$$

• 定义：完全参考构型的应力度量
• 对称性：$\mathbf{S} = \mathbf{S}^T$
• 功共轭：$\mathbf{S} : \dot{\mathbf{E}} = J\boldsymbol{\sigma} : \mathbf{D}$

4. Kirchhoff应力张量 $\boldsymbol{\tau}$： $$\boldsymbol{\tau} = J\boldsymbol{\sigma}$$

• 应用：不可压缩材料（$J = 1$）时等于Cauchy应力
• 优点：在大变形问题中数值稳定性更好

守恒方程： 连续介质必须满足基本的物理守恒定律：

1. 质量守恒（连续性方程）： - 积分形式：$\frac{d}{dt}\int_{\Omega_t} \rho \, dv = 0$ - 拉格朗日形式：$\rho_0 = J\rho$（质量密度关系） - 欧拉形式：$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\mathbf{v}) = 0$ - 物质导数形式：$\frac{D\rho}{Dt} + \rho\nabla \cdot \mathbf{v} = 0$

对于不可压缩材料：$\nabla \cdot \mathbf{v} = 0$

2. 动量守恒（运动方程）： - 积分形式：$\frac{d}{dt}\int_{\Omega_t} \rho\mathbf{v} \, dv = \int_{\Omega_t} \rho\mathbf{b} \, dv + \int_{\partial\Omega_t} \mathbf{t} \, da$ - 拉格朗日形式：$\rho_0 \frac{\partial^2 \phi}{\partial t^2} = \nabla_0 \cdot \mathbf{P} + \rho_0\mathbf{b}$ - 欧拉形式（Cauchy运动方程）：$\rho\frac{D\mathbf{v}}{Dt} = \nabla \cdot \boldsymbol{\sigma} + \rho\mathbf{b}$

其中物质导数：$\frac{D(\cdot)}{Dt} = \frac{\partial(\cdot)}{\partial t} + \mathbf{v} \cdot \nabla(\cdot)$

3. 角动量守恒： - 积分形式：$\frac{d}{dt}\int_{\Omega_t} \mathbf{x} \times \rho\mathbf{v} \, dv = \int_{\Omega_t} \mathbf{x} \times \rho\mathbf{b} \, dv + \int_{\partial\Omega_t} \mathbf{x} \times \mathbf{t} \, da$ - 局部结果：Cauchy应力张量的对称性 $\boldsymbol{\sigma} = \boldsymbol{\sigma}^T$

4. 能量守恒（热力学第一定律）： $$\rho\frac{De}{Dt} = \boldsymbol{\sigma} : \mathbf{D} - \nabla \cdot \mathbf{q} + \rho r$$ 其中：

• $e$：比内能
• $\boldsymbol{\sigma} : \mathbf{D}$：应力功率（机械功）
• $\mathbf{q}$：热流向量（Fourier定律：$\mathbf{q} = -k\nabla T$）
• $r$：体积热源

总能量形式：$\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)$

5. 熵不等式（热力学第二定律）： $$\rho\frac{Ds}{Dt} + \nabla \cdot \left(\frac{\mathbf{q}}{T}\right) - \frac{\rho r}{T} \geq 0$$ 其中 $s$ 是比熵，$T$ 是绝对温度。

边界条件和初始条件： 完整的边值问题需要指定：

1. 边界条件： - Dirichlet（本质）：$\mathbf{u} = \overline{\mathbf{u}}$ 在 $\partial\Omega_u$ 上 - Neumann（自然）：$\mathbf{t} = \boldsymbol{\sigma} \cdot \mathbf{n} = \overline{\mathbf{t}}$ 在 $\partial\Omega_t$ 上 - Robin（混合）：$\alpha\mathbf{u} + \beta\mathbf{t} = \mathbf{g}$

2. 初始条件： - 位置：$\mathbf{x}(\mathbf{X}, 0) = \mathbf{x}_0(\mathbf{X})$ - 速度：$\mathbf{v}(\mathbf{X}, 0) = \mathbf{v}_0(\mathbf{X})$

虚功原理和弱形式： 强形式的平衡方程等价于虚功原理，对任意运动学可容许的虚位移 $\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的理论基础，通过分部积分将二阶导数降为一阶，降低了对解的光滑性要求。

12.1.2 MPM的基本框架

MPM巧妙地结合了拉格朗日粒子法和欧拉网格法的优势。物质点（粒子）携带所有的历史信息和材料属性，而背景网格仅用于求解动量方程，避免了网格畸变问题。

MPM的核心思想：

1. 双重表示：材料既由粒子表示（携带历史），又通过网格求解（高效计算）
2. 信息流动：粒子→网格（收集信息）→求解→网格→粒子（更新状态）
3. 自动处理拓扑变化：无需显式追踪界面或重新网格化
4. 混合方法优势：结合了拉格朗日方法的材料追踪能力和欧拉方法的计算效率

理论基础： 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$ 代表一小块材料：

• 初始体积：$V_p^0 = \int_{\Omega_p^0} dV$
• 质量守恒：$m_p = \rho_p^0 V_p^0 = \text{const}$
• 当前体积：$V_p = J_p V_p^0$
• 密度更新：$\rho_p = m_p / V_p = \rho_p^0 / J_p$

粒子初始化： 粒子的初始分布影响数值精度：

1. 规则分布：在每个网格单元内均匀放置 $2^d$ 个粒子（$d$ 是维度）
2. 随机扰动：添加小的随机偏移避免人工各向异性
3. 自适应分布：根据应力梯度或变形程度调整粒子密度

基本变量：

1. 粒子变量（拉格朗日）： - 位置：$\mathbf{x}_p \in \mathbb{R}^d$ - 速度：$\mathbf{v}_p \in \mathbb{R}^d$ - 质量：$m_p \in \mathbb{R}^+$（常数） - 初始体积：$V_p^0 \in \mathbb{R}^+$（常数） - 当前体积：$V_p = J_p V_p^0$ - 变形梯度：$\mathbf{F}_p \in \mathbb{R}^{d \times d}$，$\det(\mathbf{F}_p) > 0$ - 应力张量：$\boldsymbol{\sigma}_p \in \mathbb{R}^{d \times d}$ - 内部变量：

   • 塑性应变：$\boldsymbol{\epsilon}_p^p$
   • 等效塑性应变：$\bar{\epsilon}_p^p$
   • 损伤变量：$d_p \in [0,1]$
   • 温度：$T_p$
   • 相场变量：$\phi_p$（多相材料）

2. 网格变量（欧拉）： - 节点位置：$\mathbf{x}_i \in \mathbb{R}^d$（固定） - 节点速度：$\mathbf{v}_i \in \mathbb{R}^d$ - 节点质量：$m_i \in \mathbb{R}^+$ - 节点动量：$\mathbf{p}_i = m_i \mathbf{v}_i$ - 节点内力：$\mathbf{f}_i^{\text{int}} \in \mathbb{R}^d$ - 节点外力：$\mathbf{f}_i^{\text{ext}} \in \mathbb{R}^d$

MPM的计算循环（单个时间步 $[t^n, t^{n+1}]$）：

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}$）。

2. 力的计算：

内力（来自应力散度）： $$\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{body}} = \sum_p w_{ip} m_p \mathbf{b}_p$
• 表面力贡献：$\mathbf{f}_i^{\text{traction}} = \int_{\partial\Omega_t} N_i \mathbf{t} \, dA$
• 总外力：$\mathbf{f}_i^{\text{ext}} = \mathbf{f}_i^{\text{body}} + \mathbf{f}_i^{\text{traction}}$

阻尼力（可选）： $$\mathbf{f}_i^{\text{damp}} = -c_d m_i \mathbf{v}_i$$ 其中 $c_d$ 是阻尼系数。

3. 网格上的动量更新：

显式欧拉（一阶精度）： $$\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}$$

4. 边界条件处理：

Dirichlet边界（位移/速度约束）：

• 固定边界：$\mathbf{v}_i = \mathbf{0}$
• 移动边界：$\mathbf{v}_i = \mathbf{v}_{\text{prescribed}}(t)$
• 滑移边界：$\mathbf{v}_i \cdot \mathbf{n} = 0$，$\mathbf{v}_i = \mathbf{v}_i - (\mathbf{v}_i \cdot \mathbf{n})\mathbf{n}$

Neumann边界（力/应力约束）： 已包含在 $\mathbf{f}_i^{\text{traction}}$ 中

接触边界（单侧约束）：

• 法向：$\mathbf{v}_n = \max(0, \mathbf{v}_i \cdot \mathbf{n})$
• 切向（库仑摩擦）： $$\mathbf{v}_t = \begin{cases} \mathbf{0} & \text{if } ||\mathbf{v}_t|| \leq \mu |\mathbf{v}_n| \ (1-\mu|\mathbf{v}_n|/||\mathbf{v}_t||)\mathbf{v}_t & \text{otherwise} \end{cases}$$ 其中 $\mu$ 是摩擦系数。

5. 网格到粒子的传输（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}$$

6. 变形梯度更新：

速度梯度计算： $$\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$$

7. 应力更新：

根据本构模型计算新的应力： $$\boldsymbol{\sigma}_p^{n+1} = \mathcal{C}(\mathbf{F}_p^{n+1}, {\text{history}_p})$$ 其中 $\mathcal{C}$ 是本构关系，history包含：

• 塑性应变历史
• 损伤演化
• 温度历史
• 相变状态

返回映射算法（塑性）：

1. 弹性预测：$\boldsymbol{\sigma}^{\text{trial}} = \mathcal{C}^{\text{elastic}}(\mathbf{F}^{n+1})$
2. 检查屈服：$f(\boldsymbol{\sigma}^{\text{trial}}) \leq 0$？
3. 塑性修正：如果 $f > 0$，投影回屈服面

算法特点与优势：

1. 自动质量守恒：粒子质量不变，总质量精确守恒
2. 动量守恒：通过网格求解，动量守恒到机器精度
3. 历史依赖：粒子携带所有历史变量，自然处理路径依赖材料
4. 大变形能力：无网格畸变，可处理任意大变形和拓扑变化
5. 多物理耦合：易于添加热传导、相变、化学反应等
6. 并行友好：P2G和G2P阶段高度并行

数值考虑：

1. 网格穿越误差：粒子穿越网格单元时可能引入误差
2. 数值断裂：粒子分离可能导致非物理断裂
3. 粒子聚集：某些区域粒子过度聚集影响精度
4. 边界处理：需要特殊技术处理复杂边界

12.1.3 粒子-网格交互

粒子与网格之间的信息传递是MPM的核心，形函数的选择直接影响数值精度、稳定性和计算效率。本节深入探讨形函数理论、实现细节和性能优化。

形函数的数学基础： 形函数 $N_i(\mathbf{x})$ 是定义在网格上的基函数，必须满足以下性质：

1. 插值性（分割单位性）：$\sum_i N_i(\mathbf{x}) = 1$ 对所有 $\mathbf{x} \in \Omega$
2. 紧支性：$N_i(\mathbf{x}) = 0$ 当 $||\mathbf{x} - \mathbf{x}_i|| > r_{\text{support}}$
3. 非负性：$N_i(\mathbf{x}) \geq 0$（保证物理意义）
4. 光滑性：至少 $C^0$ 连续，理想情况下 $C^{k-1}$（$k$ 阶精度需要）
5. 局部性：每个点最多被有限个形函数覆盖

理论基础： 从重构理论角度，任意函数 $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$ 是多项式精度阶数。

常用形函数族：

1. 线性形函数（帐篷函数/三线性插值）：

一维基函数： $$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}$$ 特点：

• 支撑域：$2^D$ 个节点（2×2×2 in 3D）
• 连续性：$C^0$（位置连续，速度不连续）
• 计算效率：最高（简单的加减乘除）
• 数值问题：网格穿越噪声、力的不连续性

2. 二次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}$$ 特点：

• 支撑域：$3^D$ 个节点（3×3×3 in 3D）
• 连续性：$C^1$（速度连续）
• 精度：$O(h^3)$ 收敛率
• 应用：流体模拟首选（平滑性好）

3. 三次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}$$ 特点：

• 支撑域：$4^D$ 个节点（4×4×4 in 3D）
• 连续性：$C^2$（加速度连续）
• 精度：$O(h^4)$ 收敛率
• 代价：计算量大，内存带宽需求高

4. 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}$$ 优势：

• 减少粒子穿透
• 更好的动量传递
• 自适应分辨率

变种：

• cpGIMP：收缩粒子域GIMP
• uGIMP：均匀GIMP（简化计算）

权重和梯度的高效计算：

给定粒子位置 $\mathbf{x}_p$ 和网格节点 $i$ 的位置 $\mathbf{x}_i$：

1. 坐标变换： 基节点索引: base = floor(x_p / h) 局部坐标: ξ = (x_p - x_i) / h 相对索引: offset = i - base

2. 权重计算优化： - 查表法：预计算 $N(\xi)$ 在细分网格上的值 - 分段多项式：直接计算，利用分支预测 - SIMD向量化：同时计算多个维度

3. 梯度计算技巧： $$\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}$$ 避免重复计算基函数值。

误差分析与收敛性：

1. 插值误差估计： 对于 $k$ 阶精度的形函数： $$||f - \Pi_h f||_{L^2(\Omega)} \leq C h^{k+1} ||f||_{H^{k+1}(\Omega)}$$ 其中 $\Pi_h$ 是插值算子。

2. 积分精度： 数值积分误差： $$\left|\int_\Omega f dV - \sum_p V_p f_p\right| \leq C h^{k+1} ||f||_{W^{k+1,1}(\Omega)}$$

3. 条件数分析： 质量矩阵条件数：$\kappa(M) = O(1)$（对角占优） 刚度矩阵条件数：$\kappa(K) = O(h^{-2})$（与FEM相同）

数值稳定性考虑：

1. 质量聚集问题： 当 $m_i = \sum_p w_{ip} m_p < \epsilon_m$ 时：

• 跳过该节点的力计算
• 或使用正则化：$\tilde{m}_i = \max(m_i, \epsilon_m)$

2. 梯度爆炸预防： 限制速度梯度： $$||\mathbf{L}_p|| \leq L_{\max} = \frac{c_{\max}}{h}$$ 其中 $c_{\max}$ 是材料波速上界。

3. 边界附近的处理： - 使用修正的形函数保证分割单位性 - 或采用虚拟节点技术

高效实现策略：

1. 数据结构优化： ``` struct GridNode { float mass; vec3 momentum; vec3 force; // 填充到缓存行边界 };

struct Particle { vec3 position; vec3 velocity; mat3 F; // 变形梯度 mat3 stress; float volume; // 按访问模式组织 }; ```

2. 并行化策略： - P2G并行：使用原子操作或图着色避免竞争 - 力计算：完全并行（只读粒子数据） - G2P并行：每个粒子独立更新

3. 内存访问优化： - 粒子排序：按空间位置排序（Z-order、Hilbert曲线） - 缓存预取：显式预取下一个粒子数据 - AoS vs SoA：根据访问模式选择

4. 计算优化： ``` // 避免重复计算 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; ```

自适应形函数选择：

根据局部物理状态动态选择形函数：

1. 应变率准则： $$\text{order} = \begin{cases} 1 & ||\mathbf{D}|| > D_{\text{thresh}} \ 2 & \text{otherwise} \end{cases}$$

2. 材料类型准则： - 固体：二次B样条（平滑应力） - 流体：三次B样条（高精度） - 颗粒：线性（计算效率）

3. 混合形函数： 在过渡区域使用加权组合： $$N^{\text{mix}} = \alpha N^{\text{low}} + (1-\alpha) N^{\text{high}}$$ 形函数选择决策树：

材料类型？
├─ 弹性固体
│  └─ 二次B样条（C¹连续性）
├─ 流体
│  ├─ 低雷诺数 → 三次B样条
│  └─ 高雷诺数 → 二次B样条
├─ 颗粒材料
│  └─ 线性 + GIMP（防穿透）
└─ 大变形塑性
   └─ 自适应（应变率相关）

12.1.4 时间积分方案

MPM中常用的时间积分方案包括：

1. 显式积分： - 前向欧拉：简单但稳定性受限 - 辛欧拉：保持能量守恒性质 - Runge-Kutta方法：高阶精度

2. 半隐式积分： - 对某些项（如弹性力）进行隐式处理 - 提高稳定性，允许更大的时间步长

CFL条件： $$\Delta t \leq C \frac{h}{c}$$ 其中 $c = \sqrt{E/\rho}$ 是材料中的声速，$C$ 是CFL数（通常取0.1-0.5）。

12.1.5 本构模型

MPM可以处理各种材料模型：

1. 线弹性： $$\boldsymbol{\sigma} = \lambda \text{tr}(\boldsymbol{\epsilon})\mathbf{I} + 2\mu\boldsymbol{\epsilon}$$

2. 超弹性（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})$。

3. 塑性模型： - von Mises屈服准则 - Drucker-Prager模型（适用于颗粒材料） - 雪的本构模型

应力计算通常涉及：

1. 计算应变或变形梯度
2. 应用本构关系获得应力
3. 考虑塑性修正（如果适用）

12.2 MLS-MPM（移动最小二乘物质点法）

12.2.1 传统MPM的局限性

传统MPM虽然在处理大变形问题上有优势，但存在一些固有问题：

1. 数值耗散：粒子-网格-粒子的传输过程引入人工粘性
2. 角动量守恒：FLIP（Fluid Implicit Particle）方法保持角动量但引入噪声
3. 计算效率：需要存储和更新大量粒子信息
4. 精度问题：低阶形函数限制了应力计算的精度

12.2.2 MLS-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)$。

12.2.3 MLS-MPM算法流程

1. 初始化： - 设置粒子位置 $\mathbf{x}_p$、速度 $\mathbf{v}_p$、质量 $m_p$ - 初始化变形梯度 $\mathbf{F}_p = \mathbf{I}$ - 计算初始 APIC 矩阵 $\mathbf{C}_p = \mathbf{0}$

2. 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)]$$

3. 网格动量更新（与传统MPM相同）： $$\mathbf{v}_i^{n+1} = \mathbf{v}_i^n + \frac{\Delta t}{m_i} \mathbf{f}_i$$

4. 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}$$

5. 粒子更新： $$\mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \Delta t \mathbf{v}_p^{n+1}$$ $$\mathbf{F}_p^{n+1} = (\mathbf{I} + \Delta t \mathbf{C}_p^{n+1}) \mathbf{F}_p^n$$

12.2.4 优化技巧

1. 惯性张量的高效计算： 对于规则网格和对称核函数： $$\mathbf{D}_p = \frac{4h^2}{3}\mathbf{I}$$（二维情况）$$\mathbf{D}_p = \frac{h^2}{2}\mathbf{I}$$（三维情况）

2. 体积守恒： 通过修正变形梯度保持体积： $$\mathbf{F}_p^{\text{corrected}} = J^{1/d} \mathbf{F}_p / \det(\mathbf{F}_p)^{1/d}$$

3. 数值稳定性增强： - 限制变形梯度的奇异值 - 使用隐式积分处理刚性材料 - 自适应时间步长

12.2.5 与APIC的关系

MLS-MPM可以看作是APIC（Affine Particle-In-Cell）方法的推广：

1. APIC：保存仿射速度场 $\mathbf{v}(\mathbf{x}) = \mathbf{v}_p + \mathbf{C}_p(\mathbf{x} - \mathbf{x}_p)$
2. MLS-MPM：通过MLS重构获得最优的 $\mathbf{C}_p$
3. 理论联系：MLS-MPM在数学上等价于使用特定权重的APIC

优势对比：

• APIC：实现简单，但精度受限
• MLS-MPM：更高精度，更好的动量守恒
• 计算成本相近

12.3 高级MPM技术与应用

12.3.1 多材料MPM

实际应用中经常需要模拟多种材料的相互作用，如流体-固体耦合、多相流等。

多材料表示方法：

1. 分离粒子法： - 不同材料使用不同的粒子集 - 在网格上进行耦合 - 适用于界面清晰的情况

2. 混合粒子法： - 单个粒子可以包含多种材料 - 使用体积分数 $\phi_k$ 表示材料 $k$ 的占比 - 适用于混合材料

界面处理：

1. 接触算法： $$\mathbf{v}_i^{k,\text{final}} = \mathbf{v}_i^k + \frac{\Delta t}{m_i^k} \mathbf{f}_i^{\text{contact}}$$ 接触力通过罚函数或拉格朗日乘子计算。

2. 混合理论： 应力通过体积加权平均： $$\boldsymbol{\sigma} = \sum_k \phi_k \boldsymbol{\sigma}_k$$ 表面张力： 对于流体界面，需要考虑表面张力： $$\mathbf{f}_{\text{surface}} = \gamma \kappa \mathbf{n}$$ 其中 $\gamma$ 是表面张力系数，$\kappa$ 是曲率，$\mathbf{n}$ 是法向量。

12.3.2 自适应MPM

为了提高计算效率和精度，可以使用自适应技术：

1. 空间自适应： - 自适应网格细化（AMR）：在需要高精度的区域使用更细的网格 - 粒子分裂与合并：根据变形程度调整粒子密度

分裂准则： $$\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}}$$

2. 时间自适应： 根据局部CFL条件调整时间步长： $$\Delta t_p = \min\left(C_{\text{CFL}} \frac{h}{||\mathbf{v}_p|| + c_p}, \Delta t_{\text{max}}\right)$$

3. 模型自适应： - 在小变形区域使用线性模型 - 在大变形区域使用非线性模型 - 动态切换本构模型

12.3.3 GPU加速技术

MPM天然适合GPU并行化，关键优化策略包括：

1. 数据结构优化： - 粒子排序：按空间位置排序减少内存访问冲突 - 网格哈希：使用空间哈希表快速查找粒子-网格映射 - 稀疏网格：只分配活跃网格节点

2. 并行化策略： - P2G阶段：每个粒子独立计算，使用原子操作累加到网格 - 网格更新：每个网格节点独立更新 - G2P阶段：每个粒子独立更新

3. 性能优化： - 共享内存使用：缓存频繁访问的网格数据 - Warp级别优化：确保线程束内的一致性 - 混合精度计算：在不影响精度的地方使用float16

典型加速比：

• CPU单线程 → GPU：100-1000倍
• 实时模拟成为可能（>30 FPS）

12.3.4 工程应用案例

1. 计算机动画： - 雪崩模拟（迪士尼《冰雪奇缘》） - 沙尘效果 - 流体-固体交互

2. 地质工程： - 滑坡预测 - 土壤液化分析 - 爆破效果模拟

3. 生物力学： - 软组织变形 - 血流动力学 - 细胞力学

4. 工业设计： - 3D打印过程模拟 - 材料成型分析 - 碰撞安全测试

12.3.5 前沿研究方向

1. 机器学习增强： - 使用神经网络学习本构模型 - 加速结构预测 - 参数自动调优

2. 多尺度MPM： - 宏观-介观-微观耦合 - 分子动力学与连续介质耦合 - 自适应尺度转换

3. 拓扑优化： - 基于MPM的形状优化 - 材料分布优化 - 多目标优化

4. 实时交互系统： - VR/AR应用 - 触觉反馈 - 实时物理编辑

本章小结

本章深入介绍了物质点法（MPM）及其现代变种MLS-MPM，这是计算机图形学中模拟大变形材料的强大工具。

关键概念：

1. MPM基础：结合拉格朗日粒子和欧拉网格的混合方法
2. 核心算法：P2G → 网格求解 → G2P → 状态更新
3. MLS-MPM：通过移动最小二乘改进精度和稳定性
4. 高级技术：多材料、自适应、GPU加速

重要公式：

• 变形梯度更新：$\mathbf{F}^{n+1} = (\mathbf{I} + \Delta t \mathbf{L}) \mathbf{F}^n$
• 内力计算：$\mathbf{f}_i^{\text{int}} = -\sum_p V_p \boldsymbol{\sigma}_p \nabla w_{ip}$
• MLS速度场：$\mathbf{v}(\mathbf{x}) = \mathbf{c}_0 + \mathbf{C}_1 (\mathbf{x} - \mathbf{x}_p)$
• CFL条件：$\Delta t \leq C h/c$

算法选择指南：

• 传统MPM：简单实现，适合教学
• MLS-MPM：生产环境推荐，更好的精度
• GIMP：需要精确接触处理时使用
• 自适应MPM：大规模场景优化

练习题

基础题

1. 形函数性质 证明二次B样条形函数满足分割单位性：$\sum_i N_i^2(\mathbf{x}) = 1$。

提示：考虑一维情况，利用B样条的递归定义。

2. 变形梯度的物理意义 给定变形梯度 $\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}$。

3. CFL稳定性分析 对于杨氏模量 $E = 10^6$ Pa、密度 $\rho = 1000$ kg/m³ 的材料，网格间距 $h = 0.01$ m，计算最大允许时间步长（取CFL数 $C = 0.3$）。

提示：声速 $c = \sqrt{E/\rho}$。

挑战题

4. 能量守恒分析 证明在无外力和无耗散的情况下，显式MPM的总能量（动能+势能）在一个时间步内的变化为 $O(\Delta t^2)$。

提示：考虑辛结构和中点法则。

5. MLS-MPM的理论分析 推导MLS-MPM中APIC矩阵 $\mathbf{C}_p$ 的更新公式，并说明其如何保持角动量守恒。

提示：从最小二乘问题出发，考虑惯性张量的作用。

6. 多材料耦合 设计一个算法处理流固耦合问题，其中固体使用Neo-Hookean模型，流体使用不可压缩条件。讨论界面条件的处理。

提示：考虑分离式求解和耦合约束。

常见陷阱与错误

1. 数值不稳定 - 错误：使用过大的时间步长导致爆炸 - 解决：严格遵守CFL条件，考虑自适应时间步

2. 质量损失 - 错误：网格节点质量过小导致除零 - 解决：设置最小质量阈值 $m_{\min} = 10^{-10} m_{\text{avg}}$

3. 应力振荡 - 错误：使用线性形函数导致力的不连续 - 解决：改用B样条或GIMP形函数

4. 内存泄漏 - 错误：动态分配网格但忘记释放 - 解决：使用稀疏数据结构和智能指针

5. 并行化错误 - 错误：P2G阶段的竞争条件 - 解决：使用原子操作或着色算法

6. 边界处理不当 - 错误：粒子穿透固定边界 - 解决：正确实现碰撞检测和响应

最佳实践检查清单

算法设计

• [ ] 选择合适的形函数（精度vs效率权衡）
• [ ] 实现稳定的时间积分（CFL条件）
• [ ] 正确处理边界条件
• [ ] 考虑能量守恒性质

数值稳定性

• [ ] 避免极小质量的除法
• [ ] 限制变形梯度的条件数
• [ ] 使用适当的浮点精度
• [ ] 实现数值阻尼（如需要）

性能优化

• [ ] 利用空间局部性（粒子排序）
• [ ] 实现稀疏网格结构
• [ ] 考虑GPU并行化
• [ ] 使用高效的数据结构

物理正确性

• [ ] 验证质量守恒
• [ ] 检查动量守恒
• [ ] 测试能量行为
• [ ] 验证本构模型实现

代码质量

• [ ] 模块化设计（分离物理和数值）
• [ ] 充分的单元测试
• [ ] 性能基准测试
• [ ] 清晰的文档和注释

调试建议

• [ ] 实现可视化调试工具
• [ ] 添加物理量监控（质量、动量、能量）
• [ ] 使用简单测试案例验证
• [ ] 逐步增加复杂度