第9章：软体与可变形物体

章节大纲

9.1 引言与动机

• 软体仿真在具身智能中的重要性
• 从刚体到连续介质的跨越
• 计算挑战与权衡

9.2 有限元方法（FEM）基础

• 9.2.1 连续介质力学回顾
• 9.2.2 弱形式与变分原理
• 9.2.3 离散化与形函数
• 9.2.4 质量矩阵与刚度矩阵
• 9.2.5 时间积分策略

9.3 Position Based Dynamics (PBD)

• 9.3.1 PBD的核心思想
• 9.3.2 约束投影算法
• 9.3.3 常见约束类型
• 9.3.4 收敛性分析
• 9.3.5 与传统方法的比较

9.4 Extended Position Based Dynamics (XPBD)

• 9.4.1 PBD的局限性
• 9.4.2 拉格朗日乘子的引入
• 9.4.3 柔度参数与物理意义
• 9.4.4 数值稳定性改进
• 9.4.5 能量守恒特性

9.5 材料模型与本构关系

• 9.5.1 线弹性模型
• 9.5.2 超弹性材料
• 9.5.3 粘弹性与塑性
• 9.5.4 各向异性材料
• 9.5.5 损伤与断裂模型

9.6 强化学习应用：软体操作任务

• 9.6.1 软体抓取的挑战
• 9.6.2 仿真精度与速度权衡
• 9.6.3 可微分软体仿真
• 9.6.4 实际案例分析

9.7 历史人物：Terzopoulos与物理动画

• 9.7.1 1987年的开创性工作
• 9.7.2 弹性模型的引入
• 9.7.3 对计算机图形学的影响
• 9.7.4 从动画到仿真的演进

9.8 高级话题：MPM与拓扑变化

• 9.8.1 Material Point Method原理
• 9.8.2 欧拉-拉格朗日耦合
• 9.8.3 大变形处理
• 9.8.4 切割与撕裂仿真

9.9 本章小结

9.10 练习题

9.11 常见陷阱与错误

9.12 最佳实践检查清单

9.1 引言与动机

软体与可变形物体的仿真是物理引擎中最具挑战性的领域之一。与刚体不同，软体具有无限多个自由度，其运动由偏微分方程（PDE）控制，需要在计算效率和物理精度之间找到合适的平衡。在具身智能领域，软体仿真对于理解和控制柔性机器人、软体抓取、布料操作等任务至关重要。

本章将深入探讨三种主流的软体仿真方法：有限元方法（FEM）、基于位置的动力学（PBD）和扩展基于位置的动力学（XPBD）。我们将从连续介质力学的基本原理出发，逐步推导各种方法的数学基础，并讨论它们在强化学习环境中的应用。

9.2 有限元方法（FEM）基础

9.2.1 连续介质力学回顾

在连续介质力学中，物体的运动由以下控制方程描述：

$$\rho \frac{\partial^2 \mathbf{u}}{\partial t^2} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f}$$ 其中，$\rho$ 是材料密度，$\mathbf{u}$ 是位移场，$\boldsymbol{\sigma}$ 是应力张量，$\mathbf{f}$ 是体力。应力与应变的关系由本构方程确定： $$\boldsymbol{\sigma} = \mathcal{C} : \boldsymbol{\varepsilon}$$ 这里 $\mathcal{C}$ 是四阶弹性张量，$\boldsymbol{\varepsilon}$ 是应变张量。对于小变形，格林应变张量可以简化为： $$\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$$

9.2.2 弱形式与变分原理

FEM的核心在于将强形式的PDE转换为弱形式。通过虚功原理，我们得到： $$\int_{\Omega} \delta \boldsymbol{\varepsilon} : \boldsymbol{\sigma} \, dV = \int_{\Omega} \delta \mathbf{u} \cdot \mathbf{f} \, dV + \int_{\Gamma} \delta \mathbf{u} \cdot \mathbf{t} \, dS$$ 其中 $\Omega$ 是物体域，$\Gamma$ 是边界，$\mathbf{t}$ 是表面力，$\delta \mathbf{u}$ 是虚位移。这个弱形式是有限元离散化的基础。

9.2.3 离散化与形函数

将连续域离散为有限个单元，位移场可以表示为： $$\mathbf{u}(\mathbf{x}) = \sum_{i=1}^{n} N_i(\mathbf{x}) \mathbf{u}_i$$ 其中 $N_i$ 是形函数，$\mathbf{u}_i$ 是节点位移。对于四面体单元，常用线性形函数： $$N_i(\xi, \eta, \zeta) = a_i + b_i\xi + c_i\eta + d_i\zeta$$ 形函数的导数用于计算应变： $$\boldsymbol{\varepsilon} = \mathbf{B} \mathbf{u}$$ 其中 $\mathbf{B}$ 是应变-位移矩阵，包含形函数的空间导数。

9.2.4 质量矩阵与刚度矩阵

将离散化的位移场代入弱形式，得到系统的运动方程： $$\mathbf{M} \ddot{\mathbf{u}} + \mathbf{C} \dot{\mathbf{u}} + \mathbf{K} \mathbf{u} = \mathbf{F}$$ 质量矩阵： $$\mathbf{M} = \int_{\Omega} \rho \mathbf{N}^T \mathbf{N} \, dV$$ 刚度矩阵： $$\mathbf{K} = \int_{\Omega} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV$$ 其中 $\mathbf{D}$ 是材料刚度矩阵。对于各向同性线弹性材料： $$\mathbf{D} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix}$$

9.2.5 时间积分策略

对于动态问题，常用的时间积分方法包括：

显式欧拉法： $$\mathbf{u}_{n+1} = \mathbf{u}_n + \Delta t \dot{\mathbf{u}}_n$$ $$\dot{\mathbf{u}}_{n+1} = \dot{\mathbf{u}}_n + \Delta t \mathbf{M}^{-1}(\mathbf{F}_n - \mathbf{K}\mathbf{u}_n)$$ 隐式欧拉法： $$(\mathbf{M} + \Delta t \mathbf{C} + \Delta t^2 \mathbf{K}) \mathbf{u}_{n+1} = \mathbf{M} \mathbf{u}_n + \Delta t \mathbf{M} \dot{\mathbf{u}}_n + \Delta t^2 \mathbf{F}_{n+1}$$ 隐式方法虽然需要求解线性系统，但允许更大的时间步长，适合刚度较大的系统。

9.3 Position Based Dynamics (PBD)

9.3.1 PBD的核心思想

Position Based Dynamics 直接操作顶点位置而非力和加速度，通过投影约束来实现物理行为。这种方法特别适合实时应用，因为它具有无条件稳定性和直观的参数调节。

PBD的基本流程：

1. 预测位置：$\mathbf{p}^* = \mathbf{x}_n + \Delta t \mathbf{v}_n + \Delta t^2 \mathbf{M}^{-1} \mathbf{F}_{ext}$
2. 约束投影：迭代求解约束 $C(\mathbf{p}) = 0$
3. 更新速度：$\mathbf{v}_{n+1} = (\mathbf{p}_{n+1} - \mathbf{x}_n) / \Delta t$
4. 更新位置：$\mathbf{x}_{n+1} = \mathbf{p}_{n+1}$

9.3.2 约束投影算法

对于约束 $C(\mathbf{p}) = 0$，使用一阶泰勒展开： $$C(\mathbf{p} + \Delta \mathbf{p}) \approx C(\mathbf{p}) + \nabla C(\mathbf{p})^T \Delta \mathbf{p} = 0$$ 结合质量加权，位置修正为： $$\Delta \mathbf{p}_i = -\frac{w_i}{\sum_j w_j |\nabla_{p_j} C|^2} \nabla_{p_i} C \cdot C(\mathbf{p})$$ 其中 $w_i = 1/m_i$ 是质量的倒数。

9.3.3 常见约束类型

距离约束： $$C(\mathbf{p}_1, \mathbf{p}_2) = |\mathbf{p}_1 - \mathbf{p}_2| - d_0$$ 梯度： $$\nabla_{p_1} C = \frac{\mathbf{p}_1 - \mathbf{p}_2}{|\mathbf{p}_1 - \mathbf{p}_2|}$$ 体积保持约束： 对于四面体： $$C = \frac{1}{6}[(\mathbf{p}_2 - \mathbf{p}_1) \times (\mathbf{p}_3 - \mathbf{p}_1)] \cdot (\mathbf{p}_4 - \mathbf{p}_1) - V_0$$ 弯曲约束： $$C = \arccos(\mathbf{n}_1 \cdot \mathbf{n}_2) - \theta_0$$ 其中 $\mathbf{n}_1, \mathbf{n}_2$ 是相邻三角形的法向量。

9.3.4 收敛性分析

PBD的收敛性取决于多个因素：

Gauss-Seidel vs Jacobi迭代：

• Gauss-Seidel：串行更新，收敛更快但难以并行化
• Jacobi：并行更新，收敛较慢但适合GPU实现

收敛速度分析：设约束误差为 $e_k = C(\mathbf{p}_k)$，对于线性约束： $$e_{k+1} = (1 - \omega) e_k$$ 其中 $\omega$ 是松弛因子。最优松弛因子： $$\omega_{opt} = \frac{2}{1 + \sqrt{1 - \rho^2}}$$ $\rho$ 是迭代矩阵的谱半径。

多约束耦合问题： 当多个约束作用于同一顶点时，投影顺序影响收敛性。使用约束图着色可以识别独立约束集，实现并行投影：

约束依赖图示例：
    C1 --- v1 --- C2
           |
          C3
           |
          v2

9.3.5 与传统方法的比较

PBD vs FEM：

| 特性 | PBD | FEM |

适用场景：

• PBD：游戏、实时仿真、交互式应用
• FEM：工程分析、科学计算、高精度仿真

9.4 Extended Position Based Dynamics (XPBD)

9.4.1 PBD的局限性

传统PBD存在以下问题：

1. 迭代次数影响物理行为
2. 时间步长影响刚度
3. 缺乏真实的物理参数对应

XPBD通过引入拉格朗日乘子和柔度参数解决这些问题。

9.4.2 拉格朗日乘子的引入

XPBD将约束力显式建模为拉格朗日乘子 $\lambda$： $$\mathbf{M} \frac{\mathbf{x}_{n+1} - 2\mathbf{x}_n + \mathbf{x}_{n-1}}{\Delta t^2} = \mathbf{F}_{ext} + \nabla C^T \lambda$$ 通过隐式积分和线性化，得到更新公式： $$\Delta \lambda = \frac{-C - \tilde{\alpha} \lambda}{\nabla C \mathbf{M}^{-1} \nabla C^T + \tilde{\alpha}}$$ 其中 $\tilde{\alpha} = \alpha / \Delta t^2$ 是正则化的柔度参数。

9.4.3 柔度参数与物理意义

柔度参数 $\alpha$ 与材料属性直接相关：

弹簧约束： $$\alpha = \frac{1}{k}$$ 其中 $k$ 是弹簧刚度。

体积约束（体积模量K）： $$\alpha = \frac{1}{K \cdot V_0}$$ 剪切约束（剪切模量G）： $$\alpha = \frac{1}{G \cdot A_0}$$ 这种对应关系使得XPBD的参数具有明确的物理意义。

9.4.4 数值稳定性改进

XPBD的隐式积分形式提供了更好的数值稳定性：

能量分析： 系统总能量： $$E = \frac{1}{2} \dot{\mathbf{x}}^T \mathbf{M} \dot{\mathbf{x}} + \sum_i \frac{1}{2\alpha_i} C_i^2$$ XPBD保证能量有界： $$E_{n+1} \leq E_n + \Delta t \mathbf{F}_{ext} \cdot \dot{\mathbf{x}}$$ 刚度矩阵条件数： XPBD的有效刚度矩阵： $$\mathbf{K}_{eff} = \sum_i \frac{1}{\alpha_i + \Delta t^2 / m_{eff}} \nabla C_i \nabla C_i^T$$ 条件数受 $\alpha$ 和 $\Delta t$ 共同控制，避免了病态系统。

9.4.5 能量守恒特性

XPBD可以通过适当的积分方案实现能量守恒：

辛积分形式： $$\begin{cases} \mathbf{p}_{n+1/2} = \mathbf{p}_n - \frac{\Delta t}{2} \nabla_q H(\mathbf{q}_n, \mathbf{p}_{n+1/2}) \\ \mathbf{q}_{n+1} = \mathbf{q}_n + \frac{\Delta t}{2} [\nabla_p H(\mathbf{q}_n, \mathbf{p}_{n+1/2}) + \nabla_p H(\mathbf{q}_{n+1}, \mathbf{p}_{n+1/2})] \\ \mathbf{p}_{n+1} = \mathbf{p}_{n+1/2} - \frac{\Delta t}{2} \nabla_q H(\mathbf{q}_{n+1}, \mathbf{p}_{n+1/2}) \end{cases}$$ 这种形式在无耗散情况下严格保持能量守恒。

9.5 材料模型与本构关系

9.5.1 线弹性模型

线弹性是最简单的材料模型，应力与应变成线性关系： $$\boldsymbol{\sigma} = 2\mu \boldsymbol{\varepsilon} + \lambda \text{tr}(\boldsymbol{\varepsilon}) \mathbf{I}$$ 其中 $\mu$ 和 $\lambda$ 是拉梅常数： $$\mu = \frac{E}{2(1+\nu)}, \quad \lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}$$ 应变能密度： $$\Psi = \frac{\lambda}{2}[\text{tr}(\boldsymbol{\varepsilon})]^2 + \mu \text{tr}(\boldsymbol{\varepsilon}^2)$$

9.5.2 超弹性材料

对于大变形，需要使用超弹性模型。常用的Neo-Hookean模型： $$\Psi = \frac{\mu}{2}(I_1 - 3) - \mu \ln J + \frac{\lambda}{2}(\ln J)^2$$ 其中 $I_1 = \text{tr}(\mathbf{F}^T\mathbf{F})$ 是第一不变量，$J = \det(\mathbf{F})$ 是变形梯度的行列式。

第一Piola-Kirchhoff应力： $$\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}} = \mu(\mathbf{F} - \mathbf{F}^{-T}) + \lambda \ln J \mathbf{F}^{-T}$$ Mooney-Rivlin模型： 更复杂的橡胶类材料： $$\Psi = C_1(I_1 - 3) + C_2(I_2 - 3) + \frac{K}{2}(J - 1)^2$$

9.5.3 粘弹性与塑性

标准线性固体模型（SLS）：

弹簧-阻尼器模型：
    k1
---/\/\/\---
           |
         ----
        |    |
        | k2 | c

        | k2 | c
        |    |

         ----
           |

本构关系： $$\sigma + \frac{\eta}{E_2} \dot{\sigma} = E_1 \varepsilon + \eta(1 + \frac{E_1}{E_2})\dot{\varepsilon}$$ Von Mises塑性准则： 屈服条件： $$f = \sqrt{\frac{3}{2}\mathbf{s}:\mathbf{s}} - \sigma_y = 0$$ 其中 $\mathbf{s} = \boldsymbol{\sigma} - \frac{1}{3}\text{tr}(\boldsymbol{\sigma})\mathbf{I}$ 是偏应力张量。

塑性流动法则： $$\dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}}$$

9.5.4 各向异性材料

对于纤维增强材料或生物组织： $$\Psi = \Psi_{iso}(\mathbf{F}) + \sum_{i=1}^{N} \Psi_{fiber}^{(i)}(I_4^{(i)})$$ 其中 $I_4^{(i)} = \mathbf{a}_0^{(i)} \cdot (\mathbf{C} \mathbf{a}_0^{(i)})$ 是纤维方向的伸长，$\mathbf{a}_0^{(i)}$ 是第i组纤维的初始方向。

横观各向同性： 刚度矩阵具有特殊结构： $$\mathbf{D} = \begin{bmatrix} E_L & \nu_{LT}E_L & \nu_{LT}E_L & 0 & 0 & 0 \\ \nu_{LT}E_L & E_T & \nu_{TT}E_T & 0 & 0 & 0 \\ \nu_{LT}E_L & \nu_{TT}E_T & E_T & 0 & 0 & 0 \\ 0 & 0 & 0 & G_{LT} & 0 & 0 \\ 0 & 0 & 0 & 0 & G_{LT} & 0 \\ 0 & 0 & 0 & 0 & 0 & G_{TT} \end{bmatrix}$$

9.5.5 损伤与断裂模型

连续损伤力学： 引入损伤变量 $D \in [0,1]$： $$\boldsymbol{\sigma}_{eff} = (1-D)\boldsymbol{\sigma}$$ 损伤演化方程： $$\dot{D} = \begin{cases} \frac{1}{\eta} \langle \frac{Y - Y_0}{Y_c - Y_0} \rangle & \text{if } Y > Y_0 \\ 0 & \text{otherwise} \end{cases}$$ 其中 $Y$ 是能量释放率。

相场断裂模型： 总能量泛函： $$\mathcal{E} = \int_\Omega g(d) \psi^+(\boldsymbol{\varepsilon}) + \psi^-(\boldsymbol{\varepsilon}) \, dV + \int_\Omega G_c \left( \frac{d^2}{2l} + \frac{l}{2}|\nabla d|^2 \right) dV$$ 其中 $d$ 是相场变量，$G_c$ 是临界能量释放率，$l$ 是长度尺度参数。

9.6 强化学习应用：软体操作任务

9.6.1 软体抓取的挑战

软体操作在强化学习中面临独特挑战：

状态空间复杂性：

• 无限维度：软体具有连续变形场
• 部分可观测：内部应力状态难以直接感知
• 非线性动力学：接触和大变形导致强非线性

奖励设计困难：

# 软体抓取的多目标奖励函数示例
reward = w1 * grasp_stability 

       + w2 * minimal_deformation
       - w3 * stress_concentration
       - w4 * energy_consumption

仿真到现实差距：

• 材料参数不确定性（杨氏模量变化±30%）
• 接触模型简化（实际摩擦系数的空间变化）
• 计算精度限制（网格分辨率vs计算时间）

9.6.2 仿真精度与速度权衡

多分辨率策略：

| 方法 | 节点数 | 时间步(ms) | 实时因子 | 精度 |

自适应精度控制：

训练阶段：
  早期探索 → 低精度PBD (10x实时)
  策略细化 → XPBD (4x实时)
  最终优化 → FEM验证 (0.5x实时)

9.6.3 可微分软体仿真

可微分仿真允许通过梯度优化直接学习控制策略：

隐式微分方法： 对于平衡方程 $\mathbf{f}(\mathbf{x}, \mathbf{u}) = 0$： $$\frac{d\mathbf{x}}{d\mathbf{u}} = -\left(\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\right)^{-1} \frac{\partial \mathbf{f}}{\partial \mathbf{u}}$$ 伴随法计算梯度： 损失函数 $L = \ell(\mathbf{x}_T)$ 的梯度： $$\frac{dL}{d\mathbf{u}} = \sum_{t=0}^{T-1} \left(\frac{\partial \mathbf{f}_t}{\partial \mathbf{u}}\right)^T \boldsymbol{\lambda}_t$$ 其中伴随变量 $\boldsymbol{\lambda}_t$ 通过反向递推求解： $$\boldsymbol{\lambda}_{t-1} = \left(\frac{\partial \mathbf{f}_t}{\partial \mathbf{x}_t}\right)^T \boldsymbol{\lambda}_t$$

9.6.4 实际案例分析

案例1：软体夹持器设计优化

目标：优化气动软体夹持器的形状和控制策略

仿真设置：

• 材料：硅橡胶（Neo-Hookean, E=1MPa, ν=0.45）
• 网格：2000个四面体单元
• 控制：3个独立气腔压力

学习结果：

• 训练时间：48小时（10万episodes）
• 成功率提升：65% → 92%
• 关键发现：非对称压力模式提高适应性

案例2：布料折叠任务

挑战：学习复杂的布料操作序列

技术选择：

• PBD用于快速原型
• XPBD用于最终策略训练
• 约束：10000个距离约束 + 5000个弯曲约束

策略架构：

观测 → CNN特征提取 → LSTM序列建模 → 动作预测
(深度图像)  (ResNet-18)   (256隐藏单元)  (末端执行器轨迹)

9.7 历史人物：Terzopoulos与物理动画

9.7.1 1987年的开创性工作

Demetri Terzopoulos在1987年发表的论文"Elastically Deformable Models"开创了基于物理的计算机动画时代。这项工作首次将连续介质力学引入计算机图形学，使得虚拟物体能够像真实材料一样变形。

核心贡献：

1. 将拉格朗日动力学应用于图形学
2. 提出了弹性能量最小化框架
3. 开发了高效的多重网格求解器

9.7.2 弹性模型的引入

Terzopoulos的弹性模型基于变分原理： $$\mathcal{L} = \int_\Omega \left[ \frac{\rho}{2} \left| \frac{\partial \mathbf{r}}{\partial t} \right|^2 - \mathcal{E}(\mathbf{r}) \right] dV$$ 其中 $\mathcal{E}$ 是弹性势能： $$\mathcal{E} = \int_\Omega \left[ \alpha |\boldsymbol{\varepsilon}|^2 + \beta |\nabla^2 \mathbf{r}|^2 \right] dV$$ 这个公式结合了拉伸能（$\alpha$项）和弯曲能（$\beta$项），为后续的物理仿真奠定了基础。

9.7.3 对计算机图形学的影响

Terzopoulos的工作产生了深远影响：

直接影响：

• SIGGRAPH社区开始关注物理仿真
• 催生了物理动画子领域
• 启发了后续的布料、流体、头发仿真研究

技术演进时间线：

1987: 弹性模型
1988: 粘弹性和塑性扩展
1992: 自适应网格细化
1995: 实时物理仿真起步
1998: 大规模并行计算
2005: GPU加速成为主流
2012: 深度学习与物理结合
2020: 可微分仿真兴起

9.7.4 从动画到仿真的演进

从纯视觉效果到物理准确仿真的转变：

早期（1987-1995）：

• 目标：视觉真实感
• 精度要求：定性正确
• 应用：电影特效、游戏

中期（1995-2010）：

• 目标：交互式仿真
• 精度要求：实时稳定
• 应用：虚拟手术、训练系统

现代（2010-至今）：

• 目标：预测性仿真
• 精度要求：定量准确
• 应用：机器人学习、数字孪生

9.8 高级话题：MPM与拓扑变化

9.8.1 Material Point Method原理

MPM结合了拉格朗日粒子和欧拉网格的优点：

基本流程：

1. 粒子携带材料信息（质量、动量、应变）
2. 将粒子信息映射到背景网格
3. 在网格上求解动量方程
4. 将网格解映射回粒子
5. 更新粒子位置和变形梯度

动量方程： $$m_i \mathbf{v}_i^{n+1} = m_i \mathbf{v}_i^n + \Delta t \sum_p \left[ -\frac{\partial \Psi_p}{\partial \mathbf{F}_p} \mathbf{F}_p^T V_p \nabla w_{ip} + m_p \mathbf{g} w_{ip} \right]$$

9.8.2 欧拉-拉格朗日耦合

传输方程： 粒子到网格（P2G）： $$m_i = \sum_p m_p w_{ip}$$ $$m_i \mathbf{v}_i = \sum_p m_p \mathbf{v}_p w_{ip}$$ 网格到粒子（G2P）： $$\mathbf{v}_p^{n+1} = \sum_i \mathbf{v}_i^{n+1} w_{ip}$$ $$\mathbf{F}_p^{n+1} = \left( \mathbf{I} + \Delta t \sum_i \mathbf{v}_i^{n+1} \nabla w_{ip} \right) \mathbf{F}_p^n$$

9.8.3 大变形处理

MPM自然处理大变形和自碰撞：

变形梯度更新： $$\mathbf{F}_{n+1} = \mathbf{F}_{n+1}^E \mathbf{F}_{n+1}^P$$ 塑性投影（Snow plasticity）： $$\mathbf{F} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^T$$ $$\boldsymbol{\Sigma}_{clamped} = \text{clamp}(\boldsymbol{\Sigma}, [1-\theta_c, 1+\theta_s])$$

9.8.4 切割与撕裂仿真

损伤追踪： 每个粒子维护损伤变量 $d_p \in [0,1]$： $$d_p^{n+1} = \min(1, d_p^n + \Delta d)$$ $$\Delta d = H(|\boldsymbol{\sigma}| - \sigma_{crit}) \cdot k_{damage} \cdot \Delta t$$ 拓扑变化处理： 当 $d_p > d_{threshold}$ 时：

1. 标记粒子为"断裂"
2. 修改权重函数以断开连接
3. 生成新的自由表面

9.9 本章小结

本章深入探讨了软体与可变形物体的仿真方法，涵盖了从传统有限元到现代基于位置动力学的各种技术。关键要点包括：

核心概念：

1. FEM提供了物理准确但计算密集的仿真框架
2. PBD通过直接位置操作实现稳定的实时仿真
3. XPBD改进了PBD的物理一致性和参数化
4. 材料模型决定了仿真的物理行为
5. MPM处理大变形和拓扑变化

关键公式汇总：

| 方法 | 核心方程 | 计算复杂度 |

实践指南：

• 选择方法时考虑精度vs速度需求
• 使用域随机化处理参数不确定性
• 实现多分辨率策略优化训练效率
• 验证仿真结果的物理合理性

9.10 练习题

基础题

练习9.1：形函数推导 对于一个二维三角形单元，顶点坐标为 $(x_1, y_1)$、$(x_2, y_2)$、$(x_3, y_3)$。推导线性形函数 $N_1$、$N_2$、$N_3$ 的表达式，并验证它们满足：

• $N_i(x_j, y_j) = \delta_{ij}$
• $\sum_{i=1}^3 N_i = 1$

提示：使用面积坐标方法

练习9.2：PBD距离约束 两个质点质量分别为 $m_1 = 2kg$ 和 $m_2 = 3kg$，当前位置为 $\mathbf{p}_1 = (1, 0, 0)$ 和 $\mathbf{p}_2 = (2.5, 0, 0)$，静止长度 $d_0 = 1m$。计算一次约束投影后的新位置。

提示：使用质量加权的位置修正公式

练习9.3：杨氏模量与泊松比转换 给定材料的杨氏模量 $E = 100 MPa$，泊松比 $\nu = 0.3$，计算： a) 拉梅常数 $\lambda$ 和 $\mu$ b) 体积模量 $K$ 和剪切模量 $G$ c) 波速 $c_p$（纵波）和 $c_s$（横波），假设密度 $\rho = 1000 kg/m^3$

提示：使用弹性常数之间的转换关系

挑战题

练习9.4：XPBD柔度参数设计 设计一个XPBD弹簧系统，要求：

• 弹簧刚度相当于 $k = 1000 N/m$
• 时间步长 $\Delta t = 0.01s$
• 系统在10Hz正弦激励下的稳态振幅放大因子不超过2

推导柔度参数 $\alpha$ 的取值范围，并分析数值阻尼的影响。

提示：考虑XPBD的频率响应特性

练习9.5：Neo-Hookean能量最小化 考虑一个Neo-Hookean材料的单轴拉伸问题，变形梯度为： $$\mathbf{F} = \begin{bmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \lambda_3 \end{bmatrix}$$ 在不可压缩约束 $\det(\mathbf{F}) = 1$ 和单轴应力条件下，推导主伸长 $\lambda_1$ 与名义应力 $P_{11}$ 的关系。

提示：利用不可压缩性得到 $\lambda_2 = \lambda_3 = 1/\sqrt{\lambda_1}$

练习9.6：MPM稳定性分析 分析MPM方法的CFL条件。给定：

• 网格尺寸 $\Delta x = 0.01m$
• 材料密度 $\rho = 1000 kg/m^3$
• 杨氏模量 $E = 1 MPa$
• 泊松比 $\nu = 0.3$

计算最大允许时间步长，并讨论APIC和CPIC方法对稳定性的改进。

提示：考虑弹性波传播速度

练习9.7：开放性思考题 设计一个混合仿真框架，结合FEM的精度和PBD的稳定性，用于机器人抓取软体食品（如豆腐）的强化学习训练。讨论：

1. 如何划分FEM和PBD区域？
2. 界面耦合如何处理？
3. 如何自适应切换方法？
4. 训练效率与仿真精度的权衡策略？

提示：考虑接触区域的重要性和计算资源分配

9.11 常见陷阱与错误

FEM相关陷阱

1. 锁定现象（Locking）

问题：不可压缩材料使用低阶单元时出现过度刚化
症状：变形远小于预期，应力分布不均匀

解决方案：

• 使用选择性积分或减缩积分
• 采用混合有限元（u-p公式）
• 使用高阶单元或B-bar方法

2. 沙漏模式（Hourglassing）

问题：减缩积分导致零能量模式
症状：网格出现非物理的锯齿状变形

预防措施：

• 添加沙漏控制（刚度或粘性）
• 使用完全积分单元
• 细化网格或使用稳定化技术

3. 负雅可比（Negative Jacobian）

错误信息：Element Jacobian determinant < 0
原因：单元过度扭曲或翻转

调试步骤：

1. 检查初始网格质量
2. 减小时间步长
3. 使用重网格化技术
4. 改进约束处理

PBD/XPBD陷阱

4. 约束顺序依赖

# 错误：结果依赖约束求解顺序
for constraint in constraints:
    project(constraint)  # 顺序影响收敛

# 正确：使用图着色并行化
parallel_sets = color_constraints(constraints)
for set in parallel_sets:
    parallel_project(set)

5. 刚度与迭代次数耦合

陷阱：增加迭代次数意外改变了材料刚度
表现：调试时行为不一致

解决：

• 迁移到XPBD（物理参数独立于迭代）
• 标准化迭代次数
• 使用收敛准则而非固定迭代

6. 质量矩阵错误

# 常见错误：忘记考虑密度变化
mass = volume * density  # 错误：使用初始体积

# 正确：更新质量或使用守恒形式
mass = volume0 * density0  # 质量守恒

材料模型陷阱

7. 能量不一致

问题：应力和能量函数不匹配
症状：系统能量异常增长或衰减

检查清单：

• 确保 $\boldsymbol{\sigma} = \frac{\partial \Psi}{\partial \boldsymbol{\varepsilon}}$
• 验证切线刚度矩阵
• 测试能量守恒性

8. 大变形下的线性假设

# 错误：大变形使用线性应变
strain = 0.5 * (grad_u + grad_u.T)  # 仅适用于小变形

# 正确：使用格林应变
F = I + grad_u
strain = 0.5 * (F.T @ F - I)  # 格林-拉格朗日应变

数值积分陷阱

9. CFL违背

症状：仿真爆炸或产生NaN
诊断：监控最大速度和加速度

# 自适应时间步
dt_cfl = cfl_number * dx / max_velocity
dt = min(dt_target, dt_cfl)

10. 能量漂移

问题：长时间仿真中总能量持续增长/衰减
原因：数值积分误差累积

缓解策略：

• 使用辛积分器
• 定期能量修正
• 监控并记录能量偏差

调试技巧

能量监控器实现：

class EnergyMonitor:
    def __init__(self):
        self.history = []

    def compute_energy(self, system):
        KE = 0.5 * np.sum(system.mass * system.velocity**2)
        PE = system.compute_potential_energy()
        return KE + PE

    def check_conservation(self, tolerance=0.01):
        if len(self.history) > 1:
            drift = abs(self.history[-1] - self.history[0])
            relative_drift = drift / abs(self.history[0])
            if relative_drift > tolerance:
                warnings.warn(f"Energy drift: {relative_drift:.2%}")

约束违背检测：

def validate_constraints(positions, constraints, tolerance=1e-6):
    violations = []
    for constraint in constraints:
        error = constraint.evaluate(positions)
        if abs(error) > tolerance:
            violations.append({
                'constraint': constraint,
                'error': error,
                'relative': error / constraint.rest_value
            })
    return violations

9.12 最佳实践检查清单

设计阶段

方法选择

• [ ] 明确精度要求和实时性约束
• [ ] 评估可用计算资源（CPU/GPU）
• [ ] 考虑并行化需求
• [ ] 确定材料参数范围
• [ ] 评估预期变形程度

网格设计

• [ ] 网格质量检查（长宽比、扭曲度）
• [ ] 分辨率与特征尺寸匹配
• [ ] 边界层网格加密
• [ ] 考虑自适应细化策略
• [ ] 验证网格收敛性

实现阶段

数据结构

• [ ] 选择合适的稀疏矩阵格式
• [ ] 实现高效的邻居搜索
• [ ] 优化内存布局（缓存友好）
• [ ] 预分配内存避免动态分配
• [ ] 实现对象池减少GC压力

算法实现

• [ ] 单元测试每个约束类型
• [ ] 验证雅可比矩阵计算
• [ ] 实现数值微分验证解析导数
• [ ] 添加断言检查不变量
• [ ] 实现回退机制处理失败

验证阶段

物理验证

• [ ] 标准测试案例（悬臂梁、拉伸试验）
• [ ] 能量守恒检查
• [ ] 动量守恒验证
• [ ] 对称性测试
• [ ] 极限情况分析

数值验证

• [ ] 网格收敛性研究
• [ ] 时间步长敏感性分析
• [ ] 条件数监控
• [ ] 残差收敛曲线
• [ ] 与解析解对比

优化阶段

性能优化

• [ ] 性能剖析识别瓶颈
• [ ] SIMD向量化关键循环
• [ ] GPU核函数优化
• [ ] 减少内存分配
• [ ] 实现多级并行

鲁棒性增强

• [ ] 添加自适应时间步长
• [ ] 实现碰撞检测后备方案
• [ ] 处理退化情况
• [ ] 添加数值阻尼选项
• [ ] 实现自动参数调节

集成阶段

RL训练集成

• [ ] 定义清晰的观测空间
• [ ] 设计可微分奖励函数
• [ ] 实现高效重置机制
• [ ] 支持批量仿真
• [ ] 提供梯度信息接口

监控与调试

• [ ] 实时可视化工具
• [ ] 性能指标仪表板
• [ ] 日志关键事件
• [ ] 支持确定性重放
• [ ] 导出调试快照

部署阶段

文档完善

• [ ] API文档完整
• [ ] 参数说明详细
• [ ] 提供使用示例
• [ ] 性能基准数据
• [ ] 已知限制说明

维护准备

• [ ] 版本控制策略
• [ ] 自动化测试套件
• [ ] 性能回归测试
• [ ] 用户反馈机制
• [ ] 更新升级路径

具体配置建议

FEM配置模板：

fem_config:
  element_type: "tet10"  # 10节点四面体
  material:
    model: "neo_hookean"
    youngs_modulus: 1.0e6  # Pa
    poisson_ratio: 0.45
  solver:
    type: "conjugate_gradient"
    preconditioner: "incomplete_cholesky"
    tolerance: 1.0e-6
    max_iterations: 1000
  time_integration:
    method: "implicit_euler"
    timestep: 0.001
    adaptive: true

PBD/XPBD配置模板：

pbd_config:
  solver_iterations: 10
  substeps: 4
  constraints:

    - type: "distance"
      stiffness: 0.9

    - type: "volume"
      stiffness: 0.95

    - type: "bending"
      stiffness: 0.1
  collision:
    margin: 0.01
    friction: 0.3
    restitution: 0.1

通过遵循这份检查清单，可以显著提高软体仿真系统的质量、性能和可维护性。记住，最佳实践需要根据具体应用场景进行调整。