第7章：混合欧拉-拉格朗日视角（下）

本章深入探讨物质点法（Material Point Method, MPM）——一种强大的混合欧拉-拉格朗日方法，它结合了拉格朗日粒子的灵活性和欧拉网格的计算效率。我们将从MPM的理论基础出发，详细分析各种本构模型，深入理解MLS-MPM算法的数学原理，并探索数值断裂和连续介质伤害力学在MPM框架下的实现。通过本章的学习，你将掌握模拟复杂材料行为所需的核心技术，包括弹塑性变形、大变形、断裂等现象的数值方法。

7.1 物质点法（MPM）理论基础

物质点法最初由Sulsky等人于1994年提出，作为粒子元胞法（PIC）在固体力学中的扩展。MPM的核心思想是使用拉格朗日粒子携带材料信息，同时利用背景欧拉网格进行动量方程求解，从而避免了纯拉格朗日方法中的网格畸变问题，同时保持了材料历史信息的准确追踪。

7.1.1 从PIC到MPM的演化

粒子元胞法（PIC）最初由Harlow和Welch于1965年提出用于流体模拟，其基本思想是将连续介质离散为拉格朗日粒子，在固定的欧拉网格上求解动量方程。这种混合方法结合了两种视角的优势：粒子追踪材料历史，网格提供高效的空间导数计算。MPM继承了PIC的核心框架，但做出了关键改进：

材料点概念：MPM中的粒子不仅携带质量和动量，还携带完整的材料状态信息，包括： - 变形梯度 $\mathbf{F}_p$ - 描述局部变形状态 - 应力张量 $\boldsymbol{\sigma}_p$ - 当前应力状态 - 材料参数（如杨氏模量 $E$、泊松比 $\nu$、密度 $\rho_0$ 等） - 塑性应变 $\boldsymbol{\epsilon}^p_p$ 等历史变量 - 损伤变量 $d_p$ 或其他内部状态变量
本构模型集成：MPM可以自然地集成各种复杂的本构模型，这是PIC难以实现的。通过在每个材料点上独立计算应力更新，MPM能够模拟： - 弹塑性材料（金属、土壤） - 粘弹性材料（聚合物、生物组织） - 损伤和断裂（脆性破坏、韧性断裂） - 各向异性材料（纤维增强复合材料、晶体） - 多相材料（饱和土、泡沫材料）
守恒性质：MPM在离散层面严格保证质量守恒，并在适当的时间积分格式下保证动量守恒。能量守恒虽然不是严格满足的，但可以通过改进的传输算法显著改善： - 传统MPM：存在显著的数值耗散 - FLIP（Fluid Implicit Particles）：减少耗散但可能引入噪声 - APIC（Affine PIC）：改善角动量守恒 - MLS-MPM：进一步提高守恒性和稳定性

PIC与MPM的关键区别

虽然MPM脱胎于PIC，但两者在应用领域和技术细节上存在重要差异：

| 特性 | PIC | MPM |

特性	PIC	MPM
主要应用	流体（等离子体、可压缩流）	固体、流固耦合
材料历史	仅追踪位置和速度	完整的变形历史
应力计算	通常基于压力	完整的应力张量
网格使用	每步重新初始化	可能保留部分信息
边界处理	相对简单	需要特殊处理

7.1.2 连续介质力学回顾

在深入MPM的数学框架之前，我们需要回顾连续介质力学的基本概念。这些概念构成了MPM理论的基础。

运动学描述

考虑一个连续体在参考构型 $\Omega_0$（$t=0$ 时的初始构型）和当前构型 $\Omega_t$（$t$ 时刻的变形构型）之间的运动。材料点的运动由映射 $\boldsymbol{\phi}: \Omega_0 \times [0,T] \rightarrow \mathbb{R}^d$ 描述：

$$\mathbf{x} = \boldsymbol{\phi}(\mathbf{X}, t)$$ 其中 $\mathbf{X} \in \Omega_0$ 是材料坐标（拉格朗日坐标），$\mathbf{x} \in \Omega_t$ 是空间坐标（欧拉坐标）。

变形梯度张量 $\mathbf{F}$ 是运动学的核心量，定义为： $$\mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \nabla_{\mathbf{X}} \boldsymbol{\phi}$$ 在分量形式下： $$F_{ij} = \frac{\partial x_i}{\partial X_j}$$ 变形梯度包含了局部变形的完整信息：

体积变化：$J = \det(\mathbf{F})$ （雅可比行列式）
局部旋转和拉伸：通过极分解 $\mathbf{F} = \mathbf{R}\mathbf{U} = \mathbf{V}\mathbf{R}$
主拉伸：$\mathbf{U}$ 或 $\mathbf{V}$ 的特征值

速度场和加速度场分别定义为： $$\mathbf{v}(\mathbf{X}, t) = \frac{\partial \boldsymbol{\phi}}{\partial t}\bigg|_{\mathbf{X}}, \quad \mathbf{a}(\mathbf{X}, t) = \frac{\partial^2 \boldsymbol{\phi}}{\partial t^2}\bigg|_{\mathbf{X}}$$ 速度梯度张量 $\mathbf{L}$ 在空间描述下定义为： $$\mathbf{L} = \nabla_{\mathbf{x}} \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)$ 是自旋张量。

守恒方程

连续介质的运动由以下守恒方程控制：

质量守恒（连续性方程）：

欧拉描述（空间形式）： $$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0$$ 或等价地： $$\frac{D\rho}{Dt} + \rho \nabla \cdot \mathbf{v} = 0$$ 拉格朗日描述（材料形式）： $$\rho J = \rho_0$$ 其中 $J = \det(\mathbf{F})$ 是雅可比行列式，表示体积比。这个关系表明材料点的质量在变形过程中保持不变。

动量守恒（运动方程）：

空间形式（使用柯西应力）： $$\rho \frac{D\mathbf{v}}{Dt} = \nabla \cdot \boldsymbol{\sigma} + \rho \mathbf{b}$$ 材料形式（使用第一Piola-Kirchhoff应力）： $$\rho_0 \frac{\partial^2 \boldsymbol{\phi}}{\partial t^2} = \nabla_{\mathbf{X}} \cdot \mathbf{P} + \rho_0 \mathbf{b}_0$$ 其中：

$\boldsymbol{\sigma}$ 是柯西应力张量（真实应力）
$\mathbf{P} = J\boldsymbol{\sigma}\mathbf{F}^{-T}$ 是第一Piola-Kirchhoff应力
$\mathbf{b}$ 是当前构型下的体力密度
$\frac{D}{Dt} = \frac{\partial}{\partial t} + \mathbf{v} \cdot \nabla$ 是材料导数

角动量守恒：

要求柯西应力张量的对称性： $$\boldsymbol{\sigma} = \boldsymbol{\sigma}^T$$ 这在存在体力偶或极性材料时需要修正。

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

$e$ 是比内能
$\mathbf{d} = \text{sym}(\mathbf{L})$ 是变形率张量
$\boldsymbol{\sigma} : \mathbf{d}$ 是应力功率
$\mathbf{q}$ 是热流向量
$r$ 是单位质量的热源

应力度量与客观性

不同的应力度量在MPM中都有应用：

柯西应力 $\boldsymbol{\sigma}$：当前构型上的真实应力
第一Piola-Kirchhoff应力 $\mathbf{P} = J\boldsymbol{\sigma}\mathbf{F}^{-T}$：功共轭于 $\mathbf{F}$
第二Piola-Kirchhoff应力 $\mathbf{S} = \mathbf{F}^{-1}\mathbf{P} = J\mathbf{F}^{-1}\boldsymbol{\sigma}\mathbf{F}^{-T}$：材料描述
Kirchhoff应力 $\boldsymbol{\tau} = J\boldsymbol{\sigma}$：在近似不可压缩材料中有用

客观性要求：本构关系必须满足框架无关性（frame indifference），即在刚体运动下保持不变。

7.1.3 MPM的数学框架

MPM的数学框架建立在弱形式（weak form）的基础上，这是有限元方法的核心思想。弱形式避免了对解的二阶导数要求，使得数值方法更加稳定和灵活。

从强形式到弱形式

动量守恒的强形式为： $$\rho \mathbf{a} = \nabla \cdot \boldsymbol{\sigma} + \rho \mathbf{b} \quad \text{in } \Omega$$ 边界条件：

Dirichlet边界条件：$\mathbf{u} = \bar{\mathbf{u}}$ on $\Gamma_u$
Neumann边界条件：$\boldsymbol{\sigma} \cdot \mathbf{n} = \mathbf{t}$ on $\Gamma_t$

通过乘以试函数 $\mathbf{w}$（满足 $\mathbf{w} = \mathbf{0}$ on $\Gamma_u$）并在域上积分，应用分部积分和散度定理，得到弱形式： $$\int_{\Omega} \rho \mathbf{a} \cdot \mathbf{w} \, d\Omega = -\int_{\Omega} \boldsymbol{\sigma} : \nabla \mathbf{w} \, d\Omega + \int_{\Omega} \rho \mathbf{b} \cdot \mathbf{w} \, d\Omega + \int_{\Gamma_t} \mathbf{t} \cdot \mathbf{w} \, d\Gamma$$ 这个弱形式是MPM的出发点，其优势在于：

降低对解的光滑性要求
自然处理Neumann边界条件
保证离散系统的对称性（对于自伴随问题）

MPM离散化策略

MPM采用独特的双重离散化策略：

材料点（拉格朗日离散）：将连续体离散为 $N_p$ 个材料点，每个材料点 $p$ 携带： - 质量 $m_p = \rho_0 V_p^0$（保持不变） - 初始体积 $V_p^0$，当前体积 $V_p = J_p V_p^0$ - 位置 $\mathbf{x}_p(t)$ - 速度 $\mathbf{v}_p(t)$ - 变形梯度 $\mathbf{F}_p(t)$ - 柯西应力 $\boldsymbol{\sigma}_p(t)$ - 内部状态变量（如塑性应变、损伤等）
背景网格（欧拉离散）：使用固定的背景网格进行空间导数计算和动量方程求解。网格节点 $i$ 具有： - 位置 $\mathbf{x}_i$（固定） - 形函数 $N_i(\mathbf{x})$ - 临时存储的质量、动量、内力等

离散化的数学表达

将连续场用材料点近似： $$\rho(\mathbf{x}, t) \approx \sum_p m_p \delta(\mathbf{x} - \mathbf{x}_p(t))$$ $$\mathbf{v}(\mathbf{x}, t) \approx \sum_p \mathbf{v}_p(t) \frac{m_p}{\rho(\mathbf{x}, t)} \delta(\mathbf{x} - \mathbf{x}_p(t))$$ 粒子到网格的映射（P2G）

网格节点的质量通过粒子质量的加权和计算： $$m_i = \sum_p m_p N_i(\mathbf{x}_p)$$ 网格节点的动量通过粒子动量的加权和计算： $$(m\mathbf{v})_i = \sum_p m_p \mathbf{v}_p N_i(\mathbf{x}_p)$$ 网格节点的速度： $$\mathbf{v}_i = \frac{(m\mathbf{v})_i}{m_i}$$ 网格节点的内力由应力散度的离散形式给出： $$\mathbf{f}_i^{int} = -\sum_p V_p \boldsymbol{\sigma}_p \nabla N_i(\mathbf{x}_p)$$ 这个公式的推导来自于弱形式中的应力项： $$-\int_{\Omega} \boldsymbol{\sigma} : \nabla \mathbf{w} \, d\Omega \approx -\sum_p V_p \boldsymbol{\sigma}_p : \nabla \mathbf{w}(\mathbf{x}_p)$$ 选择 $\mathbf{w} = N_i \mathbf{e}_j$（其中 $\mathbf{e}_j$ 是单位向量），得到节点力的第 $j$ 分量。

网格更新

在网格上求解离散化的动量方程： $$m_i \mathbf{a}_i = \mathbf{f}_i^{int} + \mathbf{f}_i^{ext}$$ 其中外力可能包括：

体力：$\mathbf{f}_i^{body} = \sum_p m_p \mathbf{b} N_i(\mathbf{x}_p)$
边界牵引力：通过边界积分计算

时间积分方案的选择：

显式积分（常用于动态问题）： $$\mathbf{v}_i^{n+1} = \mathbf{v}_i^n + \Delta t \mathbf{a}_i^n$$ CFL稳定性条件：$\Delta t \leq \alpha \frac{h}{c}$，其中 $c = \sqrt{K/\rho}$ 是声速
隐式积分（用于准静态或刚性问题）：需要求解非线性系统，通常采用Newton-Raphson迭代

网格到粒子的映射（G2P）

粒子速度通过网格速度插值更新： $$\mathbf{v}_p^{n+1} = \sum_i \mathbf{v}_i^{n+1} N_i(\mathbf{x}_p^n)$$ 粒子位置更新： $$\mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \Delta t \mathbf{v}_p^{n+1}$$ 变形梯度更新（关键步骤）： $$\mathbf{F}_p^{n+1} = \left(\mathbf{I} + \Delta t \sum_i \mathbf{v}_i^{n+1} \otimes \nabla N_i(\mathbf{x}_p^n)\right) \mathbf{F}_p^n$$ 这个更新公式来自于速度梯度的定义： $$\mathbf{L} = \nabla \mathbf{v} = \dot{\mathbf{F}}\mathbf{F}^{-1}$$ 离散化后： $$\frac{\mathbf{F}^{n+1} - \mathbf{F}^n}{\Delta t} \approx \mathbf{L}^{n+1} \mathbf{F}^n$$

7.1.4 离散化与权重函数

权重函数（形函数）的选择对MPM的精度、稳定性和计算效率至关重要。不同的权重函数在精度和计算成本之间提供不同的权衡。

线性插值函数（Tent Function）

对于规则的笛卡尔网格，多维线性插值函数通过张量积构造： $$N_i(\mathbf{x}) = \prod_{d=1}^{dim} N^1\left(\frac{x_d - x_{i,d}}{h}\right)$$ 其中一维线性基函数（帽函数）为： $$N^1(\xi) = \begin{cases} 1 - |\xi| & \text{if } |\xi| < 1 \\ 0 & \text{otherwise} \end{cases}$$ 梯度计算： $$\frac{\partial N^1}{\partial \xi} = \begin{cases} -\text{sign}(\xi) & \text{if } |\xi| < 1 \\ 0 & \text{otherwise} \end{cases}$$ 优点：计算简单，有限支撑（2×2×2网格单元）缺点：$C^0$ 连续，梯度不连续可能导致网格交叉误差（grid crossing error）

二次B样条

二次B样条提供 $C^1$ 连续性： $$N^2(\xi) = \begin{cases} \frac{3}{4} - \xi^2 & \text{if } |\xi| < \frac{1}{2} \\ \frac{1}{2}(\frac{3}{2} - |\xi|)^2 & \text{if } \frac{1}{2} \leq |\xi| < \frac{3}{2} \\ 0 & \text{otherwise} \end{cases}$$ 梯度： $$\frac{\partial N^2}{\partial \xi} = \begin{cases} -2\xi & \text{if } |\xi| < \frac{1}{2} \\ \text{sign}(\xi)(|\xi| - \frac{3}{2}) & \text{if } \frac{1}{2} \leq |\xi| < \frac{3}{2} \\ 0 & \text{otherwise} \end{cases}$$ 支撑域：3×3×3网格单元优点：梯度连续，减少网格交叉误差缺点：计算成本略高，支撑域更大

三次B样条

三次B样条提供 $C^2$ 连续性： $$N^3(\xi) = \begin{cases} \frac{1}{2}|\xi|^3 - \xi^2 + \frac{2}{3} & \text{if } |\xi| < 1 \\ -\frac{1}{6}|\xi|^3 + \xi^2 - 2|\xi| + \frac{4}{3} & \text{if } 1 \leq |\xi| < 2 \\ 0 & \text{otherwise} \end{cases}$$ 支撑域：4×4×4网格单元优点：高阶连续性，数值性质优良缺点：计算成本最高，支撑域最大

GIMP（Generalized Interpolation Material Point）

GIMP方法考虑粒子的有限尺寸，使用粒子特征函数 $\chi_p$ 与网格形函数的卷积： $$S_{ip}(\mathbf{x}_p) = \int_{\Omega_p} \chi_p(\mathbf{x}) N_i(\mathbf{x}) d\mathbf{x}$$ 对于矩形粒子域，可以解析计算。GIMP提高了大变形下的稳定性。

权重函数的性质要求

理想的权重函数应满足以下性质：

分片常数性质（Partition of Unity）： $$\sum_i N_i(\mathbf{x}) = 1 \quad \forall \mathbf{x}$$ 保证常应力/应变状态的精确表示
有限支撑（Compact Support）：权重函数只在有限区域内非零，保证计算效率
连续性要求： - 至少 $C^0$ 连续（位移连续） - $C^1$ 连续更优（应变连续） - 更高阶连续性改善数值性质
对称性： $$N(\xi) = N(-\xi)$$ 保证动量守恒
非负性： $$N_i(\mathbf{x}) \geq 0$$ 避免非物理的负质量
插值性质（对于某些方法）： $$N_i(\mathbf{x}_j) = \delta_{ij}$$ 在节点处精确插值

权重函数的选择策略

动态问题：线性插值通常足够，计算效率高
大变形问题：二次或三次B样条，减少数值误差
断裂问题：GIMP或自适应方法，处理粒子分离
准静态问题：高阶B样条，提高精度

7.2 本构模型与材料行为

本构模型（constitutive model）描述了应力与应变之间的关系，是材料力学行为的数学表征。在MPM框架中，本构模型在每个材料点上独立计算，这使得MPM能够自然地处理材料非均匀性、历史相关性和多相材料的耦合。本节将详细介绍各种常用本构模型的数学原理和实现细节。

7.2.1 弹性本构模型

弹性材料的应力-应变关系是可逆的，不依赖于加载历史。卸载后材料能完全恢复原形，不产生永久变形。我们从最基本的超弹性模型开始。

超弹性模型的基本原理

超弹性（hyperelastic）材料的应力由存储的弹性能量决定。定义应变能密度函数（Helmholtz自由能）$\Psi(\mathbf{F})$，应力通过能量对变形的导数得到： $$\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}}$$ 其中 $\mathbf{P}$ 是第一Piola-Kirchhoff应力，功共轭于变形梯度。

其他应力度量的计算：

柯西应力：$\boldsymbol{\sigma} = \frac{1}{J} \mathbf{P} \mathbf{F}^T = \frac{1}{J} \frac{\partial \Psi}{\partial \mathbf{F}} \mathbf{F}^T$
第二Piola-Kirchhoff应力：$\mathbf{S} = \mathbf{F}^{-1}\mathbf{P} = 2\frac{\partial \Psi}{\partial \mathbf{C}}$
Kirchhoff应力：$\boldsymbol{\tau} = J\boldsymbol{\sigma} = \mathbf{P}\mathbf{F}^T$

为保证客观性，应变能密度函数必须仅依赖于变形不变量。对于各向同性材料，$\Psi$ 只能是以下不变量的函数：

$I_1 = \text{tr}(\mathbf{C}) = \text{tr}(\mathbf{F}^T\mathbf{F})$
$I_2 = \frac{1}{2}[(\text{tr}(\mathbf{C}))^2 - \text{tr}(\mathbf{C}^2)]$
$I_3 = \det(\mathbf{C}) = J^2$

其中 $\mathbf{C} = \mathbf{F}^T\mathbf{F}$ 是右柯西-格林变形张量。

Neo-Hookean模型

Neo-Hookean模型是最简单且广泛使用的超弹性模型，适用于橡胶类材料的大变形： $$\Psi(\mathbf{F}) = \frac{\mu}{2}(\mathrm{tr}(\mathbf{F}^T\mathbf{F}) - d) - \mu \ln J + \frac{\lambda}{2}(\ln J)^2$$ 其中：

$\mu = \frac{E}{2(1+\nu)}$ 是剪切模量
$\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}$ 是第一拉梅常数
$E$ 是杨氏模量，$\nu$ 是泊松比
$d$ 是空间维度（2D或3D）
$J = \det(\mathbf{F})$ 是体积比

第一Piola-Kirchhoff应力： $$\mathbf{P} = \mu(\mathbf{F} - \mathbf{F}^{-T}) + \lambda \ln J \mathbf{F}^{-T}$$ 柯西应力： $$\boldsymbol{\sigma} = \frac{\mu}{J}(\mathbf{b} - \mathbf{I}) + \frac{\lambda}{J}\ln J \mathbf{I}$$ 其中 $\mathbf{b} = \mathbf{F}\mathbf{F}^T$ 是左柯西-格林变形张量。

数值实现注意事项：

当 $J \rightarrow 0$ 时，$\ln J \rightarrow -\infty$，需要数值稳定化
对于近不可压缩材料（$\nu \rightarrow 0.5$），$\lambda \rightarrow \infty$，需要特殊处理
可以引入体积模量 $K = \lambda + \frac{2\mu}{3}$ 重新参数化

固定协旋模型（Fixed Corotated Model）

对于大旋转但小应变的情况（如弯曲的梁、薄壳），固定协旋模型更为合适： $$\Psi(\mathbf{F}) = \mu \mathrm{tr}[(\mathbf{F} - \mathbf{R})^T(\mathbf{F} - \mathbf{R})] + \frac{\lambda}{2}(J-1)^2$$ 其中 $\mathbf{R}$ 是 $\mathbf{F}$ 的旋转部分，通过极分解获得： $$\mathbf{F} = \mathbf{R}\mathbf{S} = \mathbf{V}\mathbf{R}$$ 极分解的计算步骤：

计算 $\mathbf{C} = \mathbf{F}^T\mathbf{F}$
特征分解：$\mathbf{C} = \mathbf{Q}\boldsymbol{\Lambda}^2\mathbf{Q}^T$
$\mathbf{S} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^T$
$\mathbf{R} = \mathbf{F}\mathbf{S}^{-1}$

第一Piola-Kirchhoff应力： $$\mathbf{P} = 2\mu(\mathbf{F} - \mathbf{R}) + \lambda(J-1)J\mathbf{F}^{-T}$$ 该模型的优点：

在刚体旋转下不产生应力
小应变时退化为线性弹性
适合模拟固体材料

St. Venant-Kirchhoff模型

另一种常用的大变形弹性模型： $$\Psi(\mathbf{F}) = \frac{\lambda}{2}[\text{tr}(\mathbf{E})]^2 + \mu \text{tr}(\mathbf{E}^2)$$ 其中 $\mathbf{E} = \frac{1}{2}(\mathbf{C} - \mathbf{I})$ 是Green-Lagrange应变张量。

该模型在大压缩时可能不稳定，但在小应变情况下表现良好。

各向异性弹性

对于各向异性材料，需要引入材料方向。以横观各向同性材料为例，定义纤维方向 $\mathbf{a}_0$，应变能函数可写为： $$\Psi(\mathbf{F}, \mathbf{a}_0) = \Psi_{iso}(\mathbf{F}) + \Psi_{aniso}(I_4, I_5)$$ 其中 $I_4 = \mathbf{a}_0 \cdot (\mathbf{C} \mathbf{a}_0)$，$I_5 = \mathbf{a}_0 \cdot (\mathbf{C}^2 \mathbf{a}_0)$，$\mathbf{C} = \mathbf{F}^T\mathbf{F}$ 是右柯西-格林变形张量。

7.2.2 塑性与屈服准则

塑性变形是不可逆的，需要将变形梯度分解为弹性和塑性部分： $$\mathbf{F} = \mathbf{F}^e \mathbf{F}^p$$ 屈服准则

塑性流动发生在应力状态满足屈服条件时。常用的屈服准则包括：

von Mises准则（适用于金属）： $$f(\boldsymbol{\sigma}) = \sqrt{\frac{3}{2}\mathbf{s}:\mathbf{s}} - \sigma_y = 0$$ 其中 $\mathbf{s} = \boldsymbol{\sigma} - \frac{1}{3}\mathrm{tr}(\boldsymbol{\sigma})\mathbf{I}$ 是偏应力，$\sigma_y$ 是屈服应力。
Drucker-Prager准则（适用于土壤、岩石）： $$f(\boldsymbol{\sigma}) = \sqrt{J_2} + \alpha I_1 - k = 0$$ 其中 $I_1 = \mathrm{tr}(\boldsymbol{\sigma})$，$J_2 = \frac{1}{2}\mathbf{s}:\mathbf{s}$，$\alpha$ 和 $k$ 是材料参数。

塑性流动法则

塑性应变率由流动法则确定： $$\dot{\boldsymbol{\epsilon}}^p = \dot{\gamma} \frac{\partial g}{\partial \boldsymbol{\sigma}}$$ 其中 $\dot{\gamma}$ 是塑性乘子，$g$ 是塑性势函数。对于关联流动法则，$g = f$。

返回映射算法

在MPM中，通常使用返回映射算法更新应力：

弹性预测：假设增量全部为弹性 $$\boldsymbol{\sigma}^{trial} = \mathbb{C} : (\boldsymbol{\epsilon}^{n+1} - \boldsymbol{\epsilon}^p_n)$$
检查屈服： $$f^{trial} = f(\boldsymbol{\sigma}^{trial})$$
塑性修正：如果 $f^{trial} > 0$，求解： $$\boldsymbol{\sigma}^{n+1} = \boldsymbol{\sigma}^{trial} - \Delta\gamma \mathbb{C} : \frac{\partial f}{\partial \boldsymbol{\sigma}}$$ 使得 $f(\boldsymbol{\sigma}^{n+1}) = 0$。

硬化法则

材料的屈服面可能随塑性变形而演化：

各向同性硬化：屈服面均匀扩大 $$\sigma_y = \sigma_{y0} + H \epsilon^p_{eq}$$ 其中 $H$ 是硬化模量，$\epsilon^p_{eq}$ 是等效塑性应变。
运动硬化：屈服面平移 $$f = ||\mathbf{s} - \boldsymbol{\alpha}|| - \sigma_y$$ 其中 $\boldsymbol{\alpha}$ 是背应力张量。

7.2.3 粘弹性材料

粘弹性材料表现出时间相关的力学行为，结合了弹性固体和粘性流体的特性。

广义Maxwell模型

粘弹性行为可以用弹簧和阻尼器的组合来描述。广义Maxwell模型的应力演化方程为： $$\boldsymbol{\sigma} + \tau \dot{\boldsymbol{\sigma}} = \mathbb{C}_{\infty} : \boldsymbol{\epsilon} + \tau \mathbb{C}_0 : \dot{\boldsymbol{\epsilon}}$$ 其中 $\tau$ 是松弛时间，$\mathbb{C}_0$ 和 $\mathbb{C}_{\infty}$ 分别是瞬时和长期弹性张量。

Prony级数表示

复杂的粘弹性行为可以用Prony级数表示： $$G(t) = G_{\infty} + \sum_{i=1}^N G_i e^{-t/\tau_i}$$ 其中 $G(t)$ 是剪切松弛模量，$G_i$ 和 $\tau_i$ 是Prony参数。

MPM中的实现

在MPM中，每个材料点需要存储历史变量来追踪粘弹性状态。常用方法是存储内变量 $\mathbf{q}_i$： $$\boldsymbol{\sigma} = \mathbb{C}_{\infty} : \boldsymbol{\epsilon} + \sum_{i=1}^N \mathbf{q}_i$$ 内变量通过以下方程更新： $$\dot{\mathbf{q}}_i = -\frac{1}{\tau_i}\mathbf{q}_i + \frac{G_i}{\tau_i}\boldsymbol{\epsilon}$$

7.2.4 各向异性材料

各向异性材料在不同方向上具有不同的力学性质，常见于复合材料、生物组织等。

正交各向异性弹性

对于正交各向异性材料，弹性张量在材料主轴坐标系下具有特殊形式： $$\mathbb{C} = \begin{bmatrix} C_{1111} & C_{1122} & C_{1133} & 0 & 0 & 0 \\ C_{1122} & C_{2222} & C_{2233} & 0 & 0 & 0 \\ C_{1133} & C_{2233} & C_{3333} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{2323} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{1313} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{1212} \end{bmatrix}$$ 纤维增强材料

对于纤维增强复合材料，应变能可以分解为： $$\Psi = \Psi_{matrix}(\mathbf{F}) + \sum_{i=1}^{N_f} \Psi_{fiber}^{(i)}(I_4^{(i)})$$ 其中 $I_4^{(i)} = \mathbf{a}_0^{(i)} \cdot (\mathbf{C} \mathbf{a}_0^{(i)})$ 是第 $i$ 个纤维方向的不变量。

典型的纤维应变能函数： $$\Psi_{fiber}(I_4) = \begin{cases} \frac{k_1}{2k_2}[\exp(k_2(I_4-1)^2) - 1] & \text{if } I_4 > 1 \\ 0 & \text{if } I_4 \leq 1 \end{cases}$$ 这确保纤维只在拉伸时承载。

材料坐标系变换

在MPM中处理各向异性材料时，需要在全局坐标系和材料坐标系之间进行变换。设 $\mathbf{Q}$ 为旋转矩阵，将材料坐标系映射到全局坐标系，则：

应力变换： $$\boldsymbol{\sigma}_{global} = \mathbf{Q} \boldsymbol{\sigma}_{material} \mathbf{Q}^T$$
弹性张量变换： $$C_{ijkl}^{global} = Q_{im}Q_{jn}Q_{kp}Q_{lq}C_{mnpq}^{material}$$ 损伤演化的各向异性

各向异性材料的损伤演化也表现出方向性。可以引入损伤张量 $\mathbf{D}$： $$\boldsymbol{\sigma}_{effective} = (\mathbf{I} - \mathbf{D}) : \boldsymbol{\sigma}$$ 损伤演化方程： $$\dot{\mathbf{D}} = \frac{1}{\eta} \langle f(\boldsymbol{\sigma}, \mathbf{D}) \rangle \frac{\partial g}{\partial \boldsymbol{\sigma}}$$ 其中 $\eta$ 是损伤粘度参数，$\langle \cdot \rangle$ 是Macaulay括号。

7.3 MLS-MPM算法详解

移动最小二乘物质点法（Moving Least Squares MPM，MLS-MPM）是Hu等人于2018年提出的改进算法，通过在粒子-网格传输中引入局部仿射速度场，显著提高了动量和角动量守恒性，同时简化了实现。

7.3.1 APIC的回顾与局限

在介绍MLS-MPM之前，我们先回顾APIC（Affine Particle-in-Cell）方法。APIC通过在每个粒子上存储仿射速度场来改善角动量守恒： $$\mathbf{v}_p(\mathbf{x}) = \mathbf{v}_p + \mathbf{C}_p(\mathbf{x} - \mathbf{x}_p)$$ 其中 $\mathbf{C}_p$ 是速度梯度矩阵。

APIC的P2G传输： $$m_i\mathbf{v}_i = \sum_p m_p[\mathbf{v}_p + \mathbf{C}_p(\mathbf{x}_i - \mathbf{x}_p)]N_i(\mathbf{x}_p)$$ APIC的G2P传输： $$\mathbf{v}_p = \sum_i \mathbf{v}_i N_i(\mathbf{x}_p)$$ $$\mathbf{C}_p = \sum_i \mathbf{v}_i \otimes \nabla N_i(\mathbf{x}_p)$$ 虽然APIC改善了角动量守恒，但仍存在一些局限：

需要额外存储 $\mathbf{C}_p$ 矩阵
在某些情况下仍有数值耗散
实现相对复杂

7.3.2 MLS-MPM的核心思想

MLS-MPM的关键创新是将APIC的仿射速度场推广到更一般的框架，并通过移动最小二乘（MLS）插值来构建局部速度场。

移动最小二乘插值

给定散布数据点 $\{(\mathbf{x}_j, f_j)\}$，MLS插值在点 $\mathbf{x}$ 处构建局部多项式近似： $$f(\mathbf{x}) \approx \mathbf{p}^T(\mathbf{x})\mathbf{a}(\mathbf{x})$$ 其中 $\mathbf{p}(\mathbf{x})$ 是基函数向量（如 $[1, x, y, xy]^T$），$\mathbf{a}(\mathbf{x})$ 是系数向量，通过最小化加权误差获得： $$\mathbf{a}(\mathbf{x}) = \arg\min_{\mathbf{a}} \sum_j w_j(\mathbf{x})[f_j - \mathbf{p}^T(\mathbf{x}_j)\mathbf{a}]^2$$ 其中 $w_j(\mathbf{x})$ 是权重函数。

MLS-MPM中的速度重构

在MLS-MPM中，我们在每个粒子位置构建局部仿射速度场： $$\mathbf{v}(\mathbf{x}) = \mathbf{v}_p + \mathbf{C}_p(\mathbf{x} - \mathbf{x}_p)$$ 系数通过MLS拟合确定： $$[\mathbf{v}_p, \mathbf{C}_p] = \arg\min \sum_i m_i N_i(\mathbf{x}_p)||\mathbf{v}_i - \mathbf{v}_p - \mathbf{C}_p(\mathbf{x}_i - \mathbf{x}_p)||^2$$

7.3.3 MLS-MPM的数学推导

动量守恒的保证

MLS-MPM通过特殊的传输算法保证动量守恒。定义矩阵： $$\mathbf{D}_p = \sum_i N_i(\mathbf{x}_p)(\mathbf{x}_i - \mathbf{x}_p) \otimes (\mathbf{x}_i - \mathbf{x}_p)$$ 这是局部惯性张量的离散近似。

修正的P2G传输

MLS-MPM的P2G传输包含动量和应力贡献：

动量传输： $$m_i\mathbf{v}_i = \sum_p m_p[\mathbf{v}_p + \mathbf{C}_p(\mathbf{x}_i - \mathbf{x}_p)]N_i(\mathbf{x}_p)$$
应力传输： $$\mathbf{f}_i = -\sum_p V_p[\boldsymbol{\sigma}_p - m_p\mathbf{C}_p\mathbf{D}_p^{-1}]\nabla N_i(\mathbf{x}_p)$$ 注意应力项中的修正项 $m_p\mathbf{C}_p\mathbf{D}_p^{-1}$，这是MLS-MPM的关键创新。

简化的G2P传输

MLS-MPM的一个优势是G2P传输可以大幅简化： $$\mathbf{v}_p = \sum_i \mathbf{v}_i N_i(\mathbf{x}_p)$$ $$\mathbf{C}_p = 4\sum_i \mathbf{v}_i \otimes (\mathbf{x}_i - \mathbf{x}_p)N_i(\mathbf{x}_p) / h^2$$ 其中 $h$ 是网格间距，系数4来自于二维情况下的归一化。

7.3.4 算法实现细节

完整的MLS-MPM算法流程

1. 初始化粒子属性：位置 x_p, 速度 v_p, 变形梯度 F_p, 体积 V_p
2. 对每个时间步：
   a. P2G: 将粒子数据传输到网格

      - 计算网格质量: m_i = Σ_p m_p N_i(x_p)
      - 计算网格动量: (mv)_i = Σ_p m_p[v_p + C_p(x_i - x_p)]N_i(x_p)
      - 计算内力: f_i^int = -Σ_p V_p[σ_p - m_p C_p D_p^{-1}]∇N_i(x_p)

   b. 网格更新

      - 速度更新: v_i = (mv)_i / m_i
      - 施加外力: f_i = f_i^int + f_i^ext
      - 加速度: a_i = f_i / m_i
      - 更新速度: v_i^{n+1} = v_i + Δt a_i
      - 施加边界条件

   c. G2P: 将网格数据传输回粒子

      - 更新速度: v_p = Σ_i v_i^{n+1} N_i(x_p)
      - 更新C矩阵: C_p = 4Σ_i v_i^{n+1} ⊗ (x_i - x_p)N_i(x_p) / h²
      - 更新位置: x_p^{n+1} = x_p + Δt v_p
      - 更新变形梯度: F_p^{n+1} = (I + Δt C_p)F_p

   d. 本构模型更新

      - 根据F_p计算应力σ_p
      - 更新塑性变量等内部状态

数值稳定性考虑

CFL条件： $$\Delta t \leq \frac{h}{c + ||\mathbf{v}||_{max}}$$ 其中 $c = \sqrt{K/\rho}$ 是声速，$K$ 是体积模量。
变形梯度限制：对于大变形问题，需要限制变形梯度的条件数： $$\mathbf{F} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^T$$ $$\sigma_i = \max(\min(\sigma_i, \sigma_{max}), \sigma_{min})$$
粒子分布质量：监控每个网格单元的粒子数，必要时进行粒子重采样。

边界条件处理

MLS-MPM中的边界条件通常在网格上施加：

固定边界： $$\mathbf{v}_i = \mathbf{0} \quad \text{for } \mathbf{x}_i \in \Gamma_D$$
滑移边界： $$\mathbf{v}_i \cdot \mathbf{n} = 0, \quad (\boldsymbol{\sigma} \mathbf{n}) \times \mathbf{n} = \mathbf{0}$$
摩擦边界：使用库仑摩擦模型，详见第6章。

7.3.5 MLS-MPM的优势与扩展

相比传统MPM的优势

更好的守恒性：严格的动量守恒，改善的角动量守恒
更低的数值耗散：通过局部仿射重构减少信息损失
更简单的实现：统一的传输公式，无需特殊处理
更好的稳定性：对粒子分布不均匀性更鲁棒

算法扩展

高阶MLS-MPM：使用二次或更高阶多项式基函数： $$\mathbf{v}(\mathbf{x}) = \mathbf{v}_p + \mathbf{C}_p(\mathbf{x} - \mathbf{x}_p) + \mathbf{H}_p : [(\mathbf{x} - \mathbf{x}_p) \otimes (\mathbf{x} - \mathbf{x}_p)]$$
自适应MLS-MPM：根据变形程度动态调整基函数阶数和支撑域大小。
隐式MLS-MPM：将MLS框架扩展到隐式时间积分，提高大时间步长下的稳定性。

7.4 数值断裂与连续介质伤害力学

材料的断裂和损伤是物理模拟中的重要现象，特别是在计算机图形学、工程仿真和材料科学中。本节探讨如何在MPM框架下实现材料的断裂和损伤演化。

7.4.1 断裂力学基础

应力集中与断裂准则

裂纹尖端存在应力奇异性，应力强度因子（Stress Intensity Factor）描述了这种奇异性的强度： $$K_I = \lim_{r \to 0} \sqrt{2\pi r} \sigma_{yy}(r, \theta=0)$$ 其中 $K_I$ 是I型（张开型）裂纹的应力强度因子。

能量释放率

Griffith能量平衡理论认为，裂纹扩展的驱动力是弹性能释放率 $G$： $$G = -\frac{\partial \Pi}{\partial a}$$ 其中 $\Pi$ 是系统总势能，$a$ 是裂纹长度。

裂纹扩展条件： $$G \geq G_c$$ 其中 $G_c$ 是材料的断裂韧性。

J积分

Rice提出的J积分是路径无关的，适用于弹塑性材料： $$J = \int_{\Gamma} \left(W n_x - \mathbf{t} \cdot \frac{\partial \mathbf{u}}{\partial x}\right) ds$$ 其中 $W$ 是应变能密度，$\Gamma$ 是围绕裂纹尖端的任意路径。

7.4.2 MPM中的数值断裂

在MPM中实现断裂面临独特的挑战和机遇。由于材料点的拉格朗日特性，断裂可以通过粒子分离自然表示。

基于应力的断裂准则

最简单的断裂准则是最大主应力准则：

计算应力张量的特征值： $$\boldsymbol{\sigma} \mathbf{v}_i = \lambda_i \mathbf{v}_i$$
检查最大主应力： $$\sigma_1 > \sigma_{crit}$$
如果满足断裂条件，在主应力方向引入损伤。

基于应变的断裂准则

对于韧性材料，基于应变的准则更合适： $$\epsilon_{eq} = \sqrt{\frac{2}{3}\boldsymbol{\epsilon}^p : \boldsymbol{\epsilon}^p} > \epsilon_{crit}$$ 其中 $\boldsymbol{\epsilon}^p$ 是塑性应变张量。

断裂面的表示

在MPM中，断裂可以通过以下方式表示：

粒子分离法： - 检测满足断裂准则的粒子 - 在断裂面两侧创建新的粒子 - 调整粒子质量和体积以保证守恒
损伤场方法： - 引入连续的损伤变量 $d \in [0,1]$ - 修正应力：$\boldsymbol{\sigma}_{eff} = (1-d)\boldsymbol{\sigma}$ - 损伤演化遵循特定的法则

7.4.3 连续介质伤害力学（CDM）

连续介质伤害力学提供了描述材料逐渐劣化的框架，特别适合MPM实现。

各向同性损伤模型

最简单的损伤模型使用标量损伤变量 $d$： $$\boldsymbol{\sigma} = (1-d)\mathbb{C}:\boldsymbol{\epsilon}^e$$ 损伤演化方程： $$\dot{d} = \frac{1}{\eta}\langle f(Y) - Y_0 \rangle$$ 其中 $Y$ 是损伤驱动力（如应变能释放率），$Y_0$ 是损伤阈值，$\eta$ 是粘性参数。

损伤驱动力的计算

常用的损伤驱动力包括：

应变能密度： $$Y = \frac{1}{2}\boldsymbol{\epsilon}:\mathbb{C}:\boldsymbol{\epsilon}$$
等效应变： $$Y = \sqrt{\boldsymbol{\epsilon}:\boldsymbol{\epsilon}}$$
修正的von Mises应力： $$Y = \frac{\sigma_{eq}^2}{2E(1-d)^2}$$ 各向异性损伤模型

对于复杂材料，需要使用张量损伤变量 $\mathbf{D}$： $$\boldsymbol{\sigma} = (\mathbf{I} - \mathbf{D}):\mathbb{C}:(\mathbf{I} - \mathbf{D}):\boldsymbol{\epsilon}$$ 损伤张量的演化： $$\dot{\mathbf{D}} = \frac{1}{\eta}\langle f(\mathbf{Y}) \rangle \frac{\partial g}{\partial \mathbf{Y}}$$

7.4.4 MPM中的CDM实现

算法框架

在MPM中实现CDM的基本步骤：

1. 在每个材料点存储损伤变量 d_p（或损伤张量 D_p）
2. 在本构更新中：
   a. 计算试应力（假设无损伤）
   b. 计算损伤驱动力 Y
   c. 更新损伤变量：
      d^{n+1} = d^n + Δt * damage_evolution(Y, d^n)
   d. 应用损伤修正：
      σ_p = (1 - d^{n+1}) * σ_trial

3. 处理完全损伤的粒子（d ≈ 1）

损伤局部化与正则化

损伤演化容易导致应变局部化，产生网格依赖性。常用的正则化方法：

非局部损伤模型： $$\bar{Y}_p = \frac{\sum_q w_{pq} Y_q V_q}{\sum_q w_{pq} V_q}$$ 其中 $w_{pq}$ 是基于距离的权重函数。
梯度增强模型：引入损伤梯度项： $$\Psi = \Psi_0(d) + \frac{1}{2}c|\nabla d|^2$$
粘性正则化： $$\dot{d} = \frac{1}{\tau}(d_{eq} - d)$$ 其中 $\tau$ 是特征时间，$d_{eq}$ 是平衡损伤值。

裂纹路径追踪

在MPM中追踪裂纹路径的方法：

相场方法：引入相场变量 $\phi$ 表示裂纹： $$\Gamma_l = \int_{\Omega} \gamma(\phi, \nabla\phi) d\Omega$$ 其中 $\gamma(\phi, \nabla\phi) = \frac{1}{2l}\phi^2 + \frac{l}{2}|\nabla\phi|^2$。
等值面方法：使用等势集函数 $\psi$ 表示裂纹面： $$\Gamma = \{\mathbf{x} : \psi(\mathbf{x}) = 0\}$$

7.4.5 高级断裂模型

内聚力模型（Cohesive Zone Model）

内聚力模型在裂纹面上引入牵引力-位移关系： $$\mathbf{t} = \mathbf{T}(\Delta)$$ 其中 $\Delta$ 是裂纹面的位移跳跃。

典型的内聚力法则：

线性软化： $$t_n = \begin{cases} k\delta_n & \text{if } \delta_n < \delta_0 \\ k\delta_0\frac{\delta_f - \delta_n}{\delta_f - \delta_0} & \text{if } \delta_0 \leq \delta_n < \delta_f \\ 0 & \text{if } \delta_n \geq \delta_f \end{cases}$$
指数软化： $$t_n = k\delta_n \exp\left(-\frac{\delta_n}{\delta_c}\right)$$ 动态断裂

动态断裂涉及惯性效应和波传播：

裂纹扩展速度限制： $$v_{crack} \leq c_R$$ 其中 $c_R$ 是Rayleigh波速。
动态应力强度因子： $$K_I^{dyn} = k(v) K_I^{static}$$ 其中 $k(v)$ 是速度相关的修正因子。

断裂的多尺度建模

原子-连续耦合：在裂纹尖端使用原子模拟，远场使用连续介质。
自适应细化：在高应力梯度区域增加粒子密度。