第2章：拉格朗日视角（上）

在物理模拟中，我们有两种基本的观察视角：拉格朗日视角和欧拉视角。拉格朗日视角追踪每个质点在空间中的运动轨迹，就像跟随一叶扁舟顺流而下；而欧拉视角则固定在空间中的某些点，观察物质如何流经这些点。本章将深入探讨拉格朗日方法，从最简单的弹簧质点系统开始，逐步构建起布料、流体等复杂系统的模拟框架。通过本章学习，你将掌握粒子系统的核心算法、时间积分的稳定性分析，以及高效的空间数据结构。

学习目标

完成本章学习后，你将能够：

1. 理解拉格朗日力学框架：掌握质点系统的数学描述和力的计算方法
2. 实现弹簧质点系统：构建可扩展的质点-弹簧框架，支持各种约束和外力
3. 分析数值稳定性：理解不同时间积分器的稳定性条件，选择合适的积分方案
4. 掌握SPH方法：理解光滑粒子流体动力学的数学基础和实现要点
5. 优化邻居搜索：实现高效的空间数据结构，加速粒子间相互作用计算

2.1 弹簧质点系统：第一个物理模拟器

弹簧质点系统是物理引擎的"Hello World"。尽管简单，它却包含了动力学模拟的所有核心要素：状态表示、力的计算、时间积分和约束处理。

2.1.1 系统状态表示

对于包含 $N$ 个质点的系统，我们需要存储：

• 位置向量：$\mathbf{x}_i \in \mathbb{R}^3$，$i = 1, 2, ..., N$
• 速度向量：$\mathbf{v}_i \in \mathbb{R}^3$
• 质量：$m_i > 0$

系统的完整状态可以表示为： $$\mathbf{q} = [\mathbf{x}_1^T, \mathbf{x}_2^T, ..., \mathbf{x}_N^T]^T \in \mathbb{R}^{3N}$$ $$\dot{\mathbf{q}} = [\mathbf{v}_1^T, \mathbf{v}_2^T, ..., \mathbf{v}_N^T]^T \in \mathbb{R}^{3N}$$

2.1.2 力的计算

每个质点受到的总力 $\mathbf{f}_i$ 包含多个分量：

弹簧力（Hooke定律）

对于连接质点 $i$ 和 $j$ 的弹簧： $$\mathbf{f}_{ij}^{spring} = k_s \left( |\mathbf{x}_j - \mathbf{x}_i| - l_0 \right) \frac{\mathbf{x}_j - \mathbf{x}_i}{|\mathbf{x}_j - \mathbf{x}_i|}$$ 其中：

• $k_s$ 是弹簧刚度系数
• $l_0$ 是弹簧原长
• $|\mathbf{x}_j - \mathbf{x}_i|$ 是当前长度

注意这里的符号约定：当弹簧被拉伸时（$|\mathbf{x}_j - \mathbf{x}_i| > l_0$），力的方向指向对方质点。

数值稳定性考虑：当两质点重合时（$|\mathbf{x}_j - \mathbf{x}_i| \approx 0$），需要特殊处理避免除零： $$\mathbf{f}_{ij}^{spring} = \begin{cases} k_s \left( |\mathbf{d}| - l_0 \right) \frac{\mathbf{d}}{|\mathbf{d}|} & \text{if } |\mathbf{d}| > \epsilon \\ \mathbf{0} & \text{otherwise} \end{cases}$$ 其中 $\mathbf{d} = \mathbf{x}_j - \mathbf{x}_i$，$\epsilon \approx 10^{-6}$ 是数值容差。

力的雅可比矩阵（用于隐式积分）：

弹簧力对位置的导数是实现隐式积分的关键。对于弹簧力 $\mathbf{f}_{ij}^{spring}$： $$\frac{\partial \mathbf{f}_{ij}}{\partial \mathbf{x}_i} = -k_s \left[ \left(1 - \frac{l_0}{|\mathbf{d}|}\right)\mathbf{I} + \frac{l_0}{|\mathbf{d}|^3}\mathbf{d}\mathbf{d}^T \right]$$

$$\frac{\partial \mathbf{f}_{ij}}{\partial \mathbf{x}_j} = -\frac{\partial \mathbf{f}_{ij}}{\partial \mathbf{x}_i}$$ 其中 $\mathbf{I}$ 是单位矩阵。这个雅可比矩阵具有以下性质：

• 对称性：$\mathbf{J}_{ij} = \mathbf{J}_{ji}^T$
• 正定性：当弹簧被拉伸时，矩阵是正定的

非线性弹簧模型：

真实材料往往表现出非线性弹性行为。几种常见的扩展：

1. 分段线性模型： $$k(l) = \begin{cases} k_1 & \text{if } l < l_1 \\ k_2 & \text{if } l_1 \leq l < l_2 \\ k_3 & \text{if } l \geq l_2 \end{cases}$$

2. 多项式模型（Neo-Hookean类似）： $$\mathbf{f}_{ij}^{spring} = k_s \left( \frac{l}{l_0} - \frac{l_0}{l} \right) \frac{\mathbf{d}}{l}$$

3. 指数模型（用于生物组织）： $$\mathbf{f}_{ij}^{spring} = k_s \frac{l_0}{l} \left( e^{\alpha(l/l_0 - 1)} - 1 \right) \mathbf{d}$$ 张力限制：

为防止过度拉伸，可以设置最大张力： $$\mathbf{f}_{ij}^{spring} = \min\left(|\mathbf{f}_{ij}^{spring}|, f_{max}\right) \cdot \frac{\mathbf{f}_{ij}^{spring}}{|\mathbf{f}_{ij}^{spring}|}$$

阻尼力

为了模拟能量耗散，我们添加阻尼力： $$\mathbf{f}_{ij}^{damp} = -k_d \frac{(\mathbf{v}_j - \mathbf{v}_i) \cdot (\mathbf{x}_j - \mathbf{x}_i)}{|\mathbf{x}_j - \mathbf{x}_i|^2} (\mathbf{x}_j - \mathbf{x}_i)$$ 这里 $k_d$ 是阻尼系数。注意这是投影到弹簧方向的相对速度阻尼，避免了旋转运动的过度衰减。

物理解释：这种阻尼模型只衰减沿弹簧方向的相对运动，保持了系统的旋转不变性。相比于简单的速度阻尼 $\mathbf{f} = -k_d \mathbf{v}$，它不会错误地消除刚体旋转。

临界阻尼条件：对于单个弹簧-质点系统，临界阻尼系数为： $$k_d^{critical} = 2\sqrt{k_s \cdot m_{reduced}}$$ 其中 $m_{reduced} = \frac{m_i m_j}{m_i + m_j}$ 是约化质量。

外力

常见的外力包括：

• 重力：$\mathbf{f}_i^{gravity} = m_i \mathbf{g}$，其中 $\mathbf{g} = [0, -9.8, 0]^T$
• 风力：$\mathbf{f}_i^{wind} = c_w A_i (\mathbf{v}_{wind} - \mathbf{v}_i)$
• 用户交互力：鼠标拖拽等

高级外力模型：

1. 空气阻力（考虑雷诺数）： $$\mathbf{f}_{drag} = -\frac{1}{2} \rho_{air} C_d A |\mathbf{v}_{rel}| \mathbf{v}_{rel}$$ 其中 $C_d$ 是阻力系数，$A$ 是迎风面积，$\mathbf{v}_{rel}$ 是相对空气的速度。

2. 磁力/电力（用于特殊效果）： $$\mathbf{f}_{field} = q(\mathbf{E} + \mathbf{v} \times \mathbf{B})$$

3. 涡旋力（增强真实感）： $$\mathbf{f}_{vortex} = \epsilon m (\mathbf{\omega} \times \mathbf{v})$$ 其中 $\mathbf{\omega}$ 是局部涡度，$\epsilon$ 是涡旋强度参数。

2.1.3 运动方程

根据牛顿第二定律，质点 $i$ 的运动方程为： $$m_i \ddot{\mathbf{x}}_i = \sum_j \mathbf{f}_{ij} + \mathbf{f}_i^{ext}$$ 整个系统可以写成矩阵形式： $$\mathbf{M} \ddot{\mathbf{q}} = \mathbf{f}(\mathbf{q}, \dot{\mathbf{q}}, t)$$ 其中 $\mathbf{M} = \text{diag}(m_1, m_1, m_1, m_2, ..., m_N, m_N, m_N)$ 是质量矩阵。

拉格朗日形式：

从变分原理出发，系统的拉格朗日量： $$L = T - V = \frac{1}{2}\dot{\mathbf{q}}^T\mathbf{M}\dot{\mathbf{q}} - V(\mathbf{q})$$ Euler-Lagrange方程： $$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\mathbf{q}}}\right) - \frac{\partial L}{\partial \mathbf{q}} = \mathbf{0}$$ 导出： $$\mathbf{M}\ddot{\mathbf{q}} = -\frac{\partial V}{\partial \mathbf{q}} = \mathbf{f}_{conservative}$$ 哈密顿形式：

引入广义动量 $\mathbf{p} = \mathbf{M}\dot{\mathbf{q}}$，哈密顿量： $$H = \frac{1}{2}\mathbf{p}^T\mathbf{M}^{-1}\mathbf{p} + V(\mathbf{q})$$ 正则方程： $$\dot{\mathbf{q}} = \frac{\partial H}{\partial \mathbf{p}} = \mathbf{M}^{-1}\mathbf{p}$$ $$\dot{\mathbf{p}} = -\frac{\partial H}{\partial \mathbf{q}} = -\frac{\partial V}{\partial \mathbf{q}}$$ 这种形式特别适合辛积分器的设计。

约束动力学：

当系统包含约束时（如固定长度的刚性杆），需要引入约束方程： $$\mathbf{g}(\mathbf{q}) = \mathbf{0}$$ 使用拉格朗日乘子法，增广的运动方程： $$\mathbf{M}\ddot{\mathbf{q}} = \mathbf{f} + \mathbf{G}^T\boldsymbol{\lambda}$$ $$\mathbf{g}(\mathbf{q}) = \mathbf{0}$$ 其中 $\mathbf{G} = \frac{\partial \mathbf{g}}{\partial \mathbf{q}}$ 是约束雅可比矩阵。

例子：双摆系统：

考虑由两个质点和两根刚性杆组成的双摆：

• 质点质量：$m_1, m_2$
• 杆长：$l_1, l_2$
• 角度：$\theta_1, \theta_2$

拉格朗日量： $$L = \frac{1}{2}(m_1 + m_2)l_1^2\dot{\theta}_1^2 + \frac{1}{2}m_2l_2^2\dot{\theta}_2^2 + m_2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1-\theta_2)$$ $$- (m_1 + m_2)gl_1\cos\theta_1 - m_2gl_2\cos\theta_2$$ 这导出了著名的双摆混沌方程组。

2.1.4 边界条件处理

实际系统中常需要固定某些质点。处理方法包括：

1. 投影法：每步积分后将固定点的位置和速度重置

if (particle[i].is_fixed) {
    particle[i].position = fixed_position[i];
    particle[i].velocity = vec3(0);
}

优点：简单直接；缺点：破坏了系统的对称性

2. 拉格朗日乘子法：添加约束力确保 $\mathbf{x}_i = \mathbf{x}_i^{fixed}$

增广系统： $$\begin{bmatrix} \mathbf{M} & \mathbf{J}^T \\ \mathbf{J} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{q}} \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{f} \\ \mathbf{b} \end{bmatrix}$$ 其中 $\mathbf{J}$ 是约束雅可比矩阵，$\boldsymbol{\lambda}$ 是拉格朗日乘子。

3. 罚函数法：添加大刚度弹簧将质点"钉"在固定位置 $$\mathbf{f}_{penalty} = -k_{penalty}(\mathbf{x}_i - \mathbf{x}_i^{fixed}) - k_{damp}\mathbf{v}_i$$ 其中 $k_{penalty} \gg k_s$ 是惩罚刚度。

4. 消去法：直接从系统矩阵中移除固定自由度

将系统分块： $$\begin{bmatrix} \mathbf{M}_{ff} & \mathbf{M}_{fc} \\ \mathbf{M}_{cf} & \mathbf{M}_{cc} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{q}}_f \\ \ddot{\mathbf{q}}_c \end{bmatrix} = \begin{bmatrix} \mathbf{f}_f \\ \mathbf{f}_c \end{bmatrix}$$ 其中下标 $f$ 表示自由节点，$c$ 表示约束节点。只求解自由部分： $$\mathbf{M}_{ff} \ddot{\mathbf{q}}_f = \mathbf{f}_f - \mathbf{M}_{fc} \ddot{\mathbf{q}}_c$$

2.1.5 数据结构设计

高效的实现需要合理的数据结构：

struct Particle {
    vec3 position;
    vec3 velocity;
    float mass;
    bool is_fixed;
};

struct Spring {
    int particle_i, particle_j;
    float rest_length;
    float stiffness;
    float damping;
};

class MassSpringSystem {
    vector<Particle> particles;
    vector<Spring> springs;
    // 邻接表加速力计算
    vector<vector<int>> spring_neighbors;
};

内存布局优化：

1. SoA vs AoS（Structure of Arrays vs Array of Structures）：

// AoS - 缓存友好的单粒子访问
struct Particle { vec3 pos, vel; float mass; };
Particle particles[N];

// SoA - SIMD友好的批量计算
struct Particles {
    float pos_x[N], pos_y[N], pos_z[N];
    float vel_x[N], vel_y[N], vel_z[N];
    float mass[N];
};

2. 稀疏矩阵存储：对于大规模系统，使用CSR（Compressed Sparse Row）格式存储连接关系：

struct SparseConnectivity {
    vector<int> row_ptr;     // 每行的起始索引
    vector<int> col_idx;     // 列索引
    vector<Spring> springs;  // 对应的弹簧数据
};

3. 缓存对齐：确保热点数据结构对齐到缓存行边界（通常64字节）：

struct alignas(64) ParticleBlock {
    vec3 positions[16];
    vec3 velocities[16];
};

高级数据结构设计：

4. 层次化表示：

struct HierarchicalMassSpring {
    // 细节层级
    vector<MassSpringSystem> levels;
    // 层级间映射
    vector<InterLevelMapping> mappings;
    // 自适应细化准则
    function<bool(const Spring&)> refinement_criterion;
};

5. 时间相关缓存：

struct TemporalCache {
    // 上一步的力（用于高阶积分）
    vector<vec3> forces_prev;
    // 预测位置（用于碰撞检测）
    vector<vec3> positions_predicted;
    // 雅可比矩阵（可重用）
    SparseMatrix jacobian;
};

6. 并行友好的分区：

struct ParallelPartition {
    vector<int> particle_to_partition;
    vector<vector<int>> partition_particles;
    vector<vector<int>> boundary_particles;
    vector<SpinLock> partition_locks;
};

图着色优化：

为了并行计算弹簧力而不产生竞争条件，使用图着色算法：

struct GraphColoring {
    int num_colors;
    vector<int> spring_colors;
    vector<vector<int>> color_groups;

    void compute_coloring(const vector<Spring>& springs) {
        // 贪心着色算法
        // 确保同色弹簧不共享质点
    }
};

内存池管理：

template<typename T>
class MemoryPool {
    vector<T*> blocks;
    size_t block_size;
    size_t current_block;
    size_t current_offset;

public:
    T* allocate(size_t n) {
        // 从预分配的块中分配
        // 避免频繁的malloc/free
    }
};

SIMD友好的批处理：

struct SimdBatch {
    static constexpr int BATCH_SIZE = 8;

    // 8个弹簧的数据打包在一起
    __m256 rest_lengths;
    __m256 stiffnesses;
    __m256i particle_indices_i;
    __m256i particle_indices_j;

    void compute_forces(__m256* positions_x, 
                       __m256* positions_y,
                       __m256* positions_z,
                       __m256* forces_x,
                       __m256* forces_y,
                       __m256* forces_z);
};

2.1.6 能量分析与守恒

物理系统的能量分析对于验证实现正确性和选择合适的积分器至关重要。

系统总能量

弹簧质点系统的总能量包括：

1. 动能： $$T = \frac{1}{2} \sum_{i=1}^N m_i |\mathbf{v}_i|^2$$

2. 弹性势能： $$V_{elastic} = \frac{1}{2} \sum_{\text{springs}} k_s (|\mathbf{x}_j - \mathbf{x}_i| - l_0)^2$$

3. 重力势能： $$V_{gravity} = \sum_{i=1}^N m_i \mathbf{g} \cdot \mathbf{x}_i$$ 总机械能：$E = T + V_{elastic} + V_{gravity}$

能量守恒检验

理想情况下（无阻尼），系统应该保持能量守恒。实际计算中： $$\frac{dE}{dt} = -P_{damping} - P_{numerical}$$ 其中：

• $P_{damping}$ 是阻尼功率
• $P_{numerical}$ 是数值误差导致的能量变化

监测指标：

• 相对能量误差：$\epsilon_E = \frac{|E(t) - E(0)|}{E(0)}$
• 能量漂移率：$\dot{E}_{drift} = \frac{1}{t} \ln\frac{E(t)}{E(0)}$

阻尼功率的精确计算：

对于速度相关的阻尼力，瞬时功率： $$P_{damping} = -\sum_{i,j} \mathbf{f}_{ij}^{damp} \cdot (\mathbf{v}_j - \mathbf{v}_i)$$ 总耗散能量： $$E_{dissipated} = \int_0^t P_{damping} dt$$ 辛积分器与能量守恒：

辛积分器保持相空间体积（Liouville定理），对于哈密顿系统： $$\frac{\partial \dot{\mathbf{q}}}{\partial \mathbf{q}} + \frac{\partial \dot{\mathbf{p}}}{\partial \mathbf{p}} = 0$$ 修正哈密顿量： $$\tilde{H} = H + \Delta t^2 H_2 + \Delta t^4 H_4 + O(\Delta t^6)$$ 辛积分器精确保持 $\tilde{H}$，而不是原始的 $H$。

能量误差的后验修正：

1. 速度缩放法： $$\mathbf{v}_{corrected} = \sqrt{\frac{E_{target} - V}{T_{current}}} \cdot \mathbf{v}$$

2. 拉格朗日乘子法： 最小化 $|\mathbf{v}_{new} - \mathbf{v}_{old}|^2$ 满足 $E(\mathbf{v}_{new}) = E_0$

3. 投影法： 将系统投影到能量曲面 $E = E_0$ 上

稳定性的Lyapunov分析：

定义Lyapunov函数（通常选择总能量）： $$V(\mathbf{q}, \dot{\mathbf{q}}) = E = T + U$$ 系统稳定的充分条件： $$\dot{V} = \frac{\partial V}{\partial \mathbf{q}} \cdot \dot{\mathbf{q}} + \frac{\partial V}{\partial \dot{\mathbf{q}}} \cdot \ddot{\mathbf{q}} \leq 0$$ 对于有阻尼的系统，这自然满足。

频谱分析：

系统的固有频率通过线性化分析获得： $$\mathbf{M}\ddot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{0}$$ 其中 $\mathbf{K} = \frac{\partial^2 V}{\partial \mathbf{q}^2}$ 是刚度矩阵。

特征值问题： $$\det(\mathbf{K} - \omega^2\mathbf{M}) = 0$$ 最大频率 $\omega_{max}$ 决定了显式积分的稳定性条件。

相空间轨迹分析：

对于低维系统，可以绘制相空间轨迹 $(\mathbf{q}, \mathbf{p})$ 来直观检验：

• 封闭轨道：周期运动
• 螺旋轨道：耗散系统
• 发散轨道：数值不稳定

能量分配定理：

在热平衡下，每个二次自由度的平均能量： $$\langle \frac{1}{2}k_i q_i^2 \rangle = \frac{1}{2}k_B T$$ 这可用于验证系统的统计性质。

2.2 布料模拟的数学模型

布料是弹簧质点系统的经典应用。通过精心设计的弹簧网络，我们可以捕捉布料的拉伸、剪切和弯曲行为。

2.2.1 弹簧网格拓扑

对于一块矩形布料，我们通常使用规则网格，质点按行列排布。三种类型的弹簧赋予布料不同的力学特性：

结构弹簧（Structural Springs）

连接水平和垂直相邻的质点，抵抗拉伸变形：

• 水平：$(i,j) \leftrightarrow (i+1,j)$
• 垂直：$(i,j) \leftrightarrow (i,j+1)$

剪切弹簧（Shear Springs）

连接对角相邻的质点，抵抗剪切变形：

• 主对角：$(i,j) \leftrightarrow (i+1,j+1)$
• 副对角：$(i,j) \leftrightarrow (i+1,j-1)$

弯曲弹簧（Bending Springs）

跨越一个质点连接，抵抗弯曲：

• 水平：$(i,j) \leftrightarrow (i+2,j)$
• 垂直：$(i,j) \leftrightarrow (i,j+2)$

2.2.2 连续介质视角

从连续介质力学角度，布料的应变能密度可以写为： $$\Psi = \frac{1}{2} \lambda (\text{tr}(\mathbf{E}))^2 + \mu \text{tr}(\mathbf{E}^2) + \alpha |\mathbf{K}|_F^2$$ 其中：

• $\mathbf{E}$ 是Green应变张量
• $\mathbf{K}$ 是曲率张量
• $\lambda, \mu$ 是拉梅常数
• $\alpha$ 是弯曲刚度

离散化后，这些项对应不同类型的弹簧。

Green应变张量

对于二维布料参数化 $\mathbf{x}(u,v)$，变形梯度： $$\mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \begin{bmatrix} \frac{\partial \mathbf{x}}{\partial u} & \frac{\partial \mathbf{x}}{\partial v} \end{bmatrix}$$ Green应变张量： $$\mathbf{E} = \frac{1}{2}(\mathbf{F}^T\mathbf{F} - \mathbf{I})$$

离散化对应关系

1. 拉伸项 $\lambda (\text{tr}(\mathbf{E}))^2$ → 结构弹簧 - 抵抗面积变化 - 对应体积模量

2. 剪切项 $\mu \text{tr}(\mathbf{E}^2)$ → 剪切弹簧 - 抵抗形状畸变 - 对应剪切模量

3. 弯曲项 $\alpha |\mathbf{K}|_F^2$ → 弯曲弹簧 - 抵抗曲率变化 - 对应弯曲刚度

离散弯曲能量

对于四个共面质点形成的"蝴蝶"结构，二面角弯曲能量： $$E_{bend} = \frac{1}{2} k_b (\theta - \theta_0)^2$$ 其中 $\theta$ 是两个三角形之间的二面角： $$\cos\theta = \mathbf{n}_1 \cdot \mathbf{n}_2$$ 梯度计算涉及复杂的几何导数，实际中常用近似方法。

2.2.3 各向异性材料

真实布料往往具有各向异性，经纬方向的力学性质不同。我们可以为不同方向的弹簧设置不同刚度： $$k_{weft} \neq k_{warp}$$ 甚至可以引入非线性刚度模型： $$k(l) = k_0 + k_1 \left(\frac{l - l_0}{l_0}\right)^2$$

正交各向异性模型

更严格的处理使用正交各向异性本构关系： $$\boldsymbol{\sigma} = \begin{bmatrix} E_1/(1-\nu_{12}\nu_{21}) & \nu_{12}E_2/(1-\nu_{12}\nu_{21}) & 0 \\ \nu_{21}E_1/(1-\nu_{12}\nu_{21}) & E_2/(1-\nu_{12}\nu_{21}) & 0 \\ 0 & 0 & G_{12} \end{bmatrix} \boldsymbol{\epsilon}$$ 其中：

• $E_1, E_2$ 是经纬方向的杨氏模量
• $\nu_{12}, \nu_{21}$ 是泊松比
• $G_{12}$ 是剪切模量

对称性要求：$\frac{\nu_{12}}{E_1} = \frac{\nu_{21}}{E_2}$

纤维方向追踪

在大变形下，需要追踪材料主方向：

1. 初始配置：定义参考方向 $\mathbf{d}_1^0, \mathbf{d}_2^0$
2. 当前配置：$\mathbf{d}_1 = \mathbf{F}\mathbf{d}_1^0$, $\mathbf{d}_2 = \mathbf{F}\mathbf{d}_2^0$
3. 方向相关刚度： $$k_{effective} = k_1 \cos^2\phi + k_2 \sin^2\phi$$ 其中 $\phi$ 是弹簧与主方向的夹角。

高级各向异性模型

1. Kawabata模型（精确描述织物力学）：

拉伸行为： $$F = \begin{cases} a_1 \epsilon & \epsilon < \epsilon_1 \\ a_2 (\epsilon - \epsilon_1) + F_1 & \epsilon_1 \leq \epsilon < \epsilon_2 \\ a_3 (\epsilon - \epsilon_2)^2 + F_2 & \epsilon \geq \epsilon_2 \end{cases}$$ 剪切行为： $$F_s = K_s \gamma + H_s \gamma^3$$ 其中 $K_s$ 是初始剪切刚度，$H_s$ 是非线性项。

2. 超弹性模型（用于大变形）：

应变能密度函数： $$\Psi = \Psi_{iso}(\bar{I}_1, \bar{I}_2) + \Psi_{aniso}(I_4, I_6)$$ 其中：

• $\bar{I}_1, \bar{I}_2$ 是各向同性不变量
• $I_4 = \mathbf{d}_1^0 \cdot (\mathbf{C} \mathbf{d}_1^0)$ 是纤维拉伸
• $I_6 = \mathbf{d}_2^0 \cdot (\mathbf{C} \mathbf{d}_2^0)$ 是第二纤维方向拉伸

具体形式（Holzapfel-Gasser-Ogden模型）： $$\Psi_{aniso} = \frac{k_1}{2k_2} \sum_{i=4,6} \left[ e^{k_2(I_i-1)^2} - 1 \right]$$

3. 微观结构模型：

考虑纱线级别的相互作用： $$\mathbf{f}_{yarn} = \mathbf{f}_{tension} + \mathbf{f}_{bending} + \mathbf{f}_{contact}$$ 纱线间接触力： $$\mathbf{f}_{contact} = k_c \left( \frac{1}{d} - \frac{1}{d_0} \right)^2 \mathbf{n}$$ 其中 $d$ 是纱线间距离。

参数识别

从实验数据识别材料参数：

双轴拉伸测试：

1. 施加已知载荷 $\mathbf{F}_1, \mathbf{F}_2$
2. 测量变形 $\lambda_1, \lambda_2$
3. 最小化目标函数： $$J = \sum_{tests} |\mathbf{F}_{measured} - \mathbf{F}_{model}(\mathbf{p})|^2$$ Kawabata评价系统（KES）：

• 拉伸测试：测量 $E_1, E_2, \nu_{12}$
• 剪切测试：测量 $G_{12}, H_s$
• 弯曲测试：测量弯曲刚度 $B$
• 压缩测试：测量厚度方向特性

逆向工程方法： 使用优化算法从视频数据推断材料参数： $$\mathbf{p}^* = \arg\min_{\mathbf{p}} \sum_{frames} |\mathbf{x}_{observed} - \mathbf{x}_{simulated}(\mathbf{p})|^2$$

2.2.4 碰撞处理

布料模拟的一大挑战是碰撞处理，包括：

布料-刚体碰撞

使用惩罚力或约束方法： $$\mathbf{f}_{penalty} = -k_c \cdot \text{penetration} \cdot \mathbf{n}$$ 其中 $\mathbf{n}$ 是接触法向。

改进的接触模型：

1. 非线性惩罚力： $$\mathbf{f}_{penalty} = -k_c \cdot \text{penetration}^{3/2} \cdot \mathbf{n}$$ 模拟Hertz接触理论

2. 摩擦力（Coulomb模型）： $$\mathbf{f}_{friction} = \begin{cases} -k_t \Delta \mathbf{x}_t & \text{if } |\mathbf{f}_t| < \mu |\mathbf{f}_n| \\ -\mu |\mathbf{f}_n| \frac{\mathbf{v}_t}{|\mathbf{v}_t|} & \text{otherwise} \end{cases}$$

自碰撞检测

空间哈希加速潜在碰撞对的筛选，然后进行精确的边-边、点-面测试。

层次化检测流程：

1. 粗筛选：AABB包围盒测试
2. 中筛选：OBB或球体包围盒
3. 精确测试： - 点-三角形距离 - 边-边距离

代表性算法：

• 空间哈希：$O(n)$ 平均复杂度
• BVH树：$O(n\log n)$ 构建，$O(\log n)$ 查询
• 连续法向锥：利用法向信息剪枝

连续碰撞检测（CCD）

对于高速运动，需要求解： $$\min_{t \in [0, \Delta t]} \text{distance}(\mathbf{x}_i(t), \text{obstacle})$$ 三次方程求解： 对于线性轨迹，点-三角形CCD归结为求解： $$\det[\mathbf{x}_1(t)-\mathbf{x}_0(t), \mathbf{x}_2(t)-\mathbf{x}_0(t), \mathbf{x}_3(t)-\mathbf{x}_0(t)] = 0$$ 这是关于 $t$ 的三次方程，可用解析或数值方法求解。

保守估计： 使用运动包围盒的保守CCD： $$\text{AABB}_{swept} = \text{AABB}(\mathbf{x}^n) \cup \text{AABB}(\mathbf{x}^{n+1})$$

2.2.5 性能优化技巧

1. 层次化表示：使用包围盒树加速碰撞检测 - 动态BVH更新策略 - 懒惰重建机制

2. 自适应时间步长：根据最大速度动态调整 $$\Delta t = \min\left(\Delta t_{max}, \frac{CFL \cdot h}{|\mathbf{v}|_{max}}\right)$$ 其中 $h$ 是最小边长，$CFL \approx 0.5$ 是Courant数

3. GPU并行化：弹簧力计算天然并行

__global__ void compute_spring_forces(Particle* p, Spring* s, vec3* forces) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    if (tid < num_springs) {
        // 每个线程处理一个弹簧
        vec3 force = compute_single_spring_force(p, s[tid]);
        atomicAdd(&forces[s[tid].i], force);
        atomicAdd(&forces[s[tid].j], -force);
    }
}

4. 稀疏性利用：大部分质点间无直接相互作用 - 使用邻接表存储连接关系 - 缓存友好的数据布局

5. 向量化优化： - SIMD指令处理多个弹簧 - 数据对齐和预取

6. 多分辨率模拟： - LOD（Level of Detail）技术 - 远处布料使用粗网格 - 视锥体剔除

高级优化策略：

7. 异步时间积分（Asynchronous Time Integration）： 不同区域使用不同时间步长： $$\Delta t_i = \frac{\Delta t_{base}}{2^{level_i}}$$ 其中 $level_i$ 根据局部刚度和速度确定。

8. 子空间模拟： 将高维系统投影到低维子空间： $$\mathbf{q} = \mathbf{q}_0 + \mathbf{U}\mathbf{r}$$ 其中 $\mathbf{U}$ 是模态矩阵，$\mathbf{r}$ 是降维坐标。

降维动力学方程： $$\mathbf{U}^T\mathbf{M}\mathbf{U}\ddot{\mathbf{r}} = \mathbf{U}^T\mathbf{f}(\mathbf{q}_0 + \mathbf{U}\mathbf{r})$$

9. GPU内存优化：

// 纹理内存存储只读数据
texture<float4, 1> tex_positions;
texture<float4, 1> tex_velocities;

// 共享内存缓存局部数据
__shared__ float3 shared_positions[BLOCK_SIZE];

// 合并内存访问
float4 pos_vel = reinterpret_cast<float4*>(particle_data)[tid];

10. 混合精度计算：

    • 力计算使用float16
    • 位置更新使用float32
    • 能量守恒检查使用float64

11. 时空相干性利用：

    • 缓存上一帧的碰撞对
    • 增量更新空间数据结构
    • 预测-校正碰撞检测

12. 自适应网格细化（AMR）： 高应变区域动态细化：

if (strain > threshold) {
    subdivide_triangle(tri);
    redistribute_mass();
    update_connectivity();
}

2.2.6 高级布料效果

撕裂模拟

当弹簧应变超过阈值时触发撕裂： $$\epsilon = \frac{l - l_0}{l_0} > \epsilon_{tear}$$ 撕裂传播需要考虑：

1. 应力集中系数
2. 裂纹扩展方向（最大主应力方向）
3. 能量释放率

褶皱和折痕

几何褶皱：通过增加网格分辨率自然产生 程序化褶皱：叠加高频位移场 $$\mathbf{x}_{wrinkled} = \mathbf{x} + A \sin(k \mathbf{u} \cdot \mathbf{x}) \mathbf{n}$$

湿布料效果

修改材料参数模拟湿润效果：

• 增加质量：$m_{wet} = m_{dry}(1 + \alpha_{water})$
• 降低刚度：$k_{wet} = k_{dry} \beta_{soft}$
• 增加阻尼：$d_{wet} = d_{dry} \gamma_{damp}$

2.3 显式与隐式时间积分器

时间积分是物理模拟的核心。选择合适的积分器直接影响模拟的稳定性、精度和性能。

2.3.1 显式欧拉法

最简单的时间积分方法： $$\mathbf{v}^{n+1} = \mathbf{v}^n + \Delta t \mathbf{M}^{-1} \mathbf{f}(\mathbf{x}^n, \mathbf{v}^n)$$ $$\mathbf{x}^{n+1} = \mathbf{x}^n + \Delta t \mathbf{v}^{n+1}$$ 优点：

• 实现简单，每步只需一次力计算
• 无需求解线性系统
• 适合大规模并行

缺点：

• 受CFL条件限制：$\Delta t < \frac{2}{\omega_{max}}$
• 对于刚性系统（大刚度），需要极小的时间步长

2.3.2 隐式欧拉法

使用未来时刻的力： $$\mathbf{v}^{n+1} = \mathbf{v}^n + \Delta t \mathbf{M}^{-1} \mathbf{f}(\mathbf{x}^{n+1}, \mathbf{v}^{n+1})$$ $$\mathbf{x}^{n+1} = \mathbf{x}^n + \Delta t \mathbf{v}^{n+1}$$ 这导致非线性方程组。通过线性化： $$\left( \mathbf{M} - \Delta t^2 \frac{\partial \mathbf{f}}{\partial \mathbf{x}} - \Delta t \frac{\partial \mathbf{f}}{\partial \mathbf{v}} \right) \Delta \mathbf{v} = \Delta t \mathbf{f}$$ 其中 $\frac{\partial \mathbf{f}}{\partial \mathbf{x}}$ 是刚度矩阵，$\frac{\partial \mathbf{f}}{\partial \mathbf{v}}$ 是阻尼矩阵。

2.3.3 Verlet积分

利用时间对称性： $$\mathbf{x}^{n+1} = 2\mathbf{x}^n - \mathbf{x}^{n-1} + \Delta t^2 \mathbf{M}^{-1} \mathbf{f}(\mathbf{x}^n)$$ 速度通过有限差分计算： $$\mathbf{v}^n = \frac{\mathbf{x}^{n+1} - \mathbf{x}^{n-1}}{2\Delta t}$$ 特点：

• 二阶精度
• 良好的能量守恒性
• 无数值阻尼

2.3.4 稳定性分析

考虑简谐振子：$m\ddot{x} + kx = 0$，频率 $\omega = \sqrt{k/m}$。

对于显式欧拉：

• 稳定条件：$\Delta t < \frac{2}{\omega}$
• 能量逐渐增长（数值不稳定）

对于隐式欧拉：

• 无条件稳定
• 能量逐渐衰减（数值阻尼）

对于Verlet：

• 稳定条件：$\Delta t < \frac{2}{\omega}$
• 能量守恒（在稳定区域内）

2.3.5 自适应时间步长

根据系统状态动态调整步长： $$\Delta t_{new} = \Delta t_{old} \cdot \min\left(1.1, \sqrt{\frac{\epsilon_{target}}{\epsilon_{actual}}}\right)$$ 其中 $\epsilon$ 是误差估计，可通过比较不同阶方法获得。

2.3.6 刚性ODE求解器

对于刚性系统，使用专门的求解器：

1. 向后微分公式（BDF）：多步隐式方法
2. Radau IIA：隐式Runge-Kutta方法
3. IMEX方法：显式处理非刚性项，隐式处理刚性项

2.4 SPH与基于位置的流体

流体模拟是物理引擎的重要应用。本节介绍两种主流的拉格朗日流体方法。

2.4.1 光滑粒子流体动力学（SPH）

SPH的核心思想是用粒子携带的物理量通过核函数插值得到连续场。

基本插值公式

对于任意物理量 $A$，其在位置 $\mathbf{r}$ 的值为： $$A(\mathbf{r}) = \sum_j m_j \frac{A_j}{\rho_j} W(\mathbf{r} - \mathbf{r}_j, h)$$ 其中：

• $m_j$ 是粒子 $j$ 的质量
• $\rho_j$ 是密度
• $W$ 是核函数
• $h$ 是光滑长度

核函数选择

常用的核函数包括：

Poly6核（用于密度计算）： $$W_{poly6}(r, h) = \frac{315}{64\pi h^9} \begin{cases} (h^2 - r^2)^3 & \text{if } r \leq h \\ 0 & \text{otherwise} \end{cases}$$ Spiky核（用于压力梯度）： $$\nabla W_{spiky}(r, h) = -\frac{45}{\pi h^6} \frac{\mathbf{r}}{r}(h-r)^2$$ Viscosity核（用于粘性力）： $$\nabla^2 W_{viscosity}(r, h) = \frac{45}{\pi h^6}(h-r)$$

密度计算

$$\rho_i = \sum_j m_j W(\mathbf{r}_i - \mathbf{r}_j, h)$$

压力计算

使用状态方程： $$p_i = k\left(\left(\frac{\rho_i}{\rho_0}\right)^\gamma - 1\right)$$ 其中 $k$ 是刚度系数，$\gamma$ 通常取7，$\rho_0$ 是静态密度。

动量方程

$$\frac{d\mathbf{v}_i}{dt} = -\sum_j m_j \left(\frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2}\right) \nabla W_{ij} + \nu \sum_j m_j \frac{\mathbf{v}_j - \mathbf{v}_i}{\rho_j} \nabla^2 W_{ij} + \mathbf{g}$$ 第一项是压力梯度，第二项是粘性力，第三项是外力。

2.4.2 基于位置的流体（Position-Based Fluids）

PBF将流体不可压缩性作为约束，通过迭代投影保证密度恒定。

密度约束

$$C_i(\mathbf{x}_1, ..., \mathbf{x}_n) = \frac{\rho_i}{\rho_0} - 1 = 0$$

位置修正

使用拉格朗日乘子法： $$\Delta \mathbf{x}_i = -\frac{C_i}{\sum_k |\nabla_{\mathbf{x}_k} C_i|^2 + \epsilon} \sum_j \nabla_{\mathbf{x}_j} C_i$$ 其中梯度为： $$\nabla_{\mathbf{x}_j} C_i = \frac{1}{\rho_0} \begin{cases} \sum_k \nabla W_{ik} & \text{if } i = j \\ -\nabla W_{ij} & \text{if } i \neq j \end{cases}$$

迭代求解算法

1. 预测位置：x* = x + Δt·v + Δt²·fext
2. 寻找邻居粒子
3. while iter < maxIter:
   a. 计算密度ρi
   b. 计算位置修正Δxi
   c. 更新位置：x* = x* + Δx

4. 更新速度：v = (x* - x) / Δt
5. 应用涡旋约束和XSPH粘性

2.4.3 表面张力

添加表面张力可以模拟水滴等现象： $$\mathbf{f}_{surface} = -\sigma \kappa \mathbf{n}$$ 其中：

• $\sigma$ 是表面张力系数
• $\kappa = -\nabla \cdot \mathbf{n}$ 是曲率
• $\mathbf{n} = \nabla c$ 是表面法向（$c$ 是颜色场）

2.4.4 边界处理

虚拟粒子法：在边界生成镜像粒子 惩罚力法：添加排斥力防止穿透 密度修正：调整边界附近的核函数积分

2.5 快速邻居搜索算法

粒子方法的计算瓶颈在于邻居搜索。对于 $N$ 个粒子，暴力搜索需要 $O(N^2)$ 的复杂度。本节介绍将复杂度降至 $O(N)$ 的算法。

2.5.1 均匀网格（Uniform Grid）

最直观的空间划分方法。

基本原理

将空间划分为边长为 $h$（影响半径）的立方体单元： $$\text{cell_index} = \left\lfloor \frac{\mathbf{x} - \mathbf{x}_{min}}{h} \right\rfloor$$ 每个粒子只需搜索其所在单元及相邻的26个单元（3D情况）。

数据结构

struct UniformGrid {
    vector<vector<int>> cells;  // cells[i] 存储第i个单元的粒子列表
    vec3 min_bound, max_bound;
    float cell_size;
    ivec3 resolution;
};

构建算法

1. 清空所有单元
2. for each particle i:
   cell_id = compute_cell_id(particle[i].position)
   cells[cell_id].push_back(i)

时间复杂度：$O(N)$

查询算法

neighbors = []
cell_id = compute_cell_id(query_position)
for each neighbor_cell in 3x3x3 neighborhood:
    for each particle_id in cells[neighbor_cell]:
        if distance(query_position, particle[particle_id].position) < h:
            neighbors.push_back(particle_id)

平均时间复杂度：$O(1)$（假设粒子均匀分布）

2.5.2 空间哈希（Spatial Hashing）

均匀网格的内存优化版本。

哈希函数

将3D网格坐标映射到一维哈希表： $$hash(i, j, k) = ((i \cdot p_1) \oplus (j \cdot p_2) \oplus (k \cdot p_3)) \mod M$$

其中 $p_1, p_2, p_3$ 是大质数，$M$ 是哈希表大小。

冲突处理

使用链表法处理哈希冲突：

struct HashEntry {
    ivec3 cell_coord;
    vector<int> particles;
};

unordered_map<int, vector<HashEntry>> hash_table;

2.5.3 紧凑哈希（Compact Hashing）

进一步优化内存访问模式。

Z-order曲线

使用Morton编码将3D坐标映射到1D：

uint64_t morton_encode(uint x, uint y, uint z) {
    uint64_t result = 0;
    for (int i = 0; i < 21; i++) {
        result |= ((x & (1 << i)) << (2*i))
                | ((y & (1 << i)) << (2*i + 1))
                | ((z & (1 << i)) << (2*i + 2));
    }
    return result;
}

排序优化

1. 计算每个粒子的Morton码
2. 根据Morton码排序粒子
3. 相邻粒子在内存中连续，提高缓存命中率

2.5.4 性能分析与优化

理论复杂度

• 构建：$O(N)$
• 单次查询：$O(k)$，其中 $k$ 是平均邻居数
• 所有粒子查询：$O(N \cdot k)$

实际优化技巧

1. SIMD加速：使用向量指令批量计算距离
2. 并行构建：使用原子操作并行插入粒子
3. 内存预分配：避免动态内存分配
4. 自适应网格：根据粒子密度调整网格分辨率

2.5.5 GPU实现要点

GPU上的邻居搜索需要特殊考虑：

1. 避免动态内存：使用固定大小的邻居列表
2. 原子操作：使用atomicAdd构建网格
3. 共享内存：缓存邻近单元的粒子数据
4. Warp发散：平衡负载避免线程空闲

__global__ void build_grid(Particle* particles, int* cell_start, int* cell_end) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= num_particles) return;

    int cell_id = compute_cell_id(particles[i].position);
    particles[i].cell_id = cell_id;

    // 后续使用基数排序 + 前缀和构建cell_start/cell_end
}

2.5.6 高级数据结构

对于非均匀分布，考虑：

1. 八叉树：自适应空间划分
2. KD-树：平衡的空间划分
3. BVH：层次包围盒

本章小结

本章深入探讨了拉格朗日视角下的物理模拟方法。我们从最基础的弹簧质点系统出发，逐步构建了复杂的布料和流体模拟框架。

核心要点回顾：

1. 弹簧质点系统：物理引擎的基石，包含了动力学模拟的所有要素 - 力的计算：弹簧力（Hooke定律）、阻尼力、外力 - 系统方程：$\mathbf{M}\ddot{\mathbf{q}} = \mathbf{f}(\mathbf{q}, \dot{\mathbf{q}}, t)$

2. 时间积分器：稳定性与精度的权衡 - 显式方法：简单高效但受CFL条件限制 - 隐式方法：无条件稳定但需求解线性系统 - Verlet积分：良好的能量守恒性

3. SPH方法：通过核函数插值实现连续场 - 核心公式：$A(\mathbf{r}) = \sum_j m_j \frac{A_j}{\rho_j} W(\mathbf{r} - \mathbf{r}_j, h)$ - 关键在于选择合适的核函数

4. Position-Based方法：将物理规律作为约束 - 迭代求解保证约束满足 - 适合实时应用

5. 邻居搜索：$O(N^2)$ → $O(N)$ 的性能飞跃 - 均匀网格：简单有效 - 空间哈希：内存优化 - 紧凑哈希：缓存友好

练习题

基础题

练习2.1 推导二维弹簧质点系统的总动能和势能表达式。

练习2.2 证明Verlet积分器的时间反演对称性。

练习2.3 计算SPH中Poly6核函数的归一化常数。

挑战题

练习2.4 设计一个自适应时间步长算法，使其能够：

• 在系统稳定时使用较大步长
• 在剧烈运动时自动减小步长
• 保证总体误差在给定容差内

练习2.5 分析并比较以下三种邻居搜索方法在不同粒子分布下的性能：

• 均匀网格
• KD-树
• 八叉树

考虑：(a) 均匀分布 (b) 高斯分布 (c) 表面分布

练习2.6 推导SPH方法的动量守恒性，并分析数值误差来源。

练习2.7 设计一个混合时间积分方案，对不同刚度的弹簧使用不同的积分器。

常见陷阱与错误

1. 数值稳定性问题

陷阱：使用过大的时间步长导致爆炸

// 错误：固定大步长
dt = 0.1;  // 对于高刚度系统可能不稳定

解决：动态调整步长或使用隐式方法

2. 邻居搜索边界错误

陷阱：忘记处理周期边界或网格边缘

// 错误：直接取模可能导致负索引
cell_x = (int)(x / cell_size) % grid_size;

解决：正确处理负坐标和越界情况

3. SPH密度计算不准

陷阱：边界附近粒子密度偏低

// 问题：边界粒子缺少邻居
density = sum(m_j * W(r_ij))

解决：使用虚拟粒子或密度修正

4. 弹簧力计算符号错误

陷阱：力的方向计算错误导致系统不稳定

// 错误：符号可能反了
F = k * (length - rest_length) * direction

解决：仔细验证力的方向，确保满足牛顿第三定律

5. 内存访问模式

陷阱：随机访问粒子数据导致缓存未命中

// 低效：随机访问
for (auto& neighbor : neighbors) {
    force += compute_force(particles[neighbor]);
}

解决：数据排序和预取优化

最佳实践检查清单

设计阶段

• [ ] 选择合适的粒子表示（质点、刚体、可变形体）
• [ ] 确定所需的物理精度和性能要求
• [ ] 评估显式vs隐式积分的权衡
• [ ] 设计可扩展的数据结构

实现阶段

• [ ] 验证力计算的对称性和守恒性
• [ ] 实现稳定性监测（能量、动量）
• [ ] 添加碰撞检测和响应
• [ ] 优化内存布局和访问模式

调试阶段

• [ ] 单元测试每个力的计算
• [ ] 验证简单系统的解析解
• [ ] 检查数值稳定性条件
• [ ] 分析性能瓶颈

优化阶段

• [ ] 并行化力计算和邻居搜索
• [ ] 使用SIMD指令加速
• [ ] 实现自适应时间步长
• [ ] 考虑GPU加速

验证阶段

• [ ] 对比已知标准测试案例
• [ ] 检验物理守恒量
• [ ] 测试极端情况（高速、大变形）
• [ ] 评估不同参数下的鲁棒性