# A multiphase irreversible-compaction model for granular-porous materials

Nicolas Favrie, Kevin Schmidmayer, Jacques Massoni

► To cite this version:

Nicolas Favrie, Kevin Schmidmayer, Jacques Massoni. A multiphase irreversible-compaction model for granular-porous materials. Continuum Mechanics and Thermodynamics, 2022, 34 (1), pp.217-241. 10.1007/s00161-021-01054-8 . hal-03257608v2

HAL Id: hal-03257608

https://hal.science/hal-03257608v2

Submitted on 2 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A multiphase irreversible-compaction model for granular-porous materials

Nicolas Favrie · Kevin Schmidmayer · Jacques Massoni

Received: date / Accepted: date

**Abstract** An Eulerian, hyperbolic, multiphase-flow model for dynamic and irreversible compaction of porous materials is constructed. A reversible model for elastic, compressible, porous material is derived. Classical homogenization results are obtained. The irreversible model is then derived in accordance with the following basic principles. First, the entropy inequality is satisfied by the model. Second, the stress coming from the elastic energy decreases in time (the material behaves as Maxwell-type materials). The irreversible model admits an equilibrium state corresponding to a Gurson-type limit which varies with the porosity. The sound velocity at the yield limit is smaller than that of the reversible model. Such an embedded model structure ensures a thermodynamically correct formulation of the porous-material model. The usual model used in the detonation community is recovered. The model is then validated on quasi-static loading-unloading experiments with HMX. The ability of the model to capture strong shock propagation in porous material as well as its ability to deal with interface between a fluid and a porous material is demonstrated and validated on Hugoniot curve of aluminium with various porosities for a unique set of empirical parameters.

**Keywords** Multiphase-flow model · compaction · porous media · granular media · Godunov-type scheme

## 1 Introduction

When a granular bed or a porous material is subjected to mechanical loading, local deformations occur resulting in rearranged deformed grains forming a compact porous solid. Such a process is irreversible. Indeed, if the loading constraint is removed, the porous solid slightly expands but never recovers its initial volume.

In the area of powder compaction and porous material, the literature provides at least three classes of models, without or with non-trivial connections:

- The first class is composed of elastic and plastic modelling in solids [14, 29, 32, 3, 17].

- The second class considers granular media at the discrete level [25].

- The third class is composed of multiphase-flow modelling of granular media [30, 2, 7, 4, 22, 33, 12].

In the first and second classes, the wave dynamics, as well as the material compressibility, are often neglected. Whereas in the last class, the wave dynamics is taken into account but not the irreversibility of the compaction process.

5

10

15

20

25

30

<footnote>N.Favrie
Aix Marseille Univ, IUSTI, UMR CNRS7343, 5 Rue E.Fermi 13013 Marseille
E-mail: nicolas.favrie@univ-amu.fr</footnote>

<footnote>K. Schmidmayer
Same address
E-mail: kevin.schmidmayer@univ-amu.fr</footnote>

<footnote>J. Massoni
Same address
E-mail: jacques.massoni@univ-amu.fr</footnote>
In this work, the irreversible compaction of porous or granular material is addressed in the context of the multiphase-flow theory of granular materials in the presence of a compressible gas phase. We consider a special case where the porosity is closed. Therefore, there is no diffusion of the gas in the solid matrix and a one-velocity model is sufficient for the dynamics description. The one-velocity approach can also be used in the case of open porosity, where the characteristic diffusion time is much larger than the characteristic time of the dynamic process. Some success has been reached in this area with Kapila's model [22]. In particular, they introduced a configuration energy (or granular energy) as a function of the volume fraction. This is the energy necessary to maintain the contact between grains. However, such an approach is unable to model hysteresis phenomena. [16] (see also [21]) proposed a model able to deal with hysteresis phenomena. However, the model does not take into account the second phase (trapped gas). Also, the model uses a new variable called "no-load volume fraction" in terms of which the "compaction energy" and inter-granular stress are expressed. The physical sense of the "no-load volume fraction" is not direct. Another attempt appears in [12] where hysteresis phenomena are taken into account but where the model is not able to capture the pure-phase behaviours.

The aim of this paper is to propose a mathematical model of granular materials expressed in terms of well-defined physical characteristics. This model is analogous to multiphase-flow models developed for solid-fluid interactions in [13, 11, 28] and is able to recover the following limits in a unique multiphase-flow model:

- Capture the quasi-static, elastic behaviour of a porous media obtained by homogenization as in [20] or by acoustic measurement.

- Capture the static and dynamic behaviour of pure solids and pure fluids.

- Capture the plastic behaviour of the porous solids [17, 41] and granular materials [21].

- Consider the compressibility of the solid and fluid phases in the porous material and be able to capture the Hugoniot of the material with different porosity and with the same set of parameters.

- Able to capture the interface between a fluid and a porous material.

The model also verifies the following properties:

- It preserves mass, momentum and total energy.

- It is hyperbolic.

- It verifies the second law of thermodynamics.

Numerical methods to solve interface conditions in the context of compressible fluids governed by different equations of state (EOS) have been paid considerable attention during the last two decades owing to the pioneering works of [1]. Multiphase-flow models have been used in this work in order to solve interface conditions with correct thermodynamics at interfaces, but in the context of fluid-fluid interfaces only. Surface tension was addressed later in [37], evaporating interfaces in [35] and solid-fluid interfaces by [11].

In Section 2, the reversible model is derived and validated on the quasi-static limit for elastic, porous material. In Section 3, dissipative effects are introduced to take into account the plastic evolution of the material. We also present the plastic limit used for modelling granular energetic materials and solid porous materials. Hydrodynamics closure for the different materials is presented in Section 4. These EOS can be calibrated using experiments on load-unload cycles. In Section 5, we present the numerical method to solve the problem. One-dimensional examples and a parameter analysis of the model are presented in Section 6. Conclusion is addressed in Section 7.

## 2 Reversible model

The aim of this section is to build a flow model describing the reversible compaction of the granular material in contact with compressible fluid. In this section, as in the whole paper, we will consider that the shear stress can be neglected. Comments on the inclusion of shear stress are made in the conclusion. The model must take into account the compressibility of the solid and the fluid (denoted with subscripts *s* and *g*, respectively) and it must be able to deal with interfaces between pure fluids, granular materials, porous materials and pure materials.
## 2.1 Derivation of the reversible model

The mass conservation law reads

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (1) $$

where $\rho = \alpha_s \rho_s + \alpha_g \rho_g$ is the mixture density, $\alpha_k$ and $\rho_k$ are the volume fraction and density of phase $k$. The volume fractions satisfy the saturation constraint $\sum_k \alpha_k = 1$. The mass and entropy of each phase are also conserved

$$ \frac{dY_k}{dt} = 0 \quad \text{and} \quad \frac{d\eta_k}{dt} = 0, \qquad (2) $$

where $Y_k = \frac{\alpha_k \rho_k}{\rho}$ and $\eta_k$ are the mass fraction and entropy of phase $k$, and with the material derivative operator

$$ \frac{d(\cdot)}{dt} = \frac{\partial(\cdot)}{\partial t} + \mathbf{u} \cdot \nabla(\cdot). \qquad (3) $$

Note that we consider here a reversible process.

We introduce herein, as in [12], a new variable $\lambda$. This variable corresponds to $\lambda = \det(F^{-1})$ where $F = \frac{\partial \mathbf{x}}{\partial \mathbf{X}}$ is the gradient deformation tensor with $\mathbf{X}$ the Lagrangian coordinates and $\mathbf{x}$ the Eulerian coordinates. Since the deformation-tensor evolution is governed by

$$ \frac{d\mathbf{F}}{dt} - \frac{\partial \mathbf{u}}{\partial \mathbf{x}}\mathbf{F} = 0, \qquad (4) $$

one can show that this variable verifies

$$ \frac{\partial \lambda}{\partial t} + \nabla \cdot (\lambda \mathbf{u}) = 0. \qquad (5) $$

As in [13] and in [12], we can write $\alpha_s \rho_s = \alpha_{s0} \rho_{s,r} \lambda$. At the initial time, we will set $\alpha_{s0}(x) = \alpha_s(x, t=0)$ and $\lambda(x) = \rho_s(x)/\rho_{s,r}$, where $\rho_{s,r}$ is a material constant. Subscript $r$ stands for reference while subscript $0$ stands for initial conditions. One can note that in the case of pure solids ($\alpha_s = 1$), we have $\rho_s = \rho_{s,r} \lambda$. Here, we modify the energy proposed in [12] and we propose the following total energy per unit mass

$$ E = Y_s \kappa(\alpha_s) f(\xi) + \sum_k Y_k e_k(\rho_k, \eta_k) + \frac{\mathbf{u}^2}{2}, \qquad (6) $$

with $\xi = \rho_{s,r} \lambda / \rho_s$. In the case of the reversible model, this variable can be rewritten $\xi = \alpha_s / \alpha_{s0}$ and is then not independant with $\alpha_s$. It therefore seems useless because functions $\kappa$ and $f$ could be merged as a unique function and because $\alpha_s$ is already computed. However, in the irreversible case detailed in Section 3, the link would be lost and this justifies the necessity of this additional variable introduced here for consistency with the irreversible formulation. In this energy, $e_k(\rho_k, \eta_k)$ is the internal energy of each phase $k$ and $\kappa(\alpha_s)$ is a positive function which increases with the solid volume fraction. This function depends on the material and will be defined later. $f$ is a positive function of the variable $\xi$ with $\frac{\partial f}{\partial \xi}\Big|_{(\xi=1)} = 0$. One can note that this variable $\xi$ will always be equal to 1 in the case of pure solids and corresponds in the case of an incompressible matrix ($\rho_s = \rho_{s,r}$) to $\xi = \lambda$. In the incompressible case, the energy corresponds to the model proposed in [12]. The added energy (first term) can be considered as a micro-mechanical energy. One can note that this approach is reminiscent of the micromorphic approach developed in [8], neglecting the second-gradient aspects.

Using Hamilton's principle and following the section 6 of Chapter 4 in [6], the momentum equation reads

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u}^2 + p \mathbf{I}) = 0, \qquad (7) $$

where $p = \alpha_g p_g + \alpha_s p_s$ is the mixture pressure. The system is closed by the pressure-equilibrium relation (details in Appendix)

$$ p_g = p_s - \rho_{s,r} \kappa \frac{\partial f}{\partial \xi} \lambda - \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f. \qquad (8) $$
115 With such a choice, the energy conservation is automatically satisfied

$$ \frac{\partial \rho E}{\partial t} + \nabla \cdot ((\rho E + p) \mathbf{u}) = 0. \qquad (9) $$

Taking the differential of the pressure equilibrium (8), one can obtain

$$ \frac{d\alpha_s}{dt} = \frac{\rho_g c_g^2 - \rho_s c_s^2 + \rho_{s,r} \kappa \frac{\partial f}{\partial \xi} \lambda + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f}{\frac{\rho_g c_g^2}{\alpha_g} + \frac{\rho_s c_s^2}{\alpha_s} + \kappa \frac{\partial^2 f}{\partial \xi^2} \frac{(\rho_{s,r} \lambda)^2}{\alpha_s \rho_s} + \rho Y_s \frac{\partial^2 \kappa}{\partial \alpha_s^2} f + 2 \frac{\partial \kappa}{\partial \alpha_s} \frac{\partial f}{\partial \xi} \rho_{s,r} \lambda} \nabla \cdot \mathbf{u}. \qquad (10) $$

with $c_k^2 = \partial p_k / \partial \rho_k$ and $c_k$ being the sound speed of phase $k$. Let's denote

$$ \frac{\rho_s \zeta}{\alpha_s} = \kappa \frac{\partial^2 f (\rho_{s,r} \lambda)^2}{\partial \xi^2} \frac{\alpha_s \rho_s}{\alpha_s} + \rho Y_s \frac{\partial^2 \kappa}{\partial \alpha_s^2} f + 2 \frac{\partial \kappa}{\partial \alpha_s} \frac{\partial f}{\partial \xi} \rho_{s,r} \lambda. \qquad (11) $$

It then becomes

$$ \frac{d\alpha_s}{dt} = \frac{\rho_g c_g^2 - \rho_s c_s^2 + p_s - p_g}{\frac{\rho_g c_g^2}{\alpha_g} + \frac{\rho_s (c_s^2 + \zeta)}{\alpha_s}} \nabla \cdot \mathbf{u}. \qquad (12) $$

Therefore, the reversible model reads

$$ \begin{aligned} \frac{d\alpha_s}{dt} &= \frac{\rho_g c_g^2 - \rho_s c_s^2 + p_s - p_g}{\frac{\rho_g c_g^2}{\alpha_g} + \frac{\rho_s (c_s^2 + \zeta)}{\alpha_s}} \nabla \cdot \mathbf{u}, \\ \frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \mathbf{u}) &= 0, \\ \frac{\partial \lambda}{\partial t} + \nabla \cdot (\lambda \mathbf{u}) &= 0, \\ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u}^2 + p \mathbf{I}) &= 0, \\ \frac{d\eta_k}{dt} &= 0, \\ \frac{\partial \rho E}{\partial t} + \nabla \cdot ((\rho E + p) \mathbf{u}) &= 0. \end{aligned} \qquad (13) $$

## 2.2 Sound speed of the equilibrium model

120 The aim of this section is to determine the sound speed of the equilibrium model to guaranty that, in the absence of dissipation, the wave propagates correctly. Moreover, the obtained sound velocity allows to compare our model against classical results from acoustics (see Section 2.3). Reader not interested in the details can jump to Eq. (17).

Let's differentiate the pressure

$$ p = \alpha_s p_s + \alpha_g p_g = p_s - \alpha_g \kappa \rho_{s,r} \lambda \frac{\partial f}{\partial \xi} - \alpha_g \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f, \qquad (14) $$
and using the mass-conservation law we can write

125

$$
\begin{align}
\frac{dp}{dt} &= \frac{\partial p_s}{\partial \rho_s} \frac{d\rho_s}{dt} + \left( \rho_{s,r} \lambda \kappa \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \frac{d\alpha_s}{dt} - \alpha_g \left( \left( \rho Y_s \frac{\partial^2 \kappa}{\partial \alpha_s^2} f + 2 \frac{\partial \kappa}{\partial \alpha_s} \frac{\partial f}{\partial \xi} \rho_{s,r} \lambda \right) + \frac{\partial^2 f (\rho_{s,r} \lambda)^2}{\partial \xi^2} \frac{\alpha_s \rho_s}{\alpha_g} \right) \frac{d\alpha_s}{dt} \nonumber \\
&\quad + \left( \alpha_g \rho_{s,r} \lambda \kappa \frac{\partial f}{\partial \xi} + \alpha_g \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \nabla \cdot \mathbf{u} \nonumber \\
&= \left( -\rho_s c_s^2 + \alpha_g \rho_{s,r} \lambda \kappa \frac{\partial f}{\partial \xi} + \alpha_g \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \nabla \cdot \mathbf{u} \nonumber \\
&\quad + \left( \rho_{s,r} \lambda \kappa \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f - \frac{\rho_s c_s^2}{\alpha_s} - \left( \alpha_g \rho Y_s \frac{\partial^2 \kappa}{\partial \alpha_s^2} f + 2 \alpha_g \frac{\partial \kappa}{\partial \alpha_s} \frac{\partial f}{\partial \xi} \rho_{s,r} \lambda + \alpha_g \frac{\partial^2 f (\rho_{s,r} \lambda)^2}{\partial \xi^2} \frac{\alpha_s \rho_s}{\alpha_g} \kappa \right) \right) \frac{d\alpha_s}{dt} \nonumber \\
&= (-\rho_s c_s^2 + \alpha_g (p_s - p_g)) \nabla \cdot \mathbf{u} + \left( p_s - p_g - \frac{\rho_s c_s^2}{\alpha_s} - \alpha_g \frac{\rho_s \zeta}{\alpha_s} \right) \frac{d\alpha_s}{dt}. \tag{15}
\end{align}
$$

Using Eq. (12), we obtain

$$
\frac{dp}{dt} = -\rho c_r^2 \nabla \cdot \mathbf{u}, \quad (16)
$$

with

$$
\rho c_r^2 = (\rho_s c_s^2 - \alpha_g (p_s - p_g)) + \left( p_g - p_s + \frac{\rho_s c_s^2}{\alpha_s} + \alpha_g \frac{\rho_s \zeta}{\alpha_s} \right) \frac{\rho_g c_g^2 - \rho_s c_s^2 + p_s - p_g}{\frac{\rho_g c_g^2}{\alpha_g} + \frac{\rho_s (c_s^2 + \zeta)}{\alpha_s}} \quad (17)
$$

The sound velocity of the reversible model $c_r$ is analogous to Wood's sound speed for fluids [43,42], i.e. when $\kappa = 0$

$$
\rho c_r^2 = \rho_s c_s^2 + \left( \frac{\rho_s c_s^2}{\alpha_s} \right) \frac{\rho_g c_g^2 - \rho_s c_s^2}{\frac{\rho_g c_g^2}{\alpha_g} + \frac{\rho_s c_s^2}{\alpha_s}} = \frac{1}{\frac{\alpha_s}{\rho_s c_s^2} + \frac{\alpha_g}{\rho_g c_g^2}}. \quad (18)
$$

2.3 Qualitative behaviour of the sound speed

130

In the small-deformation case, [20] used a first-order homogenization procedure and proposed two bounds
for the multiphase solids. The bulk modulus *K* of a solid material with spherical voids can be expressed
under the following form

$$
K_s^* = K_s + \frac{\alpha_g}{\frac{1}{K_g - K_s} + \frac{3\alpha_s}{3K_s + 4\mu_s}}, \quad (19)
$$

and when spherical particles are considered in a fluid, the bulk modulus reads

$$
K_g^* = K_g + \frac{\alpha_s}{\frac{1}{K_s - K_g} + \frac{3\alpha_g}{3K_g}}, \quad (20)
$$

where $K_k$ and $\mu_k$ are the bulk and shear modulus of phase $k$. These results are obtained in the quasi-static configuration. This corresponds to the equilibrium model. Let's consider the material at rest $p_s = p_g = 10^5$ Pa. The parameters for aluminium and air are given in Table 1. We choose herein to use a quadratic micro-mechanical energy with $f(\xi) = \chi(\xi - 1)^2$ and the following dependency on the volume fraction

$$
\kappa(\alpha_s) = \frac{\alpha_s^{n_\kappa}}{1 + a_\kappa \alpha_g}, \quad (21)
$$

with $n_κ$ and $a_κ$ some positive parameters. Other choices are possible. At the equilibrium state $\xi \approx 1$, we thus have $\zeta = 2\kappa\chi$. The equilibrium sound speed reads

135

$$
\rho c_r^2 = \rho_s c_s^2 + \left( \frac{\rho_s c_s^2}{\alpha_s} + \alpha_g \frac{\rho_s \zeta}{\alpha_s} \right) \frac{\rho_g c_g^2 - \rho_s c_s^2}{\frac{\rho_g c_g^2}{\alpha_g} + \frac{\rho_s (c_s^2 + \zeta)}{\alpha_s}} \quad (22)
$$

140
<table><caption>Table 1: Physical parameters of materials in the linear case used for [20] limits.</caption><thead><tr><th>Material</th><th>ρ<sub>k</sub> (kg m<sup>-3</sup>)</th><th>K<sub>k</sub> (Pa)</th><th>μ<sub>k</sub> (Pa)</th></tr></thead><tbody><tr><td>Air</td><td>1</td><td>0.1089 × 10<sup>6</sup></td><td>0</td></tr><tr><td>Aluminium</td><td>2700</td><td>75.75 × 10<sup>9</sup></td><td>26 × 10<sup>9</sup></td></tr></tbody></table>

On Figure 1, we plot the different bounds of Hashin and Shtrikman and that for different values of parameter $n_κ$ with $\chi = \mu$ and $a_κ = 4000$. The sound speeds are obtained using the following definition $c_k^2 = K_k/\rho_k$. One can observe that for such parameters when $n_κ = 0$, the [20] upper limit is recovered. The exponent $n_κ$ will only change the evolution for higher porosity. This example highlights that it is possible to obtain the elasticity properties of porous materials with the reversible model. Furthermore, unlike the model of [12], one can note that the pure-phase sound speeds are also recovered.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Evolution of the sound speed as function of the porosity for different parameters. The upper and lower [20] limits are plotted in dashed lines. Results for $n_κ = 0, 1$ and $6$ are given respectively with solid, dotted and dash-dotted lines.</figcaption></figure>

### 3 Irreversible model

The aim of this section is to introduce the dissipation in the previous model. The approach is similar to the one developed in [10,9,31] for the treatment of visco-plasticity and in the compaction model presented in [12]. For the sake of clarity, we introduce a multiplicative decomposition of $\lambda$: $\lambda = \lambda^p \lambda^e$, where $\lambda^p$ is the plastic part of the deformation while $\lambda^e$ is the elastic part. Note that in the reversible case, we have $\lambda = \lambda^e$ and $\lambda^p = 1$. Let's now introduce a source term $\lambda^e$, used to introduce plasticity, on the right-hand side of the equation for $\lambda^e$

$$ \frac{\partial \lambda^e}{\partial t} + \nabla \cdot (\lambda^e \mathbf{u}) = \dot{\lambda}^e. \qquad (23) $$
One can note that the equation on $\lambda$ (5) still holds and that we are losing the link $\lambda^e = \det(F^{-1})$, which means that $\alpha_s \rho_s \neq \alpha_{s0} \rho_{s,r} \lambda^e$ and therefore $\xi \neq \alpha_s / \alpha_{s0}$. This approximation can be seen as a multiplicative decomposition from plastic theory. We also introduce a relaxation term in the volume fraction equation in order to account for the differences in the acoustic behavior of both phases. This term is related to the [2] pressure non-equilibrium terms.

$$ \frac{\partial \alpha_s}{\partial t} + \mathbf{u} \cdot \nabla \alpha_s = \dot{\alpha}_s. \qquad (24) $$

Such an approach was used in [36] to replace the equilibrium condition $p_s = p_g$ by a partial differential equation to get a pressure disequilibrium. Here we use such approach to replace our pressure equilibrium condition (Eq. 8). In the following, we follow the same procedure as in [36, 12, 10] and others.

## 3.1 Entropy inequality

Using Eq. (6), we differentiate the potential energy

$$ \begin{aligned} & \rho \frac{dY_s \kappa(\alpha_s) f(\xi) + \sum Y_k e_k(\rho_k, \eta_k)}{dt} \\ &= \rho Y_s \kappa \frac{\partial f}{\partial \xi} \left( \frac{\rho_{s,r}}{\rho} \dot{\lambda}^e + \frac{\rho_{s,r} \lambda^e}{\alpha_s \rho_s} \frac{d\alpha_s}{dt} \right) + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \frac{d\alpha_s}{dt} + \rho \sum_k Y_k \frac{\partial e_k}{\partial \rho_k} \frac{d\rho_k}{dt} + \sum_k \rho Y_k \frac{\partial e_k}{\partial \eta_k} \frac{d\eta_k}{dt} \\ &= \rho Y_s \kappa \frac{\partial f}{\partial \xi} \frac{\rho_{s,r}}{\rho} \dot{\lambda}^e - \sum_k \alpha_k p_k \nabla \cdot \mathbf{u} + \left( \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \frac{d\alpha_s}{dt} + \sum_k \left( -p_k \frac{d\alpha_k}{dt} \right) + \sum_k \rho Y_k T_k \frac{d\eta_k}{dt} \\ &= \rho_{s,r} Y_s \kappa \frac{\partial f}{\partial \xi} \dot{\lambda}^e - \sum_k \alpha_k p_k \nabla \cdot \mathbf{u} + \left( p_g - p_s + \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \frac{d\alpha_s}{dt} + \sum_k \rho Y_k T_k \frac{d\eta_k}{dt}. \end{aligned} \qquad (25) $$

Using the momentum equation, we obtain

$$ \rho_{s,r} Y_s \kappa \frac{\partial f}{\partial \xi} \dot{\lambda}^e + \left( p_g - p_s + \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \frac{d\alpha_s}{dt} + \sum_k \rho Y_k T_k \frac{d\eta_k}{dt} = 0. \qquad (26) $$

This equation should be separated into two equations, because we need the entropy evolution equation for each component. Such a separation is always *phenomenological*. There is only one constraint to respect: The entropy inequality. We use here a proposition adapted from fluid-fluid and solid-fluid interaction models in [12, 36, 28, 11, 33] and others, where all the dissipation due to the plastic deformation and the elastic energy ($\kappa f$) of the solid is put in the solid phase.

$$ \rho Y_s T_s \frac{d\eta_s}{dt} = -Y_s \kappa \frac{\partial f}{\partial \xi} \rho_{s,r} \dot{\lambda}^e - \left( p_I - p_s + \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \frac{d\alpha_s}{dt} \qquad (27) $$

and

$$ \rho Y_g T_g \frac{d\eta_g}{dt} = -(p_g - p_I) \frac{d\alpha_s}{dt}. \qquad (28) $$

Here $p_I$ is the interface pressure. This pressure can be chosen for example from the Riemann problem between solid and gas [34].

$$ p_I = \frac{Z_g}{Z_s + Z_g} \left( p_s - \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} - \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) + \frac{Z_s}{Z_s + Z_g} p_g, \qquad (29) $$

where $Z_k$ are acoustic impedances of corresponding components. One can prove that with such a choice, the entropy inequality reads

$$ \begin{aligned} & \frac{\partial \rho \sum_k Y_k \eta_k}{\partial t} + \nabla \cdot (\rho \sum_k Y_k \eta_k \mathbf{u}) \\ &= -\frac{Y_s}{T_s} \kappa \frac{\partial f}{\partial \xi} \rho_{s,r} \dot{\lambda}^e - \left( \frac{1}{T_s} \left( p_I - p_s + \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) + \frac{1}{T_g} (p_g - p_I) \right) \frac{d\alpha_s}{dt} \qquad (30) \\ &= -\frac{Y_s}{T_s} \kappa \frac{\partial f}{\partial \xi} \rho_{s,r} \dot{\lambda}^e - \frac{1}{T_s T_g} \frac{T_g Z_s + T_s Z_g}{Z_g + Z_s} \left( p_g - p_s + \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \dot{\alpha}_s. \end{aligned} $$
Thus, the entropy inequality is verified if

$$ \dot{\lambda}^e = -\frac{1}{\tau} \frac{\partial f}{\partial \xi} \quad \text{and} \quad \dot{\alpha}_s = -\mu \left( p_g - p_s + \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right), \qquad (31) $$

where $\tau$ and $\mu$ are positive. $\tau$ is the characteristic relaxation time of the viscosity of the solid material (visco-plasticity) and $\frac{1}{\mu}$ is the characteristic time needed to reach the pressure equilibrium due to wave interaction in the mixture. In the following, we will consider $\mu = \mu_0 \alpha_s \alpha_g$ with $\mu_0$ a positive constant. A discussion on the value of $\tau$ is provided in Section 3.3.

## 3.2 Hyperbolicity of the irreversible model

Finally, we propose the following irreversible model of dynamical compaction for granular material where the phase energies in the system are deduced from the previously built entropy equations

$$ \begin{aligned} & \frac{\partial \alpha_s}{\partial t} + \mathbf{u} \cdot \nabla \alpha_s = -\mu (p_g - p_s + \pi), \\ & \frac{\partial \lambda^e}{\partial t} + \nabla \cdot (\lambda^e \mathbf{u}) = \dot{\lambda}^e, \\ & \frac{\partial \alpha_s \rho_s}{\partial t} + \nabla \cdot (\alpha_s \rho_s \mathbf{u}) = 0, \\ & \frac{\partial \alpha_g \rho_g}{\partial t} + \nabla \cdot (\alpha_g \rho_g \mathbf{u}) = 0, \\ & \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u}^2 + p \mathbf{I}) = \mathbf{0}, \\ & \frac{\partial \alpha_s \rho_s e_s}{\partial t} + \nabla \cdot (\alpha_s \rho_s e_s \mathbf{u}) + \alpha_s p_s \nabla \cdot \mathbf{u} = -\rho_{s,r} Y_s \kappa \frac{\partial f}{\partial \xi} \dot{\lambda}^e - (p_I + \pi) \frac{d\alpha_s}{dt}, \\ & \frac{\partial \alpha_g \rho_g e_g}{\partial t} + \nabla \cdot (\alpha_g \rho_g e_g \mathbf{u}) + \alpha_g p_g \nabla \cdot \mathbf{u} = p_I \frac{d\alpha_s}{dt}, \end{aligned} \qquad (32) $$

where

$$ \pi = \rho_{s,r} \kappa \lambda^e \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f, \qquad (33) $$

and where the internal-energy equations are compatible with the following total-energy equation

$$ \frac{\partial \rho E}{\partial t} + \nabla \cdot ((\rho E + p) \mathbf{u}) = 0. \qquad (34) $$

Since the model is invariant by the rotation group, we only consider the one-dimensional case. The model then admits the following real characteristic velocities

$$ \omega_{1,2,3,4,5,6,7} = u \quad \text{and} \quad \omega_{8,9} = u \pm \sqrt{\sum_k Y_k c_k^2}, \qquad (35) $$

and 9 linearly independent eigenvectors. The system is therefore hyperbolic.

## 3.3 How to choose the source term $\dot{\lambda}^e$

As in [12, 10, 9], we consider a classical visco-elastic theory of Maxwell-type material. The corresponding dynamical system governing such a “plastic transformation” can easily be written. For example one can take a Perzyna-type relaxation term

$$ \dot{\lambda}^e = - \left( \frac{\langle \pi - \beta \rangle}{K_y} \right)^{n_{\lambda e}} \frac{\partial f}{\partial \xi}, \qquad (36) $$

where $\langle a \rangle = \max(a, 0)$, $K_y$ is the viscosity of the material, $n_{\lambda e}$ is a parameter to deal with non-linear viscosity and $\beta$ is, in our model, a ‘plasticity limit’ defined below.
If $\mu$ is very large, we reach a pressure-equilibrium condition $p_g = p_s - \pi$. This equilibrium condition is reached in both plastic or pure-elastic cases. First, the elastic case which was studied in Section 2 and where $\lambda^e = 0$. Second, when plasticity occurs (i.e. $\lambda^e \neq 0$), we consider the limit case where $K_y$ tends to zero. In that case, the value of $\lambda^e$ is obtained solving $\pi = \beta$ and therefore $p_g = p_s - \beta$.

When the porous or granular material is undergoing compression, elastic energy will be stored in the porous matrix. Then, the porous matrix will be deformed permanently owing to plastic deformation and the rearrangement of the solid part. During the plastic process, a part of the configuration energy ($\kappa(\alpha_s)f(\xi)$) will be transformed in internal energy ($e(\rho_k, \eta_k)$). We denote $\beta$ the value of the configuration pressure on the yield limit and it can take the form of any of the two following limits that we denote $\beta^K$ and $\beta^G$ for the Kapila [22,2,33] and Gurson [17,44,41] limits, respectively:

- In hydrodynamic models of compaction, the plastic limit is often taken as a configuration energy $B(\alpha_s)$ and this choice leads to a configuration pressure $\beta^K = \rho Y_s \frac{dB(\alpha)}{d\alpha}$ (see for example [22,2,33]) with

$$B(\alpha_s) = \frac{a}{n_p}(\alpha_g \ln(\alpha_g) + (1 + \ln(\alpha_{g0}))(\alpha_s - \alpha_{s0}) - \alpha_{g0} \ln(\alpha_{g0}))^{n_p}, \quad (37)$$

where $a$ and $n_p$ are model parameters. Therefore

$$\beta^K = a \rho Y_s \ln \left( \frac{\alpha_{g0}}{\alpha_g} \right) (\alpha_g \ln(\alpha_g) + (1 + \ln(\alpha_{g0}))(\alpha_s - \alpha_{s0}) - \alpha_{g0} \ln(\alpha_{g0}))^{n_p-1}. \quad (38)$$

[[NF:

$$g = \text{Cosh} \left( \frac{\pi_s}{a \rho Y_s} \right) - \text{Cosh} \left( \ln \left( \frac{\alpha_{g0}}{\alpha_g} \right) (\alpha_g \ln(\alpha_g) + (1 + \ln(\alpha_{g0}))(\alpha_s - \alpha_{s0}) - \alpha_{g0} \ln(\alpha_{g0}))^{n_p-1} \right) \quad (39)$$

]] This configuration pressure corresponds in the model to the difference between pressures, i.e. $\beta^K = p_s - p_g$. One can note that in the previous references, the configuration energy is a reversible energy and that the loading-unloading process will follow the configuration pressure. Herein we only consider the proposed configuration pressure as a yield limit.

- Another possible choice is to consider the well known Gurson-Tvergaard-Needleman model [17,44,41] used in solid mechanics. The plastic limit reads

$$F = \left( \frac{\frac{3}{2} \mathbf{S} : \mathbf{S}}{\sigma_Y} \right)^2 + 2q_1(1 - \alpha_s)\cosh \left( q_2 \frac{3p_h}{2\sigma_Y} \right) - (1 + q_3(1 - \alpha_s)^2) = 0, \quad (40)$$

where $\mathbf{S}$ is the deviatoric part of the Cauchy stress tensor $\sigma$ and $p_h = -\frac{1}{3}\text{Tr}(\sigma)$ is the hydrostatic pressure. The yield limit of the pure material is denoted $\sigma_Y$. $q_1$, $q_2$ and $q_3$ are material parameters depending on the pore geometry. In this paper, we only consider the hydrodynamic compaction, thus, the contribution of the deviatoric part of the stress tensor is neglected. Eq. (40) can be rewritten

$$p_h = \frac{2\sigma_Y}{3q_2}\text{arccosh}\left(\frac{1+q_3\alpha_g^2}{2q_1\alpha_g}\right). \quad (41)$$

The Gurson model is originally built to deal with plastic materials under small deformation and with empty voids. Further, it does not take into account the atmospheric pressure neither the compressibility of the matrix. Let's consider the thinking experiment of a porous material with open porosity in a high-pressure liquid. At low-pressure variations, the pressure of the solid will be equal to the liquid's and there is no reason the solid plastified with the pressure increase. Hence, it seems logical to consider the Gurson limit $\beta^G$ as the previous limit $\beta$

$$\beta^G = \frac{2\sigma_Y}{3q_2}\text{arccosh}\left(\frac{1+q_3\alpha_g^2}{2q_1\alpha_g}\right). \quad (42)$$
## 3.4 Sound velocity at the yield surface

The aim of this subsection is to compare the new model with the model of [22]. In the Kapila's model, they consider a 'configuration energy' to follow what we call the yield surface. We therefore prove in this section that the sound speed of our model exactly corresponds to the Kapila's sound speed on this surface in the equilibrium condition ($p_g = p_s - \beta$). Reader who are not interested in the calculation can directly jump to Eq. (48).

To find the sound velocity, it is sufficient to consider the case where the entropy $\eta_k$ and the mass fraction $Y_k$ are constants. Taking the differential of this equilibrium condition $p_g = p_s - \beta(\alpha_s, \rho Y_s)$, we obtain

$$ \frac{\partial p_g}{\partial \rho_g} \frac{d\rho_g}{dt} = \frac{\partial p_s}{\partial \rho_s} \frac{d\rho_s}{dt} - \frac{\partial \beta}{\partial \alpha_s} \frac{d\alpha_s}{dt} - \frac{\partial \beta}{\partial \rho Y_s} \frac{d\rho Y_s}{dt}. \quad (43) $$

With the mass-conservation law we get

$$ \frac{d\alpha_s}{dt} = \frac{\rho_g c_g^2 - \rho_s c_s^2 + \rho Y_k \frac{\partial \beta}{\partial \rho Y_s}}{\frac{\rho_s c_s^2}{\alpha_s} + \frac{\rho_s c_s^2}{\alpha_s} + \frac{\partial \beta(\alpha_s, \rho Y_s)}{\partial \alpha_s}} \nabla \cdot \mathbf{u}. \quad (44) $$

By differentiating the pressure $p = \alpha_g p_g + \alpha_s p_s = p_s - \alpha_g \beta$, using the mass-conservation law and considering $\beta = \beta(\rho Y_s, \alpha_s)$, we obtain:

$$ \begin{aligned} \frac{dp}{dt} &= \frac{\partial p_s}{\partial \rho_s} \frac{d\rho_s}{dt} + \beta \frac{d\alpha_s}{dt} - \alpha_g \frac{\partial \beta}{\partial \rho Y_s} \frac{d\rho Y_s}{dt} - \alpha_g \frac{\partial \beta}{\partial \alpha_s} \frac{d\alpha_s}{dt} \\ &= \left( \alpha_g \rho Y_s \frac{\partial \beta}{\partial \rho Y_s} - \rho_s c_s^2 \right) \nabla \cdot \mathbf{u} + \left( \beta - \alpha_g \frac{\partial \beta}{\partial \alpha_s} - \frac{\rho_s c_s^2}{\alpha_s} \right) \frac{d\alpha_s}{dt}. \end{aligned} \quad (45) $$

Replacing Eq. (44) in Eq. (45), we obtain

$$ \frac{dp}{dt} = \rho c_i^2 \nabla \cdot \mathbf{u}, \quad (46) $$

with

$$ \rho c_i^2 = \left( \rho_s c_s^2 - \alpha_g \rho Y_s \frac{\partial \beta}{\partial \rho Y_s} \right) + \left( \alpha_g \frac{\partial \beta}{\partial \alpha_s} + \frac{\rho_s c_s^2}{\alpha_s} - \beta \right) \frac{\rho_g c_g^2 - \rho_s c_s^2 + \rho Y_k \frac{\partial \beta}{\partial \rho Y_s}}{\frac{\rho_s c_s^2}{\alpha_s} + \frac{\rho_s c_s^2}{\alpha_s} + \frac{\partial \beta}{\partial \alpha_s}}, \quad (47) $$

where $c_i$ is the equilibrium sound speed for the irreversible model.

If we consider the yield limit equal to the configuration energy proposed by [22] $\beta = \beta^K$, we obtain the following equilibrium sound speed identical to the one they obtained

$$ \rho c_K^2 = (\rho_s c_s^2 - \alpha_g \beta) + \left( \alpha_g \frac{\partial \beta}{\partial \alpha_s} + \frac{\rho_s c_s^2}{\alpha_s} - \beta \right) \frac{\rho_g c_g^2 - \rho_s c_s^2 + \beta}{\frac{\rho_s c_s^2}{\alpha_s} + \frac{\rho_s c_s^2}{\alpha_s} + \frac{\partial \beta}{\partial \alpha_s}}. \quad (48) $$

## 3.5 Main characteristics of the model

Finally, we propose a hyperbolic model which is compatible with the following natural model requirements:

- The entropy inequality is verified.

- The microscopic granular energy is decreasing during the relaxation process.

- The elastic properties of a porous matrix can be recovered.

- The properties of pure phases can be recovered.

- The equilibrium sound speed following plastic curves is recovered.
<table><caption>Stiffened-gas parameters</caption><thead><tr><th>Material</th><th>γ</th><th>p<sub>∞</sub> (GPa)</th></tr></thead><tbody><tr><td>Air</td><td>1.4</td><td>0</td></tr><tr><td>Aluminium 2024</td><td>3</td><td>26.4</td></tr></tbody></table>

<table><caption>Cochran–Chan parameters</caption><thead><tr><th>Material</th><th>γ</th><th>ρ<sub>0</sub> (kg m<sup>-3</sup>)</th><th>A<sub>1</sub> (GPa)</th><th>A<sub>2</sub> (GPa)</th><th>E<sub>1</sub></th><th>E<sub>2</sub></th></tr></thead><tbody><tr><td>HMX</td><td>0.93</td><td>1905</td><td>12.87</td><td>13.42</td><td>4.1</td><td>3.1</td></tr></tbody></table>

Table 2: Stiffened-gas and Cochran–Chan parameters for the materials used in the result section 6. The parameters for the HMX are taken from [26,27]. For the aluminium, the parameters are calibrated on the Shock Wave database (http://www.ihed.ras.ru/rusbank/, [23]).

## 4 Hydrodynamic closure for the different materials

250

Herein, we consider the following equations of state for phases k (subscript is omitted for clarity):

- Stiffened gas [24]

$$ \rho e = \frac{p + \gamma p_{\infty}}{\gamma - 1}, \qquad (49) $$

where $\gamma$ and $p_{\infty}$ are some given positive constant parameters, determined by using material reference curves [24].

- Cochran–Chan [5]

255

The Cochran–Chan EOS can be written under a Mie–Gruneisen form

$$ p(\rho, e) = \rho\gamma(e - e_{CC}(\rho)) + p_{CC}(\rho), \qquad (50) $$

with

$$ e_{CC}(\rho) = \frac{A_1}{\rho_0(E_1 - 1)} \left( \left( \frac{\rho}{\rho_0} \right)^{E_1 - 1} - 1 \right) - \frac{A_2}{\rho_0(E_2 - 1)} \left( \left( \frac{\rho}{\rho_0} \right)^{E_2 - 1} - 1 \right), \qquad (51) $$

and

$$ p_{CC}(\rho) = A_1 \left( \frac{\rho}{\rho_0} \right)^{E_1} - A_2 \left( \frac{\rho}{\rho_0} \right)^{E_2}, \qquad (52) $$

and where $\gamma$, $\rho_0$, $A_1$, $A_2$, $E_1$ and $E_2$ are model parameters.

The parameters for the materials used in the result section 6 are given in Table 2.

260

## 5 Numerical solution

### 5.1 Generality

In this section, we propose a Godunov-type scheme [15] for the solution of the compaction multiphase flow. Rewriting the system (32) with

$$ \frac{\partial \alpha_s \rho_s \kappa f}{\partial t} + \nabla \cdot (\alpha_s \rho_s \kappa f \mathbf{u}) = -\pi \frac{d\alpha_s}{dt}, \qquad (53) $$
it reads

$$
\begin{align}
& \frac{\partial \alpha_s}{\partial t} + \mathbf{u} \cdot \nabla \alpha_s = -\mu (p_g - p_s + \pi), \nonumber \\
& \frac{\partial \lambda^e}{\partial t} + \nabla \cdot (\lambda^e \mathbf{u}) = \dot{\lambda}^e, \nonumber \\
& \frac{\partial \alpha_s \rho_s}{\partial t} + \nabla \cdot (\alpha_s \rho_s \mathbf{u}) = 0, \nonumber \\
& \frac{\partial \alpha_g \rho_g}{\partial t} + \nabla \cdot (\alpha_g \rho_g \mathbf{u}) = 0, \nonumber \\
& \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u}^2 + p \mathbf{I}) = \mathbf{0}, \nonumber \\
& \frac{\partial \rho E}{\partial t} + \nabla \cdot ((\rho E + p) \mathbf{u}) = 0, \nonumber \\
& \frac{\partial \alpha_s \rho_s (e_s + \kappa f)}{\partial t} + \nabla \cdot (\alpha_s \rho_s (e_s + \kappa f) \mathbf{u}) + \alpha_s p_s \nabla \cdot \mathbf{u} = -\rho_{s,r} Y_s \kappa \frac{\partial f}{\partial \xi} \dot{\lambda}^e - p_I \frac{d\alpha_s}{dt}, \nonumber \\
& \frac{\partial \alpha_g \rho_g e_g}{\partial t} + \nabla \cdot (\alpha_g \rho_g e_g \mathbf{u}) + \alpha_g p_g \nabla \cdot \mathbf{u} = p_I \frac{d\alpha_s}{dt}. \tag{54}
\end{align}
$$

The thermodynamic closure is achieved using the EOS presented in Section 4.

The numerical scheme to solve this system is an extension of [36] approach with finite relaxations and is described in the following. Furthermore, since the main interest of the paper is the model, we only propose a first-order scheme to solve the system and we invite the reader interested by the modelling to jump to the result section 6. The solution is obtained through three steps (executed at each time-step stage if a high-order method is used):

- Solution of the non-relaxed, hyperbolic equations ($\mu = \dot{\lambda}^e = 0$) (Section 5.2).

- Correction of the energy of each phase to ensure conservation of the mixture total energy (Section 5.2.5).

- Finite relaxation of the compaction variable $\lambda^e$ and volume fraction $\alpha_s$ ($(\mu, \dot{\lambda}^e) \neq 0$) ensuring a slow convergence to a mechanical equilibrium (Section 5.3).

These steps are detailed hereafter.

## 5.2 Hyperbolic step

We present here the first-order Godunov method [15]. The system to solve is

$$
\begin{align}
& \frac{\partial \alpha_s}{\partial t} + \mathbf{u} \cdot \nabla \alpha_s = 0, \nonumber \\
& \frac{\partial \lambda^e}{\partial t} + \nabla \cdot (\lambda^e \mathbf{u}) = 0, \nonumber \\
& \frac{\partial \alpha_s \rho_s}{\partial t} + \nabla \cdot (\alpha_s \rho_s \mathbf{u}) = 0, \nonumber \\
& \frac{\partial \alpha_g \rho_g}{\partial t} + \nabla \cdot (\alpha_g \rho_g \mathbf{u}) = 0, \nonumber \\
& \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u}^2 + p \mathbf{I}) = \mathbf{0}, \nonumber \\
& \frac{\partial \rho E}{\partial t} + \nabla \cdot ((\rho E + p) \mathbf{u}) = 0, \nonumber \\
& \frac{\partial \alpha_s \rho_s (e_s + \kappa f)}{\partial t} + \nabla \cdot (\alpha_s \rho_s (e_s + \kappa f) \mathbf{u}) + \alpha_s p_s \nabla \cdot \mathbf{u} = 0, \nonumber \\
& \frac{\partial \alpha_g \rho_g e_g}{\partial t} + \nabla \cdot (\alpha_g \rho_g e_g \mathbf{u}) + \alpha_g p_g \nabla \cdot \mathbf{u} = 0. \tag{55}
\end{align}
$$

The characteristic velocities expressed in one dimension are

$$ u \text{ and } u \pm \sqrt{\sum_k Y_k c_k^2}. \tag{56} $$
We use a HLLC approximate Riemann solver [39]. It preserves positivity of the volume fraction and the density, and is able to deal with strong shocks. With this solver, each wave is considered as a discontinuity and consequently jump relations are needed. The system being non-conservative, shock relations are non-conventional.

In the following, the jump relations, the Riemann solver and the Godunov scheme are presented in one dimension.

### 5.2.1 Jump relations for conservative variables

The hyperbolic system is similar to the fluid-fluid system in [36]. The main difference is on the energy. The jump relations for the mass across a discontinuity (other than interface) at velocity $D$ are

$$ \alpha_k \rho_k (u - D) = \alpha_k^0 \rho_k^0 (u^0 - D) = m_k. \qquad (57) $$

The variables with superscript “0” denote the unshocked state. With $p = \sum_k \alpha_k p_k$, $\rho = \sum_k \alpha_k \rho_k$ and $m = \sum_k m_k$ the mixture pressure, density and mass flux, respectively, the momentum jump relation reads

$$ p - p^0 + m^2 \left( \frac{1}{\rho} - \frac{1}{\rho_0} \right). \qquad (58) $$

The jump of the compaction variable is

$$ \lambda^e(u - D) = \lambda^{e,0}(u^0 - D). \qquad (59) $$

The mixture-total-energy jump relation is

$$ mE + pu = mE^0 + p^0 u^0. \qquad (60) $$

The generalized Hugoniot curve reads

$$ e - e^0 + Y_s \kappa f - Y_s^0 \kappa^0 f^0 = \frac{1}{2}(p + p^0) \left( \frac{1}{\rho^0} - \frac{1}{\rho} \right), \qquad (61) $$

with $e = \sum_k Y_k e_k$.

### 5.2.2 Jump relations for non-conservative variables

The volume fraction is continuous through the jump across a discontinuity (other than interface)

$$ \alpha_k = \alpha_k^0. \qquad (62) $$

The internal-energy equations being in non-conservative form, they are not adapted for the determination of jump relations. A path has to be defined to circumvent this difficulty. Since the results depend on the path, these shock relations are only used to predict the total-energy repartition. A correction is made afterwards to preserve the total-energy evolution of the system. The convergence of this method, different from [36], is highlighted in the result section (Section 6). The following jump relations are proposed

$$ \begin{aligned} e_s - e_s^0 + \kappa f - \kappa^0 f^0 &= \frac{1}{2}(p_s + p_s^0) \left( \frac{1}{\rho_s^0} - \frac{1}{\rho} \right), \\ e_g - e_g^0 &= \frac{1}{2}(p_g + p_g^0) \left( \frac{1}{\rho_g^0} - \frac{1}{\rho} \right). \end{aligned} \qquad (63) $$

This choice generalizes classical additive principle for heterogeneous mixtures (see [33, 28, 13, 11, 40] for details) and is compatible with the following requirements:

- Single-phase limit: When a phase vanishes, the preceding jump relations correspond to the jump relations for the total energy.

- Total-energy conservation: By multiplying each equations of (63) by $Y_k$ and by summing them, we obtain Eq. (61).

With the help of these relations, the approximate Riemann solver can be built as developed hereafter.
### 5.2.3 Interface relations

The set of interface relations are

$$ 
\begin{aligned}
u &= u^0, \\
p &= p^0.
\end{aligned} 
\qquad (64) $$

### 5.2.4 HLLC Riemann solver

Consider a cell boundary with a left state (L) and a right state (R). The left- and right-facing wave speeds are obtained following Davis' estimates [39]

$$ S_R = \max(u_L + c_L, u_R + c_R), \quad S_L = \min(u_L - c_L, u_R - c_R). \qquad (65) $$

The speed of the intermediate wave (also called contact discontinuity) is estimated under the HLL approximation [39]

$$ S_M = \frac{(\rho u^2 + p)_L - (\rho u^2 + p)_R - S_L(\rho u)_L + S_R(\rho u)_R}{(\rho u)_L - (\rho u)_R - S_L \rho_L + S_R \rho_R} = u^*. \qquad (66) $$

From these wave speeds, conservative state variables in the star region (see Figure 2) are determined

$$ (\alpha_k \rho_k)_* = (\alpha_k \rho_k)_{L,R} \frac{S_{L,R} - u_{L,R}}{S_{L,R} - u^*}, \qquad (67) $$

$$ p^* = p_R + \rho_R (u_R - S_R)(u_R - u^*), \qquad (68) $$

The compaction variables are obtained using

$$ \lambda_{L,R}^{e,*} = \lambda_{L,R}^{e} \frac{S_{L,R} - u_{L,R}}{S_{L,R} - u^*}. \qquad (69) $$

The right and left mixture total energies are

$$ E_{L,R}^* = E_{L,R} + \frac{p_{L,R} u_{L,R} - p^* u^*}{\rho_{L,R}^* (u^* - S_{L,R})}. \qquad (70) $$

The jumps of the internal energies are obtained with the help of the Hugoniot Eq. (63). Indeed, the phase pressure $p_{k,L,R}^*$ are calculated as functions of the phase densities $\rho_{k,L,R}^*$ and of the compaction variable $\lambda_{L,R}^*$. The internal energies are therefore determined using EOS relations

$$ e_{k,L,R}^* = e_{k,L,R}^* (\rho_{k,L,R}^* p_{k,L,R}^*). \qquad (71) $$

With the help of this approximate solver, it is then possible to derive a Godunov-type scheme.

### 5.2.5 Godunov-type scheme combined with the correction of the energy of each phase

The conservative part of the system, in absence of relaxation terms, is updated with the conventional Godunov scheme [15]

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x} (\mathbf{F}^*(\mathbf{U}_i^n, \mathbf{U}_{i+1}^n) - \mathbf{F}^*(\mathbf{U}_{i-1}^n, \mathbf{U}_i^n)), \qquad (72) $$

with

$$ \mathbf{U} = (\lambda^e, \alpha_s \rho_s, \alpha_g \rho_g, \rho u, \rho E)^T, \qquad (73) $$

and

$$ \mathbf{F} = (\lambda^e u, \alpha_s \rho_s u, \alpha_g \rho_g u, \rho u^2 + p, \rho(E + p)u)^T. \qquad (74) $$

The previous solver and scheme for the conservative variables is similar to the one presented in [36]. However, for the energy of each phase, we consider a different scheme. The complete scheme is therefore:

- Evolve the conservative equations following the precedent scheme (72).
<figure><img src="image_3.png" /><figcaption>Fig. 2: HLLC approximate Riemann solver. Solutions in the “star” region consist of two constant states separated from each other by an intermediate wave of speed $S_M$.</figcaption></figure>

- Evolve the non-conservative equations with the following schemes. For the volume fraction, the scheme is

$$ \alpha_{s,i}^{n+1} = \alpha_{s,i}^{n} - \frac{\Delta t}{\Delta x} \left( (\alpha_s u)_{i+\frac{1}{2}}^{*} - (\alpha_s u)_{i-\frac{1}{2}}^{*} - \alpha_{s,i}^{n} (u_{i+\frac{1}{2}}^{*} - u_{i-\frac{1}{2}}^{*}) \right). \quad (75) $$

For the internal energies, the original scheme is

$$ (\alpha_k \rho_k (e_k + \kappa f))_i^{n+1} = (\alpha_k \rho_k (e_k + \kappa f))_i^n - \Delta \hat{e}_{k,i}, \quad (76) $$

where

$$ \Delta \hat{e}_{k,i} = \frac{\Delta t}{\Delta x} \left( (\alpha_k \rho_k (e_k + \kappa f) u)_{i+\frac{1}{2}}^* - (\alpha_k \rho_k (e_k + \kappa f) u)_{i-\frac{1}{2}}^* + (\alpha_k \rho_k)_i^n (u_{i+\frac{1}{2}}^* - u_{i-\frac{1}{2}}^*) \right). \quad (77) $$

One can note that the aim here is only to obtain an evaluation of the repartition of the energy.

- Evaluate the variation of the internal energy ($\sum_k \rho Y_k e_k + Y_s \kappa f$) from the total energy

330

$$ \Delta \tilde{e}_i = \left( \rho E - \frac{\rho u^2}{2} \right)_i^{n+1} - \left( \rho E - \frac{\rho u^2}{2} \right)_i^n. \quad (78) $$

- Evolve the internal energy with the following scheme

$$ (\alpha_k \rho_k (e_k + \kappa f))_i^{n+1} = (\alpha_k \rho_k (e_k + \kappa f))_i^n - \frac{\Delta \hat{e}_{k,i}}{\sum_k \Delta \hat{e}_{k,i}} * \Delta \tilde{e}_i. \quad (79) $$

With such a scheme, one can note that the equations of the internal energies are compatible with the total energy. This step will allows to deal with finite relaxation time for the volume fraction since the updated step used in [36] is no more necessary.

335
## 5.3 Relaxation step

The system to solve reads

$$ 
\begin{aligned}
\frac{\partial \alpha_s}{\partial t} &= -\mu (p_g - p_s + \pi), \\
\frac{\partial \lambda^e}{\partial t} &= \dot{\lambda}^e, \\
\frac{\partial \alpha_k \rho_k}{\partial t} &= 0, \\
\frac{\partial \rho u}{\partial t} &= 0, \\
\frac{\partial \rho E}{\partial t} &= 0, \\
\frac{\partial \alpha_s \rho_s (e_s + \kappa f)}{\partial t} &= -Y_s \kappa \frac{\partial f}{\partial \xi} \rho_{s,r} \dot{\lambda}^e + p_I \mu (p_g - p_s + \pi), \\
\frac{\partial \alpha_g \rho_g e_g}{\partial t} &= -p_I \mu (p_g - p_s + \pi),
\end{aligned}
\quad (80) $$

This system of ordinary differential equations (ODE) can be solved using your favourite ODE solver. In the following, we use a classic, first-order, explicit, Euler scheme with time-step subdivisions. The number of subdivisions is adapted at each time step to verify the following condition $\mu |p_g - p_s + \pi| \Delta t_{Euler} < 10^{-3} \alpha_s \alpha_g$.

# 6 Results

## 6.1 Quasi-static loading-unloading cycle on granular HMX

We consider the test presented in [21]. A piston presses a HMX powder sample at 1 m s⁻¹. The characteristic times $K_y$ and $1/\mu$ are sufficiently small to consider the relaxation time effects to be negligible. Loading and unloading of the material is considered. For the reversible part of the compression, we consider $\chi = 7$ GPa, $a_\kappa = 1.8 \times 10^5$ and $n_\kappa = 0.5$. For the plastic part, we consider:

- First, the theoretical Gurson law (Eq. 42) $\beta^G$ with $2\sigma_Y/3q_2 = 50$ MPa, $q_3 = 8$ and $q_1 = 2.85$.

- Second, the plastic limit (Eq. 38) $\beta^K$ used in [33, 12] where the parameters for the yield limit are $a = 3$ kPa and $n_p = 1.04$.

The elastic parameters are fitted on the last loading-unloading curve using the linearized slope. The plastic parameters are calibrated on the plastic curve of the experimental results.

The Cochran-Chan EOS parameters of pure HMX are given in Table 2.

Numerical results are compared with experimental data of [21] on Figure 3 (a) and (b), respectively. The agreement for the Gurson limit is good but a discrepancy is observable which may come from the fact that the material is not homogeneous (binder+ HMX) and the pores are not ellipses. Moreover, this fit is unable to get the beginning of the curve because the compaction pressure is 0 Pa for $\alpha_s = 0.6875$. While the agreement for the Kapila limit is very good and catches the behaviour of the low compaction pressure.

These tests highlight the ability of the model to deal with quasi-static compaction of granular material with different yield limits. Moreover, even if the elastic part of the model is calibrated using the last loading-unloading zone, it is able to capture, for different porosities, the other elastic loading-unloading path.

## 6.2 Wave and interface dynamics

Let us consider a shock tube of 1 cm length involving a high-pressure chamber on the left, filled with air at $p = 100$ MPa, and a low-pressure chamber on the right, filled with HMX at atmospheric pressure. The air is considered as a perfect gas while the HMX is computed with the Cochran-Chan EOS (Table 2). For the yield limit, Eq. (38) for $\beta^K$ is used with similar parameters than in Section 6.1. The viscosity parameter
<figure><img src="image_4.png" /><figcaption>Fig. 3: Quasi-static loading-unloading curve for granular HMX. Numerical results (thin, purple line) using the (a) Gurson and (b) Kapila limits, compared with experimental results (thick, green line) [21].</figcaption></figure>

is set to $K_y = 65 \text{ Pa s}^{-1}$ and the relaxation constant is set to $\mu = 8 \text{ s Pa}^{-1}$. The influence of these two parameters are analysed in section 6.3. The initial discontinuity is located at coordinate $x = 6$ mm. The solid volume fraction in the left chamber is set to $\alpha_s = 10^{-5}$ and in the right chamber to $\alpha_s = 0.8$. The initial density of the solid is constant in the whole domain and is equal to $\rho_s = 1900 \text{ kg m}^{-3}$. The gas-phase density in the high-pressure chamber is set to $100 \text{ kg m}^{-3}$, while in the low-pressure chamber it is set to $1 \text{ kg m}^{-3}$. Both initial states are initially at rest.

The aim of this test case is to show the ability of the model and the numerical scheme to compute a shock propagation in a granular material (capturing an elastic-plastic shock wave) and his associated compaction, the rarefaction in the gas and the interface evolution separating a nearly pure gas from a granular mixture.

Results are presented on Figure 4 at time $t = 4.6 \text{ µs}$, on a mesh with 1000 cells and with CFL = 0.5. First, an elastic shock wave propagates in the granular media. This shock wave is followed by a plastic shock wave. The rarefaction wave in the gas is well captured. There is no oscillation of the mixture pressure and of the velocity at the interface.

We also compute the error (L2 norm) for the results obtained with different mesh sizes with a reference solution obtained using a very refined mesh ($4 \times 10^4$ cells). The curve is presented on Figure 5 and we observe a convergence order of approximately 0.6.

We now consider a (non-physical) stronger test case to highlight the ability of the model and the numerical method to deal with extremely high pressures and therefore very strong shock. The pressure and air density in the left chamber are now set to 100 GPa and $\rho_g = 1000 \text{ kg m}^{-3}$, respectively. Results are presented on Figure 6 at time $t = 0.327 \text{ µs}$ and on a mesh with 1000 cells. The shock wave propagates in the porous HMX. The compaction effect can only be seen on the magnified plot on the pressure. The interface condition $u = \text{const}$ and $p = \text{const}$ are preserved and the rarefaction waves in the gas are well captured.

## 6.3 Analysis of the thickness of the plastic shock

We consider the low-amplitude shock tube studied previously and we vary separately the viscosity and relaxation parameters. A sensor is also set at coordinate $x = 7.5$ mm.

On Figure 7, results obtained with $K_y = 325 \text{ Pa s}^{-1}$, $65 \text{ Pa s}^{-1}$, $6.5 \text{ Pa s}^{-1}$ and $0.65 \text{ Pa s}^{-1}$ are presented ($\mu = 8$). One can observe for higher viscosity coefficient, the amplitude of the first shock and the plastic-shock layer increase. For small values, the results are almost superposed. Indeed, the viscous dissipation can be neglected compared to other sources of dissipation (pressure relaxation coefficient $\mu$ and numerical dissipation).

On Figure 8, results obtained with $\mu = 5 \text{ s Pa}^{-1}$, $0.5 \text{ s Pa}^{-1}$, $0.05 \text{ s Pa}^{-1}$ and $0.005 \text{ s Pa}^{-1}$ are presented ($K_y = 65 \text{ Pa s}^{-1}$). When $\mu$ is decreasing, the elastic shock is faster and with a higher amplitude. Though, the elastic and plastic shocks are more diffused. This result is similar to the effect of the plastic viscosity. Although, the relaxation parameter $\mu$ will diffuse the shock oppositely to the viscous coefficient $K_y$ which
<figure><img src="image_5.png" /><figcaption>Fig. 4: Low-amplitude shock tube with air–HMX interface.</figcaption></figure>

has no influence on the elastic shock layer. Furthermore, if $\mu$ is sufficiently large, results are superposed. Again, the dissipation due to the relaxation parameter can be neglected compared to other sources of dissipation (numerical and viscous dissipation).

## 6.4 Hugoniot of aluminium

We consider testing the model and method on an impact test problem with various impact velocities and initial porosities for porous aluminium 2024. The domain is initially cut in two sections both filled with the porous aluminium: On the left, a velocity is imposed and directed to the right, and on the right, a null velocity. The domain is also large enough so the waves don't reach the boundaries. The pores are filled with air and we measure the shock velocity and the resulting pressure for a converged mesh. From these tests, the piston velocity is half the impact velocity (from initial left section).

The stiffened-gas EOS is used for the aluminium hydrodynamic behaviour. For its plastic behaviour, the Gurson limit is considered. Air and aluminium EOS parameters are given in Table 2. Further, elastic and plastic parameters for aluminium are given in Table 3 and its initial density is $\rho_0 = 2784 \text{ kg m}^{-3}$.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Convergence test case for the low-amplitude shock tube with air-HMX interface presented on Figure 4. The error is plotted as function of the mesh size.</figcaption></figure>

<table><thead><tr><th colspan="3">Elastic parameters</th></tr><tr><th>χ (GPa)</th><th>n<sub>κ</sub></th><th>a<sub>κ</sub></th></tr></thead><tbody><tr><td>7</td><td>0.5</td><td>1.8 × 10<sup>4</sup></td></tr></tbody></table>

<table><caption>Table 3: Elastic and plastic parameters for aluminium 2024.</caption><thead><tr><th colspan="7">Plastic parameters</th></tr><tr><th>K<sub>y</sub> (Pa s<sup>-1</sup>)</th><th>μ (s Pa<sup>-1</sup>)</th><th>n<sub>λ<sup>e</sup></sub></th><th>σ<sub>Y</sub> (MPa)</th><th>q<sub>1</sub></th><th>q<sub>2</sub></th><th>q<sub>3</sub></th></tr></thead><tbody><tr><td>65</td><td>1 × 10<sup>-5</sup></td><td>1</td><td>360</td><td>1</td><td>1</td><td>1</td></tr></tbody></table>

Note that the parameters are fitted on the pure-material Hugoniot from the Shock Wave Database (http://www.ihed.ras.ru/rusbank). One should also note that for such types of tests, the elastic and plastic parameters have little influence on the results. Indeed, the hydrodynamic part is mainly dominant. Therefore, the elastic and plastic parameters should only be chosen as reasonable.

On Figure 9, the numerical results are compared with the experimental results [23] for different porosities. The agreement is almost perfect for piston velocity lower than 3.7 km s⁻¹ (corresponding to an impact at 7.4 km s⁻¹ on a fixed target). For higher velocities, the error on the shock velocity is less than 10%, which is acceptable. The small discrepancies on the pressure come from the stiffened-gas EOS which can only fit on a small range of parameters. The use of Mie–Gruneisen-type EOS may improve the agreement for stronger shocks. This test highlights the ability of the model to capture the behaviour of the porous material with a single set of parameters for a very large range of shocks.

## 6.5 Spherical-tension test case

We consider here a spherical test case where loading-unloading traction (hydrodynamic stress tensor $\sigma = -pI$) on aluminium with two different initial porosities ($\alpha_{g0} = 10^{-4}$ and $10^{-6}$). This test is done with a unique one-dimensional cell to avoid classical localization problems during tensile test. The system can be considered as quasi-static since the strain $\varepsilon$ evolves with $\mathrm{tr}(\dot{\varepsilon}) \approx \mathrm{div}(\mathbf{u}) = \pm 0.1 \, \mathrm{s}^{-1}$: One boundary has a switching imposed velocity while the other is fixed ($\mathbf{u} = 0$). We choose the relaxation parameters ($K_y$ and $\mu$) very small and because the test case can be considered as quasi-static, they have negligible influence on the results.

The hydrodynamic tensile stress, corresponding in the classical elasticity theory to $p = -\mathrm{tr}(\sigma)/3$, is presented on Figure 10 and compared with theoretical elastic curve for linear elasticity, $p = K\mathrm{tr}(\varepsilon)$, where $K$ is the compressibility modulus (Table 1). One can observe a first elastic zone until a volume variation of $\mathrm{tr}(\varepsilon) = \int_0^t \nabla \cdot (\mathbf{u}) dt = 1.5\%$ where both materials behave as the theoretical elastic law. Then plasticity is taken under account. During this step, unloading occurs and irreversible deformation can be observed. Further, the small variation on the initial volume fraction of gas will only induce different
<figure><img src="image_7.png" /><figcaption>Fig. 6: High-amplitude shock tube with air–HMX interface.</figcaption></figure>

yield limit for the material. This property might be used with a probabilistic repartition of the volume
fraction of voids to avoid numerical localization.

6.6 Comparison with a micro-mechanical modelling for pore collapse

In this section, the aim is to compare the results of our modelling during pore collapse with the micro-mechanical modelling described in [26] where a spherical domain contains one pore. The model is recalled
<figure><img src="image_8.png" /><figcaption>Fig. 7: Influence of the viscosity parameter $K_y$ on the shock structure for the low-amplitude shock tube of Section 6.2. Results presented at the sensor location $x = 7.5$ mm for $K_y = 325$ Pa s⁻¹ (very thick, purple line), 65 Pa s⁻¹ (thick, green line), 6.5 Pa s⁻¹ (thin, blue line) and 0.65 Pa s⁻¹ (very thin, yellow line).</figcaption></figure>

hereafter

$$ 
\begin{aligned}
\frac{da}{dt} &= \xi, \\
\frac{d\xi}{dt} &= \frac{1}{\rho_s a (1 - \phi^{1/3})} \left( p_g - p_s - 4K \frac{\xi}{a} (1 - \phi) + 2Y_s \text{sign}(\xi) \ln\left(\frac{a}{b}\right) - \frac{1}{2} \rho_s \xi^2 (3 - \phi^{1/3} + \phi^{4/3}) \right), \\
\frac{db}{dt} &= \left(\frac{a}{b}\right)^2 \xi, \\
\frac{d\rho_g}{dt} &= -3 \frac{\rho_g}{a} \xi, \\
\frac{dp_g}{dt} &= p_{g0} \gamma \rho^{\gamma-1} \frac{d\rho_g}{dt},
\end{aligned}
\quad (81) $$

where $a$ is the radius of the pore, $b$ is the radius of the pore control volume, $\phi = (a/b)^3$ is the porosity, $Y$ is the plastic limit of the pure HMX and $K$ is the viscosity. One should note that we neglect the heat and mass transfer which was present in the original version.

To match with the hypothesis made in [26], the pressure in HMX is imposed during all the relaxation process to $p = 1.5$ GPa. The density is set to $\rho_s = 2150$ kg m⁻³ with a volume fraction of pore $\alpha_g = 0.017$ at atmospheric pressure. We then compute the evolution of the different pressure and volume fractions of pore. The Cochran-Chan EOS is used for the HMX with the parameters in Table 2. To be able to compare the results with those of [26], we use the following plastic limit $\beta = -\sigma_y \ln(\alpha_g)$ with $\sigma_y = 2$ GPa. The relaxation time parameter used is $\mu = 0.011$ and $K = 65$ Pa s⁻¹.

One can observe on Figure 11 that the agreement between the micro-mechanical modelling and the new compaction modelling is almost perfect with only one fitting parameter ($\mu$) and without computing the whole partial differential system.
<figure><img src="image_9.png" /><figcaption>Fig. 8: Influence of the relaxation parameter $\mu$ on the shock structure for the low-amplitude shock tube of Section 6.2. Results presented at the sensor location $x = 7.5$ mm for $\mu = 0.005 \text{ s Pa}^{-1}$ (very thick, purple line), $0.05 \text{ s Pa}^{-1}$ (thick, yellow line), $0.5 \text{ s Pa}^{-1}$ (thin, blue line) and $5 \text{ s Pa}^{-1}$ (very thin, green line).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Shock velocity and pressure as function of the piston velocity. The numerical results are represented with lines. The experimental results from [23] are represented for different porosity $\alpha_s = 0.9999, 0.911, 0.8, 0.7$ and $0.59$ denoted by $+$, $\times$, $*$, $\square$ and $\blacksquare$, respectively.</figcaption></figure>

## 7 Conclusion

A multiphase flow model for the numerical simulation of irreversible compaction of porous material has been built and validated. It is able to reproduce with a good accuracy loading-unloading cycles as well as interface dynamics separating fluids and porous material in the case of isotropic traction and compression. The model can't solve problems where shear is dominant. The shear effect will be introduced in future work using for example the hyperelastic approach developed in [11]. Dispersive effects (micro-inertia and localization) must also be introduced. In many applications, the porosity can't
<figure><img src="image_11.png" /></figure>

Fig. 10: Evolution of the hydrodynamic tensile stress as a function of the volume strain tr(ε). The linear classical elastic curve (thick, blue line) and the numerical results taking into account the void increasing during plasticity phenomenon (thin, purple and normal, green lines for $\alpha_{g0} = 10^{-4}$ and $10^{-6}$, respectively) are represented.

<figure><img src="image_12.png" /></figure>

Fig. 11: Evolution of the pressure of the gas (left) and the normalized radius of the pore (right). Numerical results of the new modelling (solid line) and the micro-mechanical modelling (dashline) are represented.

be considered as closed. The extension to multi-velocity model will be done with the principle proposed in [18] and validated against impact experiments of [19]. Further, this model will be implemented in the public version of ECOGEN [38] ; an open-source CFD platform for numerical simulation of compressible multiphase flows.

470

**Declaration of competing interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Acknowledgment**

475

This work was partially supported by A*MIDEX and ANR; grants ANR-11-LABX-0092, ANR-11-IDEX-0001-02 and ANR-ASTRID project SNIP ANR-19-ASTR-0016-01.
References

1. R. Abgrall and S. Karni. Computations of compressible multifluids. *J. Comp. Phys.*, 169(2):594-623, 2001.

2. M. R. Baer and J. W. Nunziato. A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials. *Int. J. of Multiphase Flows*, 12:861-889, 1986.

3. D. J. Benson, V. F. Nesterenko, F. Jonsdottir, and M. A. Meyers. Quasistatic and dynamic regimes of granular material deformation under impulse loading. *J. Mech. Phys. Solids*, 45(11-12):1955-1999, 1997.

4. J. B. Bzil, R. Menikoff, S. F. Son, A. K. Kapila, and D. S. Steward. Two-phase modeling of a deflagration-to detonation transition in granular materials: A critical exmination of modeling issues. *Phys. Fluids*, 11(2):378-402, 1999.

5. G. Cochran and J. Chan. Shock initiation and detonation models in one and two dimensions. *Lawrence National Laboratory Report, CID-18024*, 1979.

6. F. Dell'Isola and S. Gavrilyuk. *Variational models and methods in solid and fluid mechanics*, volume 535. Springer Science & Business Media, 2012.

7. P. Embid and M. Baer. Mathematical analysis of a two-phase model for reactive granular material. *Continuum Mech. Thermodyn.*, 4:279, 1992.

8. A. C. Eringen. Theory of micropolar elasticity. In *Microcontinuum field theories*, pages 101-248. Springer, 1999.

9. N. Favrie and S. Gavrilyuk. Dynamics of shock waves in elastic-plastic solids. In *ESAIM: Proceedings*, volume 33, pages 50-67. EDP Sciences, 2011.

10. N. Favrie and S. Gavrilyuk. Mathematical and numerical model for nonlinear viscoplasticity. *Philoso. Trans. Roy. Soc. A: Math., Phys. Eng. Sc.*, 369(1947):2864-2880, 2011.

11. N. Favrie and S. Gavrilyuk. Diffuse interface model for compressible fluid-compressible elastic-plastic solid interaction. *J. Comp. Phys.*, 231(7):2695-2723, 2012.

12. N. Favrie and S. Gavrilyuk. Dynamic compaction of granular materials. *Proc. R. Soc. A: Math., Phys. Eng. Sc.*, 469(2160):20130214, 2013.

13. N. Favrie, S. L. Gavrilyuk, and R. Saurel. Solid-fluid diffuse interface model in cases of extreme deformations. *J. Comp. Phys.*, 228(16):6037-6077, 2009.

14. D. T. Gethin, V. D. Tran, R. W. Lewis, and A. K. Griffin. An investigation of powder compaction processes. *Int. J. Powder Metall.*, 30(4):385-398, 1994.

15. S. K. Godunov. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. *Mat. Sb.*, 89(3):271-306, 1959.

16. K. A. Gonthier. Modeling and analysis of reactive compaction for granular energetic solids. *Combust. Sci. Technol.*, 175(9):1679-1709, 2003.

17. A. L. Gurson. Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media. *J. Eng. Mater. Technol.*, 99(1):2-15, 1977.

18. S. Hank, N. Favrie, and J. Massoni. Modeling hyperelasticity in non-equilibrium multiphase flows. *J. Comp. Phys.*, 330:65-91, 2017.

19. S. Hank, S. Gavrilyuk, N. Favrie, and J. Massoni. Impact simulation by an Eulerian model for interaction of multiple elastic-plastic solids and fluids. *Int. J. Impact Eng.*, 109:104-111, 2017.

20. Z. Hashin and S. Shtrikman. A variational approach to the theory of the elastic behaviour of multiphase materials. *J. Mech. Phys. Solids*, 11(2):127-140, 1963.

21. V. Jogi. Predictions for multi-scale shock heating of a granular energetic material. *MS thesis: Louisiana State University, USA*, 2003.

22. A. Kapila, R. Menikoff, J. Bdzil, S. Son, and D. Stewart. Two-phase modeling of DDT in granular materials: Reduced equations. *Phys. Fluids*, 13:3002-3024, 2001.

23. V. Kogan, P. Levashov, and I. Lomov. Shock Wave Database, 2003. http://www.ihed.ras.ru/rusbank/.

24. O. Le Métayer, J. Massoni, and R. Saurel. Elaborating equations of state of a liquid and its vapor for two-phase flow models. *Int. J. Therm. Sci.*, 43:265-276, 2004.

25. C. L. Martin, D. Bouvard, and S. Shima. Study of particle rearrangement during powder compaction by the discrete element method. *J. Mech. Phys. Solids*, 51(4):667-693, 2003.

26. J. Massoni. *Un modèle micromécanique pour l'initiation par choc et la transition vers la détonation dans les matériaux solides hautement énergétiques*. PhD thesis, Aix-Marseille 1, 1999.

27. J. Massoni, R. Saurel, G. Baudin, and G. Demol. A mechanistic model for shock initiation of solid explosives. *Phys. Fluids*, 11(3):710-736, 1999.

28. S. Ndanou, N. Favrie, and S. Gavrilyuk. Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation. *J. Comp. Phys.*, 295:523-555, 2015.

29. J. Oliver, S. Oller, and J. Cante. A plasticity model for simulation of industrial powder compaction processes. *Int. J. Solids Structures*, 33(20-22):3161-3178, 1996.

30. S. L. Passman, J. W. Nunziato, and E. K. Walsh. A theory of multiphase mixtures. In *Rational thermodynamics*, pages 286-325. Springer, 1984.

31. I. Peshkov and E. Romenski. A hyperbolic model for viscous newtonian flows. *Continuum Mech. Thermodyn.*, 28(1-2):85-104, 2016.

32. G. Ragnarsson, G. Alderborn, and C. Nyström. Pharmaceutical powder compaction technology. *Pharmaceutical powder compaction technology. New York: Marcel Dekker, Inc*, pages 77-97, 1996.

33. R. Saurel, N. Favrie, F. Petitpas, M. H. Lallemand, and S. L. Gavrilyuk. Modelling dynamic and irreversible powder compaction. *J. Fluid Mech.*, 664:348, 2010.

34. R. Saurel, S. L. Gavrilyuk, and F. Renaud. A multiphase model with internal degrees of freedom: Application to shock-bubble interaction. *J. Fluid Mech.*, 495:283-321, 2003.

35. R. Saurel, F. Petitpas, and R. Abgrall. Modelling phase transition in metastable liquids: application to cavitating and flashing flows. *J. Fluid Mech.*, 607:313-350, 2008.

36. R. Saurel, F. Petitpas, and R. Berry. Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. *J. Comp. Phys.*, 228(5):1678-1712, 2009.
37. K. Schmidmayer, F. Petitpas, E. Daniel, N. Favrie, and S. Gavrilyuk. A model and numerical method for compressible flows with capillary effects. J. Comp. Phys., 334:468–496, 2017.

38. K. Schmidmayer, F. Petitpas, S. Le Martelot, and E. Daniel. ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows. Comp. Phys. Com., 251:107093, 2020.

39. E. Toro. Riemann solvers and numerical methods for fluid dynamics. Springer Verlag, Berlin, 1997.

40. R. Trunin. Shock compression of condensed materials (laboratory studies). Physics Uspekhi, 44(4):371–396, 2001.

41. V. Tvergaard and A. Needleman. Analysis of the cup-cone fracture in a round tensile bar. Acta metallurgica, 32(1):157–169, 1984.

42. G. B. Wallis. One-dimensional two-phase flow. McGraw-Hill, 1969.

43. A. B. Wood. A textbook of sound. G. Bell and Sons LTD, London, 1930.

44. Z. L. Zhang, C. Thaulow, and J. Ødegård. A complete Gurson model approach for ductile fracture. Eng. Fracture Mech., 67(2):155–168, 2000.

## A Equilibrium condition of the reversible model

In order to find the pressure-equilibrium relation let's consider the derivation of the internal energy using the entropy equation $\frac{d\eta_k}{dt} = 0$

$$ 
\begin{aligned}
& \rho \frac{dY_s \kappa(\alpha_s) f(\xi) + \sum Y_k e_k(\rho_k, \eta_k)}{dt} \\
&= \rho Y_s \kappa \frac{\partial f \rho_{s,r} \lambda}{\partial \xi \alpha_s \rho_s} \frac{d\alpha_s}{dt} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \frac{d\alpha_s}{dt} + \rho \sum_k Y_k \frac{\partial e_k}{\partial \rho_k} \frac{d\rho_k}{dt} + \sum_k \rho Y_k \frac{\partial e_k}{\partial \eta_k} \frac{d\eta_k}{dt} \\
&= -\sum_k \alpha_k p_k \nabla \cdot \mathbf{u} + \left( \rho_{s,r} \kappa \lambda \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \frac{d\alpha_s}{dt} + \sum_k \left( -p_k \frac{d\alpha_k}{dt} \right) \\
&= -\sum_k \alpha_k p_k \nabla \cdot \mathbf{u} + \left( p_g - p_s + \rho_{s,r} \kappa \lambda \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \frac{d\alpha_s}{dt}.
\end{aligned}
\quad (82) $$

Using the momentum equation we obtain

$$ \frac{\partial \rho E}{\partial t} + \text{div} \left( \rho E \mathbf{u} + \sum_k \alpha_k p_k \right) = \left( p_g - p_s + \rho_{s,r} \kappa \lambda \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) \frac{d\alpha_s}{dt}. \quad (83) $$

The total energy is conserved if $\frac{d\alpha_s}{dt} = 0$ or if $\left( p_g - p_s + \rho_{s,r} \kappa \lambda \frac{\partial f}{\partial \xi} + \rho Y_s \frac{\partial \kappa}{\partial \alpha_s} f \right) = 0$. This condition corresponds to the classical condition in the multiphase-flow community $p_g = p_s$. Note that one could also obtain such condition using Hamilton's principle [6].