1. What Problem Does the Model Solve?
Under impact (e.g., meteorite impact) or explosive loading, porous materials (sand, rubble pile, etc.) behave completely differently from dense materials. Most of the pressure is used for pore compaction, and the pore compaction releases large amounts of heat, causing expansion. In summary, when impacting porous materials, the following characteristics are observed:
- Absorbs large amounts of energy
- Generates localized high temperatures
- Significantly weakens shock wave strength
- Ultimately affects crater size and morphology
The traditional P-α model uses pressure P directly to drive compaction, requiring iterative solving of α and P (though I have implemented an iterative-free P-α model based on the literature).
Core feature of the ε-α model: Uses volumetric strain $\epsilon_V$ to directly drive compaction, decoupling α from P — first update α, then calculate P.
2. Basic Physical Quantities
Since astronomical impacts involve high pressures and large volumetric strains, the elastic compaction phase is generally ignored, and plastic compaction (irreversible) is entered directly. For more details on the elastic compaction phase, refer to papers [1][2].
| Symbol | Meaning | Formula |
|---|---|---|
| $\phi$ | Porosity | $\phi = V_{\mathrm{V}} / V$ |
| $\alpha$ | Distension | $\alpha = 1/(1-\phi) = \rho_{\mathrm{S}} / \rho$ |
| $\rho$ | Bulk density | $\rho = \rho_{\mathrm{S}} / \alpha$ |
| $\epsilon_V$ | Volumetric strain | $\epsilon_V = \ln(V/V_0)$ |
Notes:
- No porosity: $\alpha = 1$
- 50% porosity: $\alpha = 2$
- Under compression, V decreases → $\epsilon_V < 0$
3. Basic Assumptions of the Model
- Compaction is irreversible: α remains unchanged during unloading (no rebound)
- Compaction rate controlled by strain: $\mathrm{d}\alpha/\mathrm{d}\epsilon_V$ does not directly depend on pressure
- Compaction first, then solid compression: Solid compression phase only begins after pores are fully closed
4. Core Formulas
Since pore compaction is driven by volumetric strain, we first need to understand how to calculate volumetric strain. In SPH calculations, the volumetric strain rate $\dot{\epsilon}_V$ can be computed from velocity divergence:
Cartesian coordinate system:
$$ \dot{\epsilon}_V = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} $$where $u, v, w$ are velocity components in the $x, y, z$ directions respectively.
Paper [1] gives the axisymmetric cylindrical coordinate form:
$$ \dot{\epsilon}_V = \frac{\partial u}{\partial r} + \frac{u}{r} + \frac{\partial v}{\partial z} $$where $u$ is radial velocity and $v$ is axial velocity.
Integrating $\dot{\epsilon}_V$ over time gives the volumetric strain $\epsilon_V$.
Additionally, during intense compaction, pores do work that transforms into large amounts of internal energy (heat). The temperature rise causes thermal expansion of the matrix material. Paper [2] improved the $\epsilon$-$\alpha$ model by considering this “rebound” effect of internal energy increase on volumetric strain (i.e., thermal expansion partially offsets mechanical compaction). To obtain the effective volumetric strain rate that truly drives pore compaction, the thermal expansion term must be subtracted from the total strain rate:
$$ \frac{d\epsilon_V}{dt} = \left( \frac{d\epsilon_V}{dt} \right)_{\mathrm{total}} - \frac{\Gamma_{s0}}{c_{s0}^2} \frac{dE}{dt} $$where $\left( \frac{d\epsilon_V}{dt} \right){\mathrm{total}}$ is the total strain rate calculated from velocity divergence, $\Gamma{s0}$ is the initial Grüneisen parameter of the matrix material, $c_{s0}$ is the initial bulk sound speed of the matrix material, and $E$ is the specific internal energy.
At each time step, the current phase is determined based on $\epsilon_V$, and the corresponding $\alpha$ value is calculated.
4.1 Exponential Compaction (Main Compaction)
When $\epsilon_V < \epsilon_X$:
$$ \alpha = \alpha_0 \cdot e^{\kappa \epsilon_V} $$where $\kappa$ is the compaction rate, $0 < \kappa \leqslant 1$. Smaller $\kappa$ means slower compaction (requires larger strain to close pores).
4.2 Power-Law Compaction (Late-Stage Hard Compaction)
When $\epsilon_X < \epsilon_V < \epsilon_C$:
$$ \alpha = 1 + (\alpha_X - 1) \left( \frac{-\epsilon_V}{-\epsilon_X} \right)^2 $$where:
$$ \alpha_X = \alpha_0 e^{\kappa \epsilon_X} $$4.3 Complete Compaction
When $\epsilon_V \leq \epsilon_C$:
$$ \alpha = 1 $$4.4 Determination of $\epsilon_C$ (Ensuring Continuity of Derivative)
$$ \epsilon_C = 2 \cdot \frac{1 - \alpha_X}{\kappa \alpha_X} + \epsilon_X $$In practice, $\epsilon_C$ depends only on material parameters — once $\alpha_X$, $\kappa$, and $\epsilon_X$ are determined, $\epsilon_C$ can be calculated.
4.5 Summary of Parameters
| Parameter | Physical Meaning | Typical Range |
|---|---|---|
| $\alpha_0$ | Initial distension | 1.0 ~ 2.5 |
| $\kappa$ | Compaction rate | 0.7 ~ 0.98 |
| $\epsilon_X$ | Transition strain (exponential → power-law) | -0.4 ~ -0.2 |
5. Model Validation
After implementing the core computation logic described above, particularly the thermal expansion correction due to internal energy increase (subtracting the thermal strain from the total strain rate: $d\epsilon/dt = (d\epsilon/dt){\text{total}} - (\Gamma{s0}/c_{s0}^2)(dE/dt)$), we need a standard benchmark test to verify the correctness of the implementation. For this, we choose the representative case from literature [2]: the shock Hugoniot curve calculation for porous iron with 75% initial porosity.
1. Physical Background and Why This Case?
For materials with very high porosity (>50%), intense shock compression transforms the pore space collapse into extremely high internal energy (heat). This extreme heating causes significant thermal expansion of the matrix (solid) material.
Limitations of the original $\epsilon$-$\alpha$ model: The original model assumes the matrix density is always greater than the initial density (i.e., it ignores thermal expansion). As seen in the figure, in the high-pressure region, the original model (dashed line) predicts continuously increasing density, completely failing.
Improved model: After introducing the thermal volumetric strain correction, the model can reflect that at extremely high pressures, the thermal expansion effect can even exceed the mechanical compression effect, causing the macroscopic density to decrease instead of increase.
2. Computational Results and Analysis
In this validation calculation, we use parameters consistent with the literature (Tillotson EOS for the equation of state, porosity parameters: $\kappa = 0.98, \epsilon_e = 0, \alpha_0 = 4$), and compare our SPH program results with the literature data.
Figure: Hugoniot curve for porous iron with 75% initial porosity. Horizontal axis: density (g/cc), vertical axis: pressure (GPa).
Based on the figure above, we can draw the following clear conclusions:
Captured the Hugoniot curve “inversion” phenomenon: After pressure reaches about 10 GPa, our numerical results perfectly show the curve turning left. In the 10-100 GPa range, as pressure increases, material density decreases, which is caused by extreme heating-induced thermal expansion of the matrix material.
Consistent with theoretical analysis and reference hydrocode: We overlaid the state points from our program with the theoretical analytical solution (solid line) of the improved $\epsilon$-$\alpha$ model from the literature. As can be seen, they are in good agreement. Since paper [2] does not provide detailed EOS information, I can only speculate that some errors at low pressures may be because the paper did not use the Tillotson EOS.
References
[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.