Philippe Helluy
email: helluy@math.unistra.fr

Jonathan Jung
email: jonathan.jung@unistra.fr

Interpolated pressure laws in two-fluid simulations and hyperbolicity

We consider a two-fluid compressible flow. Each fluid obeys a stiffened gas pressure law. The continuous model is well defined without considering mixture regions. However, for numerical applications it is often necessary to consider artificial mixtures, because the two-fluid interface is diffused by the numerical scheme. We show that classic pressure law interpolations lead to a non-convex hyperbolicity domain and failure of well-known numerical schemes. We propose a physically relevant pressure law interpolation construction and show that it leads to a necessary modification of the pure phase pressure laws. We also propose a numerical scheme that permits to approximate the stiffened gas model without artificial mixture.

Introduction

The numerical simulation of compressible two-fluid flows with Eulerian finite volume approximation has been widely studied. We refer for instance to [START_REF] Mulder | Computing interface motion in compressible gas dynamics[END_REF][START_REF] Karni | Multicomponent flow calculations by a consistent primitive algorithm[END_REF][START_REF] Saurel | A simple method for compressible multifluid flows[END_REF] and included references. The Eulerian approach is very appealing compared to Lagrangian front tracking methods because it generally leads to much simpler algorithms. However, one has to circumvent the pressure oscillations that appear at the two-fluid interface with standard conservative schemes. The lack of accuracy of the Godunov scheme at contact waves is a well-known issue. It is not only observed in two-fluid simulations but also in one-fluid simulations when the single fluid satisfies a complex pressure law [START_REF] Gallouët | A hybrid scheme to compute contact discontinuities in one-dimensional Euler systems[END_REF].

Another less known aspect of classic two-fluid models is that their hyperbolicity domain is generally not convex. Thus in some cases, the Godunov scheme is unstable and fails after only one time iteration.

In this paper we consider the flow of a gas and a liquid modeled by two stiffened gas equations of state. The pressure law is initially only defined in the pure phases: the mass fraction of gas ϕ can take only two values ϕ = 0 or ϕ = 1.

For numerical reasons, the pressure law is often interpolated in an artificial mixture region 0 < ϕ < 1. We show that a naive numerical interpolation always leads to the non-convexity of the hyperbolicity domain and thus to the instability of the Godunov scheme.

We then propose two alternative cures to this issue:

1. The construction of a mixture pressure law based on physical and thermodynamical arguments. We are then able to recover a convex hyperbolicity domain. But we also prove that it is necessary to modify the pressure law of the liquid phase, which cannot be a simple stiffened gas anymore. 2. For keeping the simplicity of the stiffened gas model, we propose another scheme, the Random Interface Solver (RIS) [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF]. This scheme does not diffuse the mass fraction profile and allows stable computations of the two-fluid model.

Finally we present some two-dimensional numerical results obtained with the RIS scheme.

1 Two-fluid flows and Godunov scheme

We consider a two-fluid model (air and liquid water) written as a first order system of conservation laws

∂ t W + ∂ x F(W ) = 0. (1)
The space variable is x. The time variable is t. We use the notations

∂ t = ∂ /∂t, ∂ x = ∂ /∂ x. The conservative unknowns (x,t) → W (x,t) are W = (ρ, ρu, ρe, ρϕ) T , (2)
with the density ρ, velocity u, total energy e and mass fraction of gas ϕ. The total energy e is related to the internal energy ε by

e = ε + 1 2 u 2 . ( 3
)
The flux of the conservative system is given by

F(W ) = (ρu, ρu 2 + p, (ρe + p)u, ρϕu) T . (4)
The pressure law is of the form

p = P(ρ, ε, ϕ). (5)
If at the initial time t = 0 the mass fraction ϕ takes only two values 0 (pure liquid phase) and 1 (pure gas phase) then it is also true at any later time. This property implies that theoretically it is only necessary to provide the pressure laws P(ρ, ε, 0) for the liquid and P(ρ, ε, 1) for the gas. A classic choice is the stiffened gas pressure law

P(ρ, ε, ϕ) = (γ(ϕ) -1)ρε -γ(ϕ)π(ϕ), (6)
γ(1) = γ 1 > 1, π(1) = π 1 = 0, γ(0) = γ 2 > 1, π(0) = π 2 > 0. ( 7
)
The constants γ 1 , γ 2 and π 2 are obtained from physical measurements. For instance, for air and water, we can take [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF] γ 1 = 1.4, γ 2 = 3, π 2 = 8533 × 10 5 Pa.

In the following, we consider the properties of system ( 1)- [START_REF] Colella | Glimm's Method For Gas Dynamics[END_REF]. In short, we call it the two-fluid model.

In practice, it is difficult to impose ϕ = 0 or ϕ = 1 in the numerical approximation. A widely used possibility (see for instance [START_REF] Barberon | Practical computation of axisymmetrical multifluid flows[END_REF][START_REF] Golay | Numerical schemes for low Mach wave breaking[END_REF][START_REF] Karni | Multicomponent flow calculations by a consistent primitive algorithm[END_REF][START_REF] Müller | Numerical simulation of a single bubble by compressible two-phase fluids[END_REF][START_REF] Mulder | Computing interface motion in compressible gas dynamics[END_REF][START_REF] Saurel | A simple method for compressible multifluid flows[END_REF]) is to interpolate the pressure laws parameter γ(ϕ) and π(ϕ) for ϕ ∈]0, 1[.

For the stiffened gas pressure law, the sound speed is given by

c = γ(ϕ)(p + π(ϕ)) ρ ,
we thus obtain that the system (1)-( 6) is hyperbolic if W is in the hyperbolicity domain

Ω = {(ρ, ρu, ρe, ρϕ), ρ ≥ 0, ϕ ∈ [0, 1], p + π(ϕ) ≥ 0.}
The Riemann problem for the two-fluid system consists in solving the following initial value problem

∂ t W + ∂ x F(W ) = 0, W (x, 0) = W L if x < 0, W R otherwise.
The left and right constant states W L , W R are taken into the hyperbolicity domain Ω . The solution is self-similar, denoted by

R(W L ,W R , x t ) = W (x,t).
It is made of shock waves, rarefaction waves and contact waves, separated by constant states. It is well known that the solution of the Riemann problem is generally not unique. In order to reduce the set of solutions, we can for instance consider only shock waves that satisfy the Lax characteristic criterion. Below, we will also discuss the Lax entropy criterion. An essential feature of the two-fluid system is that the Riemann problem admits a unique global solution that satisfies the Lax characteristic criterion whenever the left and right initial states are in the hyperbolicity domain Ω . This global solution is constructed from standard wave parameterization [START_REF] Lax | Hyperbolic systems of conservation laws[END_REF]. It is sometimes necessary to introduce vacuum state in the gas phase (see [START_REF] Barberon | Practical computation of axisymmetrical multifluid flows[END_REF] for details). The solution is not only theoretical. It can also be computed almost analytically in an efficient way.

The Riemann problem being uniquely solvable, it is then tempting to apply the Godunov scheme to the two-fluid model with arbitrary initial data.

We consider a sequence of time t n , n ∈ N such that the time step τ n = t n+1 -t n > 0. We consider also a space step h. We define the cell centers by x i = ih. The cell C i is the interval ]x i-1/2 , x i+1/2 [. We consider an approximation

W n i ≃ 1 h x∈C i W (x,t n )dx.
A time step of the Godunov scheme is made of two stages:

• Step1: Exact resolution starting from approximated cell averages

∂ t V + ∂ x F(V ) = 0, V (x, 0) = W n i , x ∈ C i . • Step 2: Averaging of the exact solution W n+1 i = 1 h C i V (x, τ n )dx. (8)
The time marching scheme also admits a finite volume formulation

W n+1 i = W n i - h τ n F n i+1/2 -F n i-1/2 .
The numerical fluxes are computed from exact solutions of the Riemann problem

F n i+1/2 = F(R(W n i ,W n i+1 , 0)).
The time step satisfies a CFL condition

τ n ≤ h 2λ max n ,
where λ max n is an upper bound of all the wave speeds in the solutions of the interface Riemann problems at time t n .

In our application, Step 1 of the Godunov scheme is not a problem if all the W n i are in Ω because the Riemann problem admits a unique physically relevant solution W (x,t) ∈ Ω . Surprisingly, Step 2 is much more problematic when π 1 = π 2 , because of the following result [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF]:

Theorem 1.
Consider the two-fluid system (1)-( 6) and suppose a continuous interpolation of the pressure law parameters ϕ → γ(ϕ), ϕ → π(ϕ) for ϕ ∈ [0, 1] satisfying [START_REF] Croisille | Contribution à l'étude théorique et à l'approximation par éléments finis du système hyperbolique de la dynamique des gaz multidimensionnelle et multiespèces[END_REF]. Then, the hyperbolicity set Ω is never convex.

The non-convexity of Ω is a big issue, because even if we have V (x, τ n ) ∈ Ω in the averaging formula [START_REF] Gallouët | A hybrid scheme to compute contact discontinuities in one-dimensional Euler systems[END_REF], we cannot conclude that W n+1 i ∈ Ω . In practice it is possible to construct initial data for which the Godunov scheme fails after only one iteration [START_REF] Müller | Numerical simulation of a single bubble by compressible two-phase fluids[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF]. In conclusion, the Godunov scheme applied to the two-fluid model with π 1 = π 2 is generally unstable. Remark 1. The numerical resolution of the two-fluid model with stiffened gas pressure law has been extensively studied. When the two fluids are perfect gases (when π 1 = π 2 = 0) the hyperbolicity set Ω is convex. In this case, the Godunov scheme is stable. But it often gives very inaccurate results in contact waves. In the literature, this behavior is called the pressure oscillation phenomenon. Let us emphasize that the instability of the Godunov scheme for π 2 = π 1 and the pressure oscillations in contact waves are two different and independent shortcomings of the Godunov scheme applied to two-fluid flows.

Remark 2. A very popular method for avoiding pressure oscillations has been proposed by Abgrall and Saurel in [START_REF] Saurel | A simple method for compressible multifluid flows[END_REF]. The method relies on a non-conservative numerical resolution of the transport equation

∂ t ϕ + u∂ x ϕ = 0,
and a special interpolation of the pressure law coefficients γ(ϕ), π(ϕ) that ensures that the pressure and the velocity remain numerically constant in contact waves. Even if this trick improves the Godunov scheme accuracy, it does not improve the stability. Indeed, it is possible to show that the hyperbolicity set, expressed in the non conservative variables

Ω ′ = {(ρ, ρu, ρe, ϕ), ρ ≥ 0, ϕ ∈ [0, 1], p + π(ϕ) ≥ 0}
is also generally not convex. We can also exhibit physical initial state that leads to the failure of the Abgrall-Saurel scheme after only one iteration [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF].

Remark 3. A Lax entropy W → U(W ) ∈ R ∪ {+∞} is
∂ t U(W ) + ∂ x G(W ) = 0.
If a system of conservation laws admits a Lax entropy, then from Mock's theorem [START_REF] Mock | Systems of conservation laws of mixed type[END_REF], we know that the system is hyperbolic on the domain of U and thus

Ω =DomU = W ∈ R 4 ,U(W ) < +∞ .
Because the domain of a convex function is a convex set, we deduce that the twofluid model cannot possess a global Lax entropy.

In conclusion of this section, we have two alternatives for approximating the two-fluid model in a robust and precise way:

1. We can abandon the stiffened gas pressure law [START_REF] Colella | Glimm's Method For Gas Dynamics[END_REF] and construct another pressure law that ensures the convexity of the hyperbolicity set. 2. If we keep the stiffened gas pressure law we have to construct a scheme that is stable with respect to non-convex hyperbolicity set.

In Section 2, we investigate the first possibility, while in Section 3 we consider the second.

Convex mixture

In this section, we consider a mixture of a perfect gas and a liquid satisfying the stiffened gas pressure law. From physical entropy arguments, we construct a mixture pressure law. This pressure law is naturally associated to a convex Lax entropy of the two-fluid system. Mock's theorem then ensures the convexity of the hyperbolicity set. The construction is split into several steps.

Construction of an extensive mixture entropy

The pressure law is constructed as follows. First, we consider two fluids indexed by i = 1 and i = 2 of mass M i , energy E i , occupying a volume V i . We introduce the specific heat χ i of fluid i. Then, for M i > 0, V i > 0 and E i > π i V i , the entropy function of fluid i is

S i (V i , E i , M i ) = -χ i γ i M i ln M i + χ i (γ i -1)M i lnV i + χ i M i ln(E i -π i V i ). (9)
For completely rigorous proofs, we have to define the entropies for all (V i , E i , M i ) ∈ R 3 . We thus also set

S i (V i , E i , 0) = 0, V i ≥ 0, E i ≥ π i V i , (10) and S
i (V i , E i , 0) = -∞ (11)
in all the other cases. With this definition, S i are concave upper semicontinous (in short: usc) functions and

DomS i = {(V i , E i , M i ),V i > 0, E i > π i V i , M i > 0} ∪ {(V i , E i , 0),V i ≥ 0, E i ≥ π i V i } . ( 12
)
DomS i are convex cones. In addition, the entropies S i are Positively Homegeneous function of degree 1 (in short PH1)

∀λ ≥ 0, ∀W ∈ R 3 , S(λW ) = λ S(W ).
We then define the entropy of the immiscible mixture by

S(V, E, M, M 1 ) = sup V 1 ,E 1 (S 1 (V 1 , E 1 , M 1 ) + S 2 (V -V 1 , E -E 1 , M -M 1 )) . ( 13
)
This formula is physically justified by the fact that the entropy is an additive quantity and by the second principle of thermodynamics: the mixture of the two fluids evolves until it reaches a maximum of entropy. For more details, we refer to [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Helluy | Relaxation models of phase transition flows[END_REF][START_REF] Helluy | Pressure laws and fast Legendre transform[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF] and included references. Let us also observe that we do not optimize the mixture entropy with respect to M 1 because we do not consider phase transition between the liquid and the gas.

From standard convex analysis, it is possible to prove the following result [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF]:

Theorem 2. Let S be defined by [START_REF] Helluy | OpenCL simulations of two-fluid compressible flows with a random choice method[END_REF], where S 1 and S 2 satisfy ( 9)- [START_REF] Golay | Numerical schemes for low Mach wave breaking[END_REF]. Then S is a PH1 concave and usc function. Its domain is a convex cone given by

DomS = {(V, E, M, M 1 ),V > 0, E > 0, M ≥ M 1 ≥ 0} ∪ {(V, E, 0, 0),V ≥ 0, E ≥ 0} . ( 14
)

Intensive mixture entropy and pressure law

From the extensive PH1 entropy, we can go back to intensive variables. We have the following relations

ρ = M V , τ = 1 ρ = V M , ε = E M , ϕ = M 1 M , s = S M , σ = S V .
We then define the intensive specific entropy s(τ, ε, ϕ) = S(τ, ε, 1, ϕ).

From Theorem 2, the specific entropy is a concave function. We define the pressure p, temperature T and chemical potential λ of the mixture by

T = 1 ∂ ε s , p = T ∂ τ s, λ = T ∂ ϕ s, (15)
in such a way that T ds = dε + pdτ + λ dϕ.

We can also consider the volumic entropy

σ (ρ, ρε, ρϕ) = S(1, ρε, ρ, ρϕ) = ρs( 1 ρ , ρε ρ , ρϕ ρ ).
In the same way, the volumic entropy is a concave function of (ρ, ρε, ρϕ).

It is then standard [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Harten | Convex entropies and hyperbolicity for general Euler equations[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF] to deduce that the quantity U(W ) = -σ (ρ, ρε, ρϕ), with W = (ρ, ρu, ρε + 1/2ρu 2 , ρϕ) T , is a convex Lax entropy associated to the entropy flux

G(W ) = uU(W ).
Our whole construction ensures that the two-fluid system with the pressure law given by ( 15) is necessarily hyperbolic on a convex domain and that this convex domain is simply the domain of the Lax entropy U.

Explicit pressure law

It is interesting to perform the full computations in order to see how the resulting pressure P(ρ, ε, ϕ) is different from the interpolated pressure law [START_REF] Colella | Glimm's Method For Gas Dynamics[END_REF]. The computations are not very difficult but a little bit lengthy. They are rigorously detailed in [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF].

We give the final result. Let us just mention that the same formula can be obtained by formally assuming pressure and temperature equilibrium between the two phases

p 1 = p 1 , T 1 = T 2 .
The temperatures T 1 and T 2 of the two stiffened gases are given by the relation

χ i T i = 1/∂ ε i s i (τ i , ε i ) = ε i -π i τ i .
Of course, such equilibrium assumption has no meaning when ϕ = 0 or ϕ = 1 because in this case only one phase is present in the mixture. The entropy optimization procedure is more rigorous and allows handling the cases ϕ = 0 or ϕ = 1.

We take a density ρ > 0, an internal energy ε > 0. We define the heat capacity of the mixture by

χ = χ(ϕ) = χ 1 ϕ + (1 -ϕ)χ 2 .
We also define the energy fraction of the mixture by

ζ = ζ (ϕ) = χ 1 ϕ χ .
The mixture polytropic parameter is then

γ = γ(ϕ) = ζ γ 1 + (1 -ζ )γ 2 .
We also consider the following quantities

δ = -γ 2 π 2 , r = (δ + (γ -1)ρe) 2 -4δ (γ 1 -1)ζ ρe > 0.
The volume fraction of gas is then given by

α = α(ϕ) = δ + (γ -1)ρe - √ r 2δ . 1. If 0 < ϕ < 1 then P(ρ, ε, ϕ) = (γ -1)ρε -γ(1 -α)π 2 .
2. If ϕ = 1 then P(ρ, ε, ϕ) = (γ 1 -1)ρε.

If

ϕ = 0 then P(ρ, ε, 0) = max(0, (γ 2 -1)ρε -γ 2 π 2 ). ( 16
)
The main result of this analysis is that even if the mass of gas vanishes (ϕ = 0) the remaining liquid does not always obey a pure stiffened gas law. When the energy is small enough, the liquid pressure vanishes. Intuitively, this means that the liquid undergoes a cavitation phenomenon. The liquid pressure cannot be negative anymore, while it was possible in the pressure law [START_REF] Colella | Glimm's Method For Gas Dynamics[END_REF].

The pressure law obtained with the entropy optimization procedure ensures a convex hyperbolicity domain. It thus ensures the stability of the Godunov scheme. However, the pressure law is more complex than the stiffened gas law. For instance, in a pure liquid region, the two-fluid model can degenerate to a pressureless Euler system. This system is known to lead to theoretical and numerical difficulties. In addition, we have verified in numerical experiments that the pressure oscillation phenomenon is still present at contact waves. Therefore, we will also propose in the next section a practical numerical method for solving directly the two-fluid model with a non-convex hyperbolicity domain.

A new Random Interface Solver

In this section, we return to the simple stiffened gas pressure law [START_REF] Colella | Glimm's Method For Gas Dynamics[END_REF]. We define the hyperbolicity sets in the pure phases

Ω 0 = Ω ∩ {ϕ = 0}, Ω 1 = Ω ∩ {ϕ = 1}.
These two sets are convex. We consider a numerical initial condition in the pure phases ∀i,W 0 i ∈ Ω 0 ∪ Ω 1 . We show how to construct a scheme, the Random Interface Solver (RIS) that satisfies the stability condition ∀i, ∀n ≥ 0,

W n i ∈ Ω 0 ∪ Ω 1 . (17)
From the literature and the considerations above we know two things:

1. The new scheme cannot be exactly conservative at each time step;

2. If we average the mass fraction on the cells of the initial mesh, we will introduce a numerical diffusion and certainly pressure oscillations at the interface.

In order to avoid these two pitfalls we will 1. use a random sampling strategy at the interface. It allows avoiding the diffusion of the mass fraction profile. It is not perfectly conservative, but we can prove that it is statistically conservative on long times; 2. before the random sampling, we use a Lagrangian conservative finite volume scheme at the interface. In this Lagrangian step, the mass fraction is not diffused either.

We now enter into the details of the RIS scheme. Each time step of the scheme is made of two stages: an Arbitrary Lagrangian Eulerian (ALE) step and a Projection step. The idea to combine the Glimm's scheme approach [START_REF] Glimm | Solutions in the large for nonlinear hyperbolic systems of equations[END_REF] and a Lagrangian scheme approach was first proposed by Chalons and Goatin in [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF]. See also [START_REF] Chalons | Capturing infinitely sharp discrete shock profiles with the Godunov scheme[END_REF][START_REF] Chalons | Computing material fronts with a Lagrange-Projection approach[END_REF][START_REF] Helluy | OpenCL simulations of two-fluid compressible flows with a random choice method[END_REF].

ALE stage

In the first stage, we allow the cell boundary x i+1/2 to move at velocity ξ n i+1/2 . At the end of the first stage, the cell boundary is

x n+1,- i+1/2 = x i+1/2 + τ n ξ n i+1/2 .
Integrating the conservation law (1) on the moving cells, we obtain the following finite volume approximation

h n+1,- i W n+1,- i -hW n i + τ n (F n i+1/2 -F n i-1/2 ) = 0.
The new size of cell i is given by

h n+1,- i = x n+1,- i+1/2 -x n+1,- i-1/2 = h + τ n (ξ n i+1/2 -ξ n i-1/2 ).
The numerical flux is now of the ALE form

F n i+1/2 = F(W n i+1/2 ) -ξ n i+1/2 W n i+1/2 .
The intermediate state W n i+1/2 is obtained by the resolution of a Riemann problem

W n i+1/2 = R(W n i ,W n i+1 , ξ n i+1/2 ).
In practice, R can also be an approximate Riemann solver [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF]. Finally, the interface velocity is defined by

ξ n i+1/2 = u n i+1/2 if (ϕ n i -1/2)(ϕ n i+1 -1/2) < 0, 0 else. ( 18
)
The numerical flux is thus a classic Godunov flux in the pure fluid. It is a Lagrangian numerical flux at the two-fluid interface.

Projection step

The second stage of the time step is needed for returning to the initial mesh. We have to compute on the cells C i of the initial mesh the averages of W n+1,- i , defined on the moved cells C n+1,- i =]x n+1,- i-1/2 , x n+1,- i+1/2 [. Instead of a standard integral averaging method, we rather consider a random sampling averaging process. We consider a pseudo random sequence ω n ∈ [0, 1[ and we perform the following sampling

W n+1 i =        W n+1,- i-1 , if ω n < ξ n i-1/2 τ n h , W n+1,- i , if ξ n i-1/2 τ n h ≤ ω n ≤ 1 + ξ n i+1/2 τ n h , W n+1,- i+1 , if ω n > 1 + ξ n i+1/2 τ n h . (19)
A good choice for the pseudo-random sequence ω n is the (k 1 , k 2 ) van der Corput sequence. In practice, we consider the (5, 3) van der Corput sequence. The resulting scheme has the following properties [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF]:

• it is stable in the sense of ( 17);

• it is conservative in a statistical sense;

• it is entropy dissipative in a statistical sense;

• it does not produce spurious oscillations at the two-fluid interface.

Numerical results

We can extend the scheme to higher dimensions with dimensional splitting (for more details we refer to [START_REF] Helluy | OpenCL simulations of two-fluid compressible flows with a random choice method[END_REF]). It is remarkable that the same random number can be used for one time step in the x and y directions. It is also remarkable that despite dimensional splitting, the two-dimensional scheme converges towards the right solution without oscillation. Indeed, since the work of Colella [START_REF] Colella | Glimm's Method For Gas Dynamics[END_REF], it was generally admitted that applying the dimensional splitting procedure to the Glimm's scheme leads to poor numerical results. We present in Figure 1 the results of a two-dimensional shock-droplet simulation. The initial droplet is a disk. A shock-wave coming from the right of the computational domain interacts with the droplet. The computations have been realized thanks to an OpenCL/MPI implementation of the two-dimensional RIS scheme. For this test case, we use a cluster of four AMD Radeon HD 7970 GPU. The detailed description of the test case is given in [START_REF] Helluy | OpenCL simulations of two-fluid compressible flows with a random choice method[END_REF]. We display the droplet after the interaction. We observe that we are able to capture a sharp interface. The numerical noise is moderate, despite the random nature of the scheme. Because we use a very fine mesh with 20, 000 × 5, 000 cells, we are able to zoom on small Kelvin-Helmholtz vortices (see Figures 2 and3). We present in Figure 4 the results of a two-dimensional shock-bubble simulation. The initial bubble is a disk. A shock wave coming from the left of the computational domain interacts with the bubble. We display the bubble after that it has been split by the shock wave. The computations have been realized thanks to an OpenCL/MPI implementation of the two-dimensional RIS scheme. For this test case, we use a cluster of ten NVIDIA K20 GPU. The detailed description of the test case is given in [START_REF] Helluy | OpenCL simulations of two-fluid compressible flows with a random choice method[END_REF]. We use a very fine mesh with 40, 000 × 20, 000 cells. As in the previous test case, we can zoom in order to observe small details of the split region (see Figures 5 and6).

Conclusion

We have shown that a widely used two-fluid liquid-gas model has a non-convex hyperbolicity domain. This non-convexity can lead to the failure of the Godunov scheme after only one time step. This is true even if the continuous model admits a perfectly well defined solution that satisfies the Lax characteristic criterion. We have proposed a modified equation of state for recovering a convex hyperbolicity domain. The resulting pressure law is more complicated and is not a stiffened gas equation of state anymore in the pure liquid, when the gas mass fraction ϕ = 0.

For keeping the simplicity of the stiffened gas equation, we have thus constructed a finite volume scheme, the RIS scheme, which avoids the numerical diffusion of the mass fraction. The RIS scheme is based on a Lagrangian finite volume approach coupled with a random sampling at the interface.

Let us conclude that the construction of a conservative finite volume scheme that would give accurate results at contact waves for one-fluid or two-fluid general pressure laws is still an open question.

  a strictly convex function of W associated to an entropy flux W → G(W ) such that the smooth solutions of the two-fluid model satisfy the additional conservation law

Fig. 1

 1 Fig. 1 Shock-droplet simulation. Density plot. Full view of the computational domain.

Fig. 2

 2 Fig. 2 Shock-droplet simulation. Density plot. Zoom on the droplet.

Fig. 3

 3 Fig. 3 Shock-droplet simulation. Density plot. Second zoom on the droplet.

Fig. 4

 4 Fig. 4 Shock-bubble simulation. Density plot.


Mohamed Gazibo Karimou

DEGENERATE CONVECTION-DIFFUSION EQUATION WITH A ROBIN BOUNDARY CONDITION

Keywords: 2000 Mathematics Subject Classification. Primary 35F31; Secondary 00A69 Degenerate parabolic equation, Robin boundary condition, Vanishing viscosity approximation, Entropy solution, Semigroup theory

We study a Robin boundary problem for degenerate parabolic equation. We suggest a notion of entropy solution and propose a result of existence and uniqueness. Numerical simulations illustrate some aspects of solution behavior.

Introduction.

Let Ω be an open bounded domain of R with a Lipschitz boundary ∂Ω, and η the unit normal to ∂Ω outward to Ω. The purpose of this paper is to discuss existence and uniqueness of entropy solution for the following initial boundary value problem (P )

   u t + div f (u) -∆φ(u) = 0 in Q =]0, T [×Ω, u(0, x) = u 0 (x)
in Ω, b(u) -(f (u) -∇φ(u)).η = 0 on Σ =]0, T [×∂Ω.

Here, u 0 is taking values on [0, u max ] for some u max > 0. Further, the function f is a Lipschitz continuous function. Moreover, we require that f (0) = 0 and b(0) = 0. (H1)

The diffusion term φ is a continuous function. We consider that there exist a critical value u c of the unknown u such that: φ(.) is zero on [0, u c ] with 0 ≤ u c ≤ u max and φ(.) is strictly increasing else. Then problem (P ) degenerates to hyperbolic when u takes values in the region [0, u c ] where φ is flat. We suppose that the function b is a continuous non-decreasing function on Σ. In some situation, b may be a maximal monotone graph on R (see [START_REF] Andreianov | Well-posedness of general boundary-value problems for scalar conservation laws[END_REF]). Here, we assume also that b satisfies the following hypotheses: b(u max ) ≥ |f (u max ).η|. (H3) For more than a few decades, the degenerate parabolic equation in bounded domain was studied by many authors mainly in the case of Dirichlet boundary conditions (see e.g. [START_REF] Mascia | Nonhomogeneous Dirichlet problems for degenerate parabolic-hyperbolic equation[END_REF], [START_REF] Carrillo | Entropy solutions for nonlinear degenerate problems[END_REF]). The zero-flux boundary condition is studied in [START_REF] Andreianov | Uniqueness for an elliptic-parabolic problem with Neumann boundary condition[END_REF] for non-degenerate parabolic case, in [START_REF] Bürger | On the well-posedness of entropy solution to conservation laws with a zero-flux boundary condition[END_REF] for fully degenerate hyperbolic equation and recently in [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF] for the parabolic-hyperbolic problem. Remark, that the condition b(u) -(f (u) -∇φ(u)).η = 0 on Σ includes in particular Neumann (zero-flux) condition on the boundary.

We propose an adequate entropy formulation for problem (P ) which incorporates two boundary integrals. In [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF], existence and uniqueness for the zero flux boundary condition were proved, under the assumption (H3) that reads f (u max ) = 0 in the zero-flux case b ≡ 0. In contrast to the entropy formulation in [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF], where the passage to the limit in the only boundary integral is straightforward, for our entropy inequality, we need the assumption (H2), which permits to give a sense to the boundary integral with the term b(u). Indeed, we can deduce that b(u) has a trace on the boundary as a function in Sobolev space H 1 (Ω).

The proof of existence of our entropy solution for any space dimensions ≥ 1 employs a vanishing viscosity approximation. We pass to limit in the interior of the domain Q, by using the local compactness result of Panov [START_REF] Yu | On the strong pre-compactness property for entropy solutions of a degenerate elliptic equation with discontinuous flux[END_REF], for this we suppose some relation between f and φ (see Definition 3.5). One can refer to [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF] for more details. We pay particular attention to the boundary term (here (H2) is needed).

For the uniqueness result, we use nonlinear semigroup techniques (see, e.g., [START_REF] Ph | Nonlinear evolution equations in Banach spaces[END_REF]) and Kruzhkov doubling of variables methods. The main goal is to compare two solutions of (P), and it turns out that it is simpler to compare a solution of (P ) with a regular solution (in the sense that the total flux is continuous up to the boundary) of the stationary problem associated to (P ). Then we prove that entropy solution of (P ) is an integral solution, and we refer to the uniqueness of integral solutions granted by the general theory of nonlinear semigroup. Unfortunately, we are not able to obtain regular solution to the stationary problem for any space dimensions, but only in one space dimension. Then, we can deduce the uniqueness just now when Ω is a bounded open interval of R. Notice that, for the same argument as for the zero-flux boundary condition [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF], the problem of uniqueness is still open in multiple space dimensions.

The paper is organized as follows. In the next section, we give our definition of entropy solution and state some remarks useful for the well-possedness. In section 3, we prove existence result of entropy solution. In the section 4, we prove uniqueness in the case of one space dimension. The latter part is devoted to the numerical investigation of problem (P ). We adapt the approach of finite volumes in the spirit of Vovelle ([11]) to illustrate and interpret some observations in the case where the assumptions (H2) and (H3) are absent. Thereby, we justify the importance these assumptions in this paper.

2. Notion of entropy solution. Consider the following notion.

Definition 2.1. A measurable function u taking values on [0, u max ] is called entropy solution of problem (P ) if φ(u) ∈ L 2 (0, T ; H 1 (Ω)), b(u) ∈ L 2 (0, T ; H 1 (Ω)) and the following conditions hold: ∀k ∈ [0, u max ], ∀ξ ∈ C ∞ 0 ([0, T [×R ), with ξ ≥ 0: T 0 Ω |u -k|ξ t + sign(u -k) f (u) -f (k) -∇φ(u) .∇ξ dxdt + Ω |u 0 -k|ξ(0, x)dx + T 0 ∂Ω |f (k).η(x) -b(k)| ξ(t, x)dH -1 dt - T 0 ∂Ω |b(u) -b(k)|ξ(t, x)dH -1 dt ≥ 0. ( 1
)
Here H represents the ( -1)-dimensional Hausdorff measure on ∂Ω.

Remark 1.

1. The entropy solution in the sense of Definition 2.1 is in particular a weak solution. Indeed, first take in inequality (1), k = 0 and use (H1). Next, take k = u max and use (H3). 2. Let us stress that, in particular, the boundary condition (f (u) -∇φ(u)).η = b(u) is verified literally in the weak sense as in the case of zero flux boundary condition (see [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF]). This contrasts with the properties of the Dirichlet problem (see [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF]); we expect that the boundary condition should be relaxed if assumption (H3) is dropped (see [START_REF] Andreianov | Well-posedness of general boundary-value problems for scalar conservation laws[END_REF][START_REF] Andreianov | Scalar conservation laws with nonlinear boundary conditions[END_REF] and also numerical tests of section 5). 3. The integral in the boundary term is well defined due to the hypothesis (H2).

We can use the fact that the trace of b(u)(t, .) ∈ H1 (Ω) on ∂Ω is well defined in L 2 (∂Ω) for t ∈ (0, T ) a.e.

According to the idea of J. Carrillo (cf [START_REF] Carrillo | Entropy solutions for nonlinear degenerate problems[END_REF]), we give an additional property of entropy solutions, useful for the uniqueness techniques.

Proposition 1. Let ξ ∈ C ∞ ([0, T [×R ); then for all k ∈ [φ c , u max ]; for all D ∈ R
and for all entropy solution u of (P ), we have:

T 0 Ω {|u -k|ξ t + sign(u -k)(f (u) -f (k) -∇φ(u) + D).∇ξ} dxdt + Ω |u 0 -k|ξ(0, x)dx - T 0 ∂Ω |b(u) -b(k)| ξ(t, x)dH -1 dt + T 0 ∂Ω |b(k) -(f (k) -D).η(x)| ξ(t, x)dH -1 dt ≥ lim σ→0 1 σ T 0 Q∩{-σ<φ(u)-φ(k)<σ} ∇φ(u).(∇φ(u) -D)ξ(t, x)dxdt. (2)
In general, uniqueness for evolution equation of kind (P ) appear very difficult mainly for the initial boundary values problems. In this context, the use of nonlinear semigroup techniques offers many advantages. Let us present briefly another notion of solution coming from the theory of nonlinear semigroups (see, e.g., [START_REF] Ph | Nonlinear evolution equations in Banach spaces[END_REF]). Definition 2.2. Let A be an m-accretive operator (see, e.g., [START_REF] Ph | Nonlinear evolution equations in Banach spaces[END_REF]). Suppose that

h ∈ L 1 (Q), u 0 ∈ L 1 (Ω). A measurable function v ∈ C([0, T ]; L 1 (Ω; [0, u max ])) 1 is an integral solution of the abstract evolution problem v t + A(v) h(t), v(t = 0) = u 0 , (3)
if v(0, .) = u 0 (.) and for all (u, z)

∈ A d dt ||v(t) -u|| L 1 (Ω) ≤ Ω sign 0 (v(t) -u)(h(t) -z) + {v=u} |h(t) -z| in D (0, T ).
We will see that entropy and integral solution coincide in the case Ω = (a, b) an interval of R.

3.

Existence of entropy solution. The main result of this part is the following: Theorem 3.1. Let ≥ 1. Assume that (H1), (H2) and (H3) holds. Suppose that (f, φ) is non-degenerate (in the sense of Definition 3.5 below). Then, there exists an entropy solution u for the problem (P ).

To show the existence of entropy solutions, we approximate φ(u) by φ (u ) = φ(u ) + Id(u ) for each > 0 and set b (u ) = β • φ (u ). We obtain the following regularized strictly parabolic problem (P ) with unknown u

(P )    u t + div f (u ) -∆φ (u ) = 0 in Q =]0, T [×Ω, u (0, x) = u 0 (x) in Ω, b (u ) -(f (u ) -∇φ (u )).η = 0 on Σ =]0, T [×∂Ω,
where (u 0 ) is a sequence of smooth functions that converges to u 0 a.e and respects the minimum/maximum values of u 0 .

Definition 3.2. Let u 0 be a measurable [0, u max ]-valued function. A measurable function u ∈ L 2 (0, T ; H 1 (Ω)) taking values on [0, u max ] is called weak solution of problem (P ) if : ∀θ ∈ L 2 (0, T ; H 1 (Ω)) ∩ L ∞ (Q) such that θ t ∈ L 2 (Q) and θ(T, .) = 0, one has T 0 Ω {u θ t + (f (u ) -∇φ (u )).∇θ} dxdt + Ω u 0 θ(0, x)dx - T 0 ∂Ω b (u )θdH -1 dt = 0. ( 4
)
Theorem 3.3. For u 0 ∈ [0, u max ], assume (H1), (H2) and (H3) hold. Problem (P ) admits a weak solution u which is also an entropy solution. In particular, we have 0 ≤ u ≤ u max . In addition, there exists C independent on such that

|| √ ∇u || L 2 (Q) ≤ C; (5)
||φ (u )|| L 2 (0,T ;H 1 (Ω)) ≤ C; (6)
||b n (u )|| L 1 (Σ) ≤ C and Σ u b n (u ) ≤ C. (7)
This result can be proved, e.g., using Galerkin method (see e.g. [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF]).

Lemma 3.4. Assume that the sequence (Ψ j ) j is such that:

||Ψ j || L 2 (0,T ;H 1 (Ω)) ≤ C and Ψ j -→ Ψ in L 2 (Q). Then γΨ j -→ γΨ in L 2 (Σ)
, where γ is the trace operator.

The proof uses localization to a small neighborhood of Σ. To prove existence of entropy solution, we assume that the couple (f (.), φ(.)) is non-degenerate in the sense of the following definition: Definition 3.5. (Panov [START_REF] Yu | On the strong pre-compactness property for entropy solutions of a degenerate elliptic equation with discontinuous flux[END_REF]). Let φ be zero on [0, u c ], strictly increasing on [u c , u max ] and a vector f = (f 1 , ..., f ). A couple (f (.), φ(.)) is said to be nondegenerate if, for all ξ ∈ R \{0}, the functions λ -→ i=1 ξ i f i (λ) are not affine on the non-degenerate sub intervals of [0, u c ].

Theorem 3.6. (Panov [START_REF] Yu | On the strong pre-compactness property for entropy solutions of a degenerate elliptic equation with discontinuous flux[END_REF]). Assume that (f, φ) is non degenerate in the sense of Definition 3.5. Suppose u , > 0, is a sequence such that

∃d > 1, ∀s, r ∈ R with s < r T s,r (u ) t + div f (T s,r (u )) -∇φ(T s,r (u )) is pre-compact in W -1,d Loc (Q).
Moreover, suppose u , f (u ), φ (u ) are equi-integrable locally on Q. Then, there exists subsequence (u ) that converges in L 1 Loc (Q).

Proof of Theorem 3.1. (Sketched) The proof of existence of entropy solution uses Theorem 3.6 to justify the passage to the limit in Q (for more details, see [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF]) and Lemma 3.4 for boundary integral.

4.

Uniqueness result of entropy solution in one space dimension. The main result of this section is the following theorem:

Theorem 4.1. Suppose that Ω = (a, b) is a bounded interval of R, then (P ) admits a unique entropy solution.

In order to study uniqueness in the framework of nonlinear semigroup theory, we consider for all bounded function g taking values on [0, u max ], the stationary problem (S) associated to problem (P ):

(S) u + div(f (u) -∇φ(u)) = g in Ω, b(u) -(f (u) -∇φ(u)).η = on ∂Ω.
The notion of entropy solution of (S) correspond to the time-independent entropy solution of (P ) with source term g -u. In the case where Ω = (a, b) is a bounded interval of R, we have an important result, which states that, the total flux is regular at the points a and b. This kind of regularity seem hard to obtain in multiple space dimensions for (S), and even in dimension = 1 for (P ). From now, let's define the operator A f,φ,b on L 1 associated with regular solutions of (S) by its graph: (u, z) ∈ A f,φ,b = u such that u is an entropy solution of (S), with g = u + z .

Proposition 2. For all measurable function g taking values in

[0, u max ] the problem (S) admits a solution u such that (f (u) -φ(u) y ) is continuous up the boundary, i.e., (f (u) -φ(u) y ) ∈ C([a, b]). Moreover, b(u) -(f (u) -φ(u) y ).
Proposition 3. 1. A f,φ,b is accretive in L 1 (Ω). 2. For all λ sufficiently small, R(I + λA f,φ,b ) contains L 1 (Ω; [0, u max ]). 3. D(A f,φ,b ) = L 1 (Ω; [0, u max ]).
For the proof of this proposition, we can refer to [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF]. According to the general results of [START_REF] Ph | Nonlinear evolution equations in Banach spaces[END_REF], it follows existence and uniqueness of integral solution in the sense of Definition 2.1:

Corollary 1. Let Ω = (a, b), u 0 , û0 ∈ L 1 (Ω) and h, ĥ ∈ L 1 (Q). Let v, v be integral solutions of (3) (with operator A f,φ,b
) associated with the data (u 0 , h) and (û 0 , ĥ), respectively. Then for a.e. t ∈ [0, T ).

||v(t) -v(t)|| L 1 ≤ ||u 0 -û0 || L 1 + t 0 ||h(τ ) -ĥ(τ )|| L 1 dt.
Adapted to our case, we have the following result Theorem 4.2. Let Ω = (a, b). Let v be an entropy solution of (P ) and u be an entropy solution of (S). Then

d dt ||v(t) -u|| L 1 (Ω) ≤ Ω sign(v -u)(u -g)dx in D (0, T ). (8)
In particular, v is an integral solution of (3) with h = 0.

Proof of Theorem 4.2 and Theorem 4.1. We consider v = v(t, x) an entropy solution of (P ) and u = u(y) an entropy solution of (S). Consider nonnegative function ξ = ξ(t, x, y) having the property that ξ(., ., y) ∈ C ∞ ([0, T ) × Ω) for each y ∈ Ω, ξ(t, x, .) ∈ C ∞ 0 (Ω) for each (t, x) ∈ [0, T ) × Ω. Apply the doubling of variables [START_REF] Kruzkhov | First order quasi-linear equations in several independent variables[END_REF] in the spirit of [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF], we obtain this following inequality

T 0 Ω Ω |v -u|ξ t dydxdt + Ω Ω
|v 0 -u|ξ(0, x, y)dxdy

+ T 0 Ω Ω sign(v -u) (f (v) -φ(v) x ) -(f (u) + φ(u) y ) .(ξ x + ξ y )dxdydt + T 0 x∈∂Ω Ω |b(u) -(f (u) -φ(u) y ).η(x)| ξdydσdt + T 0 Ω y∈∂Ω |b(v) -(f (v) -φ(v) x ).η(y)| ξdσdxdt + T 0 Ω Ω sign(v -u)(u -g(y))ξdydxdt ≥ T 0 x∈Ω y∈∂Ω |b(u) -b(v)| ξdσdxdt + T 0 y∈Ω x∈∂Ω |b(u) -b(v)| ξdσdydt + lim σ→0 1 σ T 0 Ω c x ×Ω c y ∩{-σ<φ(v)-φ(u)<σ} |φ(v) x -φ(u) y | 2 ξdydxdt ≥ 0. (9)
Next, following the idea of [START_REF] Andreianov | Uniqueness for an elliptic-parabolic problem with Neumann boundary condition[END_REF], we take the test function ξ(t, x, y) = θ(t)ρ n (x, y), where θ ∈ C ∞ 0 (0, T ), θ ≥ 0, ρ n (x, y) = δ n (∆) and ∆ = (1 - with the calculation detailed in [START_REF] Andreianov | Entropy formulation of degenerate parabolic equation with zero-flux boundary condition[END_REF], we deduce that

T 0 Ω Ω θsign(v -u) (f (v) -φ(v) x ) -(f (u) -φ(u) y ) . (ρ n ) x + (ρ n ) y dydxdt → 0.
Hence, we get (8) by passing to the limit in [START_REF] Kruzkhov | First order quasi-linear equations in several independent variables[END_REF] with the above choice of ξ. Thus, the entropy solution v of the problem (P ) is an integral solution of (3). This proves that v is a unique entropy solution due to Corollary 1.

5.

Role of hypotheses (H2), (H3) and some numerical illustrations. The numerical analysis of (P ) is not the aim of this paper, although we consider this alternative in a future work. We assume (H1) holds, u c = 0.6 and u max = 1. We present briefly the importance of the hypotheses (H2), (H3). We apply now the ideas developed e.g., in the work of Vovelle ([11]) to construct a monotone finite volume scheme which take into account the boundary condition. The interval [0, 1] is divided into I cells. We initialize the scheme by: ∀i ∈ {1, ..., I} :

u 0 i = 1 δx iδx (i-1)δx u 0 (x)dx, (10)
the numerical approximation solution at t = nδt in the cell number i ∈ {2, ..., I -1} is :

u n+1 i = u n i - δt δx F (u n i , u n i+1 )-F (u n i-1 , u n i )- φ(u n i+1 )-2φ(u n i )+ φ(u n i-1 ) δx (11)
with the boundary conditions taken into account via Here, F is a numerical flux which we assume monotone, consistent, Lipschitz continuous (see [START_REF] Vovelle | Convergence of finite volume monotones schemes for scalar conservation laws on bounded domains[END_REF]). In the sequel, we take u 0 (x) = 0.7 if x ∈ [ 1 2 , 1] and u 0 (x) = 0 if x ∈ [0, 1 2 [. We take δx = 0.01, δt = δx 2 5 , φ(u) = (u -0.6) + and consider a numerical solution at time t = 0.12. Initially, we remove the hypothesis (H3), by taking f (u) = u 2 2 and b(u) = φ(u), we observe numerically the loss of maximum principle (see Figure 1 ), this mean that the the solution u can be greater than u max . Our entropy formulation requires to choose b(u) in the functional space that permit to define the trace of b(u) on the boundary. In the context where assumption (H2) is not taken into account, b(u) = u and f (u) = u(1 -u)1 [0,1] ; numerically, we observe a boundary layer (see Figure 2 ) and this is confirmed by theoretical results of [START_REF] Andreianov | Well-posedness of general boundary-value problems for scalar conservation laws[END_REF]. Now, taking into account assumptions (H3), (H2), with data f (u) = u(1 -u)1 [0,1] ; b(u) = φ(u) the numerical observation shows that the boundary condition at x = 0 and x = 1 is verified literally and the numerical solution respect the maximum principle (see Figure 3).

u n+1 1 = u n 1 - δt δx F (u n 1 , u n 2 )- φ(u n 2 )-φ(u n 1 ) δx -b(u n 1 ) . (12

b

  = β • φ where β is a non-decreasing Lipschitz continuous function. (H2)

  η(y) is zero at y = a and y = b. (Here η(a) = -1 and η(b) = +1).

  1 n(b-a) )x -y + a+b 2n(b-a) . Then, ρ n ∈ D(Ω × Ω) and ρ n | Ω×∂Ω (x, y) = 0. Due to this choice, T 0 x∈Ω y∈∂Ω |b(v) -(f (v) -φ(v) x ).η(y)| ρ n θdydσdt = 0. By Proposition 2, b(u) -(f (u) -φ(u) y ).η(y) ∈ C 0 ([a, b]). Therefore we have |b(u) -(f (u) -φ(u) y ).η(x)| -→ 0 when x → y, i.e, as n -→ ∞. We conclude that lim n→∞ T 0 x∈∂Ω y∈Ω |b(u) -(f (u) -φ(u) y ).η(x)| ρ n θdydσdt = 0.

  n I ) -F (u n I-1 , u n I )+ φ(u n I )-φ(u n I-1 ) δx .(13)

Figure 1 .

 1 Figure 1.Figure 2.

Figure 2 .

 2 Figure 1.Figure 2.

Figure 3 .

 3 Figure 3.

Here, we will write L 1 (Ω; [0, umax]) for the set of all measurable functions from Ω to [0, umax].

Acknowledgments. I would like to thank Boris Andreianov for his thorough reading and helpful remarks which helped me improve this paper.


Emilien Azéma

Farhang Radjaï

Internal structure of inertial granular flows

Keywords: numbers: 45.70.-n, 83.80.Fg, 45.70.Mg, 81.05.Rm

We analyze inertial granular flows and show that, for all values of the inertial number I, the effective friction coefficient µ arises from three different parameters pertaining to the contact network and force transmission: 1) contact anisotropy, 2) force chain anisotropy and 3) friction mobilization. Our extensive 3D numerical simulations reveal that µ increases with I mainly due to an increasing contact anisotropy and partially by friction mobilization whereas the anisotropy of force chains declines as a result of the destabilizing effect of particle inertia. The contact network undergoes topological transitions, and beyond I ≃ 0.1 the force chains break into clusters immersed in a background "soup" of floating particles. We show that this transition coincides with the divergence of the size of fluidized zones characterized from the local environments of floating particles and a slower increase of µ with I.

From large-scale geological events to a variety of industrial processes involving bulk materials and powders, the flow behavior of granular materials has been a subject of intensive research since many years [1][2][3][4][5][6]. The diversity of boundary conditions and confining geometries made it difficult for a long time to extract the intrinsic rheology of dense inertial flows until a unification was achieved by analyzing several experimental and numerical data in terms of a single dimensionless inertial number I, defined as the ratio of the particle relaxation time (m/pd) 1/2 , under an average or confining stress p and for a particle of mass m and diameter d, to shear time γ-1 imposed by the flow rate γ [7,8]. The model arising from this seminal work is based on a generic dependence of the effective friction coefficient µ and packing fraction ν on I.

This empirical model of steady granular flows, in combination with continuum conservation equations, correctly predicts the velocity and stress fields in various flow geometries [7,[9][10][11]. However, it still lacks a clear particle-scale foundation. The increase of µ with I despite an increasingly lower packing fraction is a nontrivial property that reveals a genuine microstructure. A few studies reported on the microstructure of inertial flows [8,[12][13][14][15] show that, as I increases, the force chains become more sparse, the correlation length of connected particles decrease [5,16,17] , the contact lifetimes decline, and an increasing number of impulsive forces [18,19] and frictionally mobilized contacts [4,8] come into play. But a challenging issue is how to connect such particle-scale observations with the rheology.

In this Letter, we analyze inertial granular flows by means of a stress partition that readily links the µ(I) rheology to three different particle-scale mechanisms: 1) contact anisotropy, 2) force chain anisotropy and 3) fric- * Electronic address:

emilien.azema@univ-montp2.fr;franck. radjai@univ-montp2.fr tion mobilization. The data are obtained from extensive contact dynamics simulations of homogeneous shear flow with spherical particles for a broad range of I varied from low to very high values. Stress partition has been successfully applied to quasi-static deformations [20][21][22][23][24] but never to inertial flows. As we shall see below, the relative importance of local mechanisms evolves with I and therefore the flow structure undergoes qualitative transitions that underlie the evolution of the effective friction. Such transitions are consistent with the correlative evolution of the statistics of fluidized zones that we analyze by considering the clustering of floating particles.

Contact dynamics simulations were carried out with N p = 24000 spheres. The packing is sheared between two parallel rough walls by imposing a constant horizontal velocity V y on the top wall and periodic boundary conditions along the flow in the y direction and along the transversal x direction. The packing is confined by a constant compressive stress σ zz applied on the top wall along the z-direction; see Fig. 5(a). A small polydispersity (1 ± 0.1)d in particle diameters is introduced to prevent crystallization at the walls. The gravity is set to zero in order to avoid strain localization at the boundaries and to ensure uniform stress field in the bulk. In the contact dynamics method, the particles are treated as perfectly rigid so that I is the only relevant dimensionless parameter of flow [24][25][26]. Contact dissipation is modeled in terms of normal and tangential restitution coefficients e n and e t as well as a friction coefficient µ s between particles. We set µ s = 0.4 and e n = e t = 0. This choice corresponds to a highly dissipative packing. We note that the values of restitution coefficients have nearly no influence on dense granular flows as inelastic collisions occur at high frequency and dissipate the kinetic energy at time scales far shorter than those of shear and particle relaxation [7,8].

Under the action of the applied shear strain, all samples dilate from their initial high density and tend to a steady flow characterized by a linear velocity profile. We performed 17 simulations for a broad range of I varied from 6 × 10 -4 to 0.72 by keeping the same shear rate and changing the confining pressure [36]. The data presented in this Letter, are average values over the steady state with standard deviations used as error bars. The values of I are obtained from the average stress p, which fluctuates in the steady state around the average stress. The error bars on the values of I represent these fluctuations. Obviously, the fluctuations increase with I, and thus the data are to be sampled more frequently in the steady state in order to reach meaningful statistics.

The stress state being invariant along the x direction, we consider here only the restriction σ of the stress to the shear plane yz. It may be expressed as [START_REF] Moreau | Friction, Arching, Contact Dynamics[END_REF]

σ αβ = n c f c α ℓ c β , (1)
where n c is the number density of the contacts and the average is taken over the contacts c with contact force component f c α and branch vector ℓ c β joining the centers of contacting particles. Note that, the contribution of particle velocity fluctuations ( mv α v β /V ) remains very small compared to that of contact forces in all the simulations considered here. According to the Mohr-Coulomb model, the effective friction coefficient during shear is given by µ ≃ q/p, where p = (σ 1 + σ 2 )/2 is the mean stress and q = (σ 1 -σ 2 )/2 is the stress deviator, and σ 1 and σ 2 are the principal values on the shear plane [START_REF] Andreotti | Granular Media: Between Fluid and Solid[END_REF].

Fig. 1 displays µ and ν obtained from our simulations and a compilation of available published numerical and experimental data from several authors for different boundary conditions as a function of I [1-5, 29] [37]. We see that our numerical data collapse well with all other data. The effective friction coefficient increases and tends to saturate with increasing I whereas the packing fraction declines from 0.59 in the quasi-static state to 0.50 for our highest values of I. Note that 0.59 corresponds to the density of a packing of frictional spheres under continuous quasi-static shearing, as also evidenced by experiments [30].

Relying on our numerical data, we now focus on the stress partition in connection with the effective friction. Let us express the average in (1) as an integral:

σ αβ = n c f α ℓ β P ℓf n df dℓ dn, (2)
where P ℓf n is the joint probability density of forces and branch vectors ℓ = ℓn projected on the shear plane. At the lowest-order description of the microstructure, we neglect the force-fabric correlations and split P as a product P ℓf n = P ℓ (ℓ)P f (f )P n (n). Integration over f and ℓ yields:

σ αβ ≃ n c ℓ 0 Ω f α (n)n β P n (n)dn, ( 3
)
where Ω is the angular domain of integration and f (n) is the average force as a function of n and ℓ = ℓ 0 . The contact force on the shear plane can be decomposed into its normal and tangential components f n (n) and f t (n), and n is parametrized by its orientation θ. The three functions P (θ), f n (θ) and f t (θ) are πperiodic and, as shown in Fig. 3, they can be well approximated by their lowest-order Fourier expansions [20][21][22][23]31]:

P (θ) ≃ 1/π{1 + a c cos 2(θ -θ c )}, f n (θ) ≃ f n {1 + a n cos 2(θ -θ n )}, f t (θ) ≃ -f n a t sin(θ -θ t ), (4)
where a c , a n , and a t are anisotropy parameters, and θ c ≃ θ n ≃ θ t are the corresponding privileged directions on the shear plane, and nearly coincide with the major principal stress direction in the steady state. Now, introducing Eqs. (4) in the integral (2), and neglecting the cross products of the anisotropy parameters, one gets the simple relation:

µ ≃ 1 2 (a c + a n + a t ). (5)
The predicted values of µ by this equation from the anisotropy parameters are shown in Fig. 2 together with those obtained from the stress tensor as a function of I.

We see that Eq. ( 5) approximates excellently the effective friction for all values of I. This result indicates also that the expression (1) of the stress tensor holds correctly for high inertia where impulsive forces prevail. The evolution of the three anisotropies with I is plotted in Fig. 3. Interestingly, the normal force anisotropy a n decreases and tends to a constant value whereas the friction force anisotropy slightly increases with I. At the same time, the contact anisotropy a c is a sharply increasing and nonlinear function of I. This means that, since by virtue of Eq. ( 5) the three anisotropies add up to build the effective friction, the contact anisotropy is the principal microstructural cause of the increase of effective friction as a function of the inertial number. In general, a c varies oppositely with the coordination number z, which declines in our simulations from 4.5 to 1.1 as I varies from 10 -4 to 0.7. This is because the contact anisotropy is mainly a consequence of the loss of contacts in the extension direction [32]. The normal force anisotropy a n reflects the force chains, which are increasingly destabilized by particle inertia causing a n to decrease. Hence, the friction force anisotropy a t , which reflects friction mobilization ( |f t | / f n ∝ a t ), grows to re-stabilize the force chains and thus take more actively part in force transmission [23].

Since the contact anisotropy and coordination number seem to be sensibly important with respect to the effective friction, we now turn to the connectivity of the contact network in order to obtain morphological clues to the evolution of the force network. The connectivity can be described by the proportion P k of force-bearing particles having k ≥ 1 contacts and the proportion P 0 of floating particles, i.e. not participating in the force network. Note that k≥1 P k = 1 and k≥1 kP k = z. Fig. 4 displays P k for k = 1 to 8 and P 0 as a function of I. We observe several nontrivial topological transitions. For I < 0.01, the effect of inertia leads to the reduction of highly-connected particles (k > 4) in favor of low-connected particles (k < 4) while P 4 remains nearly constant. In this range, P 0 ≃ 0.07. At higher values of I, the force network is further destabilized and P 4 begins to decline whereas P 3 keeps increasing up to a peak value at I ≃ 0.1 beyond which P 3 also begins to decline whereas P 2 continues to increase. The loss of particles having 3 contacts is a dramatic change in the microstructure as multiple contacts may thereafter occur mainly in the form of linear force chains without branching, as observed in Fig. 6(d). P 0 increases slightly in this range from 0.07 to 0.2 but undergoes a sharp increase beyond I ≃ 0.1. For I > 0.3, the flow is dominated by P 2 and P 1 (corresponding to the collisions of floating particles). The sharpt transition observed at I ≃ 0.1 for P 0 coincides with the slowing down of µ(I); see Fig. 1 on linear scale. At the same time, the normal force anisotropy ceases to decline and takes a constant value a n ≃ 0.2; see Fig. 3. The fact that beyond I ≃ 0.1 the force anisotropy does not follow the contact anisotropy indicates that the force chain formation during flow is hindered by collisions, and the contacts do not live long enough to sustain percolating force chains. Indeed, at I ≃ 0.1 a fraction P * 0 = 0.2 of particles is floating and this fraction of particles disconnected from the force network is consistent with the site percolation threshold of a hexagonal compact packing [33]. This means that, beyond I ≃ 0.1, the effect of inertia may be analyzed more sensitively in terms of the evolution of floating particles rather than connected chains, whose correlations have already been a subject of several studies to characterize the transition from inertial regime to quasi-static regime [13,16].

The floating particles actually provide a complementary picture in terms of fluidized zones defined from the neighborhoods of floating particles. For each floating particle i, let P i 0 (I, R) be the proportion of floating particles for a spherical volume of radius R centered on i. We determine the size R i of the fluidized zone by requiring P i 0 (I, R i ) = P * 0 . In words, this is the size of the spherical volume in which the density of floating particles is equal to the percolation threshold P * 0 . The distribution P (R) of the sizes of fluidized zones for given value of I can be evaluated from the set {R i } and, their mean size R mean is determined by double averaging R i over all floating particles and during shear. This construction is nearly dual to the 2-point cluster function introduced by Torquato et al. [34,[START_REF] Torquato | Random heterogeneous materials[END_REF].

Figure 5 shows four snapshots of floating particles for different values of I. At low I the floating particles are mostly isolated, but they tend to cluster into fluidized zones as I increases. Fig. 6 displays the clustering length R mean as a function of I. The inset shows the statistical distributions of R i for several values of I. We see that R mean is nearly constant and equal to 2d for I < 0.1, as expected for isolated floating particles. Thereafter, R mean grows rapidly with I and reaches the system size for I ≃ 0.3. This evolution reflects the coalescence process of fluidized zones and a transition to the collisional regime at I ≃ 0.3 where P 0 > 0.6. This process is also clearly evidenced by the evolution of the size distributions of fluidized zones, which broaden as I increases up to I ≃ 0.1. For I > 0.1, the distribution develops a second peak for R equal to the system size whereas the mean size of the fluidized zones continues to increase.

To summarize, in the range I < 0.1 the force anisotropy reflects percolating force chains, which are progressively destabilized by inertial effetcs with increasing I. As a result, the force anisotropy declines towards a residual constant value a n ≃ 0.2 for I > 0.1. This residual force anisotropy is essentially induced by collisions due to shearing, which begin to affect the microstructure in the form of fluidized zones of increasing size that coalesce at I = 0.3. Hence, in the whole range I > 0.1, the microstructure may be described as composed of shortlived and impulsive force chains embedded in a "soup" of floating particles. The contact anisotropy grows due to enduring force chains and by loss of contacts for I < 0.1 and due to both impulsive force chains and shear-induced collisions of floating particles beyond I = 0.1. This geometrical anisotropy provides the main additive contribution to the effective shear friction of the flow according to Eq. ( 5). The transition to a fully collisional regime occurs at I ≃ 0.3 where the whole system is in a fluidized state. Let us also note that the value I = 0.3 may be identified with the reference value I 0 in the fitting form introduced in [9,10] and shown in Fig. 1.

The above picture reveals the highly nonlinear evolution of the microstructure with increasing inertial number that we analyzed in terms of three anisotropy parameters, which underly additively the effective friction. The dynamics of local structures may be further investigated by considering in detail the contact lifetimes, the process of contact gain and loss and the role of impulsive forces, as discussed in [18,19] regarding the jamming transition in a suspension.

FIG. 1 :

 1 FIG. 1: (Color online) Effective friction coefficient (a) and packing fraction (b) as a function of I. The data analyzed in this paper are in blue triangles. The other data are extracted from [1-5, 29]. The dashed line shows the fitting form µ = µ0 + µ 1 -µ 0 1+I 0 /I introduced in [9, 10].

FIG. 2 :FIG. 3 :

 23 FIG.2:(Color Online) Effective friction coefficient µ (black circles) together with its harmonic approximation (Eq.5) (red squares) as a function of inertial number I.

FIG. 4 :FIG. 5 :

 45 FIG. 4: (Color online) Proportion of floating particles P0 and the connectivity P k of the particles (inset) as a function of I.

4 FIG. 6 :

 46 FIG. 6: Mean size of critical floating clusters as a function of I. The inset shows the statistical distributions of R i for several values of I.