Quentin Carmouze

Richard Saurel

Emmanuel Lapebie

Riemann solver with internal reconstruction (RSIR) for compressible single-phase and non-equilibrium two-phase flows

Keywords: two-phase, dense-dilute, weakly hyperbolic, Riemann solver

. Second the two-phase flow model is considered. A locally conservative formulation is built and validated removing issues related to nonconservative terms. However, two extra major issues appear from numerical experiments: The solution appears not self-similar and multiple contact waves appear in the dispersed phase. Building HLLC-type or any other solver appears consequently challenging. The modified Linde (2002) method is thus examined for the considered flow model. Some basic properties of the equations are used, such as shock relations of the dispersed phase and jump conditions across the contact wave. Thanks to these ingredients the new Riemann solver with internal reconstruction (RSIR), modification of the Linde method, handles stationary volume fraction discontinuities, presents low dissipation for transport waves and handles shocks and expansion waves accurately. It is validated on various test problems showing method's accuracy and versatility for complex flow models. Its capabilities are illustrated on a difficult two-phase flow instability problem, unresolved before.

I.

Introduction The present contribution addresses building of a robust Riemann solver with limited dissipation for complex flow models. The reconstruction method of [START_REF] Linde | A practical, general-purpose, two-state HLL Riemann solver for hyperbolic conservation laws[END_REF] is revisited and improved in terms of accuracy and robustness. This effort is mainly motivated by the numerical approximation of a twophase non-equilibrium flow model developed by the authors that involves a series of theoretical challenges, presented in the following. In the present introduction the modelling context is recalled first, and the numerical approach is introduced secondly.

Modelling context

It is well accepted that hyperbolic models are mandatory to deal with phenomena involving wave propagation. This is the case for multiphase flows in many situations such as in particular shocks and detonations propagation in granular explosives and in fuel suspensions, as well as liquid-gas mixtures with bubbles, cavitating and flashing flows, as soon as motion is intense and governed by pressure gradients. This is thus the case of most unsteady two-phase flow situations. Wave propagation is important as it carries pressure, density and velocity disturbances. Sound propagation is also very important as it determines critical (choked) flow conditions and associated mass flow rates. It has also fundamental importance on sonic conditions of detonation waves when the two-phase mixture is exothermically reacting [START_REF] Petitpas | Modelling detonation waves in condensed energetic materials: Multiphase CJ conditions and multidimensional computations[END_REF]. Hyperbolicity is also related to the causality principle, meaning that initial and boundary conditions are responsible of time evolution of the solution. When dealing with first-order partial differential equations it means that the Riemann problem must have a solution, and the Riemann problem is correctly posed only if the equations are hyperbolic. However, only a few two-phase flow models are hyperbolic in the whole range of parameters. The [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] model (BN) seemed the only formulation able to deal with such requirement. However, in the dilute limit at least, the acoustic properties of this model seemed inconsistent (Lhuillier et al., 2013). Indeed, with this model, the dispersed phase sound speed corresponds to the one of the pure phase, while this phase is not continuous and unable to propagate sound in reality, at least at a scale larger than particle's one. When the phase is not continuous (dispersed drops in a gas, dispersed bubbles in a liquid), the associated sound speed should vanish, such effect being absent in the formulation. In the low particle's concentration limit, the [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF] model is preferred. This model corresponds to the Euler equations with source terms for the gas phase and pressureless gas dynamic equations for the particle phase (see also [START_REF] Zeldovich | Gravitational instability: An approximate theory for large density perturbations[END_REF]. This model is thermodynamically consistent and hyperbolic as well, except that the particle phase equations are weakly hyperbolic. In this model, contrarily to the BN model, sound doesn't propagate in the particles phase, this behaviour being more physical in this limit. However, the Marble model has a limited range of validity as the volume of the dispersed phase is neglected, this assumption having sense only for low (less than per cent) condensed phase volume fraction. Recently, the gap between these two models has been filled (Saurel et al., 2017a). Modifying the volume fraction equation in the BN model resulted in a flow model where sound propagates only in the carrier phase. The model is hyperbolic with same 4 wave speeds as Marble's one and is thermodynamically consistent in the stiff pressure relaxation limit. Moreover, in the stiff velocity relaxation limit, the [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF] model, important for diffuse interface computations [START_REF] Saurel | Diffuse interface capturing methods for compressible two-phase flow[END_REF], is recovered. Saurel et al. (2017a) model has been solved in the same reference with a Godunov type scheme based on [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves and obstacles[END_REF] solver. However, as well known, this solver is quite diffusive for stationary discontinuities and linearly degenerate fields, such as in the present context, volume fraction discontinuities and contact waves. This is precisely the motivation of the present work, focused on the building of a Riemann solver with enhanced accuracy.

Riemann solver with internal reconstruction (RSIR)

In the quest of Riemann solver with low dissipation for this flow model several issues appear:

-The flow model, as most two-phase flow models, presents non-conservative terms; -Numerical experiments of typical initial value problems (IVPs) achieved with the Godunov-Rusanov method exhibit non self-similar solutions. Such behavior appears as a combination of non-conservative terms, acting as a drag force (in differential form), and stiff pressure relaxation, mandatory for this specific flow model. -Governing equations of the dispersed phase are hyperbolic degenerate, as a single eigenvalue is responsible for characteristic waves propagation. Therefore, it is impossible to determine a basis of eigenvectors and associated Riemann invariants. Moreover, the solution can be multivalued, as for the Marble's model [START_REF] Saurel | Two-phase flows -Second-order schemes and boundary conditions[END_REF]. It means that multiple volume fraction waves may be present in the solution, rendering the analysis and design of any Riemann solver intricate. Several attempts for the building of approximate Riemann solver were done by the authors for this flow model on the basis of, -characteristic relations for the carrier phase and jump conditions for the dispersed one, -HLLC-type approximation based on a local conservative formulation, that will be presented later. None of these attempts yielded efficient solver. The authors consequently move to another type of solver, based on internal reconstruction of intermediate states, computed from a simple and robust intercell state, such as [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves and obstacles[END_REF] or HLL (Harten et al., 1983). This research direction has been investigated by [START_REF] Linde | A practical, general-purpose, two-state HLL Riemann solver for hyperbolic conservation laws[END_REF], [START_REF] Miyoshi | A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics[END_REF] and many others, mainly in the frame of magnetohydrodynamics equations that also involve many waves in the Riemann problem. The aim is to build two intermediate states instead of one. Doing so, the method should maintain stationary discontinuities and reduce numerical diffusion during simple transport. In the present work the Linde method is revisited, and extra physics is embedded to enhance robustness and accuracy. This is done to the price of generality loss, in the sense that the method becomes model dependent, as most Riemann solvers. The paper is organized as follows. The [START_REF] Linde | A practical, general-purpose, two-state HLL Riemann solver for hyperbolic conservation laws[END_REF] method is recalled in Section II in its basic version and computational examples are shown with the Euler equations. It sometimes works perfectly, but oscillations appear depending on the initial conditions. This observation motivates insertion of extra physics in the solver, resulting in significant improvements, yielding robust and accurate solutions. Similar accuracy and robustness as the HLLC solver are observed. In Section III the internal reconstruction method is extended to the two-phase flow model of interest. As before extra physics is inserted in the closure relations. The (trivial) Rankine-Hugoniot relations of the dispersed phase are used as well as jump conditions across contact wave of the dispersed phase. Thanks to these ingredients the flow solver becomes very efficient. One-dimensional computational examples are shown in the same section. The method is then embedded in the DALPHADT unstructured meshes code. It is used to compute fingering instability occurring during explosive dispersion of particle clouds. Such instability seems misunderstood and not reproduced by existing flow models. Intensive experimental and numerical studies were done in this area recently, as for example in [START_REF] Rodriguez | Solid-particle jet formation under shock-wave acceleration[END_REF][START_REF] Mcgrath | Modeling compressible multiphase flows with dispersed particles in both dense and dilute regimes[END_REF], [START_REF] Osnes | Numerical simulation of particle jet formation induced by shock wave acceleration in a Hele-Shaw cell[END_REF] and [START_REF] Xue | Dual hierarchical particle jetting of a particle ring undergoing radial explosion[END_REF]. Thanks to the new model and present numerical method, the formation stage of this instability seems correctly predicted at least qualitatively. Computational examples are shown in Section IV. Conclusions are given in Section V.

II-Riemann solver with internal reconstruction (RSIR) for the Euler equations

The Riemann solver with internal reconstruction is a modification of the [START_REF] Linde | A practical, general-purpose, two-state HLL Riemann solver for hyperbolic conservation laws[END_REF] solver. The original Linde solver is recalled in the frame of the Euler equations and modifications are addressed next. The Euler equations of compressible fluids consist in a system of conservation laws,

U F 0 t x       (II.1) with   T U , u, E     and   T F u, u ² p, ( E p)u      
where  denotes the density, u the velocity and

1 E( e u ²) 2  
the total energy. The pressure p is given by a convex equation of state (EOS), as a function for example of internal energy e and density  : p p( , e)   . The stiffened gas (SG) EOS will be used frequently in the present contribution as, p( , e) ( 1) e p         , where  and p  are typical constants for a given fluid.

This system is strictly hyperbolic with wave's speeds

1 u   , 2 u c    et 3 u c    . The sound speed is defined by s p c      
where s denotes the entropy.

The various reconstruction methods considered in the present paper are based on HLL (Harten et al., 1983) approximate solution, or its simplified version by [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves and obstacles[END_REF]. In the HLL solver, the extreme waves only are used, and are estimated, following Davis (1988)

by,   L L L R R S min u c , u c    and   R L L R R S max u c , u c    , (II.2)
where subscripts L and R denote the left and right states in the initial data of the Riemann problem.

      * * * R L HLL R M R M L L S S U S S U S S U      .
(II.4) Relation (II.4) can be expressed as,

* * * HLL R R L L U U U     , with R M R R L S S S S     and M L L R L S S S S     .
The contact wave speed is given by,

            R L L L R R L R M L L L R R R p p u S u u S u S S u S u              , Relation (

II.1 Linde reconstruction

Linde reconstruction is based on the following relation,

  * * R L R L U U U U     (II.5
) where  represents a viscosity parameter, 0 1    . When  is taken equal to zero, the HLL approximation is recovered. When 1   the reconstruction tends to the HLLC representation but is not equivalent, as interface conditions are ignored in this approach. It is worth to note that when 1   , isolated density discontinuities are preserved, an important property in CFD solvers. Relation (II.5) is then combined to Relation (II.4) resulting in: Accuracy of the conventional Linde method is consequently highly dependent of initial data and viscosity parameter β. It is worth to mention that in Linde ( 2002) a method is given to determine parameter vector for β to improve the solution. But it is important to note that as soon as β is strictly less than 1, the solver loses property of isolated stationary discontinuity preservation. Modification is thus addressed in the next paragraph.

    * * L HLL R R L * * R HLL L R L U U U U U U U U                (II.
    * * R R R R R * * L L L L L F F S U U F F S U U            . (II.7) Solution sampling is achieved through, L L * L L M * * R R M R R F if S 0 F if S 0 and S 0 F F if S 0 and S 0 F if S 0              

II.2 New reconstruction method (RSIR)

The reconstruction derived hereafter is based on two ingredients:

-Quasi-isentropic variations across right-and left-facing waves; -Insertion of interface conditions across the contact wave. Although not strictly correct, in particular across strong shocks, thermodynamic evolutions through right-and left-facing waves are approximated as isentropic. Isentropic evolutions are themselves approximated through sound speed definition as,

* 2 R R R * R R p p c      , * 2 L L L * L L p p c      . (II.8)
Across the contact wave SM, the pressure interface condition reads,

* * * L R p p p   . Thus (II.8) becomes,   * 2 * R R R R p p c      ,   * 2 * L L L L p p c      .
Taking the difference of these two relations the following one is obtained,

2 2 * * L L L R R L R L 2 2 2 R R R c c p p c c c                .
For simplicity reasons another approximation is now done with respect to the square sound speeds, that are assumed closed to a certain average:

2 2 2 L R c c c   . Consequently, * * L R R L R L 2 p p c          with 2 2 2 R L c c c 2   . (II.9)
This relation corresponds to a modification of the first relation of System (II.5).

To maintain flexibility of the reconstruction method, parameter β is reintroduced as,

* * L R R L R L 2 p p c                 . (II.10)
Indeed, parameter β is convenient to control numerical viscosity. From (II.10) similar relation is deduced for momentum and energy jumps across the contact wave.

As * * L R M u u S   , (II.10) implies the following relation, * * L R R L R L M 2 p p ( u) ( u) S c                 (II.11)
Then, for the energy jump, using once more the pressure interface condition,

* * * L R p p p 
 , it appears, provided that the fluid is governed by ideal gas or SG EOS,

    * * L R e e    . Consequently,       * * * * 2 *2 * L R R L R L 2 R R L R p p 1 1 E e u e u 2 2 c                                   . Alternatively, it reads,     2 * * L R M R L 2 R L p p S E E c 2                 . (II.12)
Relations (II.10-II.12) summarize as,

* * R ,k L,k k U U    , k=1,..,3. (II.13) with, L R R L 2 p p c               and T 2 M M S 1,S , 2         .
Combining (II.4) and (II.13) the various intermediate states are rebuilt as, The RSIR method is consequently efficient, at least provided that the fluid is governed by the ideal gas or SG EOS, as considered in the derivation of (II.12). Considering a more sophisticated but convex EOS, such as for example NASG (Le Metayer and Saurel, 2016), Relation (II.12) is thus replaced by,

* * L,k HLL,k R k * * R,k HLL,k L k U U U U              , k=1,.
          * * * * * * * * * * 2 R R R L L L R L M R L 1 E E e p , e p , S 2              (II.15)
The RSIR solver now reads,

* * L,k HLL,k R k * * R ,k HLL,k L k U U U U              , k=1,..,3, (II.16) with,       L R R L 2 L R R L M 2 * * * * * * * * 2 R R R L L L R L M p p c p p S c 1 e p , e p , S 2                                                       
. 2017a) is now examined. As already mentioned, this system is weakly hyperbolic. Moreover, it will be shown that solutions are not selfsimilar. For such a complex flow model, present authors tried but failed in the derivation of HLLCtype Riemann solver.

This

III -Extension to dense-dilute two-phase flow model

The two-phase flow model of Saurel et al. (2017a), considering a dispersed phase 1 in a carrier fluid 2 is recalled hereafter. Pressure and velocity relaxation terms only are considered as interaction effects:

    1 1 1 2 u p p t x          ,        1 1 u 0 t x         ,       2 1 1 1 I 2 1 u p u p u u t x x                ,             1 1 1 I I 2 1 I 1 2 E p u E u p u u u μp p p t x x                   , (III.1)     2 2 u 0 t x         ,       2 2 2 2 I 2 1 u p u p u u t x x                ,             2 2 1 I I 2 1 I 1 2 E p u E u p u u u μp p p t x x                    .
Same notations as for the Euler equations are used. Additionally k  represent the volume fraction of phase k, index I is related to interfacial variables and relaxation parameters are denoted by  and  , with respect to velocity and pressure relaxation. Appropriate relations are given for example in Saurel et al. (2017a). It is important to mention that this flow model is valid only in the stiff pressure relaxation limit (    ). Appropriate pressure relaxation solvers are given for example in [START_REF] Lallemand | On the quest for a hyperbolic effective-field model of disperse flows[END_REF]. This is not equivalent to strict pressure equilibrium models that are non-hyperbolic, or conditionally hyperbolic. Also, this is not a restrictive assumption for most two-phase flow applications, except possibly extreme situations, such as hot spots ignition in condensed energetic materials (Saurel et al., 2017b), where pressure relaxation is the driving effect for hot spot appearance. The present model is hyperbolic with wave speeds,

1 4 1 u    , 5 2 u   , 6 2 2 u c    and 7 2 2
u c    . These waves speeds are the same as [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF] model but are significantly different to those of the [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] model. Combination of the equations of System (III.1) result in the following mixture entropy equation, that guarantees non-negative evolutions,

1 1 1 2 2 2 1 1 1 1 2 2 2 2 I 1 I 2 2 1 1 2 s s s u s u (u u ) (u u ) (u u ) t x T T                          . (III.2)
Indeed, admissible estimates for the interfacial velocity are,  the following Rankine-Hugoniot system is obtained (Saurel et al., 2017a):

0 1 1    , 0 1 1    , 0 1 1 e e  , 0 1 u u  , (III.3) 0 0 2 2 2 (u ) (u )        , 0 0 0 0 2 2 2 2 2 u (u ) p u (u ) p          , 0 0 0 0 0 2 2 2 2 2 2 2 E (u ) u p E (u ) u p         
. These relations will be used in the RSIR derivation. Numerical resolution of System (III.1) has been done in Saurel et al. (2017a) with the help of [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves and obstacles[END_REF] solver. This solver being quite diffusive, the aim is now to build an improved solver. System (III.1) involves however three main difficulties:

-It is non-conservative; -The eigenvalue 1 u is multiple. Therefore System (III.1) admits multivalued solutions [START_REF] Forestier | Multivalued solutions to some non-linear and non-strictly hyperbolic systems[END_REF][START_REF] Saurel | Two-phase flows -Second-order schemes and boundary conditions[END_REF][START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF].

-Solutions are not self-similar, as will be shown later with the help of numerical experiments. These issues are addressed gradually in the following.

III.1 Local conservative formulation and Rusanov-type solvers

In Saurel et al. (2017a) a Rusanov-type method was derived to determine qualitative solutions of the new flow model and validations against both exact solutions and experimental data. This method is recalled hereafter, and an improved version based on a local conservative formulation is built. The aim is to show that non-conservative terms are treated correctly through the local conservative formulation.

a) Basic Rusanov version System (III.1) is considered in non-conservative form and in the absence of relaxation terms as,

  F U U U H U, 0 t x x                , (III.4) with,                                       T 1 1 1 1 2 2 2 T 2 2 1 1 2 1 2 1 2 T 1 1 1 1 I I I I U , , u , E , , u , E F u , u , u p E p u , u , u p , E p u u u H 0, 0, p , p , 0, p , p x x x x                                            Let us denote by,   k k k L R S Max ,    ,
the maximum wave speed separating two states L and R. The Rusanov flux reads,

  * R L R L 1 F F F S U U 2         ,
and the Godunov scheme associated to system (III.4) reads,

  n 1 n * * i i 1 1 i i i 2 2 t U U F F tH x           , where i
H is an approximation of non-conservative terms. Following Saurel et al. (2017a), based on [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] method for a slightly different flow model, approximation of these terms read;

-For the momemtum equation,

* * 1 1 1,i 1,i n 2 2 i,u i H p x        with 1,i 1,i 1 * 1 1,i 2 2        .
-For the energy equations,

    * * 1 1 1,i 1,i n 2 2 i,E i u u H p x        with       * 1,i 1 1,i 1 1,i 2 u u u 2        .
These formulas are built to respect mechanical equilibrium condition. Another version is examined hereafter.

b) Local conservative formulation

The interfacial pressure I p appears in the presence of non-conservative terms, such as

  k I p x
. As I p has been assumed equal to the dispersed phase pressure, the following assumption is made:

1,L 1,L 1,R I 1,R 1,L 1,R p if p p if             . ( III
  1 1 u 0 t x        , 1 1 ( ) ( u) 0 t x         ,       2 I 1 1 u p p u 0 t x            ,       I 1 1 E p p u E 0 t x            , (III.6)   1 1 2 u 0 t x        , 2 2 ( ) ( u) 0 t x         ,       2 I 2 2 u p p u 0 t x            ,         I 1 2 2 E p u u p E 0 t x             . In compact form it reads,   U U 0 t x       , with obvious definition for   U  .
The associated Rusanov flux is immediate,

  * R L R L 1 S U U 2            .
(III.7)

From *  , the * F flux of formulation (III.4) is deduced as,

* * k k I * k * k 0 0 F p ( )                   , k=1,2. (III.8)  * k are determined form Rusanov state as,   k,R k,L * k k,R k,L ( ) ( ) 1 2 S                    . (III.9)
The fluxes are inserted in the same Godunov scheme as before,

  n 1 n * * i i 1 1 i i i 2 2 t U U F F tH x           , except that i
H are now given by,

* * 1 1 1,i 1,i n 2 2 i,u i H p x        * * 1 1 1 1 i i n 2 2 i,E i ( ) ( ) H p x          .
In these expressions Length of the domain has been increased to 5 m to avoid interaction with the boundaries, and the initial discontinuity is placed at 2.5 m. Stiff pressure relaxation is used at any time. 1 000 computational cells are used with CFL=0.5. Computations are done with the MUSCL method and Minmod limiter. Solutions for the carrier phase (2) appear self-similar but solutions of the dilute phase 1 are not, regarding volume fraction and velocity. Volume fraction of the dilute phase keeps increasing without converging towards a constant state. The same tendency is observed in the velocity of the dilute phase. These observations are mesh and solver independent.

Origin of this interesting behavior is examined hereafter with the help of the [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] model. The same shock tube test problem is rerun with this model with the help of the HLLC-type solver of [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF]. In the computations of Figure III.4 pressure relaxation is absent. These difficulties are omitted when dealing with single state Riemann solver, such as HLL or Rusanov.

From this robust solution basis, we now address internal reconstruction to reduce numerical smearing.

III.2 Riemann solver with internal reconstruction (RSIR) for the two-phase model

The present flow solver is not based on variations across the various waves but only on rebuilding two intermediate states to preserve isolated volume fraction discontinuities and reduce artificial smearing during transport. In this direction, the intermediate wave speed is based on the phase 1 contact wave:

        * HLL 1 M1 * HLL 1 U u S U    , (III.10)
with * HLL U given by (III.11) in the context of System (III.6) and same wave speed estimates as before. As detailed with the Euler equations the method to solve System (III.1) proceeds in two steps:

-Determine average state with HLL based on (III.6); -Rebuild the solution. Thanks to the local conservative formulation the first step is immediate. System (III.6) is expressed as,

  U U 0 t x       ,
with I p given by (III.5). The average state is obtained as,

* R L L L R R HLL L R S U S U U S S        . (III.11)
The next step consists in the approximation of states L* and R*. Knowledge of the intermediate wave speed M1 S from (III.10) enables following decomposition of the average state,

* * * HLL R R L L U U U     (III.12) with, R M1 R R L S S S S     and M1 L L R L S S S S     .
As System (III.12) involves two unknown variable vectors, another set of relations is needed.

Reconstruction

For the dispersed phase, the Rankine-Hugoniot relations (III.3) imply, * 1R 1R

   and * 1L 1L

   . The difference of these two relations reads,

* * 1R 1L 1R 1L        . Parameter  is introduced to control numerical diffusion,   * * 1R 1L 1R 1L         .
As phase 1 density has no jump across left-and right-facing waves, similar relation is obtained:

          * * 1R 1L 1R 1L
        .

Regarding momentum jump across the intermediate wave, the same relation is used as in the context of the Euler equations,

          * * M1 1R 1L 1R 1L u u S         .
The energy jump relation is based on the interface condition related to the momentum equation of phase 1:

            * * * * * * * * * * 1R 1R M1 1R 1R I 1L 1L M1 1L 1L I 1R 1L u u S p p u u S p p            .
Velocity jump conditions across left-and right-facing waves of the dispersed phase are introduced as,

* 1R 1R u u  and * 1L 1L u u  . Momentum jump across the contact wave consequently reads,             * * * * * * 1R 1R M1 1R 1R I 1L 1L M1 1L 1L I 1R 1L u u S p p u u S p p            .
Assuming fluid 1 governed by the stiffened gas EOS, it becomes,

                * * * * * * 1R 1 1 1, 1L 1 1 1, 1L 1L M1 1R 1R M1 1L 1R 1R 1R * * 1R 1L I 1 e p 1 e p u u S u u S p                                 
The internal energy jump thus reads,

                    * * * * 1L M1 1R M1 I 1 1, * * 1L 1R 1R 1L 1R 1L 1L 1R 1 1 1 u S u S p p e e u u 1 1 1                       .
The total energy jump for phase 1 follows,

                        2 * * * * I 1 ,1 * * M1 1R 1L 1R 1L 1R 1L 1 * * 1L M1 1R M1 1L 1R 1L 1R 1 1 p p S E E 1 2 u S u S u u 1 1                               .
For the first phase, these relations summarize as,

* * R ,k L,k k U U    , k 1,.., 4  , (III.13) with,                               1R 1L 1R 1L M1 1R 1L 2 * * 1L M1 1R M1 I 1 ,1 M1 1R 1L 1L 1R 1R 1L 1L 1R 1 1 1 S u S u S p p S u u 1 2 1 1                                                                     
 For the dispersed phase, combining (III.12) and (III.13) the intermediate states are computed as,

* * R ,k HLL,k L k * * L,k HLL,k R k U U U U              , k 1,.., 4  . (III.14)
Reconstruction is now addressed for the second phase. For the sake of simplicity, let us consider that the volume fraction of the second phase is considered additionally. To preserve the saturation constraint (

1 2 1     ) the jump relation across the intermediate wave reads,   * * 2R 2L 2R 2L         .
Another assumption is now introduced. The density of the carrier phase is assumed uniform in the two intermediate states. In other words,

  * 2,HLL 2 * 2,HLL     .
This assumption is needed as the flow solver behaves incorrectly when the product of two discontinuous functions are present (volume fraction and density).

The assumption made at this level is similar to prolongated formulations used in Ghost Fluid Methods [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], immersed boundary methods and diffuse interface methods [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF][START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Massoni | Proposition de méthodes et modèles eulériens pour les problèmes à interfaces entre fluides compressibles en présence de transfert de chaleur: Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF][START_REF] Saurel | Diffuse interface capturing methods for compressible two-phase flow[END_REF]. In diffuse interface methods, numerical diffusion of volume fractions k  and apparent densities k k   results in automatic prolongment of phase's density k  . In this frame, apparent densities are rebuilt as,

  * * 2 2R 2R     ,   * * 2 2L 2L     . Consequently,       * * 2R 2L 2 2R 2L          .
Momentum jump across the intermediate wave is rebuilt as,

          * * * * 2 M 2R 2L M 2 2R 2L 2R 2L u u S S                
 where the contact wave speed of the second phase is computed as,

        * HLL 2 M 2 * HLL 2 U u S U    .
It remains to determine total energy jump of the second phase. At this level no distinction is made between the left and right velocities of that phase:

                                The internal energy jump thus reads,                  * * * * M 2 M 2 M1 2R 2L I 2 ,2 2L 2R * * 2R 2L 2 S S S p p e e 1                  .
The total energy jump for phase 2 follows,

                  2 * * * * M2 M1 I 2 ,2 * * M2 M2 2R 2L 2R 2L 2R 2L 2 2 S S p p S E S 2 1 1                                 .
The carrier phase is thus rebuilt as,

* * R ,k L,k k U U    , k 5,..,8  , (III.15) with,           2 2R 2L 2 M 2 2 M 2 M1 I 2 ,2 M 2 2 M 2 2 2 1 S S S p p S S 2 1 1                                                Last,
    * * R R R R R * * L L L L L S U U S U U                . (III.16)
Fluxes of System (III.1) are computed as,

* * k k k F k L, R      , (III.17) with         T * * * * k I 1k 1k 1k 1k p 0, 0, , u , 0, 1 , u         .
 * k are determined form HLL state as,

k,R k,L L k,L R k,R * k L R ( ) ( ) S S S S             . (III.18)
The associated Godunov-type scheme including non-conservative terms reads,

  n 1 n * * i i 1 1 i i i 2 2 t U U F F tH x           , with i
H given by,

* * 1 1 1,i 1,i n 2 2 i,u i H p x            * * 1 1 1 1 i i n 2 2 i,E i H p x          .
In these expressions

* 1 k,i 2   and   * 1 1 i 2  
 are given by (III.18) and (III.16) respectively. Gas expansion during explosion fragments the liquid layer to a cloud of droplets that form highly dynamical particle jets. The same observation is reported when a granular bed is used instead of a liquid layer. Dispersion is consequently not spherical or cylindrical in the sense that one-dimensional computations result in significant errors in the predictions of materials presence. 2018). Impulsive motion of a particle ring by a gas flow induces well defined particles fingers flowing to the center direction, oppositely to the gas flow. At later times, here at 8 ms, short wavelength instabilities also appear at the external surface. As time evolves, external surface instabilities grow and become dominant, as shown in the third picture at time 57 ms.

III.3 Examples and validations

Internal jets are thus observed at early times, followed by external ones at late times.

As reported by [START_REF] Rodriguez | Solid-particle jet formation under shock-wave acceleration[END_REF] and [START_REF] Xue | Dual hierarchical particle jetting of a particle ring undergoing radial explosion[END_REF], instabilities appear first at the inner interface and second at the outer one. Shape of these fingers is singular, in the sense that they do not qualitatively compare to the Richtmyer-Meshov instabilities or Rayleigh-Taylor ones, nor any other known instability. Indeed, mushroom type shape is observed with these two instabilities, while fingers are observed in the present context With the help of the new model and present RIRS solver attempt is done to reproduce at least qualitatively these instabilities. Former attempts by the authors were based on [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF] The mesh is made of triangles and contains about 15 cells in the radial direction of the particles ring. Although not precisely identified from the present numerical experiments, the formation mechanism of this fingering instability appears closely related to non-conservative terms. They play the role of a differential drag force, acting intensively at cloud boundaries and vanishing in the wake, when volume fraction gradients disappear. At the modeling level, non-conservative terms present similar form as capillary ones (Brackbill et al., 1992, Perigaud and[START_REF] Perigaud | A compressible flow model with capillary effects[END_REF] except that curvature effects are absent in the present two-phase formulation. Another major difference is that cloud boundaries are obviously highly permeable in the present context, while interfaces are not permeable in conventional hydrodynamic instabilities, except those considering flames and phase transition, where low permeability is present.

V. Conclusion

A Riemann solver with internal reconstruction (RSIR) has been built as an extension of [START_REF] Linde | A practical, general-purpose, two-state HLL Riemann solver for hyperbolic conservation laws[END_REF] solver. It has been first developed for the Euler equations and shown to provide similar accuracy and robustness as the HLLC solver, while being not as systematic as this one. The method has been secondly extended to a two-phase non-equilibrium model developed by the authors. This model presented serious difficulties as it is weakly hyperbolic and valid only in the presence of stiff pressure relaxation, rendering solutions not self-similar. Thanks to the RSIR approach, a low dissipation solver has been developed. It has been validated against solution obtained with more conventional, but dissipative solvers. The new method has been applied in the last section to a difficult problem of fingering instability in granular media and has shown possible explanation of the formation mechanism. Extra work is still needed in this special two-phase flow topic to achieve understanding of this instability.

Regarding the Riemann solver, it seems flexible for many applications where most of the physics is governed by the two extreme waves and an intermediate one. Moreover, a parameter is present to control dissipation when flow conditions are particularly severe. In all present computations the solver appeared robust with parameter (β=1), corresponding to the minimum dissipation.

  , an extra relation is needed.[START_REF] Linde | A practical, general-purpose, two-state HLL Riemann solver for hyperbolic conservation laws[END_REF] postulate is considered first.

Figure

  Figure II.1 -Schematic representation of the two intermediate states * L U and * R U rebuilt from * HLL U .

U

  are computed the various fluxes are computed through the Rankine-Hugoniot relations,

.

  Typical solutions obtained by this solver embedded in a Godunov-type code are shown hereafter. Let us first consider transport of a density discontinuity in a uniform velocity and pressure flow with a gas governed by the ideal gas equation of state (EOS) with ꙋ=1.4. Corresponding results are shown in Figure II.2.

Figure

  Figure II.1 -Computed results with the original Linde solver (β=1) for the transport of a density discontinuity in a uniform pressure and velocity flow. MUSCL extension of the Godunov scheme is used with Minmod limiter. 100 computational cells are used with CFL=0.5. Results are shown at time t=6ms. The density discontinuity is correctly transported, and mechanical equilibrium is maintained.

Figure II. 3 -

 3 Figure II.3 -Computed results with the original Linde solver (β=1) for a shock tube test case. MUSCL extension of the Godunov scheme is used with Minmod limiter. 100 computational cells are used with CFL=0.5. Results are shown at time t=300µs. The various waves and states are computed correctly.

Figure

  Figure II.4 -Computed results with the original Linde solver (β=1) for a double expansion test case. MUSCL extension of the Godunov scheme is used with Minmod limiter. 100 computational cells are used with CFL=0.5. Results are shown at time t=850µs. Oscillations appear. With β=0.9 same behavior of the solution is observed.

Figure II. 5 -

 5 Figure II.5 -Computed results with the original Linde solver (β=0.5) for the double expansion test case of Figure II.4. MUSCL extension of the Godunov scheme is used with Minmod limiter. 100 computational cells are used with CFL=0.5. Results are shown at time t=850µs. Oscillations decrease but are still present.

Figure

  Figure II.6 -Computed results with the original Linde solver (β=1) for a double shock test case. MUSCL extension of the Godunov scheme is used with Minmod limiter. 100 computational cells are used with CFL=0.5. Results are shown at time t=850µs. Oscillations appear as before. They decrease when β tends to zero.

  these intermediate states fluxes are computed with the help of (II.7) as before. The Riemann solver thus consists in (II.3) to compute the HLL state and (II.14) to rebuild the two intermediate states. It is thus particularly simple and contains parameter β to control numerical diffusion. Comparison with the HLLC solver is now addressed on various test problems. In all tests that follow, β is taken equal to 1 to minimize numerical diffusion. The first test consists in the computation of a density discontinuity at rest, as shown in Figure II.7.

Figure II. 7 -

 7 Figure II.7 -Comparison of HLLC and new solver (β=1) embedded in the Godunov-MUSCL method with Minmod limiter for the computation of a contact discontinuity at rest. 100 computational cells are used with CFL=0.5. Results are shown at time t=6ms. For the sake of clarity only 50 symbols out of 100 are plotted for both computations. Both computed results are merged. The double-expansion and double-shock tests of Figures II.4 and II.6 are rerun with the new method and HLLC. Results are shown in Figures II.8 and II.9 showing oscillation free solutions.

Figure

  Figure II.8 -Comparison of HLLC and new solver (β=1) embedded in the Godunov-MUSCL method with Minmod limiter for the computation of a double-expansion test. 100 computational cells are used with CFL=0.5. Results are shown at time t=850µs. For the sake of clarity only 50 symbols out of 100 are plotted for both computations. Both computed results are merged.

Figure

  Figure II.9 -Comparison between HLLC and the new solver (β=1) embedded in the Godunov-MUSCL method with Minmod limiter for the computation of a double-shock test. 100 computational cells are used with CFL=0.5. Results are shown at time t=850µs. For the sake of clarity only 50 symbols out of 100 are plotted for both computations. Both computed results are merged.

  the covolume, the present method requires modification. Combining (II.4) and (II.10) star densities * R  and * L  are determined. With the help of (II.8) the star pressures * R p and * L p are estimated. The associated internal energies are computed with the help of the EOS,  , k=L, R.

  RSIR method is compared to computations based on HLLC solver for the transport of a density discontinuity in Figure II.10 and on shock tube computations in Figure II.11. NASG EOS parameters used in these computations are those of liquid water: ꙋ=4.4, p∞=6.10 8 Pa and b=5.10 -5 m 3 /kg.

Figure

  Figure II.10 -Comparison between HLLC and the new solver (β=1) embedded in the Godunov-MUSCL method with Minmod limiter for the transport of a density discontinuity in a uniform pressure and velocity flow using the NASG EOS. The initial discontinuity is located at x=0.2m. 100 computational cells are used with CFL=0.5. Results are shown at time t=6ms. For the sake of clarity only 50 symbols out of 100 are plotted for both computations. Both computed results are merged.

Figure

  Figure II.11 -Comparison between HLLC and the new solver (β=1) embedded in the Godunov-MUSCL method with Minmod limiter for the computation of a shock tube test with liquid water and NASG EOS. 100 computational cells are used with CFL=0.5. Results are shown at time t=75µs. For the sake of clarity only 50 symbols out of 100 are plotted for both computations. Both computed results are merged.

  in the righthand side of the momentum and energy equations. However, assuming I 1 p p

  equal to the pressure of phase 1 when this phase is present. Possible situations are schematized in Figure III.1.

Figure

  Figure III.1 -Schematic representation of the estimate for I p . As a consequence of jump conditions (III.3) both volume fraction and pressure of phase 1 are invariant across the extreme waves 2

  given by (III.9) and (III.7) respectively. The overall scheme is consequently quite different of the basic version presented before. The two Rusanov solvers are now compared on an arbitrary shock tube test problem, in the absence of velocity relaxation, but with stiff pressure relaxation. Corresponding results are shown in Figure III.2.

Figure

  Figure III.2 -Comparison of the two Rusanov solvers on a two-phase shock tube test. Phase 1 corresponds to the dispersed fluid, considered here as liquid water, with SG EOS parameters ( 1 4.4   and 8 1 p 610 Pa  ). Phase 2 represents the carrier phase, here air considered as ideal gas ( 2 1.4   ). Stiff pressure relaxation is used at any time. 100 computational cells are used with CFL=0.5. Computations are done with the MUSCL reconstruction method and Minmod limiter. Results are shown at t=300µs. For the sake of clarity only 50 symbols out of 100 are plotted for both computations. Both results are perfectly merged validating the approach based on local conservative formulation.

Figure III. 3 -

 3 Figure III.3 -Shock tube test problem of Figure III.2 considered at various times. Length of the domain has been increased to 5 m to avoid interaction with the boundaries, and the initial discontinuity is placed at 2.5 m. Stiff pressurerelaxation is used at any time. 1 000 computational cells are used with CFL=0.5. Computations are done with the MUSCL method and Minmod limiter. Solutions for the carrier phase (2) appear self-similar but solutions of the dilute phase 1 are not, regarding volume fraction and velocity. Volume fraction of the dilute phase keeps increasing without converging towards a constant state. The same tendency is observed in the velocity of the dilute phase. These observations are mesh and solver independent.

Figure

  Figure III.4 -Two-phase shock tube problem of Figure III.2 computed with the[START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] model in the absence of pressure relaxation. Length of the domain has been increased to 5m to avoid waves interaction with the boundaries, and the initial discontinuity is placed at 2.5m. 1 000 computational cells are used with CFL=0.5. Computations are done with the MUSCL method and Minmod limiter. The solution is now self-similar and consists in two decoupled shock tube solutions, as well known.

Figure

  Figure III.5 -Two-phase shock tube problem of Figure III.2 computed with the[START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] model in the presence of stiff pressure relaxation. Length of the domain has been increased to 5m to avoid waves interaction with the boundaries, and the initial discontinuity is placed at 2.5m. 1 000 computational cells are used with CFL=0.5. Computations are done with the MUSCL method and Minmod limiter. The solution is not self-similar. This is not surprising as source terms related to pressure relaxation are present.

Figure III. 6 -

 6 Figure III.6 -Schematic representation of the multivalued phase 1 velocity. Differences in the right and left star velocities of the dispersed phase may result in particle agglomeration as time evolve. This effect is present in the computations of Fig. III.3 and III.5.

  relations (III.14) and (III.15) are used to compute * . Local conservative fluxes are computed as,

  Validations of the flow solver and comparisons with the former Rusanov method are addressed first. Second computational examples are shown showing method's capabilities. a) Validation of the RSIR solver A volume fraction discontinuity at rest is considered to check method capability to maintain such stationary wave. Corresponding results are shown in Figure III.7.

Figure

  Figure III.7 -Results obtained by the new solver (β=1) embedded in the Godunov-MUSCL method with Minmod limiter for the computation of a contact discontinuity at rest. 100 computational cells are used with CFL=0.5. Results are shown at time t=6ms. The discontinuity is well preserved and spurious pressure and velocity oscillations are absent.

Figure III. 8 -

 8 Figure III.8 -Results obtained by the new solver (β=1) embedded in the Godunov-MUSCL method with Minmod limiter for the computation of volume fraction discontinuity transport in a uniform pressure and velocity flow. 100 computational cells are used with CFL=0.5. Results are shown at time t=6ms.

Figure

  Figure III.9 -Comparison of the results obtained by the new solver (β=1) and the local conservative solver of Rusanov, both embedded in the Godunov-MUSCL method with Minmod limiter for the computation of the shock tube presented in Figure III.2. 100 computational cells are used with CFL=0.5. Results are shown at time t=300µs. Stiff pressure relaxation is used.With the RSIR method, main improvements appear in the volume fraction and phase 1 density computations. Indeed, the RSIR solution lies between computed results with the Rusanov method with 100 and 500 cells.

Figure

  Figure IV.1 -A cylindrical explosive charge is initially surrounded by a liquid layer. When the charge explodes the liquid layer transforms to a cloud of droplets forming highly dynamical particle jets.

  Experimental and numerical studies of this phenomenon have been achieved by[START_REF] Zhang | Explosive dispersal of solid particles[END_REF][START_REF] Milne | Dynamic fragmentation of blast mitigants[END_REF],[START_REF] Frost | Jet formation during explosive particle dispersal[END_REF],[START_REF] Parrish | Dynamic jet formation from mitigation materials[END_REF] to cite a few. Simplified situations have been considered in[START_REF] Rodriguez | Solid-particle jet formation under shock-wave acceleration[END_REF] and[START_REF] Xue | Dual hierarchical particle jetting of a particle ring undergoing radial explosion[END_REF]. The explosive is replaced by a shock tube and the matter to disperse is placed between two plates, in a Hele-Shaw cell. Other simplified situations have been considered for example in McGrath et al. (2018), Osnes et al. (2018), Carmouze et al. (2018) to study possible clustering effects due to aerodynamic forces. It seems that formation and size selection mechanism is still unidentified. In the following a configuration like the one considered in Rodriguez et al. (2013) with a Hele-Shaw cell and a particle ring is studied. Typical results reported in Rodriguez thesis are shown in Figure IV.2 at times 5 ms, 8 ms and 57 ms after rupture of the shock tube diaphragm, inducing shock wave and gas flow through a ring of flour particles.

Figure

  Figure IV.2 -Typical interfacial instabilities reported in Rodriguez thesis and papers, as well as in[START_REF] Xue | Dual hierarchical particle jetting of a particle ring undergoing radial explosion[END_REF]. Impulsive motion of a particle ring by a gas flow induces well defined particles fingers flowing to the center direction, oppositely to the gas flow. At later times, here at 8 ms, short wavelength instabilities also appear at the external surface. As time evolves, external surface instabilities grow and become dominant, as shown in the third picture at time 57 ms.Internal jets are thus observed at early times, followed by external ones at late times.

  model and fail to reproduce any, inner or outer, instability. In the present attempt computations are based on flow model III.1 extended to 2D and resolved numerically in the DALPHADT code on triangular cells. Compared to the Marble model, the present one has a fundamental difference. Non-conservative terms are present in the momentum and energy equations. These non-conservative terms are often called 'nozzeling terms' in reference to the Euler equations with variable cross section. We prefer howether to interpret them as 'differential drag force'. As shown for example in Figure III.9, phase 1 is set to motion as a result of these terms, in the absence of conventional drag force. Indeed, drag parameter  is set to zero in these computations. As shown later, it seems that non-conservative terms are responsible for appearance of these instabilities and size selection. Considered computational domain and initial data are reported in Figure IV.3.

Figure IV. 3 -

 3 Figure IV.3 -Computational domain and initial data. The portion represents 1/16 of the complete disc with angle θ = π/8. The domain denoted A corresponds to the high-pressure chamber, filled with gas only. The domain denoted B represents the initial particles ring of 4cm width. The domain denoted C corresponds to the low-pressure chamber at atmospheric conditions.

Figure IV. 4 -

 4 Figure IV.4 -Representation of the mesh size used in the particle jetting computations. In the entire domain 361 222 computational cells are used.

Acknowledgements. of this work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir French Government program managed by the French National Research Agency (ANR). We also acknowledge funding from ANR through project SUBSUPERJET ANR-14-CE22-0014.

qualitatively both fingering instabilities issued of the inner and outer surfaces. In these computations drag force is considered with constant particles diameter (1 mm) and air kinematic viscosity (18.10 -6 Pa.s). Drag force is modelled through [START_REF] Clift | Motion of entrained particles in gas streams[END_REF]

Drag force D F is inserted in the right-hand side of the momentum equation of phase 1 and its opposite in phase 2. The power of this force D 1 F .u is inserted similarly in associated energy equations. In the absence of drag effects internal particles jets selection is observed as well but quantitative differences appear at later times. It is worth to mention that particle volume fraction is initially uniform in the cloud. No arbitrary wavelength is introduced in the initial data.