SPH基础(五): Runge-Kutta自适应时间积分

1. 引言

在流体动力学乃至更广泛的科学计算领域中，光滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH）是一种强大的无网格拉格朗日方法。SPH模拟的核心之一是对描述系统演化的常微分方程组（ODEs）进行时间积分，以更新每个粒子的物理状态（如位置、速度、内能等）。

传统的时间积分方案（如简单的欧拉法或固定步长的龙格-库塔法）虽然实现简单，但在处理复杂动态过程时面临效率与稳定性的两难困境：

步长过大：可能导致数值不稳定，模拟结果迅速发散，功亏一篑。
步长过小：虽然能保证稳定性，但会极大地增加计算成本，尤其是在模拟过程相对平稳、变化缓慢的阶段，造成了不必要的资源浪费。

为了在保证计算精度的同时最大化效率，自适应步长（Adaptive Time-Stepping） 的积分方法应运而生。本文将详细介绍自适应步长的核心思想，重点讲解龙格-库塔（Runge-Kutta）方法家族中的 RK23 和 RK45 算法，并阐述如何将其与 SPH 的物理计算流程相结合。

2. 龙格-库塔（Runge-Kutta）方法简介

2.1 广义龙格-库塔方法

龙格-库塔方法是一类用于求解常微分方程的、应用广泛的显式和隐式迭代法。其核心思想是通过在当前时间步内评估多个“中间”斜率，并用这些斜率的加权平均值来更新解，从而获得比简单欧拉法更高的精度。

对于一个形如 $\frac{dy}{dt} = f(t, y)$ 的初值问题，一个 $s$ 级的显式RK方法可以表示为：

$$ \begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + c_2 h, y_n + h a_{21} k_1) \\ k_3 &= f(t_n + c_3 h, y_n + h (a_{31} k_1 + a_{32} k_2)) \\ &\vdots \\ k_s &= f(t_n + c_s h, y_n + h \sum_{j=1}^{s-1} a_{sj} k_j) \end{aligned} $$

最终的解通过这些中间斜率的加权和来计算：

$$ y_{n+1} = y_n + h \sum_{i=1}^{s} b_i k_i $$

其中 $h$ 是时间步长，系数 $a_{ij}$, $c_i$, 和 $b_i$ 是预先确定的常数，它们的选择决定了方法的精度阶数和稳定性。例如，经典的四阶RK方法（RK4）就是这个家族的一个特例。

2.2 嵌入式RK方法：自适应步长的关键

固定步长的RK方法虽然精度高，但无法感知模拟过程的动态变化。自适应步长的精髓在于“在积分的同时估计误差”。嵌入式龙格-库塔方法（Embedded Runge-Kutta Methods） 正是为此而生。

这类方法（也称 RKF 或 Fehlberg 方法）通过一组精心设计的系数，在一次计算中同时得到两个不同阶数的解：

一个 $p$ 阶精度的解 $y_{n+1}$。
一个 $p-1$ 阶精度的嵌入解 $\hat{y}_{n+1}$。

这两个解共享大部分（甚至全部）的 $k_i$ 计算，因此额外开销很小。它们的差值则可以作为局部截断误差 $E_{n+1}$ 的一个可靠估计：

$$ E_{n+1} = \| y_{n+1} - \hat{y}_{n+1} \| $$

通过将这个误差估计 $E_{n+1}$ 与用户设定的容忍度 tol 进行比较，我们就可以动态地调整时间步长 $h$：

若 $E_{n+1} \le \text{tol}$，则接受当前步（通常使用更高阶的解 $y_{n+1}$），并可尝试在下一步增大大步长。
若 $E_{n+1} > \text{tol}$，则拒绝当前步，缩小步长并重新计算。

RK23 和 RK45 都是这个家族中的杰出代表。

3. 两种经典的自适应RK算法

3.1 RK23 (Bogacki-Shampine) 方法

RK23 方法是一种广泛应用的低阶嵌入式方法，它同时计算一个三阶解和一个二阶嵌入解。

特点：
- 它需要进行3次函数求值（计算 $k_1, k_2, k_3$）来得到一个三阶精度的解。
- 一个二阶精度的解可以通过这些 $k_i$ 的不同线性组合得到，用于误差估计。
- 它具有 FSAL (First Same As Last) 特性：一个步长计算结束时所用的 $k_3$（在新的 $y_{n+1}$ 处的值），可以作为下一个步长的 $k_1$，从而每步实际只需要额外计算2次函数求值，非常高效。
适用场景：
- 对精度要求不是特别苛刻，但希望有自适应步长能力的场景。
- 当函数 $f(t, y)$ 的计算成本较高时，其较少的函数求值次数是一个优势。

其步长调整和误差控制逻辑与高阶方法完全相同，只是系数和阶数不同。

3.2 RK45 (Dormand-Prince 5(4)) 方法

RK45 是自适应积分方法中的“黄金标准”，也是 MATLAB ode45 的默认选择。它通过一次计算得到一个五阶解和一个四阶嵌入解。

核心思想：
- 它需要进行 7 次函数求值（计算 $k_1$ 到 $k_7$）。
- 使用这些 $k_i$ 的线性组合，分别构造出五阶解 $y_{n+1}^{(5)}$ (用于更新状态) 和四阶解 $y_{n+1}^{(4)}$ (用于误差估计)。
- Dormand-Prince 系数经过特别优化，使得误差估计 $| y_{n+1}^{(5)} - y_{n+1}^{(4)} |$ 相对于步长 $h$ 更加平滑和精确。
- 同样具备 FSAL 特性，计算 $k_7$ 的函数值可以复用于下一步的 $k_1$，使得每个成功步长的平均函数求值次数约为6次。
误差控制策略 (模仿 ode45)：为了使误差控制更具鲁棒性，我们不使用固定的绝对误差，而是结合相对容忍度 (RelTol) 和绝对容忍度 (AbsTol)。对于状态向量 y 的每个分量 y_i，容忍度阈值 Tol_i 定义为：
$$ \text{Tol}_i = \text{RelTol} \times |y_i| + \text{AbsTol} $$
这个策略的优点是：
- 当解的数值很大时，误差控制主要由相对容忍度决定。
- 当解趋近于零时，由绝对容忍度托底，防止步长被无限压缩。
最终，我们计算一个归一化的误差率 err_rate：
$$ \text{err\_rate} = \sqrt{ \frac{1}{N} \sum_{i=1}^{N} \left( \frac{E_{n+1, i}}{\text{Tol}_i} \right)^2 } $$
其中 $E_{n+1, i}$ 是第 $i$ 个分量的误差估计。
步长调整决策：
- 如果 err_rate <= 1.0: 接受当前步。使用更高阶的解 $y_{n+1}^{(5)}$ 更新状态，并计算下一个建议步长 $h_{\text{new}}$。
- 如果 err_rate > 1.0: 拒绝当前步。状态回退，使用一个更小的步长 $h_{\text{new}}$ 重新计算。
步长调整的经典公式为：
$$ h_{\text{new}} = h_{\text{old}} \times \text{safe} \times \left( \frac{1.0}{\text{err\_rate}} \right)^{p} $$
- safe: 安全因子，通常取 0.9。
- 指数 p: 对于RK45，通常取 0.2 (即 $1/5$)。

4. SPH右端项（RHS）的计算

在前面的讨论中，我们反复提到函数 $f(t, y)$，它代表了系统状态量的时间导数，即常微分方程的右端项（Right-Hand Side, RHS）。在 SPH 模拟中，这个函数 compute_derivatives 的任务就是根据当前所有粒子的状态，计算出它们各自的时间导数。

对于复杂的物理过程，尤其是涉及固体力学的弹塑性、损伤和断裂时，RHS 的计算远不止一个简单的压力梯度。下面是构成 SPH 中 RHS 的核心部分：

速度导数（加速度 $\mathbf{a}$）: 这是动量方程的右端项。加速度由作用在粒子上的所有力的总和除以质量得到。
$$ \frac{d\mathbf{v}_i}{dt} = \mathbf{a}_i = \frac{1}{m_i} \sum_j \mathbf{F}_{ij} $$
力 $\mathbf{F}_{ij}$ 包括：
- 压力梯度力: 由压强 $P$ 产生。
- 粘性力: 人工粘性，用于处理冲击波。
- 应力散度力: 这是固体力学中的关键项。总应力张量 $\boldsymbol{\sigma}$ 可以分解为各向同性的压力 $P$ 和偏应力张量 $\mathbf{S}$。由偏应力引起的力代表了材料的抗剪切和形变能力。 $$ \boldsymbol{\sigma} = -P\mathbf{I} + \mathbf{S} $$ 因此，加速度的计算需要准确的应力信息。
应力导数（应力率 $\dot{\mathbf{S}}$）: 应力本身不是一个守恒量，它会随着材料的变形而演化。为了计算应力的变化，我们需要本构模型（Constitutive Model）。对于弹塑性材料，通常使用客观应力率（如 Jaumann 率）来描述偏应力张量的时间导数：
$$ \frac{d\mathbf{S}_i}{dt} = \text{JaumannRate}(\mathbf{S}_i, \dot{\boldsymbol{\epsilon}}_i, \boldsymbol{\Omega}_i) $$
其中 $\dot{\boldsymbol{\epsilon}}$ 是应变率张量，$\boldsymbol{\Omega}$ 是自旋张量，均由速度梯度计算得到。
屈服与塑性: 当应力达到材料的**屈服强度（Yield Strength）**时，材料进入塑性流动状态。屈服模型（Yield Model），如 von Mises 或 Drucker-Prager，定义了这个边界。对于岩石等材料，屈服强度还依赖于压力。一旦屈服，就需要通过“径向返回算法”将应力拉回到屈服面上，这个过程是非线性的，并决定了塑性变形的能量耗散。
损伤导数（损伤率 $\dot{D}$）: 为了模拟材料的开裂和失效，我们引入一个内部状态变量——损伤 $D$ (从0到1)。损伤模型（Damage Model），如 Grady-Kipp 模型，描述了损伤如何随着拉伸或应变而累积。
$$ \frac{dD_i}{dt} = g(\boldsymbol{\sigma}_i, \dot{\boldsymbol{\epsilon}}_i, D_i, \dots) $$

总结：在每一个时间积分步中，compute_derivatives 函数的计算流程大致如下：

根据当前密度 $\rho$ 和内能 $u$，通过状态方程计算压力 $P$。
利用屈服模型判断当前应力状态。
通过本构模型计算弹性试探应力，如果屈服则进行塑性修正，得到最终的偏应力 $\mathbf{S}$。
利用损伤模型更新损伤变量 $D$。
最后，将压力 $P$ 和偏应力 $\mathbf{S}$ 代入动量方程，计算出最终的加速度 $\mathbf{a}$。
同时，计算内能变化率 $\dot{u}$、密度变化率 $\dot{\rho}$ 等其他变量的导数。

这些导数共同构成了时间积分器所需要的 RHS 向量。

5. 应用算例简述

在实现复杂的自适应积分器后，用一些已知解或具有守恒律的简单问题进行验证是至关重要的一步。

5.1 太阳系模拟

这是一个经典的 N 体问题。每个行星（粒子）的 RHS 就是其他所有天体对其施加的万有引力之和。

$$ \frac{d\mathbf{v}_i}{dt} = \mathbf{a}_i = \sum_{j \ne i} \frac{G m_j (\mathbf{r}_j - \mathbf{r}_i)}{\| \mathbf{r}_j - \mathbf{r}_i \|^3} $$

验证重点：
- 长期能量守恒和角动量守恒：是衡量积分器好坏的关键指标。
- 轨道精度：能否准确再现行星的椭圆轨道。
自适应步长优势：对于具有高偏心率轨道的天体（如彗星），它在靠近太阳时（速度快，引力变化剧烈）会自动采用小步长，而在远离时（速度慢，引力平缓）则采用大步长，兼顾了精度和效率。

5.2 单摆模拟

一个简单的单摆系统由以下一阶方程组描述：

$$ \begin{cases} \frac{d\theta}{dt} = \omega \\ \frac{d\omega}{dt} = -\frac{g}{L} \sin(\theta) \end{cases} $$

验证重点：
- 周期稳定性：对于小角度摆动，周期应接近 $2\pi\sqrt{L/g}$。
- 能量守恒：在无阻尼情况下，总能量（动能+势能）应保持不变。
自适应步长优势：当摆锤经过最低点（速度最快）时，步长会自然减小；在最高点（速度为零）时，步长会增大。这展示了积分器对系统动态的灵敏响应。

6. 结合SPH的计算伪代码

现在，我们将上述 RK45 误差控制逻辑整合到 SPH 的主循环中。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
// --- 全局参数 ---
double RelTol = 1e-6; // 相对容忍度
double AbsTol = 1e-9; // 绝对容忍度
double h_min = 1e-8, h_max = 1e-2; // 步长上下限
double safe_factor = 0.9;
double max_increase = 5.0, min_decrease = 0.2;

// --- 主循环 ---
double t = 0.0;
double h = initial_dt; // 初始步长
// y_n 包含了所有粒子的位置、速度、内能、应力等状态量
StateVector y_n = get_initial_conditions(); 

// 预计算下一步的 k1 (利用 FSAL 特性)
RHSVector k1 = compute_derivatives(t, y_n);

while (t < T_max) {
    bool step_accepted = false;

    while (!step_accepted) {
        if (h < h_min) {
            error("Timestep smaller than h_min");
            break;
        }

        // 1. RK45 核心计算：利用已有的 k1 计算 k2, ..., k7
        //    并得到 y_next_4 (四阶解) 和 y_next_5 (五阶解)
        //    (此处省略繁杂的系数计算，但会返回 k_next 用于 FSAL)
        auto [y_next_4, y_next_5, k_next] = perform_rk45_step(t, y_n, h, k1);

        // 2. 计算误差率 err_rate (模仿 MATLAB)
        double err_rate = calculate_error_rate(y_n, y_next_4, y_next_5, RelTol, AbsTol);

        // 3. 决策与步长调整
        double h_new = h * safe_factor * pow(err_rate, -0.2);

        if (err_rate <= 1.0) {
            // --- 接受步长 ---
            step_accepted = true;
            t += h;
            y_n = y_next_5; // 更新状态为更高阶的解
            k1 = k_next;    // FSAL: 下一步的 k1 已经算好

            // 限制步长增幅
            h = min({h * max_increase, h_new, h_max});

        } else {
            // --- 拒绝步长 ---
            // 状态 y_n, t, k1 保持不变

            // 限制步长降幅
            h = max({h * min_decrease, h_new, h_min});
        }
    }

    // (可选) 在每个成功的时间步后，更新邻域、输出数据等
    UpdateAndSaveData(t);
}

// `compute_derivatives` 函数实现了第4节描述的RHS计算逻辑。
// `perform_rk45_step` 和 `calculate_error_rate` 分别实现了第3节的算法核心和误差控制。