Thierry Gallouët

Philippe Helluy

Jean-Marc Hérard

Julien Nussbaum

Thierry Gallou Ët

Jean-Marc H Érard

Hyperbolic relaxation models for granular flows

Keywords: 1991 Mathematics Subject Classification. 76M12, 65M12 two-phase flow, hyperbolicity, relaxation, finite volume, entropy

HAL is

Introduction

We are interested in the numerical modeling of a two-phase (granular-gas) flow with two velocities and two pressures p 1 and p 2 . In one space dimension, the model is made up of seven non-homogeneous partial differential equations: two mass balance laws, two momentum balance laws, two energy balance laws and one volume fraction evolution equation. It is similar to the initial model proposed by Baer-Nunziato [START_REF] Baer | A two phase mixture theory for the deflagration to detonation transition (ddt) in reactive granular materials[END_REF]. The main feature of this model is that the left hand side of the equations is hyperbolic. This property is very important because it ensures the mathematical stability of the model.

However, in many industrial applications it is not realistic to admit two independent pressures. Generally, an algebraic relation between the two pressures is assumed. An example (among many others) of such a modeling in the framework of internal ballistics is given by Gough in [START_REF] Gough | Modeling of two-phase flows in guns[END_REF]. For a general presentation of twophase flow models, we refer to the book of Gidaspow [START_REF] Gidaspow | Continuum and kinetic theory descriptions[END_REF]. The relation between the two pressures is classically of the form p 2 = p 1 + R where R is the granular stress. In the general case, the granular stress depends on all the thermodynamic variables of the two phases.

Because of the pressure relation, the system is now overdetermined. The volume fraction equation can be eliminated and a six-equation model is obtained. Unfortunately, the new model has a reduced hyperbolicity domain. The worst situation corresponds to a vanishing granular stress R = 0. In this case, the model is almost never hyperbolic.

In the case of a vanishing granular stress, several authors have proposed to relax the algebraic relation p 2 = p 1 by adding an adequate source term to the volume fraction evolution [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF], [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF], [START_REF] Hurisse | Couplage interfacial instationnaire de modèles diphasiques[END_REF], etc. An important parameter of the source term is the characteristic equilibrium time. When the equilibration time tends to zero, the six equation model is recovered. When the equilibrium time is positive the stability of the model is improved.

In this paper, we extend the relaxation approach to non-vanishing granular stress R > 0. With a positive granular stress, the hyperbolicity domain of the six-equation model is slightly extended. But it is generally not possible to remove all the elliptic regions. Therefore, we apply a relaxation source term to the volume fraction evolution equation of the seven-equation model. This source term takes into account the granular stress (see (20)). When the relaxation parameter τ p tends to 0, we recover the equilibrium six-equation model. When the relaxation parameter τ p > 0, the stability of the model is enhanced.

An important aspect of the model is that the granular stress cannot be chosen arbitrarily once the pressure law of the solid phase is fixed. Indeed, it has to satisfy some thermodynamical relations in order that an entropy dissipation equation can be established. We illustrate this fact when the equation of state of the solid phase is a stiffened gas law. We propose a very simple but useful expression (see ( 14)) for the granular stress, which is mathematically and physically coherent.

Then we propose a numerical method to solve the relaxed system. For that purpose, we use a splitting algorithm. Each time step of the algorithm consists in:

• evolving the seven equation model without the source term;

• solving the relaxed pressure equilibrium with granular stress;

• solving the other source terms.

In the second stage the relaxed pressure equilibrium implies to solve an update for the volume fraction, keeping constant the conserved variables. Under some monotony hypothesis on the granular stress, and when the pressure laws of the two phases are stiffened gas laws, we are able to prove the existence and uniqueness of the new volume fraction in the interval ]0, 1[. The relaxed equilibrium becomes an exact equilibrium when the relaxation parameter τ p = 0. Thus our method can also be used to solve the equilibrium sixequation model on coarse meshes, which is important in industrial applications.

Finally, we propose some numerical experiments. For academic test cases, we demonstrate some behavior of the relaxed approach in the case of a non-stable (elliptic) case. We then compare the results of our new approach with the standard Gough model [START_REF] Nussbaum | Numerical simulations of gas-particle flows with combustion[END_REF] in the case of a simplified internal ballistics problem.

Notations and model

We consider a two-phase flow of a granular solid mixed with a compressible gas. The solid is denoted by the index (2) and the gas by the index [START_REF] Baer | A two phase mixture theory for the deflagration to detonation transition (ddt) in reactive granular materials[END_REF]. For more generality, the solid is supposed to be compressible. The unknowns are, for each phase k = 1, 2, the partial density ρ k , the velocity u k , the internal energy e k . The volume fractions α k satisfy α 1 + α 2 = 1. The solid volume fraction α 2 is also called the porosity in the context of granular flows. The pressure of each phase is given by an equation of state of the form (1)

p k = p k (ρ k , e k ).
We note

α k ρ k = m k .
The balance of mass, momentum and energy read

(2)

m k,t + (m k u k ) x = ±M, (m k u k ) t + (m k u 2 k + α k p k ) x -p I α k,x = ±Q, (m k E k ) t + ((m k E k + α k p k )u k ) x + p I α k,t = ±S, α k,t + v I α k,x = ±P, where (3) E k = e k + u 2 k 2 .
The right hand side terms M , Q, P , S are internal exchange source terms that will be discussed later. Here, the sign ± = + if k = 1 and ± = -if k = 2. For the moment, we suppose that there are no external force and energy source (this explains the ± signs in the source terms). The quantities p I and v I are respectively the interface pressure and the interface velocity. In this paper, we take the special choice of Baer-Nunziato ( 4)

p I = p 1 , v I = u 2 .
which enjoys good properties (see [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF], [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: reduced equations[END_REF]):

• the left hand side of the system is hyperbolic (the proof is recalled in Section 7); • this choice ensures that the non-conservative products are well defined. This is due to the fact that the volume fraction only jumps in linearly degenerated fields. In a linearly degenerated field, the jump relations are simply provided by the Riemann invariants of this field. See [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF]; • in the applications, the Baer-Nunziato model is particularly adapted to granular flows.

Using the volume fraction equation, the time derivative α k,t can be replaced by space derivatives (5)

m k,t + (m k u k ) x = ±M, (m k u k ) t + (m k u 2 k ) x + α k p k,x + (p k -p I )α k,x = ±Q, (m k E k ) t + (m k E k u k ) x + α k (p k u k ) x + (p k u k -p I v I )α k,x = ±S ∓ p I P, α k,t + v I α k,x = ±P.
The equations can then be written under the form of a first order non-conservative and non-homogeneous system (6)

W t + F (W ) x + A(W )L(W ) x = S(W ), with (7)
W = (m 1 , m 1 u 1 , m 1 E 1 , m 2 , m 2 u 2 , m 2 E 2 , α 1 ) T , F (W ) = m 1 u 1 , m 1 u 2 1 , m 1 E 1 u 1 , m 2 u 2 , m 2 u 2 2 , m 2 E 2 u 2 , 0 T , L(W ) = (p 1 , p 1 u 1 , p 2 , p 2 u 2 , α 1 ) T , A(W )L(W ) x = (0, α 1 p 1,x , α 1 (p 1 u 1 ) x + p 1 (u 1 -u 2 )α 1,x , 0, α 2 p 2,x + (p 2 -p 1 )α 2,x , α 2 (p 2 u 2 ) x + u 2 (p 2 -p 1 )α 2,x , v 2 α 1,x ) T , S(W ) = (M, Q, S -p 1 P, -M, -Q, -S + p 1 P, P ) T .
Let us note that this writing is not unique. We have chosen a formulation in which the non-conservative terms vanish when the pressures and velocities are constant, i.e. [START_REF] Hurisse | Couplage interfacial instationnaire de modèles diphasiques[END_REF] p 1 = p 2 = p 0 = Cst and u 1 = u 2 = u 0 = Cst.

Entropy dissipation

In this section, we establish an entropy dissipation equation. This equation is very important because it permits to select the source terms that are compatible with the second principle of thermodynamics.

For that purpose we first rewrite the system as follows (9)

m k (u k,t + u k u k,x ) + (α k p k ) x -p 1 α k,x = ±Q ∓ u k M, m k (e k,t + u k e k,x ) + u k 2 (m k u k ) t + (m k u 2 k + α k p k ) x + (α k p k ) x + 1 2 m k u k (u k,t + u k u k,x ) + α k p k u k,x + p 1 α k,t = ±S ∓ e k M.
The last equation also reads [START_REF] Nussbaum | Modélisation et simulation numérique d'un écoulement diphasique de la balistique intérieure[END_REF]

m k (e k,t + u k e k,x ) + u k 2 [p 1 α k,x ± Q + (α k p k ) x ] + 1 2 u k (-(α k p k ) x + p 1 α k,x ± Q ∓ u k M ) + α k p k u k,x + p 1 α k,t = ±S ∓ e k M. and (11)
m k (e k,t + u k e k,x ) + α k p k u k,x + p 1 (u k -u 2 ) α k,x = ±S ∓ e k M ∓ u k Q ± 1 2 u 2 k M ∓ p 1 P
Finally, we obtain (12)

α k,t + u 2 α k,x = ±P, α k (ρ k,t + u k ρ k,x ) + ρ k (u k -u 2 )α k,x + m k u k,x = ±M ∓ ρ k P, m k (u k,t + u k u k,x ) + (α k p k ) x -p 1 α k,x = ±Q ∓ u k M, m k (e k,t + u k e k,x ) + α k p k u k,x + p 1 (u k -u 2 ) α k,x = ±S ∓ e k M ∓ u k Q ± 1 2 u 2 k M ∓ p 1 P
Now, we introduce entropies s k for the two phases. The entropy of the phase k = 1, 2 satisfies the following first order partial differential equation (13)

T k ds k = de k - p k ρ 2 k dρ k - R k m k dα k = T k ∂s k ∂e k de k + T k ∂s k ∂ρ k dρ k + T k ∂s k ∂α k dα k .
The temperature of phase k is denoted by T k . The granular stress in the phase k is denoted by R k . Physically, the term -R k /m k dα k represents the work of the granular stress due to a change of volume dα k . For sake of simplicity and in order to be more realistic, we assume that the "granular" stress R 1 vanishes in the gas phase, thus R 1 = 0. Without ambiguity, we can also denote the granular stress in the solid phase by R = R 2 . Generally, the granular stress R depends on (ρ 2 , e 2 , α 2 ). When the solid phase is modeled by a stiffened gas equation of state (with a polytropic parameter γ 2 ), a very nice choice is

(14) R = R(ρ 2 , α 2 ) = λ (ρ 2 α 2 ) γ2 .
This choice is discussed in Section 3. The chemical potential of phase k is noted µ k and is defined by

(15) µ k = e k + p k ρ k -T k s k
We multiply the last equation in (??) by 1/T k , the second by -p k /ρ k /T k , the first by -R k /T k and take the sum (16)

m k (s k,t + u k s k,x ) + p 1 -p k T k (u k -u 2 )α k,x = 1 T k ±P (-R k + p k -p 1 ) ± M ( u 2 k 2 - p k ρ k -e k ) ± Q(-u k ) ± S , (m k s k ) t + (m k u k s k ) x = 1 T k ±P (p k -R k -p 1 ) ± M (T k s k + u 2 k 2 - p k ρ k -e k ) ± Q(-u k ) ± S .
Adding now the two entropy equations leads to the entropy dissipation partial differential equation that we sum up in the following proposition.

Proposition 1. Consider a smooth solution of the system (2) and two entropy functions s 1 and s 2 satisfying (13). Then, the smooth solution satisfies the following entropy dissipation equation

(17) ( m k s k ) t + ( m k u k s k ) x = P T 2 (p 1 + R -p 2 ) + M u 2 1 2T 1 - u 2 2 2T 2 - µ 1 T 1 + µ 2 T 2 + Q u 2 T 2 - u 1 T 1 + S 1 T 1 - 1 T 2 .
Remark 1. According to the second principle of thermodynamics the right hand side of (17) has to be non-negative. But each term in the formula (17) has not a clear physical meaning. It is often more convenient to rewrite the source term in a different way. For example, we can set

(18) Q = Q 0 + u 1 M, S = S 0 + u 1 Q 0 + u 2 1 2 M + µ 1 M.
In this way, the dissipation rate becomes

(19) P T 2 (p 1 + R -p 2 ) + M T 2 (µ 2 -µ 1 - (u 2 -u 1 ) 2 2 ) + Q 0 T 2 (u 2 -u 1 ) + S 0 T 1 T 2 (T 2 -T 1 ).
It is > 0 if each term in the sum is > 0. The source S 0 can then be interpreted as the heat flux (it is > 0 when T 2 > T 1 , i.e. when the phase 2 heats the phase 1). The source Q 0 is the drag force. Finally, M is the mass transfer due to chemical reaction. When u 1 = u 2 , we recover that the chemical reaction tends to create the phase with the smallest chemical potential.

Remark 2. Generally, the equations (13) satisfied by the entropies s k have not a unique solution, once the pressure laws are given. For example, if s k is a solution, -s k is also a solution. A supplementary condition has thus to be given in order to fix the sign of the entropy dissipation rate. In the case of conservative systems the entropies are supposed to satisfy some convexity property. For a non-conservative system, it is not possible to apply the Godunov-Mock theorem and it is difficult to extend naturally the convexity approach. We propose here only to forbid the change s → -s by imposing that the temperature remains > 0. It implies

(20) 1 T k = ∂s k ∂e k > 0.
In this paper, we will concentrate on the pressure relaxation source term. We will assume the following form, which ensures a positive entropy dissipation

(21) P = 1 τ p α 1 α 2 (p 1 + R -p 2 ) , τ p > 0,
where τ p is the relaxation parameter. An instantaneous relaxation corresponds to the limit τ p → 0.

Remark 3. Let p ref be a reference pressure. We can define a characteristic time for the pressure equilibrium by

(22) t ref = τ p p ref .
The knowledge of this characteristic time is important for a proper modeling.

3. Application to stiffened gas laws: admissible granular stress 3.1. Admissible granular stress. In many papers, the granular stress is supposed to depend only on the solid volume fraction α 2 . This hypothesis is reasonable when the solid phase is incompressible. However, this choice is not compatible with the existence of an entropy satisfying (13) in the case of a compressible phase. The choice of the granular stress expression cannot be arbitrary. That is why in this section, we compute a simple expression of it when the pressure law of the solid phase is a stiffened gas equation of state.

Let us note

(23) Θ(1/ρ 2 , e, α 2 ) = R 2 (1/ρ 2 , e, α 2 ) α 2 ρ 2 .
We omit now the subscript k = 2 because we concentrate only on the solid phase. We have to find an entropy s, a temperature T and a function Θ (containing the granular modeling) such that (24)

T ds = de -p ρ 2 dρ -Θdα.

We note

(25) τ = 1/ρ, T = 1/ϕ.
Then, ϕ = ϕ(τ, e, α) is an integrating factor for the form (26) de + pdτ -Θdα, which reads (27) ds = ϕde + ϕp(τ, e)dτ -ϕΘdα.

In order to construct a practical and simple model, we suppose that the granular stress only depends on the density and the volume fraction of the solid phase. This leads to the choice Θ = Θ(τ, α). The differential form is closed if (28)

ϕ α = -Θϕ e , pϕ α = -Θϕ τ -ϕΘ τ , pϕ e + p e ϕ = ϕ τ .
The general case corresponds to ϕ e = 0, ϕ α = 0 and ϕ = 0. We then have necessarily

(29) Θ τ Θ = -p e .
3.2. Practical example. Now we propose some computations when the pressure law is a stiffened gas EOS (30) p(ρ, e) = (γ -1)ρe -γπ.

The parameter γ must be > 1. The parameter π has the dimension of a pressure and can be arbitrary. But in practice, for a solid phase, it is positive and large compared to a characteristic pressure of the low. In the case of a stiffened gas equation, we thus find

(31) Θ τ Θ = - γ -1 τ ⇒ Θ(1/ρ, α) = θ(α)ρ γ-1 , which leads to (32) R(τ, α) = αθ(α)ρ γ .
In this paper, we will perform numerical experiments with a very simple particular choice

(33) R(τ, α) = λρ γ α γ .
The value of λ can be adjusted to experiments. With this choice, the parameter R 0 = λρ γ2 2 has the dimension of a pressure. It represents the maximal stress corresponding to the maximal compaction α 2 = 1. This model is not so different from classical approaches (as described for example in [START_REF] Goldshtein | Mechanics of colisional motion of granular materials. part 3: Self similar shock wave propagation[END_REF]): usually, the granular stress vanishes under some critical porosity α c (dilute case) and increases with α when α > α c (packed case). In our approach, the parameter γ allows to ensure that the granular stress is small when α is small. Actually, the higher γ is, the faster the granular stress tends to zero when α tends to zero.

Of course, it would be also possible to consider the most general case where the granular stress also depends on the internal energy (34) Θ = Θ(α, τ, e).

However, we will see that the choice (33) is very interesting for the modeling and the numerics because it ensures a maximum principle on the volume fraction during the pressure equilibrium resolution. Is also possible to compute the whole thermodynamic underlying model. The expressions for the associated entropy and temperature are given in Section 8.

Finite volume approach

For the numerical implementation, we consider a finite volume discretisation and a splitting approach. The convection step is solved by a standard Rusanov scheme already described in many works as [START_REF] Nussbaum | Numerical simulations of gas-particle flows with combustion[END_REF]. In the second stage, the source terms are applied. We concentrate on the pressure relaxation source term in the next section.

For the finite volume scheme, we consider a step h, a time step ∆t. The cells are intervals ]x i-1/2 , x i+1/2 [. For simplicity, we consider a regular mesh x i = ih (but this of course is not mandatory). The vector W is approximated in each cell at time t n by (35)

W n i ≃ W (x i , t n ). The numerical scheme for the convective term is a standard Rusanov scheme for non-conservative systems, which reads

(36) h(W n+1/2 i -W n i ) + ∆t(F n i+1/2 -F n i-1/2 ) + ∆tA(W n i ) L n i+1 -L n i-1 2 = 0.
The conservative numerical flux is given by ( 37)

F (W L , W R ) = F (W L ) + F (W R ) 2 -ζ W R -W L 2 , ζ = max (ρ(B(Y L )), ρ(B(Y R ))) .
where ρ(B) denotes the spectral radius of the matrix B. The convection matrix B(Y ) in the primitive variables Y is given in Section 7.

Our particular choice of the non-conservative terms L ensures that constant velocity-pressure states will be maintained by the Rusanov scheme.

Relaxation algorithm

In this section, we address now the numerical approximation of the pressure relaxation source term. As usual, we use a fractional step method in order to separate the convection step, the pressure equilibrium step and the step involving the remaining source terms . The standard source terms step is solved by a simple explicit first-order Euler method.

Thus, we concentrate only on the description of the pressure equilibrium step, which can be formally written (38)

α k,t = ±P, m k,t = u k,t = 0, (m k e k ) t + p 1 α k,t = 0.
In order to simplify the notations, we denote now by a 0 superscript the physical values in a given cell i at the end of the advection step. These values are computed from the vector W n+1/2 i given by the Rusanov scheme (36). The updated values at time n + 1 are noted without any superscript.

Owing to mass and momentum conservation we have m k = m 0 k and u k = u 0 k . We have now to compute (α 1 , p 1 , p 2 ) in order to pursue the computation. The system is (39)

p 2 -R -p 1 = τ p α 2,t , m 1 e 1 + m 2 e 2 = m 0 1 e 0 1 + m 0 2 e 0 2 , (m 1 e 1 -m 0 1 e 0 1 ) + p 1 (α 1 -α 0 1 ) = 0.
We assume now that the two phases obey stiffened gas equations of state

(40) p k (ρ k , e k ) = (γ k -1)ρ k e k -γ k π k , γ k > 1, k = 1, 2.
It is also physically reasonable to suppose that (41)

π 2 > π 1 ,
because phase (2) is the solid phase. We have also to state some assumptions for the granular stress. We suppose that

(42) R = R(ρ 2 , α 2 ) = ρ γ2 2 α 2 θ(α 2 )
(as already discussed in Section 3). Because of the stiffened gas law, we have

(43) m k e k = α k p k + γ k π k γ k -1 .
Thus we have to solve for (α 1 , p 1 , p 2 ) the following three-equation system, at each time step and in each cell (44)

p 2 -α 2 ρ γ2 2 θ(α 2 ) -p 1 = τ p α 2 (1 -α 2 ) α 2,t , α 2 p 2 + π 2 γ 2 -1 -α 0 2 p 0 2 + π 2 γ 2 -1 + (p 1 + π 2 )(α 2 -α 0 2 ) = 0, α 1 p 1 + π 1 γ 1 -1 -α 0 1 p 0 1 + π 1 γ 1 -1 + (p 1 + π 1 )(α 1 -α 0 1 ) = 0.
We have (45)

p 2 -α 2 ρ γ2 2 θ(α 2 ) -p 1 = τ p α 2 (1 -α 2 ) α 2,t , (α 2 + (γ 2 -1)(α 2 -α 0 2 ))(p 2 + π 2 ) -α 0 2 (p 0 2 + π 2 )- (γ 2 -1)(α 2 ρ γ2 2 θ(α 2 ) + τ p α 2,t )(α 2 -α 0 2 ) = 0, (α 1 + (γ 1 -1)(α 1 -α 0 1 ))(p 1 + π 1 ) -α 0 1 (p 0 1 + π 1 ) = 0. We then note (46) A 1 = α 0 1 (p 0 1 + π 1 ), A 2 = α 0 2 (p 0 2 + π 2 )
. For a stiffened gas law, the sound speed c is given by the formula

(47) c = γ(p + π) ρ .
It implies that the two quantities A 1 and A 2 are > 0. We suppose the following conditions on the granular constraint (48)

θ continuous on [0, 1], θ(α) 0, θ(α) = o(α γ2-2 ) when α → 0.
Let us remark that these conditions are fulfilled by the granular stress that we proposed in (33) where

(49) θ(α) = λα γ2-1 .
Of course, we also suppose that the initial volume fraction 0 < α 0 2 < 1. After the elimination of p 1 and p 2 , the system can be rewritten (50)

G(α 2 ) = (π 2 -π 1 )(α 1 + (γ 1 -1)(α 1 -α 0 1 ))(α 2 + (γ 2 -1)(α 2 -α 0 2 )) +(α 2-γ2 2 m γ2 2 θ(α 2 ) + τ p (1 -α 2 ) α 2,t -A 2 )(α 1 + (γ 1 -1)(α 1 -α 0 1 )) +A 1 (α 2 + (γ 2 -1)(α 2 -α 0 2 )) = 0. We use an implicit first-order time discretisation of α 2,t (51) α 2,t = α 2 -α 0 2 ∆t .
(recall that ∆t is the time step in the convection step (36)). The implicit approach is natural because τ p may be arbitrary small and thus the source term (21) may be stiff.

We first compute G at the left point of the interval [0, 1] (52)

G(0) = -(π 2 -π 1 )(γ 2 -1)α 0 2 (1 + (γ 1 -1)(1 -α 0 1 )) -(A 2 + τ p α 0 2 ∆t )(1 + (γ 1 -1)(1 -α 0 1 )) -A 1 (γ 2 -1)α 0 2 < 0,
because of the hypothesis (48) and A 1 , A 2 > 0. For the computation at the right point, we introduce (53)

β 1 = γ 1 -1 γ 1 α 0 1 .
We have

(54) α 1 + (γ 1 -1)(α 1 -α 0 1 ) > 0 ⇔ α 1 > β 1 ⇔ α 2 < 1 -β 1 . We compute (55) G(1 -β 1 ) = A 1 (1 -β 1 + (γ 2 -1)(α 0 1 -β 1 )) = A 1 γ 1 -1 ((γ 1 -1)(1 -β 1 ) + (γ 2 -1)β 1 ) > 0.
Thus we have existence of a solution α 2 to G(α 2 ) = 0 in the interval ]0, 1 -β 1 [. We have to check that this solution leads to a correct pressures p 1 and p 2 i.e. that p k + π k > 0. But we have (56)

p 1 + π 1 = A 1 α 1 + (γ 1 -1)(α 1 -α 0 1 )
. This quantity is > 0 if the solution satisfies α 2 < 1 -β 1 . Finally we also have (57)

p 2 + π 2 = p 1 + π 1 + α 2 ρ γ2 2 θ(α 2 ) + τ p α 2 -α 0 2 τ α 2 (1 -α 2 ) + π 2 -π 1 .
When α 2 -α 0 2 ≥ 0 the previous quantity is obviously > 0. In the case α 2 -α 0 2 < 0, we use the second equation of (44) and we find (58)

α 2 p 2 + π 2 γ 2 -1 = α 0 2 p 0 2 + π 2 γ 2 -1 -(p 1 + π 1 + π 2 -π 1 )(α 2 -α 0 2 ) > 0.
And the algorithm can continue.

Remark 4. Formula (56) shows that we have to discard any solution that is not in the interval ]0, 1 -β 1 [. Lets us notice that in many cases, we can find another solution to G(α 2 ) = 0 in ]1 -β 1 , 1[. Now, we will give a sufficient condition on the granular stress in such a way that the solution α 2 is unique in the interval ]0, 1 -β 1 [. From now on we set (59)

f (α) = α 2-γ2 θ(α) + τ p ∆tm γ2 2 α -α 0 2 1 -α .
and we assume that (60) f is convex.

We then have

(61) G ′′ (α 2 ) = -2γ 1 γ 2 (π 2 -π 1 )-2γ 1 m γ2 2 f ′ (α 2 )+(α 1 +(γ 1 -1)(α 1 -α 0 1 ))m γ2 2 f ′′ (α 2 ).
The function f is convex and satisfies f (α) > f (0) for 0 < α < 1. It implies that f is also increasing. Then the two first terms in (61) are < 0 and the last is > 0. On the other hand, it is sufficient that G is concave to prove the uniqueness. But (62)

G ′′ (α 2 ) -2γ 1 γ 2 (π 2 -π 1 )- 2γ 1 m γ2 2 f ′ (α 2 ) + γ 1 m γ2 2 (1 -α 2 )f ′′ (α 2 ) 2γ 1 m γ2 2 1 2 (1 -α 2 )f ′′ (α 2 ) -f ′ (α 2 ) - γ 2 m γ2 2 (π 2 -π 1 )
A sufficient condition to obtain uniqueness is then

(63) 1 2 (1 -α 2 )f ′′ (α 2 ) -f ′ (α 2 ) - γ 2 m γ2 2 (π 2 -π 1 ) 0.
Example 1. As suggested in (33), we take

(64) θ(α 2 ) = λα γ2-1 2 , λ > 0.
Thus we have

(65) f (α 2 ) = λα 2 + τ p ∆tm γ2 2 α 2 -α 0 2 1 -α 2 ,
wich is indeed convex. The above inequality (63) becomes

(66) -λ - γ 2 m γ2 2 (π 2 -π 1 ) 0
and it is obviously satisfied, independantly of τ p .

We sum up the previous computations in the following proposition, which is useful for the implementation of the algorithm. Proposition 2. Let the granular stress be defined by ( 64) and (32). Let (67) 0 < α 0 1 < 1, p 0 k + π k > 0, k = 1, 2, then, the algebraic system (39), (51) admits a unique solution (α 1 , p 1 , p 2 ) that complies with

(68) 0 < α 1 < 1, p k + π k > 0, k = 1, 2.
Moreover we also have

(69) α 1 > γ 1 -1 γ 1 α 0 1 .
Finally, the solution can be computed by the Newton's method by solving equation (50) for α 2 . A safe choice for the initialisation of the Newton's method is α 2 = 0 (since G is concave).

Numerical results

6.1. Academical test cases. First, we consider two one-dimensional test cases in order to evaluate the influence of the granular stress on the system stability. We take τ p = 0, which corresponds to instantaneous pressure equilibrium. It is known that generally, the equilibrium system is not hyperbolic (the computations are recalled in Section 7). The numerical parameters are taken from [START_REF] Coquel | A numerical method using upwind schemes for the resolution of two-phase flows[END_REF] (and also studied in [START_REF] Hurisse | Couplage interfacial instationnaire de modèles diphasiques[END_REF]). We consider a simple Riemann problem in the interval [-1/2, 1/2]. The two phases are supposed to satisfy perfect gas equations of state, with γ 1 = 1.0924 and γ 2 = 1.0182. The initial condition is made of two constant states jumping at x = 0.

We plot the solution at time t = 0.0008. The CFL number is fixed to 0.9. The initial data are We perform our algorithm with a granular stress R = 0. With a mesh of 1, 000 cells we observe that the solution is rather smooth but a small oscillation starts to develop in the center of the computational domain. The volume fraction α 1 , the velocities and pressures are plotted on Figures 1, 2 and 3.

(70) (L) (R
The same computation is made with 10, 000 cells. We observe on Figures 4, 5 and 6 that instabilities arise, due to the non-hyperbolic behavior of the model.

We have also performed a computation on a 100, 000 cells mesh. The oscillations clearly increase as can be seen on Figure 7 for the volume fraction.

We perform then another computation on the finer mesh with a granular stress given by (33). For the numerics, we have chosen λ = 500. We observe a very slight damping of the oscillations on Figure 8 (to be carefully compared with Figure 4). We have also plotted the two pressures on Figure 10 in order to show the difference of pressures imposed by the granular stress.

It clearly arises that the magnitude of the granular stress is not sufficient to recover an hyperbolic regime. This is confirmed on Figure 11 where we compare We observe on Figure 11 that the imaginary part slightly decreases, owing to the introduction of the granular stress; however, it does not vanish. We now evaluate the influence of the positive time scale parameter τ p , and on its ability to stabilize the model. While setting r ref = 10 7 , we compare in Figure 12 the pressures obtained with different relaxation coefficients (corresponding with true time scales τp p ref = 0, 10 -7 , 10 -6 , 10 -5 and 10 -4 ). On this rather coarse mesh of 1,000 cells, we observe that the pressures are very similar for small enough relaxation time scales. When focusing on a finer mesh of 50,000 cells, the stability of the approximations increases when the relaxation time scales are larger (Figure 13). We also observe on this test cases that even with high values of τ p , the difference between both pressures p 1 and p 2 is indeed rather small. 6.2. Simplified combustion chamber. We consider now a more realistic case taken from [START_REF] Nussbaum | Modélisation et simulation numérique d'un écoulement diphasique de la balistique intérieure[END_REF]. We are interested in the modeling of a simplified gun. The gun is modeled by a one-dimensional tube filled with a solid phase (the powder grains) The momentum source term is given by ( 73)

Q = M u 2 -D D = Cα 1 α 2 ρ 2 (u 1 -u 2 ) |u 1 -u 2 | (drag force) C = 3 4r
(simplified shape factor) The energy source terms are ( 74)

S 1 = -u 2 D + M Q ex S 2 = u 2 D
Q ex = 37.3839 × 10 6 J/kg (chemical combustion energy)

Let us remark that they do not cancel when summed up. This is only apparently a violation of the total energy conservation. Actually, we can rewrite the model in order to have opposite source terms. The rewriting is based on a translation of the internal energy in the pressure laws. , where e ′ k is the translated internal energy of phase k, and e 0 k is a reference energy for phase k, and if we define the translated total energies

(76) E ′ k = e ′ k + e 0 k + u 2 k 2 .
the energy balance equations can be rewritten

(77) (m k E ′ k ) t + ((m k E ′ k + α k p k )u k ) x + p I α k,t = ±S ∓ M e 0 k . Now the total translated energy m 1 E ′ 1 + m 2 E ′
2 is no longer conserved since the ± terms do not cancel. The term M (e 0 2 -e 0 1 ) can be identified to the chemical reaction heat.

The other parameters of the computations are (78) When the gas pressure p 1 at the shot base is greater than the resistive pressure p r , the acceleration of the projectile is given by

γ 1 = 1.4 γ 2 = 3 π 2 = 2.
m p dv dt = (p 1 -p r ) π diam 2 4 .
The algorithm to move the right boundary is based on an Arbitrary Lagrangian Eulerian (ALE) approach described in [START_REF] Nussbaum | Modélisation et simulation numérique d'un écoulement diphasique de la balistique intérieure[END_REF]. As the domain enlarges, the number of computational cells increases.

We compare our new compressible model with the classical Gough model described for example in [START_REF] Nussbaum | Numerical simulations of gas-particle flows with combustion[END_REF].

We obtain the following results for the projectile velocity at the exit time On Figure 14, we compare the pressure evolution at the breech and the shot base of the projectile. We observe a good qualitative agreement between the Gough model and the relaxation model.

Finally, we plot some quantities in the tube at the final time. The porosity, the velocities and the pressures are given in Figures 15,16, 17. We also plot on Figure 18 the density ρ 2 of the solid phase at the final time in order to check that the variations of the powder density are small when compared with the initial density ρ 2 = 1600 kg/m 3 .

Conclusion

In this paper, we have adapted the pressure relaxation method described in [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] and [START_REF] Hurisse | Couplage interfacial instationnaire de modèles diphasiques[END_REF] to the case of a non-vanishing granular stress.

Starting from the two-velocity, two-pressure multiphase model of Baer-Nunziato, we have proposed a relaxation source term in the governing equation of the void fraction that is compatible with the second principle of thermodynamics. In this study, we have shown that

• the source term increases the entropy of the phase mixture;

• the granular stress cannot have an arbitrary form. It is related to the fact that the differential form satisfied by the entropy is closed. When the relaxation time tends to zero, we have then proposed a numerical method based on the underlying two-pressure model to approximate the singlepressure model. In the pressure relaxation step, the void fraction is updated in order to equilibrate the jump of pressures with the granular stress. We have proved existence and uniqueness of the equilibrium void fraction under some hypothesis on the granular stress. Those hypothesis are satisfied by physically reasonable models.

Eventually, we have proposed some numerical experiments in order to validate our approach. In an ideal test case, we have checked that when the mesh is refined, the instability of the one-pressure model is (fortunately) not suppressed. We also checked that the introduction of the granular stress slightly improves the whole stability. We finally performed more realistic simulations where we were able to reproduce correct quantitative features of a simplified gun.

The whole approach is thus very promising and must now be extended to more sophisticated granular pressure laws, equations of state and geometries.

7. Appendix I: hyperbolicity 7.1. Relaxed system. For the sake of completness, we recall the proof of hyperbolicity of the convection part of the equations. It is convenient to study it in the variables

(79) Y = (α 1 , ρ 1 , u 1 , s 1 , ρ 2 , u 2 , s 2 ) T .
In this set of variables the system becomes (80)

Y t + B(Y )Y x = 0, with (81)
c k = ∂p(ρ k , s k ) ∂ρ k , k = 1, 2 B(Y ) =            u 2 ρ1(u1-u2) α1 u 1 ρ 1 c 2 1 ρ1 u 1 p1,s 1 ρ1 u 1 u 2 ρ 2 p1-p2 m2 c 2 2 ρ2 u 2 p2,s 2 ρ2 u 2           
The characteristic polynomial is

(82) P (λ) = (u 2 -λ) 2 (u 1 -λ)(u 1 -c 1 -λ)(u 1 + c 1 -λ)(u 2 -c 2 -λ)(u 2 + c 2 -λ)
We can then state the following proposition

Proposition 3. If |u 1 -u 2 | = c k , k = 1, 2 then, the system (2) is hyperbolic. If |u 1 -u 2 | = c k for k = 1
or 2 then the system is resonant.

7.2. Equilibrium system. We also study the hyperbolicity behavior of the equilibrium system, which correspond to the limit τ p = 0. When the granular stress vanishes, the results are given in several papers. We thus only detail the case R > 0 with a granular stress satisfying (32). The computations given below have been used to draw Figure 11. We note, for any quantity z,

D k z = z t + u k z x . (83)
At equilibrium, we can remove the transport equations in α k and replace them by the pressure relation (84)

p 2 = p 1 + α 2 ρ γ2 2 θ(α 2 ) = p 1 + ρ γ2 2 g(α 2 ). We note (85) h = g -1 .
Of course, when the granular stress vanishes, the inverse of g is not defined and the computations must be carried out in another way. At equilibrium, we have (86)

α 2 = h p 2 -p 1 ρ γ2 2 ⇒ dα 2 = δ c 2 2 -γ 2 p 2 -p 1 ρ 2 dρ 2 + p 2,s2 ds 2 -c 2 1 dρ 1 -p 1,s1 ds 1 with δ = h ′ p2-p1 ρ γ 2 2 ρ γ2 2 > 0.
Example 2. We can consider

(87) θ(α) = λα γ2-1 .
We then have

(88) δ = α 1-1/γ2 2 λγ 2 ρ γ2 2 . It is natural to introduce (89) a 2 2 = γ 2 ρ 1 γ 1 ρ 2 c 2 1 + γ 2 π 2 -π 1 ρ 2 > 0,
In such a way that we have also

(90) dα 2 = δ a 2 2 dρ 2 + p 2,s2 ds 2 -c 2 1 dρ 1 -p 1,
P = -δ a 2 2 D 2 ρ 2 + p 2,s2 D 2 s 2 -c 2 1 D 2 ρ 1 -p 1,s1 D 2 s 1
We then rewrite the equilibrium system in the variables

(92) Z = (ρ 1 , u 1 , s 1 , ρ 2 , u 2 , s 2 ) T ,
In this variables, the system is

(93) Z t + C(Z)Z x = 0.
For the sake of completeness, we give some details of the computations (94)

ρ 1,t + u 1 ρ 1,x - ρ 1 α 1 (u 1 -u 2 )δ a 2 2 ρ 2,x + p 2,s2 s 2,x -c 2 1 ρ 1,x -p 1,s1 s 1,x + ρ 1 u 1,x - ρ 1 α 1 δ a 2 2 D 2 ρ 2 + p 2,s2 D 2 s 2 -c 2 1 D 2 ρ 1 -p 1,s1 D 2 s 1 = 0, ρ 2,t + u 2 ρ 2,x + ρ 2 α 2 δ a 2 2 D 2 ρ 2 + p 2,s2 D 2 s 2 -c 2 1 D 2 ρ 1 -p 1,s1 D 2 s 1 + ρ 2 u 2,x = 0 (95) ρ 1,t + u 1 ρ 1,x - ρ 1 α 1 u 1 δ a 2 2 ρ 2,x + p 2,s2 s 2,x -c 2 1 ρ 1,x -p 1,s1 s 1,x + ρ 1 u 1,x - ρ 1 α 1 δ a 2 2 ρ 2,t + p 2,s2 s 2,t -c 2 1 ρ 1,t -p 1,s1 s 1,t = 0, (96)
u 1,t + u 1 u 1,x + 1 ρ 1 p 1,x = 0, u 2,t + u 2 u 2,x + 1 ρ 2 p 2,x + p 2 -p 1 m 2 δ a 2 2 ρ 2,x + p 2,s2 s 2,x -c 2 1 ρ 1,x -p 1,s1 s 1,x = 0 (97) s k,t + u k s k,x = 0 (98) (1 + ρ 1 c 2 1 δ α 1 )ρ 1,t - ρ 1 a 2 2 δ α 1 ρ 2,t + u 1 ρ 1,x + ρ 1 α 1 δ -a 2 2 u 1 ρ 2,x + c 2 1 u 1 ρ 1,x + (u 2 -u 1 )p 2,s2 s 2,x + ρ 1 u 1,x = 0 (1 + ρ 2 a 2 2 δ α 2 )ρ 2,t - ρ 2 c 2 1 δ α 2 ρ 1,t + u 2 ρ 2,x + ρ 2 α 2 δ a 2 2 u 2 ρ 2,x -c 2 1 u 2 ρ 1,x -p 1,s1 (u 2 -u 1 )s 1,x + ρ 2 u 2,x = 0 (99) (1 + ρ 1 c 2 1 δ α 1
)ρ It is not easy to compute the eigenvalues analytically. It is also difficult to give a practical sufficient condition on all the parameters in order to prove that the eigenvalues are all real. In the case δ = 0, corresponding to an infinite granular stress, the characteristic polynomial is

(102) P (λ) = (u 2 -λ)(u 1 -λ)(u 1 -c 1 -λ)(u 1 + c 1 -λ)(u 2 -c 2 -λ)(u 2 + c 2 -λ)
We recover the same eigenvalues as in (82). With a small λ, which corresponds to a big δ as can be seen by formula (88), we observe numerically that the system is elliptic when u 1 = u 2 . When λ increases, δ decreases and we recover an hyperbolic behavior.

Numerical application: we take We modify λ to λ = 500, the pressure p 2 is now p 2 = 0.20352544 × 10 It is also possible to compute an entropy associated to the choice (33). For this, we postulate the following form of the entropy of the solid phase (as in Section 3, we omit the subscript) (106) s = K(α)U ((e -πτ )τ γ-1 ).

This choice is justified by the fact that when K is constant, then we recover the general entropy of a stiffened gas. We can verify that our entropy and our stiffened gas law are compatible. Without the subscripts, the equation ( 13) reads (107)

T ds = de + p dτ -Θ dα.

The temperature is given by ( 108)

s e = 1 T = K (α) τ γ-1 U ′ (e -πτ ) τ γ-1
In a similar way, we can deduce a relationship between p and T U (e -πτ ) τ γ-1 U ′ ((e -πτ ) τ γ-1 ) ρ γ-1

But Θ has also to be of the form (32). It implies that We choose now the sign of the constants A and B in such way that the temperature is positive and that the function U is concave. It implies that A and B are > 0. The positivity of T and the concavity of U are required for obtaining real sound speeds in the pure phases (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]).

Figure 1 .

 1 Figure 1. Void fraction, 1, 000 cells, τ p = 0, R = 0.

Figure 2 .

 2 Figure 2. Velocities u 1 and u 2 , 1,000 cells, τ p = 0, R = 0.

Figure 3 . 6 i=1

 36 Figure 3. Pressures, 1,000 cells, τ p = 0, R = 0.

Figure 4 .

 4 Figure 4. Void fraction, 10,000 cells, τ p = 0, R = 0.

Figure 5 .

 5 Figure 5. Velocities u 1 and u 2 , 10,000 cells, τ p = 0, R = 0.

Figure 6 .

 6 Figure 6. Pressures, 10,000 cells, τ p = 0, R = 0.

Figure 7 .

 7 Figure 7. Void fraction, 100,000 cells, τ p = 0, R = 0 (the two colors only correspond to the two processors used in the MPI computation).

Figure 8 .

 8 Figure 8. Void fraction, 10,000 cells, τ p = 0, R = 500m γ2 2 .

Figure 9 .

 9 Figure 9. Velocities u 1 and u 2 , 10000 cells, with granular stress .

Figure 10 .

 10 Figure 10. Pressures, 10,000 cells, τ p = 0, R = 500m γ2 2 .

Figure 11 .

 11 Figure 11. Imaginary parts, 1,000 cells, τ p = 0, R = 500m γ2 2 or R = 0. stress.

Figure 12 .

 12 Figure 12. Pressure p 1 , 1,000 cells, with 0 ≤ τ p ≤ 1000.

Figure 13 .Remark 5 .

 135 Figure 13. Pressure, 50, 000 cells, τ p = 0, 2, 10, R = 0.

  1333 × 10 9 P a ρ 2 = 1600kg/m 3 (initial solid density) m p = 30 kg (projectile mass) p r = 10 8 Pa (resistive pressure) p 0 = 10 5 Pa (initial pressure) T 0 = 294 K (initial temperature) ρ 0 = 0.8713 kg/m 3 (initial gas density) α 2,0 = 0.5709 (initial porosity) diam = 132mm (diameter of the gun) length = 762mm (length of the tube) mpow = 9.5255kg (powder mass) λ = 0.03 (granular parameter) mmol = 21.3g/mol (molar mass of the gas)

Figure 14 .

 14 Figure 14. Pressure evolution at the breech and the shot base wrt time. Comparison between the Gough and the relaxation model.

Figure 15 .

 15 Figure 15. Porosity at the final time. Relaxation model with granular stress.

Figure 16 .

 16 Figure 16. Velocities at the final time. Relaxation model with granular stress.

Figure 17 .

 17 Figure 17. Pressures at the final time. Relaxation model with granular stress.

Figure 18 .

 18 Figure 18. Density of the solid phase at the final time. Relaxation model with granular stress.

π 1 = 8 u 2 = -u 1

 1821 π 2 = 0 α 1 = 0.25 p 1 = 0.2 × 10 8 λ = 0.01 ⇒ p 2 = 0.20000007050881 × 10

  ) (γ -1) e τ γ-2 -γ π τ γ-1 U ′ (e -πτ ) τ γ-1 K (α) τ γ-1 U ′ ((e -πτ ) τ γ-1 ) = (γ -1) e τ -γ πand we indeed recover the stiffened gas equation of state.We try now to find an expression for the function U (x). From (107) we can write (111) s α = -

U

  (x) = A exp(-Bx).

  s1 ds 1 It gives another expression of the source term P at equilibrium

	(91)

  -a 2 2 u 1 ρ 2,x + p 2,s2 (u 2 -u 1 )s 2,x + ρ 1 u 1,x = 0 u 2 ρ 1,x -p 1,s1 (u 2 -u 1 )s 1,x + ρ 2 u 2,x = 0

	1,t -)ρ 2,t -δ -c 2 δ (1 + + ρ 1 α 1 ρ 2 a 2 2 δ α 2 + ρ 2 α 2 u 1,t + u 1 u 1,x + c 2 1 ρ 1 ρ 1,x + ρ 1 a 2 2 δ α 1 ρ 2 c 2 1 δ α 2 p 1,s1 ρ 2,t + (1 + ρ 1,t + (1 + ρ 1 s 1,x = 0, 1 (100) u 2,t + u 2 u 2,x + 1 ρ 2 (c 2 2 + p 2 -p 1 α 2 δa 2 2 )ρ 2,x + ρ 1 c 2 1 δ α 1 ρ 2 a 2 2 δ α 2 p 2,s2 )u 1 ρ 1,x )u 2 ρ 2,x ρ 2 (1 +	p 2 -p 1 α 2	δ)s 2,x
	+	p 2 -p 1 m 2	δ -c 2 1 ρ 1,x -p 1,s1 s 1,x = 0
	Finally, setting
	(101)		∆ = α 1 α 2 + δ(α 1 ρ 2 a 2 2 + α 2 ρ 1 c 2 1 ),
	we find


Valerio Lucarini
email: valerio.lucarini@zmaw.de

Davide Faranda

Jeroen Wouters

Universal behavior of extreme value statistics for selected observables of dynamical systems

The main results of the extreme value theory developed for the investigation of the observables of dynamical systems rely, up to now, on the Gnedenko approach. In this framework, extremes are basically identified with the block maxima of the time series of the chosen observable, in the limit of infinitely long blocks. It has been proved that, assuming suitable mixing conditions for the underlying dynamical systems, the extremes of a specific class of observables are distributed according to the so called Generalized Extreme Value (GEV) distribution. Direct calculations show that in the case of quasi-periodic dynamics the block maxima are not distributed according to the GEV distribution. In this paper we show that, in order to obtain a universal behaviour of the extremes, the requirement of a mixing dynamics can be relaxed if the Pareto approach is used, based upon considering the exceedances over a given threshold. Requiring that the invariant measure locally scales with a well defined exponent -the local dimension -, we show that the limiting distribution for the exceedances of the observables previously studied with the Gnedenko approach is a Generalized Pareto distribution where the parameters depends only on the local dimensions and the value of the threshold. This result allows to extend the extreme value theory for dynamical systems to the case of regular motions. We also provide connections with the results obtained with the Gnedenko approach. In order to provide further support to our findings, we present the results of numerical experiments carried out considering the well-known Chirikov standard map.

I. INTRODUCTION

Extreme value Theory was originally introduced by Fisher and Tippett [START_REF] Fisher | Proceedings of the Cambridge philosophical society[END_REF] and formalised by Gnedenko [START_REF] Gnedenko | [END_REF], who showed that the distribution of the maxima of a sample of independent identically distributed (i.i.d) stochastic variables converges under very general conditions to a member of the so-called Generalised Extreme Value (GEV) distribution. The attention of the scientific community to the problem of understanding extreme values theory is growing, also because this theory is crucial in a wide class of applications for defining risk factors such as those related to instabilities in the financial markets and to natural hazards related to seismic, climatic and hydrological extreme events. Even if the probability of extreme events decreases with their magnitude, the damage that they may bring increases rapidly with the magnitude as does the cost of protection against them. From a theoretical point of view, extreme values of observables are related to large fluctuations of the corresponding underlying system. An extensive account of recent results and relevant applications is given in [3].

The traditional ('Gnedenko') approach for the statistical inference of extremes is related to the original results by Gnedenko [START_REF] Gnedenko | [END_REF]: we partition the experimental time series into bins of fixed length, we extract the maximum of each bin, and fit the selected data to the GEV distribution family using methods such as maximum likelihood estimation (MLE) or L-moments. See [4] for a detailed account of this methodology. The selection of just one maximum in a fixed period may lead to the loss of relevant information on the large fluctuations of the system, especially when there are many large values in a given period [5]. This problem can be taken care of by considering several of the largest order statistics instead of just the largest one. For such maxima distributions we expect convergence to the Generalized Pareto Distribution (GPD) introduced by Pickands III [6] and Balkema and De Haan [7] to model the exceedances over a given threshold. We call this the 'Pareto' approach. Also this approach has been widely adopted for studying empirically natural extreme phenomena such as those related to waves, winds, temperatures, earthquakes and floods [8][START_REF] Bayliss | Peaks-over-threshold flood database[END_REF][START_REF] Pisarenko | [END_REF].

Both the Gnedenko and Pareto approaches were originally designed to study extreme values for series of i.i.d variables. In this case it is well known that a strong connections exists between the two methodologies, as we have that if block maxima obey the GEV distribution, then exceedances over some high threshold will have an associated GPD. Moreover, the shape parameter of the GPD and that of the corresponding GEV distribution are identical [START_REF] Leadbetter | Extremes and related properties of random sequences and processes[END_REF]. As a result, several practical methods (e.g. Hill's and Pickands' estimators) developed for estimating the shape parameter of the GEV distribution of the extremes of a given time series are actually based upon comparing the GPD fits at various thresholds [START_REF] Embrechts | Modelling Extremal Events[END_REF][START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]. In practical terms, it appears that, while the Gnedenko and Pareto approaches provide equivalent information in the asymptotic limit of infinitely long time series, the GPD statistics is more robust when realistic, finite time series are considered (see, e.g., [START_REF] Ding | [END_REF]).

In recent years, especially under the influence of the rapid development of numerical modelling in the geophysical sciences and of its applications for the investigation of the socio-economic impacts of extreme events, it has become of great relevance to understand whether it is possible to apply the extreme value theory on the time series of observables of deterministic dynamical systems. Carefully devised numerical experiments on climate models of various degrees of complexity have shown that the speed of convergence (if any) of the statistical properties of the extremes definitely depends on the chosen climatic variable of interest [4,[15][16][17].

Several papers have addressed this issue at a more general level. A first important result is that when a dynamical system has a regular (periodic of quasi-periodic) behaviour, we do not expect, in general, to find convergence to GEV distributions for the extremes of any observable. These results have been presented by Balakrishnan et al. [18], and more recently, by Nicolis et al. [19] and by Haiman [20].

A different mathematical approach to extreme value theory in dynamical systems was proposed in the landmark paper by Collet [21], which has paved the way for the recent results obtained in the last few years [22][START_REF] Freitas | Probability Theory and Related Fields[END_REF][START_REF] Freitas | [END_REF][START_REF] Gupta | Lozi-like maps, and Lorenz-like maps[END_REF]. The starting point of all of these investigations has been to associate to the stationary stochastic process given by the dynamical system, a new stationary independent sequence which obeys one of the classical three extreme value laws introduced by Gnedenko [START_REF] Gnedenko | [END_REF]. The assumptions which are necessary to observe a GEV distribution in dynamical systems rely on the choice of suitable observables (specific functions of the distance between the orbit and the initial condition, chosen to be on the attractor) and the fulfillment of particular mixing conditions that guarantee the independence of subsequent maxima. Recent studies have shown that the resulting parameters of the GEV distributions can be expressed as simple functions of the local (around the initial condition) dimension of the attractor, and detailed numerical investigations have clarified the conditions under which convergence to the theoretical GEV distributions can be satisfactorily achieved when considering finite time series [START_REF] Faranda | [END_REF][START_REF] Faranda | [END_REF][28].

In this paper, we wish to attempt a unification of these two lines of work by using the Pareto rather than the Gnedenko approach. We choose the same class of observables presented in [22][START_REF] Freitas | Probability Theory and Related Fields[END_REF][START_REF] Freitas | [END_REF][START_REF] Gupta | Lozi-like maps, and Lorenz-like maps[END_REF][START_REF] Faranda | [END_REF][START_REF] Faranda | [END_REF][28] and show that, assuming only that the local measures scales with the local dimension [START_REF] Bandt | Physics and theoretical computer science: from numbers and languages to (quantum) cryptography security[END_REF], it is possible to obtain by direct integration a GPD for the threshold exceedences when considering a generic orbit of a dynamical systems, without requiring any special mixing properties. The parameters will depend only on the choice of the threshold and, more importantly, on the local dimension. Note that Castillo and Hadi [5] had already pointed out that in the case of periodic or quasi-periodic motion the Gnedenko approach to the evaluation of the extreme value statistics is inefficient, basically because in the limit of very large blocks, we tend to observe always the same maximum in all bins. To support our analytical results we provide numerical experiments that we carry out considering the classic Chirikov standard map [START_REF] Chirikov | [END_REF]. This paper is organised as follows. In section 2 we recapitulate the extreme value theory for dynamical systems obtained using the Gnedenko approach. In section 3 we present our general results obtained using the Pareto approach. In section 4 we provide support to our investigation by examining the results of the numerical simulations performed on the standard map. In section 5 we present our final remarks and future scientific perspectives.

II. GNEDENKO APPROACH: GENERALIZED EXTREME VALUE DISTRIBUTIONS IN DYNAMICAL SYSTEMS

Gnedenko [START_REF] Gnedenko | [END_REF] studied the convergence of maxima of i.i.d. variables

X 0 , X 1 , ..X m-1 with cumulative distribution function (cdf) F (x) = P {a m (M m -b m ) ≤ x}
where a m and b m are normalizing sequences and M m = max{X 0 , X 1 , ..., X m-1 }. Under general hypothesis on the nature of the parent distribution of data, Gnedenko [START_REF] Gnedenko | [END_REF] showed that the asymptotic distribution of maxima, up to an affine change of variable, belongs to a single family of generalized distribution called GEV distribution whose cdf can be written as:

F GEV (x; µ, α, κ) = e -t(x) (1)
where

t(x) = 1 + κ( x-µ α ) -1/κ if κ = 0 e -(x-µ)/α if κ = 0 . (2)
This expression holds for 1 + κ(x-µ)/α > 0, using µ ∈ R (location parameter) and α > 0 (scale parameter) as scaling constants in place of b m , and a m [31], in particular, in Faranda et al. [START_REF] Faranda | [END_REF] we have shown that µ = b m and α = 1/a m , where κ ∈ R is the shape parameter (also called the tail index). When κ → 0, the distribution corresponds to a Gumbel type (Type 1 distribution). When the index is positive, it corresponds to a Fréchet (Type 2 distribution); when the index is negative, it corresponds to a Weibull (Type 3 distribution).

Let us consider a dynamical system (Ω, B, ν, f ), where Ω is the invariant set in some manifold, usually R d , B is the Borel σ-algebra, f : Ω → Ω is a measurable map and ν an f -invariant Borel measure. In order to adapt the extreme value theory to dynamical systems, following [22][START_REF] Freitas | Probability Theory and Related Fields[END_REF][START_REF] Freitas | [END_REF][START_REF] Gupta | Lozi-like maps, and Lorenz-like maps[END_REF], we consider the stationary stochastic process X 0 , X 1 , ... given by:

X m (x) = g(dist(f m (x), ζ)) ∀m ∈ N (3)
where 'dist' is a distance in the ambient space Ω, ζ is a given point and g is an observable function. The partial maximum in the Gnedenko approach is defined as:

M m = max{X 0 , ..., X m-1 }. (4)
Defining r = dist(x, ζ), we consider the three classes of observables g i , i = 1, 2, 3:

g 1 (r) = -log(r) (5) g 2 (r) = r -β (6) g 3 (r) = C -r β ( 7
)
where C is a constant and β > 0 ∈ R. Using the observable g i we obtain convergence of the statistics of the block maxima of their time series obtained by evolving the dynamical system to the Type i distribution if one can prove two sufficient conditions called D 2 and D ′ , which basically imply a sort of independence of the series of extremes resulting from the mixing of the underlying dynamics [22]. The conditions cannot be simply related to the usual concepts of strong or weak mixing, but are indeed not obeyed by observables of systems featuring a regular dynamics or power-law decay of correlations.

A connection also exists between the existence of extreme value laws and the statistics of first return and hitting times, which provide information on how fast the point starting from the initial condition ζ comes back to a neighborhood of ζ, as shown by Freitas et al. [START_REF] Freitas | Probability Theory and Related Fields[END_REF] and Freitas et al. [32]. In particular, they proved that for dynamical systems possessing an invariant measure ν, the existence of an exponential hitting time statistics on balls around ν-almost any point ζ implies the existence of extreme value laws for one of the observables of type g i , i = 1, 2, 3 described above. The converse is also true, namely if we have an extreme value law which applies to the observables of type g i , i = 1, 2, 3 achieving a maximum at ζ, then we have exponential hitting time statistics to balls with center ζ. Recently these results have been generalized to local returns around balls centered at periodic points [START_REF] Freitas | [END_REF].

In Faranda et al. [START_REF] Faranda | [END_REF][START_REF] Faranda | [END_REF]28] we analised both from an analytical and numerical point of view the extreme value distribution in a wide class of low dimensional maps. We divided the time series of length k of the g i observables into n bins each containing the same number m of observations, and selected the maximum (or the minimum) value in each of them [33]. We showed that at leading order (the formulas are asymptotically correct for m, k → ∞), the GEV parameters in mixing maps can be written in terms of m (or equivalently n) and the local dimension of the attractor D. We have:

• g 1 -type observable:

α = 1 D µ ∼ 1 D ln(k/n) κ = 0 (8)
• g 2 -type observable:

α ∼ n -β D µ ∼ n -β D κ = β D (9)
• g 3 -type observable:

α ∼ n β D µ = C κ = - β D (10)
Moreover, we clearly showed that other kind of distributions not belonging to the GEV family are observed for quasi-periodic and periodic motions.

III. THE PARETO APPROACH: GENERALIZED PARETO DISTRIBUTIONS IN DYNAMICAL SYSTEMS

We define an exceedance as z = X -T , which measures by how much X exceeds the threshold T . As discussed above, under the same conditions under which the block maxima of the i.i.d. stochastic variables X obey the GEV statistics, the exceedances z are asymptotically distributed according to the Generalised Pareto Distribution [START_REF] Leadbetter | Extremes and related properties of random sequences and processes[END_REF]:

F GP D (z; ξ, σ) =    1 -1 + ξz σ -1/ξ for ξ = 0, 1 -exp -z σ for ξ = 0, ( 11
)
where the range of z is 0 ≤ z < ∞ if ξ ≤ 0 and 0 ≤ z ≤ σ/ξ if ξ > 0. We consider the same set up described in the previous section and take into account an observables g = g(dist(x, ζ)) = g(r), such that g achieves a maximum g max for r = 0 (finite or infinite) and is monotonically decreasing. We study the exceedance above a threshold T defined as T = g(r * ). We obtain an exceedence every time the distance between the orbit of the dynamical system and ζ is smaller than r * . Therefore, we define the exceedances z = g(r) -T . By the Bayes' theorem, we have that P (r < g -1 (z + T )|r < g -1 (T )) = P (r < g -1 (z + T ))/P (r < g -1 (T )). In terms of invariant measure of the system, we have that the probability H g,T (z) of observing an exceedance of at least z given that an exceedence occurs is given by:

H g,T (z) ≡ ν(B g -1 (z+T ) (ζ)) ν(B g -1 (T ) (ζ)) . (12)
Obviously, the value of the previous expression is 1 if z = 0. In agreement with the conditions given on g, the expression contained in Eq. ( 12) monotonically decreases with z and vanishes when the radius is given by g -1 (g max ). Note that the corresponding cdf is given by F g,T (z) = 1-H g,T (z). In order to address the problem of extremes, we have to consider small radii. At this regard we will invoke, and assume, the existence of the following limit

lim r→0 log ν(B r (ζ)) log r = D(ζ), for ζ chosen ν -a.e., (13)
where D(ζ) is the local dimension of the attractor [START_REF] Bandt | Physics and theoretical computer science: from numbers and languages to (quantum) cryptography security[END_REF]. Therefore, we rewrite the following expression for the tail probability of exceedance:

H g,T (z) ∼ g -1 (z + T ) g -1 (T ) D . ( 14
)
where we have dropped the ζ dependence of D to simplify the notation. By substituting g with specific observable we are considering, we obtain explicitly the corresponding extreme value distribution law.

By choosing an observable of the form given by either g 1 , g 2 , or g 3 , we derive as extreme value distribution law one member of the Generalised Pareto Distribution family given in Eq. [START_REF] Leadbetter | Extremes and related properties of random sequences and processes[END_REF]. Results are detailed below:

• g 1 -type observable:

σ = 1 D ξ = 0; (15)
• g 2 -type observable:

σ = T β D ξ = β D ; (16)
• g 3 -type observable:

σ = (C -T )β D ξ = - β D . (17)
The previous expressions show that there is a simple algebraic link between the parameters of the GPD and the local dimension of the attractor around the point ζ. This implies that the statistics of extremes provides us with a new algorithmic tool for estimating the local fine structure of the attractor. These results show that it is possible to derive general properties for the extreme values of the observables g 1 , g 2 , or g 3 independently on the qualitative properties of the underlying dynamics, be the system periodic, quasi-periodic, or chaotic. Therefore, by taking the Pareto instead of the Gnedenko approach, we are able to overcome the mixing conditions (or the requirements on the properties of the hitting time statistics) needed to derive a general extreme value theory for dynamical systems, as proposed in [21][22][START_REF] Freitas | Probability Theory and Related Fields[END_REF][START_REF] Freitas | [END_REF][START_REF] Gupta | Lozi-like maps, and Lorenz-like maps[END_REF]. In [START_REF] Faranda | [END_REF][START_REF] Faranda | [END_REF][28] we had proposed that the reason why a link between the extreme value theory and the local properties of the invariant measure in the vicinity of the point ζ can be explained by the fact that selecting the extremes of the observables g 1 , g 2 , or g 3 amounts to performing a zoom around ζ. In the case of the Gnedenko approach, such a picture is accurate only if the dynamics is mixing (time and spatial selection criteria are equivalent). Instead, in the case of the Pareto approach, this is literally what we are doing when writing Eq. ( 14), as we are remapping the radius of the ball in a monotonic fashion though the inverse of the g-functions.

A. Relationship between the Gnedenko and Pareto approaches

The relation between GEV and GPD parameters have been already discussed in literature in case of i.i.d variables [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF][START_REF] Ding | [END_REF]34,35]. Coles [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF] and Katz et al. [34] have proven that the cdf of the GEV defined as F GEV (z; µ, α, κ) can be asymptotically written as that of GPD under a high enough threshold as follows: (18) where κ = ξ, σ = α + ξ(T -µ), and T = µ + σ ξ (λ -ξ -1), with ln(α) = ln(σ) + ξ ln(λ). In the present case, we have to compare Eqs. ( 8)- [START_REF] Pisarenko | [END_REF] for GEV with Eqs. ( 15)-( 17) for GPD, keeping in mind that the GEV results hold only under the mixing conditions discussed before. While it is immediate to check that κ = ξ, the other relationships are valid in the limit of large n, as expected.

F GEV (z; µ, α, κ) ∼ F GP D (z; T, σ, ξ) = = 1 -1 -ξ z-T σ 1/ξ

IV. NUMERICAL INVESTIGATION

The standard map [36] is an area-preserving chaotic map defined on the bidimensional torus, and it is one of the most widely-studied examples of dynamical chaos in physics. The corresponding mechanical system is usually called a kicked rotator. It is defined as:

y t+1 = y t -K 2π sin(2πx t ) mod 1 x t+1 = x t + y t + 1 mod 1 ( 19
)
The dynamics of the map given in Eq. ( 19) can be regular or chaotic. For K << 1 the motion follows quasi periodic orbits for all initial conditions, whereas if K >> 1 the motion turns to be chaotic and irregular. An interesting behavior is achieved when K ∼ 1: in this case we have coexistence of regular and chaotic motions depending on the chosen initial conditions [START_REF] Ott | Chaos in Dynamical Systems[END_REF].

We perform for various values of K ranging from K = 10 -4 up to K = 10 2 an ensemble of 200 simulations, each characterised by a different initial condition ζ randomly taken on the bidimensional torus, and we compute for each orbit the observables g i , i = 1, 2, 3. In each case, the map is iterated until obtaining a statistics consisting 10 4 exceedances, where the threshold T = 7 • 10 -3 and β = 3. We have carefully checked that all the results are indeed robust with respect to the choice of the threshold and of the value of β. For each orbit, we fit the statistics of the 10 4 exceedances values of the observables to a GPD distribution, using a MLE estimation [5] implemented in the MATLAB c function gpdfit [START_REF] Martinez | Computational statistics handbook with MATLAB[END_REF]. The results are shown in Fig. 1 for the inferred values of ξ and σ and should be compared with Eqs. ( 15)- (17). When K ≪ 1, we obtain that the estimates of ξ and σ are compatible with a dimension D = 1 for all the initial conditions: we have that the ensemble spread is negligible. Similarly, for K ≫ 1, the estimates for ξ and σ agree remarkably well with having a local dimension D = 2 for all the initial conditions. In the transition regime, which occurs for K ≃ 1, the ensemble spread is much higher, because the scaling properties of the measure is different among the various initial conditions. As expected, the ensemble averages of the parameters change monotonically from the value pertaining to the regular regime to that pertaining to the chaotic regime with increasing values of K. Basically, this measures the fact that the so-called regular islands shrink with K. Note that in the case of the observable g 1 , the estimate of the ξ is robust in all regimes, even if, as expected, in the transition between low and high values of K the ensemble spread is larger. These results can also be compared with the analysis presented in Faranda et al. [START_REF] Faranda | [END_REF], where we used the Gnedenko approach. In that case, the values obtained in the regular regions were inconsistent with the GEV findings, the very reason being that the dynamics was indeed not mixing. Here, it is clear that the statistics can be computed in all cases, and we have a powerful method for discriminating regular from chaotic behaviors through the analysis of the inferred local dimension.

V. CONCLUSIONS

The growing attention of the scientific community in understanding the behavior of extreme values have led, in the recent past, to the development of an extreme value theory for dynamical systems. In this framework, it has been shown that the statistics of extreme value can be linked to the statistics of return in a neighborhood of a certain initial conditions by choosing special observables that depend on the distance between the iterated trajectory and its starting point. Until now, rigorous results have been obtained assuming the existence of an invariant measure for the dynamical systems and the fulfillment of independence requirement on the series of maxima achieved by imposing D ′ and D 2 mixing conditions, or, alternatively, assuming an exponential hitting time statistics [22][START_REF] Freitas | Probability Theory and Related Fields[END_REF][START_REF] Freitas | [END_REF][START_REF] Gupta | Lozi-like maps, and Lorenz-like maps[END_REF]. The parameters of the GEV distribution obtained choosing as observables the function g i , i = 1, 2, 3 defined above depend on the local dimension of the attractor D and numerical algorithms to perform statistical inference can be set up for mixing systems having both absolutely continuous and singular invariant measures [START_REF] Faranda | [END_REF]28,[START_REF] Vitolo | Chaos: An Interdisciplinary[END_REF]40]. Instead, when considering systems with regular dynamics, the statistics of the block maxima of any observable does not converge to the GEV family [18,19].

Taking a complementary point of view, in this paper we have studied the statistics of exceedances for the same class of observables and derived the limiting distributions assuming only the existence of an invariant measure and the possibility to define a local dimension D around the point ζ of interest. To prove that the limiting distribution is a GPD we did not use any further conditions. In particular no assumptions on the mixing nature of the maxima sequence have been made. This means that a GPD limiting distribution holds for the statistics of exceedance of every kind of dynamical systems and it depends only on the threshold value and on the local dimension once we choose the observables g i , i = 1, 2, 3. Other functions can converge to the limiting behaviour of the GPD family if they asymptotically behave like the g i 's (compare the discussion in [START_REF] Freitas | Probability Theory and Related Fields[END_REF]). Nonetheless, this requires, analytically, to perform separately the limit for the threshold T going to g max and that for the radius of the ball going to zero. In practical terms, this requires, potentially, much stricter selection criteria for the exceedances when finite time series are considered.

We note that, as the parameter ξ is inversely proportional to D, one can expect that each time we analyse systems of intermediate or high dimensionality, the distributions for g 2 and g 3 observables will be virtually indistinguishable from what obtained considering the g 1 observable: the ξ = 0 is in this sense an attracting member of the GPD family. This may also explain, at least qualitatively, why the Gumbel (k = 0 for the GEV family) distribution is so efficient in describing the extremes of a large variety of natural phenomena [3].

The universality of this approach allows to resolve the debate on whether there exists or not a general way to obtain information about extreme values for quasi-periodic motions raised in Balakrishnan et al. [18] and Nicolis et al. [19]. This is due to the fact that considering several of the largest order statistics instead of just the largest one we can study orbits where numerous exceedance are observed in a given block, as it happens for systems with periodic or quasi-periodic behaviors. As an example, one may consider a system with multiple commensurable frequencies: choosing a block length larger or equal to the smallest common period, we select always the same value for any considered observable. On the other hand for mixing maps we find, as expected, an asymptotic equivalence of the results obtained via the Gnedenko and via the Pareto approach.

The Pareto approach provides a way to reconstruct local properties of invariant measures: once a threshold is chosen and a suitable exceedance statistics is recorded, we can compute the local dimension for different initial conditions taken on the attractor. This is true also in the opposite direction: if the knowledge of the exact value for the local dimension is available, once we chose a small enough radius (threshold), it is possible to compute a priori the properties of extremes without doing any further computations. In fact, the expression for the parameters Eqs. ( 15)-( 17) do not contain any dependence on the properties of the dynamics except the local dimension.

Besides the analytical results, we have proved that Pareto approach is easily accessible for numerical investigations. The algorithm used to perform numerical simulations is versatile and computationally accessible: unlike the GEV algorithm that requires a very high number of iterations to obtain unbiased statistics, using the Pareto approach we can fix a priori a value for the threshold and the number of maxima necessary to construct the statistics. With the simulations carried out on the standard maps, we obtain meaningful results with a much smaller statistics with respect to what observed when considering the Gnedenko approach.

We hope that the present contribution may provide a tool that is not only useful for the analysis of extreme events itself, but also for characterising the dynamical structure of attractors by giving a robust way to compute the local dimensions, with the new possibility of embracing also the case of quasi-periodic motions.

Future investigations will include the systematic study of the impact on the extreme value statistics of adding stochastic noise to regular and chaotic deterministic dynamical systems, and the use of the Ruelle response theory [41,42] to study the modulation of the statistics of extremes due to changes in the internal or external pa-rameters of the system, especially in view of potentially relevant applications in geophysics such as in the case of climate studies [43,44].

FIG. 1 .

 1 FIG.1. GPD parameters for the observables gi, i = 1, 2, 3 computed over orbits of the standard map, for various values of the constant K. For each value of K, results refer to an esnsemble of 200 randomly chosen initial conditions ζ. The notation p(gi) indicates the parameter p computed using the extreme value statistics of the observable gi. a) ξ(g1) VS K, b) ξ(g2) VS K, c) ξ(g3) VS K, d) σ(g1) VS K, e) σ(g2) VS K, f ) σ(g3) VS K. Black solid lines: ensemble-average value. Black dotted lines: ensemble spread evaluated as one standard deviation of the ensemble. Green lines: theoretical values for regular orbits. Red lines: theoretical values for chaotic orbits.

ACKNOWLEDGMENTS

The authors acknowledge various useful exchanges with R. Blender and K. Fraedrich, and the financial support of the EU-ERC project NAMASTE-Thermodynamics of the Climate System.


Valerio Lucarini

Tobias Kuna

Jeroen Wouters

Davide Faranda

come