## Abstract

An Eulerian hyperbolic multiphase flow model for dynamic and irreversible compaction of granular materials is constructed. The reversible model is first constructed on the basis of the classical Hertz theory. The irreversible model is then derived in accordance with the following two basic principles. First, the entropy inequality is satisfied by the model. Second, the corresponding ‘intergranular stress’ coming from elastic energy owing to contact between grains decreases in time (the granular media behave as Maxwell-type materials). The irreversible model admits an equilibrium state corresponding to von Mises-type yield limit. The yield limit depends on the volume fraction of the solid. The sound velocity at the yield surface is smaller than that in the reversible model. The last one is smaller than the sound velocity in the irreversible model. Such an embedded model structure assures a thermodynamically correct formulation of the model of granular materials. The model is validated on quasi-static experiments on loading–unloading cycles. The experimentally observed hysteresis phenomena were numerically confirmed with a good accuracy by the proposed model.

## 1. Introduction

When a granular bed is subjected to mechanical loading, local deformations occur resulting in rearranged deformed grains that form 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, the literature provides at least three types of models, without or with non-trivial connections:

— The first class is composed of elastic and plastic modelling in solids [1–3].

— The second class considers granular media at the discrete level [4,5].

— The third class is composed of multiphase flow modelling of granular media [6–10].

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.

In this work, the irreversible compaction of powders 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. So, there is no diffusion of the gas in the solid matrix. Hence, 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 the Kapila *et al.* [9] model. In particular, they introduced a ‘configuration energy’ as a function of the volume fraction. This is the energy that is necessary to maintain the contact between grains. However, such an approach is unable to model hysteresis phenomena. Gonthier [11] (see also Jogi [12]) proposed a model capable 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 intergranular stress are expressed. The physical sense of the ‘no-load volume fraction’ is not direct.

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 interaction in Favrie *et al.* [13] and Favrie & Gavrilyuk [14].

The model is able to solve problems involving material interfaces separating a pure fluid and a granular mixture. Numerical methods to solve interface conditions in the context of compressible fluids governed by different equations of state (EOSs) have been paid considerable attention during the last two decades owing to the pioneering works of Karni [15] and Abgrall [16]. Multiphase flow models have been used in this work [17] in order to solve interface conditions with correct thermodynamics at interfaces, but in the context of fluid–fluid interfaces only. Capillary interfaces were addressed later by Perigaud & Saurel [18], evaporating interfaces by Saurel *et al.* [19] and solid–fluid interfaces by Favrie *et al.* [13].

In §2, we start from the Hertz theory to introduce a state variable which is necessary for the description of the granular media deformation. In §3, we present a reversible single-pressure, single-velocity multiphase flow model. In §4, we introduce dissipative effects in the model and formulate an irreversible compaction model. The von Mises-type yield condition is formulated in terms of the ‘configuration energy’. In §5, we present EOSs for modelling granular energetic materials. These EOSs are easy to calibrate using experiments on load–unload cycles. In §6, we present the numerical method to solve the problem. One- and two-dimensional examples are presented in §7. Conclusion is addressed in §8.

## 2. Euristic considerations

A granular media is different from a classical fluid–fluid mixture in several aspects, one of them being the presence of intergranular contacts. Owing to grain contacts, the medium presents resistance to compression linked to efforts that appear in the grains around each contact point. At the macroscopic scale, these efforts can be summarized with an additional EOS that expresses the resistance to the compression. An additional elastic energy is stored in the solid independently of the hydrodynamic energy owing to the pressure in the gas surrounding the particles.

To understand the problem, we consider a simplified case where the compressibility of the particles is neglected. We consider the one-dimensional case shown in figure 1. The spherical particles of the same diameter *L*_{0} are placed in the vacuum. A compression loading is applied. Consider the case where the particle behaves elastically (each particle will recover its initial spherical shape when the loading is suppressed).

Following Hertz, the particle displacement is [20]
2.1Here, *ν* is Poisson's ratio, *E* is Young's modulus, *L* is the actual distance between particles’ centres after the compression loading of amplitude *p* (figure 1). We will give now a simple interpretation of such an expression. Let us denote *X*_{k} the initial coordinate of the *k*-th particles and *x*_{k} its Eulerian coordinate. Then we have *L*=*x*_{k+1}−*x*_{k} and *L*_{0}=*X*_{k+1}−*X*_{k}. We introduce *λ*=*L*_{0}/*L* as a discrete presentation of the gradient of the continuous Lagrangian coordinate *X* with respect to the Eulerian coordinate *x*. As the Lagrangian coordinate *X* is conserved along trajectories
its gradient ∂*X*/∂*x* (that will also be denoted by *λ*) will satisfy the conservation law
Expression (2.1) can now be rewritten in the form:
2.2A simplified dynamic system can easily be formulated to study the wave propagation in such a system
This model is completed with the EOS *p*=*χ*(*λ*−1)^{3/2}, where the constant *χ* is defined by (2.2). With this EOS, the model is hyperbolic. The characteristic speeds *κ* can be found in explicit form:
It is interesting to note that if *λ*=1, the ‘sound velocity’ in such a system vanishes. Such a phenomenon is often called ‘sonic vacuum’ [21]. This sound velocity is an analogue of the velocity of transverse waves in solids.

## 3. Reversible model

The aim of this section is to develop a flow model describing the reversible compaction of the granular materials in contact with a compressible fluid. The model must take into account the compressibility of the solid (the corresponding variables are denoted with index *s*) and the fluid (the corresponding variables are denoted with index *g*), and be able to deal with interfaces between pure fluids and granular materials.

Let us consider the state variable denoted now *λ*_{s} instead of *λ* introduced previously. This variable defines the deformation of the elastic solid particles. This variable verifies the inequality *λ*_{s}≥1. It increases in shocks and decreases in rarefaction waves. *λ*_{s} evolves according to the following conservative equation:
This equation is complemented by the conservation of the mass of the mixture
where *ρ*=*α*_{s}*ρ*_{s}+*α*_{g}*ρ*_{g} is the mixture density, *α*_{k} is the volume fraction of the *k*-th component and *ρ*_{k} are the corresponding densities (*k*=*s*, *g*). Also, the mass of each phase is conserved
where *Y* _{k}=*α*_{k}*ρ*_{k}/*ρ* is the mass fraction. The entropy *s*_{k} of each phase is also conserved (we consider a reversible process)
The total energy per unit mass will be taken in the form:
where *e*_{s}(*ρ*_{s},*s*_{s}), *e*_{g}(*ρ*_{g},*s*_{g}) are the EOSs of the solid and the gas, respectively, satisfying the Gibbs identity
*T*_{k} being the temperatures of each component, *τ*_{k} being the specific phase volumes and *f*(*λ*_{s}) being the specific elastic energy owing to the contacts between grains. The momentum must be conserved in the model
with the total pressure defined as
The system is closed by the pressure equilibrium relation *p*=*p*_{s}=*p*_{g}. The energy equation is then automatically satisfied
To show this, it is sufficient to remark that in the case where d*Y* _{s}=0, we have
Then the energy equation is a consequence of this identity, the mass, momentum equation and the equation for *λ*_{s}.

### (a) Hyperbolicity of the equilibrium model

The equilibrium relation *p*=*p*_{s}(*ρ*_{s},*s*_{s})=*p*_{g}(*ρ*_{g},*s*_{g}), after eliminating the volume fraction of one of the components, can be written in the form: *p*=*p*(*ρ*,*Y* _{s},*s*_{s},*s*_{g}).

Its differential is
where the Wood sound velocity *c*_{W} is defined as
with the squared sound velocities of each phase are given by .

Hence,
Finally, the model can be written in the following non-conservative vector form: ∂**W**/∂*t*+*A*(∂**W**/∂*x*)=0 with **W**^{T}=(*ρ*,*λ*_{s},*u*,*s*_{s},*s*_{g},*Y* _{s}) and the matrix **A** defined as
The waves speeds are
When the contacts between grains are present, the sound speed must be greater than that in the granular material without contacts. So, we must have
We will suppose that
3.1Moreover,
In particular, this implies that the sound speed of granular materials with contacts is greater than the Wood sound speed: for all . The value corresponds to the maximum admissible elastic compression after which a plastic deformation occurs. An example of such a function is
3.2Here, the exponent *n*_{e} (the index ‘e’ means elastic) should be larger than 2. The constant *χ* and the exponent *n*_{e} (which is not necessarily a whole number) should be calibrated on a specific granular material when its behaviour is purely elastic. Other EOSs can also be proposed.

## 4. Irreversible model

The aim of this section is to introduce dissipation in the model explained previously. The approach is similar to the one developed in Favrie & Gavrilyuk [22, 23] for the treatment of elasto-plasticity.

Let us introduce a relaxation term in the equation for *λ*_{s}
4.1and the relaxation equation for *α*_{s}
The last equation will replace the equilibrium condition: *p*_{g}=*p*_{s}.

### (a) Entropy inequality

Using the mass, momentum, energy conservation law and Gibbs identity, one can find the following balance equation:
This equation can be separated into two parts
and
Here, *p*_{I} is the interface pressure (average pressure at the grain surface). That pressure can be chosen, for example, in the form coming from the solution of the Riemann problem between solid and gas [24]:
where *Z*_{k} are acoustic impedances of corresponding components. The entropy inequality
is verified if
where μ is a positive constant.

Consider the part of the stress related to grain interactions
The term *Π*_{s}=*λ*_{s}d*f*(*λ*_{s})/d*λ*_{s} has the dimension of the energy per unit mass. In particular, if , the inequalities (3.1) will imply that this energy is decreasing. The choice of will now be done in the same way as in classical viscoplasticity theory of Maxwell-type materials.

### (b) Equilibrium configuration

When the granular material is undergoing compression, some elastic energy will be stored in the grains. Then, the grains will be deformed permanently owing to plastic deformations and their rearrangement. Finally, in the equilibrium configuration, a part of the energy which is not dissipated will be stored as the energy called ‘configuration energy’ [6]. This is the energy which is necessary to maintain contact between grains. Usually, it is taken as a function of the volume fraction of the solid, so that the total energy of the granular mixture is not additive: *E*=*e*+*Y* _{s}*B*(*α*_{s})+1/2*u*^{2}, *e*=*Y* _{g}*e*_{g}+*Y* _{s}*e*_{s}. The EOS *B*(*α*_{s}) is determined from experiments on quasi-static compression of powders. The mixture volume is measured, and the solid volume fraction is deduced as a function of the applied constraint. This type of experiments is described, for example, in Kuo *et al.* [25], Elban & Chiarito [26] and more recently in Jogi [12], where successive load–unload cycles have been studied.

A general approach for constructing such a function could be described as follows. Let *g*(*α*_{s}) be a smooth function. Consider the function
4.2with some *n*_{p}>1 (index ‘p’ means plastic). The value *α*_{0} corresponds to the solid volume fraction when the granular constraint is zero. This volume fraction depends on the powder material, on its morphology, its granulometry, etc. Parameters *a* and *n*_{p} are characteristics of a given powder, and more precisely to its response during quasi-static loading. Such a ‘configuration energy’ can be seen as a sort of ‘granular material yield strength’.

It was possible to fit corresponding experimental curves by the following function [10]:
4.3This function has the following properties:
and
To avoid (for numerical reasons) the infinite derivative at *α*_{s}=1, we consider a modified function , where *m* is a little bit larger than 1 (for example, *m*=1.01). This function is monotonic with respect to *α*_{s} and is bounded with its first derivative at *α*_{s}=1. The behaviour of this function is close to that given by (4.3). The corresponding dynamical system governing such a ‘plastic transformation’ can easily be written. For example, one can take
Here, *τ* is a relaxation time. In the following, we will suppose that * τ* is very small. Finally, the equilibrium configuration is determined by the following conditions:
If , we have a one-to-one correspondence between variables

*λ*

_{s}and

*α*

_{s}varying in the intervals , (

*α*

_{0},1), respectively.

### (c) Hyperbolicity of the irreversible model

Finally, we proposed the following irreversible model of dynamical compaction of granular materials generalizing that in Saurel *et al.* [10]:
With *E*=*Y* _{g}*e*_{g}+*Y* _{s}*e*_{s}+*Y* _{s}*f*(*λ*_{s})+(1/2)*u*^{2} and *P*=(*αp*)_{g}+(*αp*)_{s}+(*αρΠ*)_{s}, *Π*_{s}=*λ*_{s}d*f*/d*λ*_{s}
Each energy *e*_{k} verifies the Gibbs identity
It can be proved that the system is hyperbolic, with the eigenvalues given by
Here, is the frozen sound speed. Comparing the sound speeds of the non-equilibrium model and the equilibrium elastic model, we have
Finally, we have proposed a hyperbolic model which is compatible with the following natural model requirements:

— the entropy inequality is verified;

— the granular pressure

*β*_{s}=(*αρΠ*)_{s}decreases during the relaxation step (Maxwell-type model).

### (d) Sound velocity at the yield surface

Consider now the equilibrium *p*_{1}=*p*_{2} at the yield surface *B*(*α*_{s})=*f*(*λ*_{s}). For the differential of the total pressure, we have
The differential d*α*_{s} can be expressed now from the equilibrium condition *p*_{s}=*p*_{g}. To find the sound velocity, it is sufficient to consider the case, where the entropies *s*_{k} and the mass fractions *Y* _{k} are constant. In this case, the pressure will be the only function of the density. Indeed, we have
But the equality of pressures *p*_{1}=*p*_{2} implies .

Hence,
where *c*^{2}_{p} is the sound velocity of ‘plastic’ waves
The velocity of ‘plastic waves’ *c*^{2}_{p} is smaller than the equilibrium sound speed if
A sufficient condition for that is *α*_{s}(d*B*/d*α*_{s})<*λ*_{s}d*f*/d*λ*_{s} along the curve *B*(*α*_{s})=*f*(*λ*_{s}) because
To resume, we constructed two mathematical models describing elastic and non-equilibrium elastic–plastic behaviour of materials, respectively. The last model admits, as an attractor, an equilibrium state belonging to the von Mises-type yield surface. Sufficient conditions on the EOSs are proposed to guarantee that the sound velocity at the yield surface was smaller than the sound velocity of the elastic equilibrium model. Finally, sound velocity of the elastic equilibrium model is smaller than the sound velocity of the non-equilibrium elastic–plastic model. This structure assures a thermodynamically correct formulation of the model of granular materials.

## 5. Equation of state

To be able to solve the governing equations, we need to know four EOSs: the bulk energies *e*_{k}(*ρ*_{k},*p*_{k}), the elastic energy equation *f*(*λ*_{s}) related to the contacts between grains (presenting a reversible energy), and finally the ‘compaction energy’ equation *B*(*α*_{s}) as a function of the volume fraction of the solid (presenting an equilibrium energy attained in a non-equilibrium process).

### (a) Hydrodynamic part of the equation of state

For the sake of simplicity, each material is assumed to obey the ‘stiffened gas’ EOS. These EOSs are suitable both for gas and condensed solids. It reads
5.1Parameters *γ*_{k} and are specific to a given material. They are determined from a given experimentally established Hugoniot curve, as detailed for example, in Le Metayer *et al.* [27]). Data for some materials are given in table 1.

### (b) Equation of state related to the grains interaction

We should distinguish two different EOSs corresponding to the reversible and to the irreversible process. The behaviour of the stress corresponding to the elastic interaction energy
as a function of the volume fraction *α*_{s} is schematically shown in figure 2. The thin lines correspond to the reversible case. The bold line corresponds to the irreversible case, where *f*(*λ*_{s})=*B*(*α*_{s}).

The elastic granular EOS is determined by using experimental data about quasi-static compression of powders. Such an equation is taken in the form (3.2). The compaction energy corresponding to the complete equilibrium is taken in the form (4.3). By using Jogi [12] experimental data, based on granular HMX, with particles size about 100 μm, we have obtained the fitting curve (dashed line) shown in figure 3.

In the next section, we perform numerical validation of the model. In particular, we will study the compaction hysteresis phenomenon. Appropriate numerical schemes have been derived in Saurel *et al.* [28] in a simplified situation of fluid–fluid model. The same algorithm can be used for the present application, except some details that have to be precise.

## 6. Numerical method

The non-equilibrium model (4.2) is reminiscent of that studied in Saurel *et al.* [28]. However, the model is free now of the empirical (but reasonable) ‘switch procedure’ used in that paper. The numerical method described in this reference proceeds in three main steps:

— Hyperbolic step: solve the non-equilibrium system without relaxation terms.

— Relaxation step: determine the equilibrium pressure and associated volume fraction. In this step, the granular elastic stress

*β*_{s}will be relaxed too.— Internal energies reset: using the equilibrium volume fraction obtained previously, determine the corrected equilibrium pressure and the internal energy of each phase on the base of the conservative mixture energy equation.

The third step of this numerical method (internal energy reset) is practically unchanged. We describe below only the first and the second step.

### (a) Hyperbolic step

In the present context of granular materials, the formulation solved in the hyperbolic step is the following:
6.1The only difference compared with the formulation used in Saurel *et al.* [28] for fluid mixtures relies in the mixture total energy definition *E*=*e*+(1/2)*u*^{2} that additionally involves granular energy *f*(*λ*_{s}).

### (b) Relaxation step

The second step (pressure relaxation and plastic relaxation) is however quite different because the mechanical equilibrium conditions are different. The relaxation step consists in solving the following system of ordinary differential equations:
where *p*_{I}=*Z*_{2}*p*_{1}+*Z*_{1}*p*_{2}/*Z*_{1}+*Z*_{2}.

We will solve this system in the limit and . In this case, an approximate solution can be given in algebraic form. The internal energy equations can be rewritten as
Integrating these equations along the shock (between 0 and Δ*t*), this system becomes
with
Here, the superscript 0 corresponds to the end of the elastic step. The average interfacial pressures can be estimated as , that is the pressure in the equilibrium state.

Two cases have to be considered here:

The first case: . Then, no plasticity in the grains occurs (∂

*λ*_{s}/∂*t*=0) and we have: .Thus, the equation is reduced to ,

*k*=*g*,*s*. This case is similar to the one presented in Saurel*et al.*[28].The second case is . Then plastic deformations occur and we have an equilibrium state at the yield surface given by

*f*_{s}=*f*(*λ*_{s})=*B*(*α*_{s}).

As , the system of algebraic equations for variables *p*, *v*_{s}, *v*_{g} is
6.2Here, the volume fractions should be replaced by the expression
To illustrate this procedure, let us consider the example of materials governed by the stiffened gas EOS
System (6.1) becomes
Let us recall that the state 0 corresponds to the state at the end of the hyperbolic step where, *a priori*, the pressures *p*^{0}_{g} and *p*^{0}_{s} are different.

The specific volumes of each phase can be found In the last equation, we have to replace . The pressure is found by the Newton method from the equation

## 7. Test problem and validation

The model is capable to describe irreversible wave propagation in granular materials. It can also be applied for the problems involving interfaces between granular materials and pure fluids. We will now demonstrate it on various test problems.

### (a) Compaction of a powder sample

Let us consider an HMX powder sample of 1.6 cm long made of the same powder as that studied by Jogi [12] with an initial solid volume fraction of 0.7. The sample is pressed by a piston moving to the left at the imposed velocity of 1 m s^{−1}. The right boundary of the sample is fixed. The situation is depicted in figure 4.

The computational domain involves 100 cells. The method described in the preceding section is used with an extended version on moving meshes. The corresponding solution is shown at times 26.7, 53.5 and 80.2 μs in figure 5.

The propagation mechanism shown in figure 5 is very fast compared with the piston motion. Multiple reflections occur between the two boundaries during slow piston motion rendering the field variables quasi-uniform: they are functions only of time.

Consider the case where the powder sample is pressed up to the solid volume fraction of 0.80 with a piston moving with the velocity of 1 m s^{−1}. Then the piston is removed (−1 m s^{−1}). After some time interval, the piston stops. The flow variables being uniform in space, their time evolution is shown in figure 6.

The first graph in figure 6 shows the piston velocity evolution. The piston trajectory is shown in the second graph. The mixture pressure increases during compression, then decreases rapidly when the piston has negative velocity. The compaction pressure rapidly increases during 0.2 ms when reversible compaction occurs. Then the compaction pressure increases slowly when the irreversible compaction occurs. During expansion the compaction pressure rapidly decreases. The solid volume fraction increases during compression and decreases during expansion. The total variation of the volume fraction is much more important during the compression than during the expansion. Thus, the granular media stay compact. The most characteristic graphs are the two last ones, where the total pressure *α*_{s}*p*_{s}+*α*_{g}*p*_{g}+*α*_{s}*ρ*_{s}*Π*_{s} and the compaction pressure *α*_{s}*ρ*_{s}*Π*_{s} are shown as a function of the volume fraction. These last graphs clearly show model's ability for the hysteresis phenomenon prediction.

We now address the same experiments as in Jogi [12] with three successive loading–unloading cycles. In these examples, the volume fraction has been recorded at the end of each unloading stage. For the first cycle, the powder is pressed up to the solid volume fraction of 0.75, and then the constraint is relaxed. Then the same sample is compressed again up to the solid volume fraction of 0.809 and the constraint relaxed again. For the third cycle, the mixture is pressed up to the solid volume fraction of 0.938 and the constraint relaxed again. The numerical results are shown in figure 7 and compared to Jogi [12] experimental data.

We now examine the model's ability to deal with large amplitude wave propagation as well as interfaces separating pure fluids and granular mixtures.

### (b) Waves and interfaces dynamics

In all following tests, the granular material corresponds to HMX powder with EOS parameters given in table 1 and figure 3.

We consider a ‘shock tube’ of 1 m length consisting of a high-pressure chamber at the left, filled with a gas at the initial pressure of (0.1 GPa) and a low-pressure chamber at the right, filled with the same granular material as studied in figures 5–7 at atmospheric pressure. The gas is governed by the ideal gas EOS with the air polytropic exponent. An initial volume fraction discontinuity is thus present in the tube at abscissa 0.65 m. In order to prevent the divisions by zero in the flow model, the solid volume fraction in the left chamber is set to 10^{−6}. In the right chamber, it is equal to 0.68. The initial density of solid is constant in the whole domain and is equal to 1903 kg m^{−3}. That of the gas phase in the high-pressure chamber is set to 100 kg m^{−3}, while in the right chamber it is set to 1 kg m^{−3}. Both initial states are initially at rest. The aim of this test problem is to show the model ability to compute right-facing shock and associated bed compaction, the left-facing expansion wave and the interface that separates the nearly pure gas from the mixture. The results are shown in figure 8 at time 0.10 ms. The mesh involves 1000 cells.

On the first graph of figure 8 the left-facing expansion wave and the right-facing shock are clearly visible. At the interface location, visible at the fifth and sixth graphs, the total pressure and velocity are perfectly continuous. The mixture pressure (*p*=*α*_{1}*p*_{1}+*α*_{2}*p*_{2}) has a jump at the interface owing to the granular pressure effect.

The same test is repeated with increased initial gas pressure in the left chamber (10 GPa) in order to simulate conditions close to detonation products contact with a granular medium. The corresponding results are shown in figure 9 at time 47 μs.

This test shows the model's ability to deal with interfaces separating fluids and granular mixtures in severe conditions, where the fluid behaviour of granular materials is recovered as a limit. The compaction pressure in the granular bed becomes negligible.

### (c) Two-dimensional illustrations

The model and the numerical method are able to deal with complex multidimensional phenomena appearing during the impact and penetration of a projectile in a granular bed. Such an example is shown in figure 10, where a copper projectile impacts onto a gas–granular mixture interface. The solid projectile is made of copper and has the initial velocity of 100 m s^{−1}. A 600×100 cells mesh is used. This type of computation involves three materials: the air, the granular powder and the copper projectile. The model and the method extension to an arbitrary number of phases are not straightforward. It will be explained in detail in further publications. In particular, we do not suppose that the copper projectile is a rigid solid: it behaves as a fluid.

In figure 11, at the left, we show the flow behaviour when only the hydrodynamic behaviour occurs, while at the right the interaction between grains is taken into account. At the first instant, the behaviour is almost the same. When the reflected waves interact with the ‘rarefaction region’ (free of particles), the solutions behaviour is different. In particular, the amount of ejected granular material is different. Indeed, more granular material is ejected in the ‘fluid’ case, with a formation of the cumulative jet in the middle.

### (d) Shock-induced powder compaction

We now examined the model's validity against experimental data for shock propagation in powder samples. Corresponding experiments have been done by Sandusky & Liddiard [29]. In the following test case, the granular material corresponds to HMX powder with the following parameters: *α*_{s0}=0.73, *a*/*n*_{p}=2.5×10^{4}, *n*_{p}=2, *n*_{e}=2.8, *χ*=10^{6} for the granular EOS. HMX parameters are given in table 1. We assume that the pores are filled with air with the density *ρ*=1.2 kg m^{−3} and the polytropic exponant *γ*=1.4. At initial time, the HMX–air mixture is at rest, and a piston impacts the granular bed at various velocities. In figure 12, we compare experimental measurements (cross symbols) and numerical results (lines). The shock wave speed, the total pressure and the per cent theoretical maximum density after compaction for various piston velocities are presented. The results are in excellent agreement with experiments.

## 8. Conclusion

A multiphase flow model for the irreversible compaction of powders 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 granular mixtures. Extension of this model is needed to include the friction effect between grains, the effects of grain scales, chemical reactions, high-frequency thermal fluctuations at the intergranular contact surfaces, etc.

## Funding statement

This work was partially supported by A*MIDEX and ANR, the grant ANR-11-LABX-0092 and ANR-11-IDEX-0001-02.

## Acknowledgements

We thank the referee for important remarks.

- Received April 2, 2013.
- Accepted September 5, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.