基于应变的孔隙度模型（ε-α 模型）

GASPHiA 技术细节

一、模型要解决什么问题？

在冲击（如陨石撞击）或爆炸加载下，多孔材料（砂土、rubble pile等）的行为与密实材料完全不同。绝大多数压力用于孔隙的压实过程，而孔隙压实又会大量放热，导致膨胀。总之冲击多孔隙度材质时，有如下特点：

传统 P-α 模型直接用压力 P 驱动压实，需要迭代求解 α 和 P（不过我根据论文实现了无需迭代的 P-α 模型）。

ε-α 模型的核心特点：改用体积应变 \(\epsilon_V\) 直接驱动压实，α 与 P 解耦，先更新 α，再算 P。

二、基本物理量

由于天文撞击中，压力一般都很大，体积应变也很大，所以大家一般都默认不考虑弹性压实阶段，直接进入塑形压实（不可逆），详细关于弹性压实阶段的信息可以参考论文[1][2].

符号含义公式
\(\phi\)孔隙度\(\phi = V_{\mathrm{V}} / V\)
\(\alpha\)膨胀系数\(\alpha = 1/(1-\phi) = \rho_{\mathrm{S}} / \rho\)
\(\rho\)体密度\(\rho = \rho_{\mathrm{S}} / \alpha\)
\(\epsilon_V\)体积应变\(\epsilon_V = \ln(V/V_0)\)

说明：

三、模型的基本假设

四、核心公式

因为孔隙压实是依靠体积应变驱动的，首先需要了解如何计算体积应变。在 SPH 计算中，体积应变率 \(\dot{\epsilon}_V\) 可以通过速度散度来计算：

笛卡尔坐标系下：

\[ \dot{\epsilon}_V = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} \]

其中 \(u, v, w\) 分别是 \(x, y, z\) 方向的速度分量。

原始论文[1]给出的是轴对称柱坐标系下的格式：

\[ \dot{\epsilon}_V = \frac{\partial u}{\partial r} + \frac{u}{r} + \frac{\partial v}{\partial z} \]

其中 \(u\) 是径向速度，\(v\) 是轴向速度。

对 \(\dot{\epsilon}_V\) 逐步积分即可得到体积应变 \(\epsilon_V\)。

此外，孔隙在强烈的压实过程中会做功并转化为大量内能（热量），温度的上升反过来会引起基质材料的热膨胀。文献[2]在改进 \(\epsilon\)-\(\alpha\) 模型时，考虑了这种内能增加对体积应变造成的"反弹"效应（即热膨胀抵消了部分机械压实）。为了获取真实驱动孔隙压实的有效体积应变率，需要从总应变率中扣除热膨胀项：

\[ \frac{d\epsilon_V}{dt} = \left( \frac{d\epsilon_V}{dt} \right)_{\mathrm{total}} - \frac{\Gamma_{s0}}{c_{s0}^2} \frac{dE}{dt} \]

其中，\(\left( \frac{d\epsilon_V}{dt} \right)_{\mathrm{total}}\) 为由速度散度计算得到的应变率，\(\Gamma_{s0}\) 为基质材料的初始 Grüneisen 参数，\(c_{s0}\) 为基质材料的初始体声速，\(E\) 为单位质量的内能。

在每个时间步中，根据当前的 \(\epsilon_V\) 判断所处的压实阶段，进而计算对应的 \(\alpha\) 值。

4.1 指数压实段（主要压实）

当 \(\epsilon_V < \epsilon_X\) 时：

\[ \alpha = \alpha_0 \cdot e^{\kappa \epsilon_V} \]

其中 \(\kappa\) 是压实速率，\(0 < \kappa \leqslant 1\)。\(\kappa\) 越小，压实越"慢"（需要更大应变才能闭合）。

4.2 幂律压实段（后期硬压实）

当 \(\epsilon_X < \epsilon_V < \epsilon_C\) 时：

\[ \alpha = 1 + (\alpha_X - 1) \left( \frac{-\epsilon_V}{-\epsilon_X} \right)^2 \]

其中：

\[ \alpha_X = \alpha_0 e^{\kappa \epsilon_X} \]

4.3 完全压实段

当 \(\epsilon_V \leq \epsilon_C\) 时：

\[ \alpha = 1 \]

4.4 \(\epsilon_C\) 的确定（保证导数连续）

\[ \epsilon_C = 2 \cdot \frac{1 - \alpha_X}{\kappa \alpha_X} + \epsilon_X \]

实际计算时，\(\epsilon_C\) 其实只与材料参数有关，确定 \(\alpha_X\)、\(\kappa\)、\(\epsilon_X\) 就可以得到 \(\epsilon_C\)。

4.5 参数汇总

参数物理含义典型范围
\(\alpha_0\)初始膨胀系数1.0 ~ 2.5
\(\kappa\)压实速率0.7 ~ 0.98
\(\epsilon_X\)过渡应变（指数→幂律）-0.4 ~ -0.2

运行方式

测试目录：


GASPHiA-Tests/Collins2011_EpsilonAlphaPorousIronHugoniot

执行以下命令启动：


./run_all.sh --source-dir /path/to/GASPHiA --cuda 0

---

五、模型验证

完成了前文所述的核心计算逻辑，特别是针对内能增加导致的热膨胀修正（即从总应变率中扣除热应变项：\(d\epsilon/dt = (d\epsilon/dt)_{\text{total}} - (\Gamma_{s0}/c_{s0}^2)(dE/dt)\)）之后，我们需要一个标准的基准测试来验证程序实现的正确性。为此，我们选取了文献[2]中的案例：初始孔隙度为 75% 的多孔铁的冲击 Hugoniot 曲线计算。

1. 物理背景与为什么选择此案例？

对于极高孔隙率（>50%）的材料，在遭受强冲击压实时，孔隙空间的塌陷会转化为极高的内能（热量）。这种极端的加热效应会导致基质（固体）材料发生显著的热膨胀。

原始 \(\epsilon\)-\(\alpha\) 模型的局限：原始模型假设基质密度永远大于初始密度（即忽略了热膨胀的贡献）。从图中可以看出，在高压区，原始模型（虚线）预测的密度持续增大，完全失效。

改进后的模型：在引入了热体积应变修正后，模型能够反映出在极高压力下，热膨胀效应甚至会超过机械压缩效应，导致宏观密度不增反降。

2. 计算结果与图表分析

在此次验证计算中，我们采用了与文献一致的参数（状态方程采用 Tillotson EOS，孔隙参数设为 \(\kappa = 0.98, \epsilon_e = 0, \alpha_0 = 4\)），并将我们的 SPH 程序计算结果与文献数据进行了对比。

Hugoniot曲线对比

图：初始孔隙度为 75% 的多孔铁的 Hugoniot 曲线。横轴为密度 (g/cc)，纵轴为压力 (GPa)。

通过观察上图，我们可以得出以下清晰的结论：

参考资料

[1] Wünnemann K, Collins GS, Melosh HJ. A strain-based porosity model for use in hydrocode simulations of impacts and implications for transient crater growth in porous targets. Icarus. 2006 Feb 1;180(2):514-27.

[2] Collins GS, Melosh HJ, Wünnemann K. Improvements to the ɛ-α porous compaction model for simulating impacts into high-porosity solar system objects. International Journal of Impact Engineering. 2011 Jun 1;38(6):434-9.