Vincent Perrier

The Chapman-Jouguet closure for the Riemann problem with vaporization

. We check that the Chapman-Jouguet theory can be applied in our context, and that it is consistent with the entropy growth criterion. As the characteristic Lax criterion does not hold for this solution, an additional relation, the kinetic closure is necessary. The common closure, i.e. the Chapman-Jouguet closure is proved to be uncorrect in general in that context.

Introduction

We are interested in the study of some problems arising in the modelling of phase transition in compressible fluids. A widely used model for phase transition [START_REF] Benzoni-Gavage | Stability of subsonic planar phase boundaries in a van der Waals fluid[END_REF][START_REF] Bedjaoui | Diffusive-dispersive traveling waves and kinetic relations. I. Nonconvex hyperbolic conservation laws[END_REF][START_REF] Slemrod | Dynamic phase transitions in a van der Waals fluid[END_REF][START_REF] Tao | A vanishing viscosity approach on the dynamics of phase transitions in van der Waals fluids[END_REF] is the Van-Der-Waals'one. Nevertheless the very physical meaning of this model is questionable, because the resulting system of partial derivative equations is not hyperbolic. Moreover according to [START_REF] Thomson | Shock waves and phase changes in a large heat capacity fluid emerging from a tube[END_REF], the shock structure found in [START_REF] Slemrod | Dynamic phase transitions in a van der Waals fluid[END_REF] does not seem to match with what experiments show. Another approach consists in modelling each phase by an equation of state, and in coupling them by optimizing the entropy [START_REF] Helluy | Finite volume simulation of cavitating flows[END_REF][START_REF] Olivier | Élaboration de lois d'état d'un liquide et de sa Vapeur pour les Modèles d'écoulements Diphasiques[END_REF] to get a mixture equation of state. As explained in [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF], the convexity of internal energy is necessary to ensure the local thermodynamic equilibrium. Whether the mixture equation of state is convex or not will be discussed in a first part.

T temperature τ specific volume ρ density s specific entropy µ chemical potential h specific enthalpy f specific free energy P pressure ε specific internal energy y mass fraction α volume fraction Table 1: Thermodynamic notations We will then concentrate on the Euler system of partial derivative equations, which models flow dynamics without viscosity nor thermal conduction. A fundamental step to approximate the solutions of the Euler system with a Godunov' method is to solve the Riemann problem, i.e. the Cauchy problem where the initial condition is composed of two differents constant states. If the solution of the Riemann problem for the Euler system is easy to solve for convex equation of state, it becomes much harder when the equation of state suffers from loss of derivative and from local non-convexity [START_REF] Voss | Exact Riemann solution for the Euler equations with Nonconvex and Nonsmooth Equation of state[END_REF][START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF][START_REF] Wendroff | The Riemann problem for materials with nonconvex equations of state. I. Isentropic flow[END_REF][START_REF]The Riemann problem for materials with nonconvex equations of state. II. General flow[END_REF][START_REF] Thomson | Shock waves and phase changes in a large heat capacity fluid emerging from a tube[END_REF][START_REF] Thomson | Wave splitting in a fluid of large heat capacity[END_REF][START_REF] Thomson | Negative shock waves[END_REF][START_REF] Tai | The Riemann problem for general systems of conservation laws[END_REF], because the common entropy growth criterion fails to ensure the existence and uniqueness of the solution. Based on the experiments of [START_REF] Roberto | Evaporation waves in superheated dodecane[END_REF], we propose to take into account out of thermodynamic equilibrium states, or metastable states. The Chapman-Jouguet theory will be used [START_REF] Courant | Supersonic Flow and Shock Waves[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], and the compatibility between the model of equation of state and that theory will be discussed. The entropy growth criterion will be shown to hold. Finally we will show that the closure usually used to close the problem [START_REF]Modelling evaporation fronts with reactive Riemann solvers[END_REF] i.e. the Chapman-Jouguet closure leads to a solution that does not depend continuously in general on its initial data, and this solution will therefore be rejected.

1 Thermodynamic preliminaries

Thermodynamic with phase transition

We suppose that we have two phases of the same fluid. For the sake of simplicity, we will consider that one phase is liquid (subscript l) and the other one is the vapor (subscript v). For each phase, we will take the notations of Table 1. If we consider that we have a mixture of the two phases, then the total specific quantities are defined by

τ tot = y l τ l + y v τ v , (1a)
ε tot = y l ε l + y v ε v , (1b)
s tot = y l s l + y v s v . (1c)
To find the thermodynamic equilibrium, the total entropy (1c) must be optimized. Of course, the optimization must be consistent with the following constraints:

conservation of total energy y l ε l + y v ε v = cste, (2a) conservation of mass

y l + y v = 1. (2b)
Moreover, we suppose that the two phases are locally non miscible, which means that

y l τ l + y v τ v = cste. (2c)
We choose to optimize (1c) with the variables τ, ε, y for each phase. The first and second principle of thermodynamic impose that for each phase ds = dε T + P T dτ.

Then the differential of s tot must belong to the set spanned by the gradients of the constraints; if we denote by λ 1 , λ 2 , λ 3 the Lagrange multipliers associated to the constraints (2a),(2b),(2c), we find (with the notations of 1)

s l = λ 1 ε l + λ 2 + λ 3 τ l , (3a)
y l T l = λ 1 y l , (3b)
y l P l T l = λ 3 y l , (3c)
s v = λ 1 ε v + λ 2 + λ 3 τ v , (3d)
y v T v = λ 1 y v , (3e)
y v P v T v = λ 3 y v . (3f)
If we suppose that both phases coexist, then equations (3b) and (3e) give

T l = T v =: T,
and equations (3c) and (3f) give

P l = P v =: P,
and finally, (3a) and (3d) lead to

µ l (P, T ) = µ v (P, T ). ( 4
)
It is well known that there exist a temperature T c and a pressure P c such that there is no more difference between the liquid and the gas above these temperature and pressure (the fluid is said to be supercritical ). If we suppose that under these values we have The thermodynamic plane (P, τ ) is divided into four parts : two parts where the single phases are stable, one part where the mixture is stable, and one part where the fluid is supercritical then we can apply the implicit functions theorem to (4) to prove that P is locally a function of T P = P sat (T ).

τ l (P, T ) = τ v (P, T ),
The limit of thermodynamic stability of the mixture is given, in the (τ, P ) plane by the functions τ v (P sat (T ), T ) and τ l (P sat (T ), T ). The set of all the physical states lying between these curves is called the saturation dome. When there is no ambiguity, we will also denote by τ v the function T → τ v (P sat (T ), T ) and by τ l the function T → τ l (P sat (T ), T ). Thus, the thermodynamic plane can be divided into four parts (see Figure 1) : two parts in which one of the two pure phases is stable, one part where a mixture of saturated liquid and vapor is stable, and where the specific thermodynamic variables are defined as in [START_REF] Bedjaoui | Diffusive-dispersive traveling waves and kinetic relations. I. Nonconvex hyperbolic conservation laws[END_REF], and one part with the supercritical fluid. Moreover, if we differentiate the equality µ v (P, T ) = µ l (P, T ) with respect to (P, T ), we find the Clausius-Clapeyron relation

dP sat dT (T ) = s v (T ) -s l (T ) τ v (T ) -τ l (T ) . ( 5
)
For most of the solid-liquid phase transition (except for special cases like bismuth or water) and for all the liquid-gas phase transition, we have, for all

T < T c τ v (T ) -τ l (T ) > 0.
This means that in the Figure 1, the phase 1 is the liquid, and the phase 2 is a gas.

For any phase transition, the entropy of the most compact constituent is lower than the entropy of the other constituent. Thus in our case we have for all T < T c s v (T ) -s l (T ) > 0, which induces dP sat dT > 0.

Adimensioned thermodynamic coefficients

We adopted the notations of the Table 1. As in [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF], we define three adimensioned parameters as

Γ = - τ T ∂T ∂τ s , γ = - τ P ∂P ∂τ s , g = P τ T 2 ∂T ∂s τ . ( 6
)
γ is called the adiabatic coefficient, and Γ is the Grüneisen coefficient. As done in [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF], we propose to use these three thermodynamic coefficients to express all the thermodynamic quantities. One can then show that the following identities hold

ds = P τ T 2 1 g dT + P T Γ g dτ, (7a) ds
= - τ T Γ γg -Γ 2 dP + P τ T 2 γ γg -Γ 2 dT, ( 7b
) ds = τ T 1 Γ dP + P T γ Γ dτ, ( 7c
) dτ = - τ P g γg -Γ 2 dP + τ T Γ γg -Γ 2 dT, ( 7d
) dh = τ Γ + 1 Γ dP + P γ Γ dτ, ( 7e
) dε = τ 1 Γ dP + P (γ -Γ) Γ dτ. ( 7f
)
Thermodynamic stability requires that ε be convex as a function of τ and s, which leads to

g ≥ 0, γ ≥ 0, γg -Γ 2 ≥ 0. ( 8
)

The fundamental derivative

In [START_REF] Thomson | A fundamental derivative in gasdynamics[END_REF], the fundamental derivative was defined as

G = - τ 2 ∂ 3 ε ∂τ 3 s ∂ 2 ε ∂τ 2 s
The sign of G determines whether the Hugoniot curve and the isentropes are convex or not in the (τ, P ) plane. We will suppose in the following that G is positive, so that no undercompressive discontinuity nor expansion fans can exist (see [START_REF] Thomson | A fundamental derivative in gasdynamics[END_REF][START_REF] Wendroff | The Riemann problem for materials with nonconvex equations of state. I. Isentropic flow[END_REF][START_REF]The Riemann problem for materials with nonconvex equations of state. II. General flow[END_REF]).

1.4 Mixture equation of state In the following, we will denote with a subscript m all the variables relative to the mixture equation of state. The mixture equation of state is naturally parameterized by y, the mass fraction of the gas, and T , the temperature. Nevertheless, in the next sections, the parameters that will intervene are mostly τ and s. They are linked by the transformation

Φ : y T → yτ v (T ) + (1 -y)τ l (T ) ys v (T ) + (1 -y)s l (T ) = τ s . (9)
Theorem 1. For all points in the saturation dome, Φ is a local diffeomorphism provided the equations of state of the liquid and of the gas are both convex.

Proof. To show that Φ is a local diffeomorphism, it is sufficient to show that its Jacobian does not vanish. Differentiation of ( 9) and using the Clausius-Clapeyron relation leads to

det(DΦ) = (τ v -τ l ) y ds v dT - dP sat dT dτ v dT + (1 -y) ds l dT - dP sat dT dτ l dT .
(10) We supposed that τ v -τ l > 0 (except at the critical point), so that there remains to show that

y ds v dT - dP sat dT dτ v dT + (1 -y) ds l dT - dP sat dT dτ l dT , (11)
Using the equations (3a),(3b) for P = P sat (T ) leads to

ds b dT - dP sat dT dτ b dT = γg -Γ 2 g P τ dτ b dT 2 + P τ T 2 1 g > 0, (13)
for b = g or l, which is positive provided each pure phase equation of state is convex. Therefore [START_REF] Tai | The Riemann problem for general systems of conservation laws[END_REF] is positive, because it is a convex combination of two terms like [START_REF] Roberto | Evaporation waves in superheated dodecane[END_REF]. As a consequence, det (DΦ) > 0.

Thanks for the parameterization (9), we can calculate the adimensioned coefficients defined by [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] .

Along an isotherm we have ds = (s v -s l )dy, and dτ = (τ v -τ l )dy, so that ∂s ∂τ T = dP dT sat .

Therefore, we find

∂T ∂s τ = dT dP sat Γ m T τ , which induces γ m g m = Γ 2 m .
As γ m ≥ 0, this means that g m ≥ 0. Therefore, we proved that g m ≥ 0, γ ≥ 0, and that γ m g m -Γ 2 m = 0, so that the convexity of energy is ensured.

Adimensioned coefficients near a phase transition boundary

In that section, we keep on denoting by the subscript m the thermodynamic coefficients of the mixture equation of state, the coefficients with no subscript being the one of the pure phase.

In [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF] (p.121), the following identity is proved

γ -γ m γ m = (γg -Γ 2 ) T τ ds b dP sat 2 > 0,
with b = g or l. This identity proves that isentropes are stiffer in the pure phases than in the mixture. In the same manner it is proved that

Γ m Γ = γ m -ξ γ -ξ , (14)
with ξ = -τ P dP dτ sat . As in [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF], we suppose that the isentropes can be parameterized by τ , so that

γ m -ξ γ -ξ > 0, (15)
so that Γ is positive too, because Γ m > 0.

Retrograde and Regular behavior

In [START_REF] Thomson | Shock waves and phase changes in a large heat capacity fluid emerging from a tube[END_REF], the retrogradicity r was introduced, to study the behavior of isentropes near a phase transition boundary

r = ∂T ∂τ P ds b dP .
Thanks for (7d) and as Γ is positive near a phase transition boundary, ∂T ∂τ P is positive, so that the sign of r is the same as the sign of ds b dP . Figure 2: The saturation dome in the (S, T ) plane. On the left, the fluid is regular: all the isentropes (drawn as arrows) cross the saturation dome from the pure phase to the mixture. On the right, the fluid is retrograde: the isentropes are crossing the liquid saturation curve from the pure phase to the mixture, whereas it is the contrary on the gas side.

We suppose now that a fluid undergoes a rarefaction isentrope : this is the only regular transformation that a fluid can undergo. In the (S, T ) plane, this transformation is drawn as a vertical line. As the transformation is undercompressive, the temperature decreases (at least near the phase transition boundary, because Γ > 0 and Γ m > 0). If r > 0 then the isentrope crosses the saturation curve from the pure phase to to the mixture phase (as on both of the sides of the left Figure of In that case, the fluid is said to be regular. If r is negative, then the isentrope crosses the saturation curve from the mixture to the pure phase, as on the gas side of the right Figure of Figure 2.

In [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF] (p.121), others expressions of r are given

r = Γ m Γ γg -Γ 2 γ m ds b dT = γ -γ m γ m ξ ξ -γ .
Experiments show that the liquid saturation curve is always regular. The gas saturation curve can be either regular or retrograde.

Validity domain of a model

To compute quickly a solution of the Riemann Problem for fluid flows, simplified equation of state (perfect gas of stiffened gas for example) are often preferred to tabulated ones. Nevertheless such equations of state have often only a narrow range of validity, out of which they do not have a physical behavior (negative energy, non convexity).

If we want to use a simplified EOS for both liquid and gas, we have to care not only about the physical behavior of the two EOS, but also about the mixture EOS computed. If we look at the properties needed in Section 1, we see that the property dP sat dT > 0 is fundamental to ensure the local convexity of energy. Nevertheless, it is not always true as we show now on examples.

Two perfect gas

This model was proposed by [START_REF] Jaouen | Étude mathématique et numérique de la stabilité pour des modèles hydrodynamiques avec transition de phase[END_REF][START_REF] Helluy | Finite volume simulation of cavitating flows[END_REF]. The two phases are modeled with a perfect gas equation of state. To complete the equation of state, we suppose moreover that C v = 1 for each fluid. We denote by Γ i the Grüneisen coefficient of the phase i. Then we have

ε i (P, τ ) = P τ Γ i , (16a)
s i (P, T ) = log   T Γ i T P Γi   , (16b)
µ i (P, T ) = (Γ i + 1)T -T log   T Γ i T P Γi   . (16c)
The equation µ 1 (P, T ) = µ 2 (P, T ) can be explicitly solved to get P = βT , with

β = exp(1) Γ Γ2 2 Γ Γ1 1 1 Γ 1 -Γ 2
. We see here that the condition dP sat dT > 0 always holds. The limits of the saturation dome are given by the equations

T = ε = P sat (T )τ i (T ) Γ i which gives τ i (T ) = Γ i β . Thus, T → τ i (T ) is a constant function.
In particular, the critical point does not exist. If we decide for example that Γ 1 < Γ 2 then we get the projections of the phase diagram in the (P, τ ) plane and in the (S, T ) plane that is drawn on Figure 3. The mixture equation of state can be explicitly calculated:

         P (τ, ε) = Γ 2 ε τ if τ ≤ τ 2 P (τ, ε) = Γ 2 ε τ 2 = Γ 1 ε τ 1 if τ 2 ≤ τ ≤ τ 1 P (τ, ε) = Γ 1 ε τ if τ 1 ≤ τ
Nevertheless, we remark that the most heavy phase is described by the lowest adiabatic coefficient, which is in contradiction with what is described for example in [START_REF] Ya | Physics of shock waves and high temperature hydrodynamic phenomena[END_REF] (chapter XI). Thus, the model with two perfect gas is a good mathematical model because the mixture equation of state can be explicitly calculated, but it cannot give a good account for the physic.

τ 1 τ 2
Figure 3: Shape of the saturation dome for two perfect gas. We note that the fluid is always retrograde.

Two stiffened gas

We model the two phases of a fluid with the Stiffened gas equation of state, for which we have (see [START_REF] Olivier | Élaboration de lois d'état d'un liquide et de sa Vapeur pour les Modèles d'écoulements Diphasiques[END_REF])

ε(P, τ ) = P + γP ∞ γ -1 τ + q (17a) s(P, T ) = C v log T γ (P + P ∞ ) γ-1 + q ′ (17b) G(P, T ) = (γC v -q ′ )T -C v T log T γ (P + P ∞ ) γ-1 + q (17c)
For this equation of state, the adimensioned coefficients are given by

γ = γ 1 + P ∞ P Γ = γ -1 g = (γ -1)P P + P ∞ γ > 0 and g > 0 are ensured if γ > 1.
In [START_REF] Olivier | Élaboration de lois d'état d'un liquide et de sa Vapeur pour les Modèles d'écoulements Diphasiques[END_REF], the coefficients q, q ′ , C v , γ, P ∞ were calculated for the gaseous and liquid phases to fit with the saturation curves near T = 298K. These coefficients are in Table 2. The resulting function P sat (T ) was drawn in Figure 4. For the two stiffened gas model, we cannot be sure that the functions T → s v (T ) -s l (T ) and T → τ v (T ) -τ l (T ) simultaneously vanish. Therefore, the critical point does not really exist. As we saw in the Section 1.4, we need that dP sat dT > 0 to ensure the convexity of the mixture equation of state. Thus, the model is valid only when τ v (T ) -τ l (T ) and s v (T ) -s l (T ) are both positive. In our example, with the coefficients of Table 2, the function T → P sat (T ) are drawn on Figure 4. We can see that the limit is near T = 970 K, for which dP sat dT vanishes.

Phase 2. For temperatures below 970 K, T → P sat (T ) increases. For T ≈ 970 K, the function T → s v (T ) -s l (T ) vanishes and its sign changes, whereas the function T → τ v (T ) -τ l (T ) does not vanish. As a consequence, T → P sat (T ) does not increase any more and the equation of state is no more valid. 2. For low temperatures (T ≤ 715 K), the fundamental derivative is positive whereas for 715 ≤ T ≤ 970 K the fundamental derivative is negative.

γ P ∞ C v q q ′ Gas 1.

Positivity of the fundamental derivative

In Section 1.4, we found some criteria to ensure the convexity of the mixture equation of state. Nevertheless, we did not manage to find a simple criterion that can ensure too the positivity of the fundamental derivative of the mixture equation of state. Actually, as it is shown numerically in the Figure 5 for the thermodynamic coefficients of the Table 2, the fundamental derivative of the mixture equation of state can be either positive or negative, even if both equation of states have a positive fundamental derivative. In that example, the equation of state can be considered as valid if T ≤ 715K, where G > 0.

2 Reminds on the Chapman-Jouguet theory (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] pp 142-160)

Based on the experiments of [START_REF] Roberto | Evaporation waves in superheated dodecane[END_REF][START_REF] Thomson | Shock waves and phase changes in a large heat capacity fluid emerging from a tube[END_REF][START_REF] Thomson | Wave splitting in a fluid of large heat capacity[END_REF], we propose to take into account out of thermodynamic equilibrium states for solving the Riemann problem. By "out of equilibrium states" we mean metastable states, or overheated states, i.e. pure fluids that have a pressure P and a specific volume τ that lie in the saturation dome. Existence of such states is due to some phenomena such as capillarity for example. In [START_REF] Roberto | Evaporation waves in superheated dodecane[END_REF][START_REF] Thomson | Shock waves and phase changes in a large heat capacity fluid emerging from a tube[END_REF][START_REF] Thomson | Wave splitting in a fluid of large heat capacity[END_REF] it was observed that phase transition waves were self-similar waves, so that Rankine-Hugoniot relations hold across them: F (U ) -σU = 0, where σ is the velocity of the discontinuity. These relations can be put in the following form (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for example)

Ṁ = u 1 -u 0 τ 1 -τ 0 , (18a)
Ṁ 2 = - P 1 -P 0 τ 1 -τ 0 , (18b)
ε 1 -ε 0 + 1 2 (P 1 + P 0 )(τ 1 -τ 0 ) = 0, ( 18c
)
where Ṁ is the flow rate across the wave: Ṁ = ρ(u -σ). The interest of writing the Rankine-Hugoniot relations as in [START_REF] Thomson | Negative shock waves[END_REF] is that the two last equations are purely thermodynamic. The equation (18b) describes the Rayleigh line in the (τ, P ) plane. Equation (18c) describes the Crussard curve. The very difference with classical shock relations is that the set of the downstream states is not described with the same equation of state as the upstream one. For that sort of wave, we can use the Chapman-Jouguet theory. Let us remind the main points of that theory (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] or [START_REF] Courant | Supersonic Flow and Shock Waves[END_REF] for the details and the proofs)

Property 1 (Position of the initial state and the Crussard curve). Suppose that the equation of state (τ, s) → P (τ, s) has the following properties

∂P ∂τ s < 0 and ∂P ∂s τ > 0, ( 19a
)
and that the reaction is exothermic

ε 1 (τ 0 , p 0 ) < ε 0 (τ 0 , p 0 ); (19b)
then the point A 0 corresponding to the upstream state is under the Crussard curve.

In this first property, note that (19a) is always true provided γ and Γ are both positive (thanks for (7c)). We will find in Section 3.1 a condition to ensure the exothermic property (19b).

Property 2. Suppose moreover that

∂ 2 P ∂τ 2 s > 0;
then the Crussard curve is convex. Hence, the Rayleigh line (18b) and the Crussard curve are crossing in zero or two points.

Note that the Property 2 supposes that the pressure can be differentiated twice, which is not the case in our application, because of the local loss of derivative due to phase transition. Nevertheless, if the properties 1 and 2 hold detonations deflagrations

τ 0 P 0 P τ A 0 Figure 6:
The Crussard curve related to an initial point (τ 0 , P 0 ). The curve is cut into three parts : the upper one is the detonation branch, the lower one is the deflagration branch, and the middle part does not match with the negative slope of the Rayleigh line.

for the equation of state of the downstream states, then the Crussard curve can be schematically drawn as in Figure 6. The Crussard curve is cut into two parts : the upper part is called the detonation branch, and the lower one is the deflagration branch. In the middle part of the curve, P 1 -P 0 τ 1 -τ 0 > 0. This does not match with the negative slope of the Rayleigh line (18b). Each part of the Crussard curve is itself cut into two parts, separated by the tangential point of the Rayleigh line with the Crussard curve (the existence of such tangential point can be shown under some assumption on the asymptotic behavior of equation of state). Both branches are schematically drawn on Figure 7.

|v 0 | > c 0 |v 1 | < c 1 , for a weak detonation : |v 0 | > c 0 |v 1 | > c 1 , for a weak deflagration : |v 0 | < c 0 |v 1 | < c 1 , for a strong deflagration : |v 0 | < c 0 |v 1 | > c 1 .
This last property is very important because it can allow to know the structure of the half Riemann Problem with a combustion wave provided we know which "family" the combustion wave belongs to. In our case, we are interested in waves which transform a heavy phase into a lighter one. Therefore, we expect that τ will increase, so that we will concentrate only on the deflagration branch of the Crussard curve. Thanks for the Property 3, we can give the structure of the Riemann Problem in the case of strong and weak deflagrations. In both cases, from the Property 2, the deflagration wave is always subsonic relative to the liquid; for example, if the liquid is on the left, we have Ṁ > 0, and u 0 -c 0 < σ < u 0 . For the position of the wave relative to the fields of 1, we have

• if the wave is a strong deflagration then σ < u 1 -c 1 • if the wave is a weak deflagration then u 1 -c 1 < σ < u 1 .
In [START_REF] Courant | Supersonic Flow and Shock Waves[END_REF] (p. 230), it is shown that under the assumption that a wave is a deflagration, and that across that wave, the mass fraction of gas always increases, then the reaction is a weak deflagration. We will suppose that we are always in that case in the following.

The structure of the Riemann problem with weak deflagrations is drawn on Figure 8. The problem for deflagrations is that the Lax characteristic criterion is not satisfied (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] p. 154), and the Riemann problem cannot be solved only with the classical relations across the sonic wave and the Rankine-Hugoniot relations across the subsonic wave. 3 Application to the solution of the Riemann problem with vaporization

Useful verifications for the use of CJ theory

In that section, we check whether the inequalities needed for applying the Chapman-Jouguet theory hold.

Theorem 3.

1. If both equations of state are convex, and if dP sat dT > 0, then inequalities (19a) hold.

2. If (P, τ ) lie in the saturation dome and under the condition

γ Γ - T P dP sat dT > 0 ( 20
)
the inequality (19b) holds.

Proof. According to the identity (7c), it is sufficient to have γ and Γ > 0. This is supposed for pure fluids, and this is ensured for mixture equation of state if dP sat dT > 0, so that (19a) holds.

Let us now check if the inequality (19b) is ensured. We suppose that (P, τ ) lie in the saturation dome, so that the corresponding equilibrium downstream state is a mixture :

ε 1 (τ, P ) = ε m (τ, P ) = y l ε l (P, τ l (P )) + (1 -y l )ε v (P, τ v (P )),
and we want to know if ε 1 (P, τ ) -ε 0 (P, τ ) < 0, the state 0 being of course described by the liquid equation of state. For that, we denote by δε(y l ) = y l ε l (P, τ l (P )) + (1 -y l )ε v (P, τ v (P )) -ε l P, y l τ l (P ) + (1 -y l )τ v (P ) , and we immediately see that δε(1) = 0. It remains to show that δε is an increasing function dδε dy l (y l ) = ε l (P, τ l (P )) -ε v (P, τ v (P ))

-(τ l (P ) -τ v (P )) ∂ε l ∂τ P (P, y l τ l (P ) + (1 -y l )τ v (P )).

Integration of the identity dε + P dτ = T ds across the saturation dome leads to ε l (P, τ l (P )) -ε v (P, τ v (P )) + P (τ l (P ) -τ v (P )) = T (s l (P ) -s v (P )), so that dδε dy l (y l ) = -P (τ l (P ) -τ v (P )) + T (s l (P ) -s v (P ))

-(τ l (P ) -τ v (P )) ∂ε l ∂τ P P, y l τ l (P ) + (1 -y l )τ v (P )

, which can be cast into the following form, thanks for ( 5) and (7f) dδε dy l (y l ) = P (τ v -τ l ) γ Γ -T P dP sat dT .

As P > 0, τ v -τ l >, and as [START_REF] Wendroff | The Riemann problem for materials with nonconvex equations of state. I. Isentropic flow[END_REF], δε increases, so that δε ≤ δε(1) = 0. Thus, (19b) holds.

Remark 1

Supposing that (τ, P ) is always in the saturation dome is not a strong assumption. Indeed, as the upstream state is few compressible, its specific volume cannot increase a lot across a sonic wave, and it is likely that a metastable liquid with a specific volume equal to the one of a gas at equilibrium cannot exist, except just near the critical point.

Remark 2

The condition (20) holds at least in the two following framework 1. The terms γ Γ and T P dP sat dT can easily be compared near the saturation curve. Indeed, if we use (7b) with saturated variables, we find

T 2 (γg -Γ 2 ) P τ ds l dT = γ Γ - T P dP sat dT , (21)
so that we have dδε dT

(1) = (τ v -τ l ) T 2 τ γg -Γ 2 Γ ds l dT .
Thus, if the liquid saturation curve is regular (which is always the case) then ds l dT > 0, so that the condition (20) is ensured.

2. For simple model, as perfect gas or stiffened gas, we have

γ Γ = γ 1 + P ∞ P γ -1 ,
so that γ Γ does not depend on the specific volume. Thus, equality [START_REF]The Riemann problem for materials with nonconvex equations of state. II. General flow[END_REF] holds for any τ , so that the condition (20) always holds.

Remark 3

The same calculations can be made for liquefaction. Then we find, near the vapor saturation curve

dδε dT (1) = (τ l -τ v ) T 2 τ γg -Γ 2 Γ ds v dT .
Thus, if the fluid is regular then ds v dT < 0, and as τ l -τ v < 0 then dδε dT (1) > 0, so that locally we have ε m -ε v < 0, and the Chapman-Jouguet theory can be used. If the fluid is retrograde, then we find that locally ε m -ε v > 0 and the Chapman-Jouguet theory may be used, but by exchanging the upstream and the downstream state. Note that the Hugoniot curves might enter the saturation dome only in the retrograde case.

The Chapman-Jouguet theory relies also a lot on the convexity properties of the Crussard curve (see Property 2 of Section 2), which is ensured if the fundamental derivative G is positive. Nevertheless, even if we suppose that the liquid and the gas equation of state have a positive fundamental derivative, the mixture equation of state can have a negative fundamental derivative, as it was shown numerically in Section 1.8. This non positivity of the fundamental derivative can lead to a wrong behavior of the Crussard curve as shown in Figure 9: the CJ points do not exist any more and all the undercompressive downstream states are strong deflagrations. If the sign of the fundamental derivative changed many times along the Crussard curve, we could expect to observe several CJ points. From now, we suppose that G > 0.

Entropy growth criterion

As the particles are crossing the front from the liquid area to a mixture or pure phase area, we have to check whether the entropy growth criterion is ensured, i.e. if the entropy of the downstream state (gas or mixture) is greater than the entropy of the upstream state (liquid). We first prove the following theorem.

Theorem 4. Let s = s 0 a liquid isentrope that crosses the liquid saturation curve. To any metastable point (τ 0 , P 0 ) on that isentrope, we map the point 0 0.01 0.02 0.03 Specific volume (kg.m-3) To any point (τ 0 , P 0 ) on a given isentrope s = s 0 , we associate the point on the Crussard curve (τ P , P 0 ). The liquid saturation curve is drawn on dashed lines. When τ 0 is on the liquid saturation curve, we have τ P = τ 0 , so that s P = s 0 . Thus, to show that s P ≥ s 0 , we only have to prove that entropy of the point P grows when τ 0 increases.

(τ P , P 0 ), point of constant pressure deflagration (see Figure 10). If we suppose that ∀(τ 0 , P 0 ) γ(τ 0 , P 0 ) > γ(τ P , P 0 ) [START_REF] Ya | Physics of shock waves and high temperature hydrodynamic phenomena[END_REF] then s(τ P , P 0 ) > s 0 .

Proof. To any point (τ 0 , P 0 ) on this isentrope, we associate the point (τ P , P 0 ), point of constant pressure deflagration (see Figure 10). For more conveniency in the notations, we suppose that (τ P , P 0 ) is a mixture state (i.e. all the linked quantities have m). τ P is defined by the implicit equation

ε m (τ P , P 0 ) -ε l (τ 0 , P 0 ) + P 0 (τ P -τ 0 ) = 0 ( 23
)
Differentiation of (23) with respect to τ P is equal to γ m Γ m , which never vanishes, so that according to the implicit function theorem, τ P is a C 1 function of τ 0 and P 0 . Moreover, we can calculate its derivative with respect to τ 0 and P 0 :

       ∂τ P ∂τ 0 P0 = γ l γ m Γ m Γ l ∂τ P ∂P 0 τ0 = τ 0 γ m P 0 Γ m (Γ l + 1) Γ l - τ P τ 0 (Γ m + 1)
Besides, as we supposed that the points (P 0 , τ 0 ) belong to the same isentrope, P 0 is actually a function of τ 0 with dP 0 dτ 0 = -γ l P 0 τ 0 , so that τ P is a function of the only variable τ 0 and

dτ P dτ 0 = ∂τ P ∂τ 0 P0 + dP 0 dτ 0 ∂τ P ∂P 0 τ0 = γ l γ m -Γ m + τ P τ 0 (Γ m + 1)
. Now, we calculate the entropy variation of the point τ P when the point (P 0 , τ 0 ) follows the isentrope s = s 0 ds dτ 0 = dτ P dτ 0

∂s ∂τ P + dP 0 dτ 0 ∂s ∂P τ = γ m P 0 T Γ m γ l Γ m γ m τ P τ 0 -1 + τ P τ 0 γ l γ m -1 .
According to the hypothesis [START_REF] Ya | Physics of shock waves and high temperature hydrodynamic phenomena[END_REF], γ l > γ m . Moreover as we have τ P -τ 0 > 0, s is an increasing function of τ 0 . Furthermore, in the limit of no overheating, we have lim

τ0→τ inf τ P (τ 0 , P 0 ) = τ inf ,
where τ inf is the crossing point of the isentrope s = s 0 with the saturation curve. Thus lim τ0→τ inf s(τ P , P 0 ) = s 0 . As a conclusion ∀τ 0 ≥ τ inf s(τ P , P 0 ) ≥ s 0 , which ends the proof.

Remark (About the hypothesis ( 22))

1. We know that near the saturation curve, we have γ l > γ m . For actual data, we have γ l ≫ γ m . So that we can suppose that any γ is greater than any γ m .

2. γ l > γ v just means that the liquid phase is very much less compressible than the gas phase (see e.g. [START_REF] Ya | Physics of shock waves and high temperature hydrodynamic phenomena[END_REF] chapter XI).

Given an initial point, we know that the entropy grows from the constant pressure point to the Chapman-Jouguet point, so that if the entropy growth criterion holds for the constant pressure deflagration point, it holds for all the downstream states between the constant pressure deflagration point and the Chapman-Jouguet point, i.e. Corollary 1. Under the same hypothesis as the Theorem 4, the entropy growth criterion holds for all the weak deflagration points.

Behavior of the Crussard curve near the gas saturation curve

In Section 1.5, the behavior of the isentropes near the saturation curves was studied. The difference of the differential behavior of the pure phase and the mixture equation of state induced kinks in isentropes. Now, we want to study the behavior of the Crussard curve when it crosses the vapor saturation curve. It is more difficult than the study of the isentrope, because the Crussard curve does not depend only on the local variables, but also on the starting point (τ 0 , P 0 ).

General study

We denote by C the point in which the Crussard curve crosses the saturation curve, and by

ζ = - τ P dP dτ |C ,
the adimensioned slope of the Crussard curve.

The first thing we will prove for the behavior of the Crussard curve near the saturation curve, is that it can be parameterized by τ , under some conditions Theorem 5. If all the equations of state are convex and if Γ > 0, then ζ > 0. With the same hypothesis, the Crussard curve can be parameterized by τ , even near the saturation curve.

Proof. As proved in [12, page 101], we have

ζ = γ Γ - ∆P 2P 1 Γ + ∆τ 2τ (24)
Across a deflagration wave, we have ∆τ > 0 and ∆P < 0. Moreover, we proved that Γ m > 0, and we suppose that Γ > 0. The conditions γ > 0 and γ m > 0 were already supposed to ensure the convexity of the specific energy. Then ζ > 0.

If we combine (24) with the identities near the saturation boundary of Section 1.5, we get

ξ -ζ Γ 1 + Γ ∆τ 2τ = ξ -ζ m Γ m 1 + Γ m ∆τ 2τ (25)
Across a deflagration, we have ∆τ > 0. Then

ξ -ζ ξ -ζ m > 0
which means that the Crussard curve, near a boundary, can be parameterized by τ . As ζ > 0, the Crussard curve is a diffeomorphism of τ in each side of the saturation curve. As the Crussard curve can locally parameterize the Crussard curve near a boundary, we conclude that the Crussard curve is a homeomorphism of τ .

Remark As the Crussard curve is a decreasing homeomorphism in τ , the point of constant pressure deflagration is uniquely defined.

To be more precise on the relative behavior of the isentropes, the Crussard curve, and the saturation curve, we will prove that Theorem 6. The relative behavior of the isentropes and the Crussard curves, which gives the nature of the deflagration on each side of the saturation curve follows the alternative

• if γ > ξ then -either ξ ≥ ζ,
then the deflagration is weak on both sides of the saturation curve, -or ξ < ζ, then point C cannot be a weak deflagration simultaneously on both sides of the saturation curve.

• if γ < ξ then either ξ ≤ ζ, then the deflagration is strong on both of the sides of the saturation curve or ξ > ζ, then the deflgration cannot be simultaneously strong on the mixture side, and weak on the pure phase side.

Proof. Equation (25) can be rewritten as

ξ -ζ ξ -ζ m = 1 Γ m + ∆τ 2τ 1 Γ + ∆τ 2τ
so that the discontinuity in the slope of the Crussard curve is directly linked with the sign of Γ m -Γ (we recall that Γ m > 0 and that Γ has the same sign as Γ m near the saturation curves). Equation ( 14) induces a separation into the following cases In that case, as shown on the Figure 11, the relative behavior of the isentrope and the Crussard curve show that in both sides of the saturation curve, the downstream state is a weak deflagration (see Figure 11). In that case, we have

• γ > ξ If γ > ξ then we also have γ m > ξ. As γ m < γ, we have γ m -ξ γ -ξ ≤ 1, so that Γ m ≤ Γ. Therefore ξ -ζ ξ -ζ m ≥ 1 γ ζ γ m ζ m P τ C
ζ ≤ ζ m ≤ ξ ≤ γ m ≤ γ. Suppose now that ξ -ζ ≤ 0. Then we have ξ ≤ ζ m ≤ ζ.
The nature of the deflagration is given by the relative position of the slope of the Crussard curve and the Rayleigh line, so that on the point saturation curve, there are three cases (see Figure 12):

-If the Rayleigh line has a lower slope than both of the slopes of the Crussard curve, then the two parts match with strong deflagrations. Thus, we have γ m ≤ ζ m and γ ≤ ζ (see Figure 12, case (a)).

-If the slope of the Rayleigh line is between the slopes of the Crussard curve, then the mixture Crussard curve matches with strong deflagrations whereas the pure phase Crussard curve matches with weak deflagrations. In that case, we have γ m ≤ ζ m and γ ≥ ζ (see Figure 12, case (b)). In that case we have

ξ ≤ ζ m ≤ γ m ≤ γ ≤ ζ.
-If the Rayleigh line has a greater slope than both of the slopes of the Crussard curve then the point C is a weak deflagration with respect to the pure and the mixture Crussard curve. Therefore, we have γ m ≥ ζ m and γ ≥ ζ (see Figure 12, case (c)).

• γ < ξ

If γ < ξ then we also have γ m < ξ. As we know that γ m ≤ γ, we have γ m -ξ γ -ξ ≥ 1, so that Γ m ≥ Γ (thanks for equation ( 14)). Therefore

ξ -ζ ξ -ζ m ≤ 1
We suppose first that ξ -ζ ≤ 0. Then we immediately have ζ ≤ ζ m . Thus, we have γ m ≤ ζ m and γ ≤ ζ, so that point C matches on both sides of the Crussard curve with strong deflagrations (see Figure 13). In that case we have

ζ m ≤ ζ ≤ ξ ≤ γ ≤ γ m .
We suppose now that ξ -ζ ≥ 0. Then we have ζ m ≤ ζ ≤ ξ. The nature of the deflagration is given by the relative position of the slope of the Crussard curve and the Rayleigh line, so that three cases may happen (see Figure 14):

-If the Rayleigh line has a lower slope than the slopes on both sides of the Crussard curve in C, then the point C is a weak deflagration with respect to the mixture and the pure phase Crussard curve. We have then γ ≤ ζ and γ m ≤ ζ m (see Figure 14, case (a)).

-If the slope of the Rayleigh line is between the slopes on each side of the Crussard curve, then point C is a strong deflagration for the mixture Crussard curve, and a weak deflagration for the pure phase Crussard curve (see Figure 14, case (b)). In that case we have

γ m ≤ ζ m ≤ ζ ≤ γ ≤ ξ.
-If the Rayleigh line has a lower slope than both of the slopes of the Crussard curve then the point C is a strong deflagration with respect to the pure and the mixture Crussard curve. Therefore, we have γ m ≥ ζ m and γ ≥ ζ (see Figure 14, case (c)).

This ends the proof.

Ill-posedness of the Chapman-Jouguet closure

The following result comes immediately from the Theorem 6

Corollary 2. Let (P 0 , τ 0 , u 0 ) be an initial state of liquid at thermodynamic equilibrium, such that the isentrope C s coming from this point enters the saturation dome. If (P ⋆ 0 , τ ⋆ 0 ) is a point in C s , we build (P ⋆ , τ ⋆ ) in the following way • if (P ⋆ 0 , τ ⋆ 0 ) is not in the saturation dome then P ⋆ = P ⋆ 0 , and τ ⋆ = τ ⋆ 0 , As a consequence, the pure phase side matches with a weak deflagrtion too, so that there exist another Chapman-Jouguet point.

• if (P ⋆ 0 , τ ⋆ 0 ) is in the saturation dome, then it is linked with (P ⋆ , τ ⋆ ) with a Chapman-Jouguet deflagration.

If (P ⋆ , τ ⋆ ) can reach the pure gas phase, then the curve (P ⋆ , τ ⋆ ) is discontinuous.

Proof. We suppose that the set described is continuous. As (P ⋆ , τ ⋆ ) can reach the saturation dome, it crosses the gas saturation curve on a point (P c , τ c ). As it is a Chapman-Jouguet point, we have γ m = ζ m , so that we are in case 4 of Figure 14. The case 4.(c) is excluded because the Rayleigh line is tangential with the Crussard curve, so that the slope of the Rayleigh line is greater than the slope of the Crussard curve in the pure phase side. Thus, the point ⋆ matches with a weak deflagration with respect to the pure phase Crussard curve. As a consequence, the Crussard curve has another Chapman-Jouguet point that lies in the pure phase area (see Figure 15), so that the curve C CJ has already a branch in the pure gas area.

Eventually, we can state the following theorem Theorem 7. With the same hypothesis of Corollary 2, if we model the vaporization wave by a Chapman-Jouguet deflagration, then the resulting solution of the Riemann problem is ill-posed in the L 1 sense: the solution does not depend continuously on the initial state.

Proof. We fix a point for x < 0 in the liquid area for which the conditions of Corollary 2 hold, and we suppose that on the right, there is some gas. The Figure 16: In red, the wave curve of the liquid side. In blue, different wave curves for the gas side, depending on the initial state. In the case (I), the gas wave curve intersects the liquid one in one mixture point. In the case (II), the gas wave curve intersects the liquid one in the two branches: one mixture and one pure gas point. In the case (III), the gas wave curve intersects the liquid one only on the pure phase branch.

Riemann problem is composed (from the left to the right) of a sonic wave, a vaporization wave (if the intermediate state is metastable), a contact discontinuity (across which P and u are constant), and a sonic wave in the gas side. As it is done usually [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], to solve the Rieman problem, we intersect the wave curve of the downstream state (sonic wave and maybe followed by a CJ vaporization) of the left side with the wave curve of the sonic wave of the right side, in the plane (P, u). Corollary 2 says that the wave curve of the left side is composed of (at least) two branches (see Figure 16). So that the gas wave curve intersects the liquid wave curve either in one mixture point (case (I)), or in two points (case (II)), or in one pure gas point (case (III)). Existence of case (I) and case (III) implies that we must jump from the mixture to the gas branch of the liquid wave curve. But jumping from one branch to the other means that we change a lot the vaporizes state (so the L 1 loc norm too), but by changing few the initial state.

We finish this paper by drawing the curve C CJ described in Corollary 2 for the models of equation of state of Section 1.

Example 1 : two perfect gas equation of state As an example, we take the model with two perfect gas. As we said before, this model enables to make all the calculation, because the mixture equation of state is explicit.

• Mixture CJ-point

In the case when the downstream state is a mixture, the equation of the Crussard curve is the following

P τ 2 Γ 2 - P 0 τ 0 Γ 2 + 1 2 (P + P 0 )(τ -τ 0 ) = 0,
which gives an expression of τ as a function of P : τ = τ 0 -2 (P τ 2 -P 0 τ 0 ) Γ 2 (P + P 0 ) .

The CJ-point is such that dτ dP (P CJ ) = τ -τ 0 P -P 0 , so that we find the following equation for

P CJ P P 0 2 -2 τ τ 0 P P 0 + 1 = 0,
whose undercompressive solution is

P CJ = P 0   τ 0 τ 2 - τ 0 τ 2 2 -1   .
τ CJ is then given by

τ CJ = τ 2         τ 0 τ 2 + 2 τ 0 τ 2 -1 Γ 2   τ 0 τ 2 + 1 - τ 0 τ 2 2 -1           .
Of course, this point can be chosen only when the mixture is stable, that means when τ CJ ≤ τ 1 .

• Vapor CJ-point

The equation of the Crussard curve is then As P CJ is a linear function of P 0 , we see that τ CJ is a linear function of τ 0 . The CJ-point of pure vapor can be chosen only when τ CJ ≥ τ 1 .

τ P Γ 1 - τ 0 P 0 Γ 2 + τ 1 τ 0 τ 2 τ 2 τ CJ
The two functions τ CJ , for vapor and mixture equation of state are drawn on Figure 17, highlighting the fact that they cannot be linked continuously.

Example 2 : model with two stiffened gas As already mentioned, the mixture equation of state cannot be computed when we deal with the two stiffened gas model. Therefore we can only show a numerical computation as an illustration. We chose the model of dodecane for which coefficients lie in Table 2. We begin on the point P 0 = 900000 Pa with a specific volume of τ 0 = 0.0025 kg.m -3 . We compute all the states 0 ⋆ that can be linked with that initial point via an isentrope. If the state 0 ⋆ is overheated (i.e. lie in the saturation dome), then we compute the Chapman-Jouguet point(s) corresponding to a mixture downstream state and/or to a pure vapor downstream state. Numerical results are on Figure 18.

Conclusion

In a first part we studied the model with two convex equation of state. In particular, we gave a necessary and sufficient condition for the convexity of the mixture equation of state resulting from an entropy optimization:

dP sat dT > 0
Then we proposed to take into account metastable states in the solution of the Riemann problem. For that, we used the Chapman-Jouguet theory. We first proved that this theory can be applied. We emphasized the link between the overheat or overcooled of the metastable state and the retrograde and regular behavior of the fluid. In a particular case, when γ/Γ does not depend on τ , the condition of regular behavior of the fluid is necessary and sufficient to ensure that the energy of a metastable liquid is lower than the energy of a mixture at thermodynamic equilibrium with the same pressure and specific volume.

For the entropy growth condition, we proved that it is ensured provided γ l > γ m and γ l > γ v .

The problem with the deflagration waves is that the Lax characteristic criterion is not ensured, so that the problem is under-determinated. The only thing that we can state with no more hypothesis is that the set of all the downstream states lies in an area limited on the top by the set of all the constant pressure deflagrations, which is continuous, and on below by the set of all the Chapman-Jouguet points, which was proved to be discontinuous thanks for a detailed study of the behavior of the Crussard curve near the saturation curve. As the set of all the Chapman-Jouguet points is discontinuous, the use of the Chapman-Jouguet closure as in [START_REF]Modelling evaporation fronts with reactive Riemann solvers[END_REF] for solving the Riemann problem leads to a solution that does not depend continuously on its initial data in general. A first step to find a right kinetic closure would be for example to study travelling waves for relaxation model as given in [START_REF] Caro | Modélisation et simulation numérique des transitions de phase liquide-vapeur[END_REF]. As we know that liquid-vapor phase transition is governed by a competition between relaxation phenomena and thermal conduction, it would be more relevant (but much harder) to study travelling waves with relaxation model and thermal conductivity.

1 Figure 1 :

 11 Figure 1: The thermodynamic plane (P, τ ) is divided into four parts : two parts where the single phases are stable, one part where the mixture is stable, and one part where the fluid is supercritical

Figure 2

 2 and of the Liquid side of the right Figure of Figure 2).

Figure 4 :

 4 Figure 4: Numerical computation of the behavior of P sat (T ) for two phases of stiffened-gas with the coefficients of Table2. For temperatures below 970 K, T → P sat (T ) increases. For T ≈ 970 K, the function T → s v (T ) -s l (T ) vanishes and its sign changes, whereas the function T → τ v (T ) -τ l (T ) does not vanish. As a consequence, T → P sat (T ) does not increase any more and the equation of state is no more valid.

Figure 5 :

 5 Figure 5: In this Figure, we draw all the T → G for y ∈ [0; 1]. The equation of state is the one of Table2. For low temperatures (T ≤ 715 K), the fundamental derivative is positive whereas for 715 ≤ T ≤ 970 K the fundamental derivative is negative.

Property 3 .

 3 Along the Crussard curve, the velocity |v| = |u -σ| has a local minimum on the CJ-detonation point, and has a local maximum on the CJ-deflagration point. More precisely, we have for a strong detonation :

Figure 7 :

 7 Figure 7: Zoom on the Crussard curve; on the left side, the detonation branch (P ≥ P 0 ) is cut into two parts by the Chapman-Jouguet point. The upper part is the part of the strong detonations, and the lower part is called the part of the weak detonations. On the right side, the deflagration branch (P ≤ P 0 ) is cut into two parts by the Chapman-Jouguet point. The upper part is the part of the weak deflagrations, and the lower part is called the part of the strong deflagrations.

Figure 9 :Figure 10 :

 910 Figure9: Wrong behavior of the mixture Crussard curve if the condition G > 0 is violated. We notice that the Crussard curve is concave which induces no existence of any CJ point. For a more complicated couple of equation of state, we could expect to observe two or three tangential points if the sign of G changed two or three times along the Crussard curve.

Figure 11 :

 11 Figure 11: Qualitative relative behavior of the isentrope and of the Crussard curve when they cross the vapor saturation curve in the case 1. Arrows represent half tangent of the Crussard curve (ζ) and of the isentrope (γ)

Figure 12 :Figure 13 :

 1213 Figure 12: Qualitative relative behavior of the isentrope and of the Crussard curve when they cross the vapor saturation curve, case 2. Arrows represent half tangent of the Crussard curve (ζ) and of the isentrope (γ).

Figure 14 :

 14 Figure 14: Qualitative relative behavior of the isentrope and of the Crussard curve when they cross the vapor saturation curve, case 4. Arrows represent half tangent of the Crussard curve (ζ) and of the isentrope (γ).

Figure 15 :

 15 Figure 15: Qualitative behavior of the Rayleigh line and the Crussard curve when the Crussard curve crosses the saturation curve on a mixture Chapman-Jouguet point. As ζ m ≤ ζ, and as the Rayleigh line is tangential to the mixture Crussard curve, we are in the case (4.c) of the proof of Theorem 6. As a consequence, the pure phase side matches with a weak deflagrtion too, so that there exist another Chapman-Jouguet point.

Figure 17 : 1 + 1 γ 1 - 1 P

 1711 Figure 17: Qualitative behavior of the mixture and the vapor Chapman-Jouguet points for the model with two perfect gas. The horizontal dashed line represents the vapor saturation curve. The increasing line (dashed line, then solid line) is the set of the Chapman-Jouguet points for the vapor equation of state. The other function (solid line, then dashed line) is the set of the mixture Chapman-Jouguet points. The solid lines of the curves correspond to the part in which they match with the equation of state used.

Figure 18 :

 18 Figure18: Set of all the Chapman-Jouguet points that can be reached from a given point. The isentrope is drawn in green, and is nearly vertical. In blue, the set of all the Chapman Jouguet points was drawn, which shows a jump between the mixture Chapman-Jouguet points and the pure vapor Chapman-Jouguet point.

Table 2 :

 2 Thermodynamic coefficients for the liquid and gaseous phase of dodecane.

	025	0	1956.45 -237547. -24485.
	Liquid 2.35 4.10 8 Pa 1077.7	-755269	0.

  There remains one indeterminate. The supplementary relation needed is often called "the kinetic closure". Structure of the half Riemann problem : the state 0 is linked with the state 0 ⋆ by a forerunner sonic wave (rarefaction wave or shock). Then the state 0 ⋆ and the state ⋆ are linked with a deflagration wave. Eventually, there is a contact discontinuity.

	contact surface
	vaporization state ⋆
	state 0 ⋆ sonic wave
	state 0
	Figure 8:

† perrier@math.u-bordeaux1.fr, Université Bordeaux 1, Mathématiques Appliquées de Bordeaux and CEntre des Lasers Intenses et Applications, 33405 Talence Cedex

(P + P 0 )(τ -τ 0 ) = 0.As in[START_REF]Modelling evaporation fronts with reactive Riemann solvers[END_REF], we first calculate the point of constant specific volume detonation, i.e. the downstream state such as τ = τ 0 : P τ = Γ 1 Γ 2 P 0 . If we take the calculations of[START_REF]Modelling evaporation fronts with reactive Riemann solvers[END_REF], we getP CJ = Γ 1 P 0 Γ 2 1 -1 -Γ 2 Γ 1 1 + Γ 2 Γ 1 + 2Γ 2 Γ 1 (γ 1 + 1).

Acknowledgements

This work was partially supported by the CEA through the contract number 123-C-BEFI. The author would like to thank Siegfried Müller from Institut für Geometrie und Praktische Mathematik (Aachen, Germany) for profitable discussions on the Riemann problem with phase transition at thermodynamic equilibrium.


Thierry Alboussière

Bound of dissipation on a plane Couette dynamo

come

Bound of dissipation on a plane Couette dynamo

Introduction

Natural dynamos exist whenever the conditions of their existence is possible, i.e. when a sufficiently large magnetic Reynolds is reached. One can imagine that this is an additional route for mechanical energy dissipation and that it is more likely for natural systems to take it than not to take it. It has been argued sometimes that a state of maximal dissipation rate should be reached and this idea has been used as a closure assumption for turbulence (e.g. Malkus [START_REF] Malkus | Borders of disorder: in turbulent channel flow[END_REF]). It is not necessary to make such an assumption, as explained clearly by Howard [START_REF] Howard | Bounds on flow quantities[END_REF] and yet one can draw useful information from the determination of rigorous lower and upper bounds on energy dissipation in turbulent flows. These bounds are obtained in a larger functional space than the solutions to the Navier-Stokes equation, hence they are not necessarily attained. However, in a number of cases, turbulent flows lead to a dissipation rate of the same magnitude as this upper bound [START_REF] Plasting | Improved upper bound on the energy dissipation rate in plane couette flow: the full solution to busse's problem and the Constantin-Doering-Hopf problem with one-dimensional background field[END_REF].

The variational approach to turbulence was introduced by Malkus [START_REF] Malkus | The heat transport and spectrum of thermal turbulence[END_REF], Busse [START_REF] Busse | On howard's upper bound for heat transport by turbulent convection[END_REF] and Howard [START_REF] Howard | Heat transport by turbulent convection[END_REF]. More recently, this approach was reformulated by Doering and Constantin [START_REF] Doering | Energy dissipation in shear-driven turbulence[END_REF] and expressed in a simpler way, using the concept of a background function (not necessarily the mean flow of turbulence) following Hopf [START_REF] Hopf | Ein allgemeiner endlichkeissatz der hydrodynamik[END_REF]. The objective is to bound energy dissipation under the constraint of horizontally averaged energy balance (for statistically plane invariant configurations). In a series of papers, the approach was improved by optimizing the spectral Lagrange parameter [START_REF] Nicodemus | Improved variational principle for bounds on energy dissipation in turbulent shear flows[END_REF] and the background function [START_REF] Plasting | Improved upper bound on the energy dissipation rate in plane couette flow: the full solution to busse's problem and the Constantin-Doering-Hopf problem with one-dimensional background field[END_REF]. The final bound is better than the nearly rigorous bound of Busse. Our objective here is not to exhaust these possibilities of optimization but rather to apply the general principle to a new configuration of a dynamo problem.

The idea of applying a variational approach to a magnetohydrodynamic flow has been applied already to a Couette and Poiseuille flow subjected to an applied transverse magnetic field [START_REF] Alexakis | Bounds on dissipation in magnetohydrodynamic Couette and Hartmann shear flows[END_REF]. This paper has been a source of motivation for the present work with two significant variations. First, there is now no imposed magnetic field and secondly, magnetic boundary conditions are different. In the paper by Alexakis et al., magnetic disturbances are constrained to vanish at the boundaries, which does not correspond to a physically plausible situation. It is assumed here that the fluid A plane Couette flow configuration is considered and our objective is to find an upper bound to the total energy dissipation when a prescribed velocity is applied. There is no applied magnetic field of external origin and the problem may look like a purely hydrodynamical one. However, the flow may support dynamo action (see [START_REF] Willis | A Taylor-Couette dynamo[END_REF][START_REF] Rogachevskii | Electromotive force and large-scale magnetic dynamo in a turbulent flow with a mean shear[END_REF]) and the amount of dissipated energy must then take into account Joule dissipation. Upper bounds will be obtained as a function of two dimensionless parameters, the Reynolds number and the magnetic Reynolds number. In addition, the effect of electrically conducting walls of finite thickness will be investigated.

Section 2 provides details on the flow configuration, notation, dimensionless variables and equations. In section 3, the principle of decomposition with background flow is presented. The horizontally averaged energy balance is obtained in section 4 and the expression for the total energy dissipation (Joule plus viscous) is given in section 5. Energy dissipation bounds are obtained in section 6 and 7 respectively, when velocity fluctuations and magnetic fluctuations are considered respectively in addition to the background flow. An improved bound is determined numerically in section 8. Section 9 is devoted to a discussion of the bounds obtained and their relevance to experimental and geophysical configurations and section 10 to the directions in which the variational approach could be extended to more relevant models of dynamo.

Plane Couette flow configuration

The dimensionless Navier-Stokes and induction equations can be written

∂u ∂t + u • ∇u = -∇p + j × B + Re -1 ∇ 2 u, (1)
∂B ∂t + u • ∇B = B • ∇u + R -1 m ∇ 2 B. (2)
In the equations above, the dimensional length and velocity scales are chosen to be half the distance, H, and differential velocity, U , between the plates (see Fig. 1) while the quantities H/U , ρU 2 , U √ ρµ, ρ/µU/H are taken as dimensional scales for time, pressure, magnetic field and electric current density respectively. The dimensionless parameters are the Reynolds number Re = U H/ν and the magnetic Reynolds number Rm = µσU H. The symbols ρ, µ, ν and σ denote density, vacuum magnetic permeability, kinematic viscosity and electrical conductivity respectively.

In these dimensionless formulation, the electric current density j is simply the curl of the magnetic field B.

Background velocity decomposition

A frame of reference (x, y, z) is chosen with x aligned with the direction of the imposed shear velocity and z in the direction perpendicular to the plates. Its origin is half-way between the plates.

The velocity field u is written as the sum of a steady parallel profile Φ(z)ex satisfying the boundary conditions Φ(±1) = ±1 and of a field v with homogeneous boundary conditions v(x, y, ±1, t) = 0, i.e.

u(x, y, z, t) = Φ(z)ex + v(x, y, z, t)

It must be stressed here that Φ -the so-called background velocity profile -needs not be the average value of the full velocity field u, neither in the temporal nor in the horizontal space average sense.

Substituting this decomposition into the Navier-Stokes and induction equations leads to evolution equations for v and B for any particular choice of a background profile:

∂v ∂t + v • ∇v + Φ ∂v ∂x + vzΦ ′ ex = -∇p + j × B + Re -1 Φ ′′ ex + ∇ 2 u , (4)
∂B ∂t + v • ∇B + Φ ∂B ∂x = B • ∇v + BzΦ ′ ex + R -1 m ∇ 2 B. (5)
where Φ ′ and Φ ′′ denote the first and second derivative of the background velocity profile Φ respectively.

Energy balance of fluctuations

The dot product of equations ( 4) and ( 5) with v and B respectively are integrated over V a large rectangular volume -L < x < L, -L < y < L and -1 < z < 1. They are then integrated over a long period of time T . The following space and time averaged equations are obtained:

Φ ′ vxvz = (j × B) • v -Re -1 Φ ′ ∂vx ∂z -Re -1 ∇v : ∇v , (6) 0
= Φ ′ BxBz + j • (v × B) -R -1 m j 2 , (7)
where the Poynting flux going out of the control volume has been assumed to be zero (there is no external source or sink) and where the following global space-time average is introduced for any quantity f :

< f >= limT →∞ limL→∞ 1 L 2 T T 0 1 -1 L 0 L 0 f dxdydzdt. (8)
The work of Lorentz forces (j × B) • v and electromotive work j • (v × B) are equal and opposite, so that the sum of ( 6) and ( 7) leads to the following energy balance:

-Φ ′ vxvz + Φ ′ BxBz -Re -1 Φ ′ ∂vx ∂z -Re -1 ∇v : ∇v -R -1 m j 2 = 0. ( 9
)

Total energy dissipation

Let us take ρνU 2 H as a dimensional scale for energy dissipation per cubic metre. The average dimensionless energy dissipation due to viscous effects is ∇u : ∇u and the average dimensionless Joule dissipation is P -1 m j 2 . When u is expressed using the decomposition (3), the total dissipation

D = ∇u : ∇u + P -1 m j 2 becomes D = ∇v : ∇v + Φ ′2 + 2 Φ ′ ∂vx ∂z + P -1 m j 2 . ( 10
)
Combining the energy constraint [START_REF] Nicodemus | Improved variational principle for bounds on energy dissipation in turbulent shear flows[END_REF], so as to remove the linear term Φ ′ ∂vx/∂z , gives the following expression for the dissipation:

D = Φ ′2 -∇v : ∇v + 2Re Φ ′ vxvz -P -1 m j 2 -2Rm Φ ′ BxBz . (11)
Dissipation is bounded by the background dissipation as long as the so-called spectral conditions are satisfied, i.e. ∇v : ∇v + 2Re Φ ′ vxvz > 0 and j 2 -2Rm Φ ′ BxBz > 0 for all admissible vector fields v and B.

Velocity spectral condition

Equation [START_REF] Willis | A Taylor-Couette dynamo[END_REF] shows that the spectral conditions can be treated independently for velocity and magnetic disturbances. The velocity spectral condition consists in ensuring that ∇v : ∇v +2Re Φ ′ vxvz remains positive for all velocity fields v.

In order to satisfy the spectral constraint at large Reynolds numbers, it is convenient to choose a background velocity profile which is uniform in most of the fluid, with a linear profile on each side so as to recover correct boundary conditions (see fig. 2). The thickness δ of these linear parts (virtual boundary layers) is a free parameter. The 'production term' 2Re Φ ′ vxvz is confined to these regions where it can hopefully be balanced by viscous 'dissipation' ∇v : ∇v for a sufficiently small thickness δ.

This problem is solved as follows (Doering and Gibbons [START_REF] Doering | Applied Analysis of the Navier-Stokes Equations[END_REF], Alexakis et al. [START_REF] Alexakis | Bounds on dissipation in magnetohydrodynamic Couette and Hartmann shear flows[END_REF]). Using the condition v = 0 at the lower wall, velocity at another position z in the fluid can be bounded as follows using Schwartz relationship:

vx = z -1+e ∂vx ∂z ′ dz ′ ≤ z -1+e dz ′ z -1+e ∂vx ∂z ′ 2 dz ′ ≤ √ z + 1 -e 1-e -1+e ∂vx ∂z ′ 2 dz ′ . ( 12
)
Combining with the corresponding equation for vz and using Young relationship leads to:

vxvz ≤ 1 2 (z + 1 -e) 1-e -1+e ∂vx ∂z ′ 2 + ∂vz ∂z ′ 2 dz ′ ≤ 1 2 (z + 1 -e)
1-e -1+e

∇v : ∇vdz.

Integrating vxvzΦ ′ over the lower boundary layer yields:

-1+e+δ -1+e vxvzΦ ′ dz ≤ 1 2δ
1-e -1+e

∇v : ∇vdz

-1+e+δ -1+e z + 1 -edz ≤ δ 4
1-e -1+e

∇v : ∇vdz

Averaging over the x and y direction and taking into account the upper boundary layer leads to:

Φ ′ vxvz ≤ δ 2 ∇v : ∇v . (15)
The velocity spectral condition is then satisfied as soon as:

δ ≤ Re -1 . (16)
In terms of energy dissipation [START_REF] Willis | A Taylor-Couette dynamo[END_REF], assuming magnetic effects to be absent, this background profile with sufficiently small boundary layers ( 16) provides an upper bound:

D ≤ 2Re (17)
7 Magnetic spectral condition

In the previous section, it has been possible to find a condition of δ fulfilling the spectral requirement for velocity perturbations. The corresponding task for magnetic perturbations cannot be exactly similar as the property v = 0 on the walls has to be changed into magnetic physical boundary conditions. We are considering that the magnetic field inside the fluid and boundary domains should match a potential field outside decaying to zero at infinity (or possibly a constant if a non-zero mean electric current density is allowed in the electrically conducting domain).

We cannot easily find a simple expression for the boundary condition for Bx or Bz in the physical space, however we obtain below a boundary condition for the x, y average product BxBz. This will then been used to bound Φ ′ BxBz in terms of j 2 .

Let us denote the horizontal average over the (x, y) plane with an overbar:

f (z) = limL→∞ 1 L 2 L 0 L 0 f (x, y, z) dxdy. ( 18
)
Outside the domain, the curl of the magnetic field is zero. So is the product of its y component with Bz averaged over x and y: Bz(∂zBx -∂xBz) = 0. ( 19) Using integration by part, divergence-free condition for B, another integration by parts and the z component of ∇ × B = 0 leads to:

Bz(∂zBx -∂xBz) = ∂zBxBz -Bx∂zBz, = ∂zBxBz + Bx∂yBy, = ∂zBxBz -By∂yBx, = ∂zBxBz + By∂xBy, = ∂zBxBz = 0, (20)
As the product BxBz is vanishing at infinite z, we conclude that BxBz is zero at the boundary of the conduction domain, i.e. when z = ±1. Strictly speaking, this is enough to see that it will be possible to bound BxBz in terms of the mean square of the gradient of B using Poincaré's theorem, and then in terms of the mean joule dissipation. A convenient way of obtaining such a bound is to repeat the calculations leading to (20) while retaining a non-zero electric current density for -1 ≤ z ≤ 1. This leads to:

∂zBxBz = jyBz -jzBy = (j × B)x, (21)
which is the expression of the x component of Lorentz forces in terms of the magnetic stress tensor, averaged on constant z planes. The product BxBz can be evaluated inside the electrically conducting domain using the boundary condition BxBz = 0 at z = ±1.

BxBz(z) = z -1 ∂zBxBzdz ′ , = z -1 (j × B)x, (22)
which can be bounded as follows:

BxBz(z) ≤ z -1 dz ′ z -1 j 2 B 2 dz ′ , ≤ √ z + 1 j 2 B 2 . ( 23
)
As shown in appendix A, the average square magnetic field in the electrically conducting domain is bounded as follows:

B 2 ≤ 4 π 2 j 2 . ( 24
)
Hence, equation ( 23) can be written:

BxBz(z) ≤ 2 π √ z + 1 j 2 . ( 25
)
The magnetic boundary condition is rather loose compared to the velocity boundary condition which leads here to a bound proportional to the square-root of z + 1 for the product BxBz rather than a linear dependence for the velocity product vxvz.

Using equation (25) it is now possible to bound BxBzΦ ′ :

BxBzΦ ′ ≤ 1 δ -1+e+δ -1+e BxBz(z) dz + 1 δ 1-e 1-e-δ BxBz(z) dz, ≤ 4 δπ j 2 e+δ e √ udu. ≤ 8 3δπ j 2 (e + δ) 3/2 -e 3/2 (26)
The magnetic spectral condition j 2 -2Rm Φ ′ BxBz ≥ 0 is satisfied as soon as:

16 Rm 3δπ (e + δ) 3/2 -e 3/2 ≤ 1. ( 27
)
This relationship is then used to express δ is terms of viscous dissipation for the background flow, hence to determine an upper bound for global dissipation. This bound is plotted on Fig. 2 for various values of the wall thickness e. The limiting case of vanishing wall thickness (or rather of an electrically insulating wall) can be easily derived from (27) and corresponds to:

16 Rm δ 1/2 3π ≤ 1, (28)
and since the dissipation of the background profile is 2/δ, we have: 17) for the velocity spectral condition (blue line), and (28) for the magnetic spectral condition (red lines), for different wall conductances e = 0, 10 -3 , 10 -5 and 10 -7 , plotted here for P m = 10 -5

D ≤ 512 R 2 m 9π 2 . ( 29
)
One must take the higher upper bound derived from the velocity and magnetic spectral conditions in Fig. 2 as both conditions must be satisfied. It is straightforward from equations ( 17) and ( 27) that the magnetic spectral condition will take over the velocity spectral condition when the magnetic Reynolds number is exceeding 9π 2 /(256 Pm) ≃ 0.347 P -1 m . It can also be derived from (27) that the critical magnetic Reynolds number at which our upper bound estimate diverges towards infinity is Rm = π/(8 √ e).

Numerical determination of upper bounds

The velocity and magnetic spectral conditions appearing in [START_REF] Willis | A Taylor-Couette dynamo[END_REF]) are tackled numerically in this section. The reason for this is not to achieve a better (i.e. lower) upper bound, although this is the case. The main reason is to check that our analysis in sections 6 and 7 provides results that are of the correct magnitude and not gross overestimates of dissipation upper bounds. For a given background profile, the numerical analysis provides Reynolds and magnetic Reynolds numbers such that the spectral conditions are satisfied for all admissible velocity and magnetic fluctuations, and this cannot be improved as one particular velocity disturbance and one particular magnetic disturbance make the spectral conditions just zero. The method is sketched here: the spectral condition can be treated as a problem of energy stability for the background flow [START_REF] Doering | Applied Analysis of the Navier-Stokes Equations[END_REF]. This is changed into an eigenvalue problem and the game consists in finding a Reynold or magnetic Reynolds number such that the maximum eigenvalue is zero. The velocity and magnetic fluctuations are Fourier transformed in the x and y directions, while they are expanded using Chebyshev collocation polynomials in the z direction. The optimal magnetic disturbance is always found in the limit of zero wave numbers kx and ky, whereas the optimal velocity disturbance is found for kx = 0 and ky finite and increasing with the Reynolds number (see Fig. ( 4)).

For this numerical analysis, the background function is no longer piecewise linear. In order to avoid numerical difficulties due to discontinuities in the velocity gradient, a hyperbolic sine function class is chosen Φ = sinh(zδ). In the limit of large dissipation, the background velocity is zero everywhere except in thin boundary layers where the velocity profile is an exponential function. Similarly to what has been done in sections 6 and 7, the background functions contain a single free parameter, i.e, the typical thickness δ of the boundary layers. Due to the new profile of the background function, the associated dissipation takes a different form that can be derived analytically:

1 -e + sinh 2 1-e δ 2δ sinh 1-e δ 2 .
(30)

For a given choice of δ, hence dissipation, a lower bound of Reynolds and magnetic Reynolds number is sought. This is repeated for various values of δ and the resulting curve can be read as an upper bound of dissipation for each value of Reynolds or magnetic Reynolds number (see Fig. 4 and Fig. 3 respectively). Both spectral conditions can be plotted on the same figure once a magnetic Prandtl number is specified.

It is clear that this numerical procedure leads to improved upper bounds compared to those obtained analytically, while similar trends are obtained.

Discussion

For the simple Couette flow configuration under investigation, and for electrically insulating boundaries, it has been found that a classical 'hydrodynamical' upper bound for dissipation holds when Rm is smaller than around P -1 m : the dimensional upper bound of dissipation per unit mass is of order U 3 /H. For larger magnetic Reynolds numbers though, another scaling is followed and the bound on dissipation becomes proportional to µ 2 σ 2 ν U 4 . This new bound is higher than the 'hydrodynamical' bound and is independent on the half-distance H between the plates.

When the boundary plates have the same electrical conductivity of the fluid and a finite thickness e, the resulting upper bound is changed dramatically. The upper bound corresponding to the magnetic part of the spectral condition increases suddenly and diverges to infinity for a value of the magnetic Reynold number depending on the thickness of the walls. This critical magnetic Reynolds number scales as e -1/2 . This value can be reached much before the change of regime between 'hydrodynamical' and (insulating) 'magnetic' branches of upper bound estimates (Rm ∼ P -1 m ). It should not be too surprising that an dissipation can be unbounded at some finite magnetic Reynolds number. Such a situation can be easily simulated with a Bullard dynamo model. This is a solid dynamo with no possibility of correction of the flow on the dynamo mechanism. Hence magnetic energy (and dissipation) grows without limit above the linear instability threshold. As soon as electrical conductivity is restricted to the fluid (e = 0), this behaviour is no longer possible. For any finite magnetic Reynolds number, energy dissipation is certainly finite. Yet the scaling law for this dissipation is distinctly above the a priori 'hydrodynamical' estimate.

It is appropriate to discuss here how energy dissipation was affected by dynamo action in those three experimental setups where a self-sustained liquid dynamo was observed. They are the Riga (Latvia) [START_REF] Gailitis | The riga dynamo experiment[END_REF], Karlsruhe (Germany) [START_REF] Müller | A two-scale hydromagnetic dynamo experiment[END_REF] and Cadarache (France) [START_REF] Monchaux | Generation of a magnetic field by dynamo action in a turbulent flow of liquid sodium[END_REF] dynamos. From these references, a table can be shown (table 9), where the estimated extra dissipation when dynamo action is present is given (in percentage) for a certain magnetic Reynolds number (expressed in percentage above the threshold). the Karlsruhe dynamo experiment can be clearly distinguished by the large amount of Riga [START_REF] Gailitis | The riga dynamo experiment[END_REF] Karlsruhe [START_REF] Müller | A two-scale hydromagnetic dynamo experiment[END_REF] Cadarache [ The authors of the dynamo experiments themselves specify that the extra dissipation is obtained by difference of the actual dissipation and the extrapolated curve of dissipation below threshold. Even though it has been used to estimate Ohmic dissipation by some authors, it should be stressed here that there is no theoretical reason why this extra dissipation should be the Ohmic contribution to dissipation. For instance, in the present work, we use the fact that there is an exact cancellation between the work of Lorentz forces and the electromotive work, but we have no access to either of these terms.

We have concentrated mainly on the low magnetic Prandtl number case, because our primary interest lies in the Earth's core dynamics and in the relevant liquid metal experiments. However the results we have obtained are equally valid for large magnetic Prandtl numbers. It is however anticipated that our dissipation bound will be a gross overestimate at low values of the hydrodynamic Reynolds number. In that case, the flow is the simple uniform shear flow between the moving walls and such a flow is not going to sustain dynamo easily. This perspective should not be ruled out as strong Lorentz forces might change the flow in such a way that it could drive a dynamo, however this transition if it exists will be severely subcritical.

In the derivation of velocity and magnetic spectral conditions, there are strong similarities. The only difference is due to boundary conditions for velocity disturbances and magnetic field disturbances. This mere difference in the boundary conditions leads to fundamentally different upper bounds for dissipation. Magnetic boundary conditions are sometimes not well treated in general textbooks (see, for instance, the energy stability of hydromagnetic flows in the otherwise excellent book by Joseph [START_REF] Joseph | Stability of fluid motions I[END_REF], volume II). It is confirmed again in the present work that one should pay great attention to the relevant physical boundary conditions.

Perspectives

It is certainly useful to derive such upper bounds for dissipation in the field of magnetofluid dynamics, and more particularly regarding dynamos. When complex situations are considered (magnetic field, global rotation, ...) there is no reliable heuristic approach that can provide good estimates for dissipation. Ordinary hydrodynamic turbulence can be roughly tackled by Kolmogorov's theory but it is not safe to extend it to other types of turbulence. When the 'Doering-Constantin-Hopf' background profile method can provided an upper bound, this is a solid reliable result. The bound can be rather constraining for instance in the case of a Couette flow with electrically insulating boundaries, at low Pm: it is not an obvious result that dissipation must remain similar to hydrodynamic dissipation up until Rm ∼ P -1 m . There are actually a lot of other configurations for which the method can be applied.

Upper bounds for dissipation can be calculated for experimental setups involving liquid metal flows. In Grenoble, LGIT laboratory, we have the DTS setup involving 40 liters of sodium in an imposed magnetic field [START_REF] Nataf | Experimental study of super-rotation in a magnetostrophic spherical Couette flow[END_REF][START_REF] Nataf | Rapidly rotating spherical couette flow in a dipolar magnetic field: an experimental study of the mean axisymmetric flow[END_REF][START_REF] Schmitt | Rotating spherical couette flow in a dipolar magnetic field: Experimental evidence of magneto-inertial waves[END_REF]. Upper bounds for dissipation would provide some information about the importance of the size of the experimental setup. This is crucial when planning a larger setup that might sustain dynamo action.

Thermal convection can also be taken into account to derive dissipation bounds. This is the case for the flow of liquid iron in the Earth's core. Upper bounds of dissipation may prove useful in terms of the thermal budget of the history of the Earth. Although the configuration of the core of the Earth is not similar to our simple Couette flow considered here, it could be very relevant to examine the role of a thin electrically conducting layer at the bottom of the mantle, as it could make a large change in the upper bound for dissipation.

Fourier series in the x and y directions, i.e. p = Σp k (z) exp(ikxx + ikyy) and a generic Fourier component satisfies:

p k + λ(p ′′ k -k 2 p k ) = 0, (37)
where k = k 2 x + k 2 y . It does not matter which component is considered, as all lead eventually to the same maximal value for λ. Let us consider p k for instance. Equation (37) admits solutions of the form:

p k = A cos z 1 λ -k 2 + B sin z 1 λ -k 2 , (38)
which given the boundary conditions p k = 0 when z = ±1, leads to the following maximal value for λ: λ < 1

π 2 4 + k 2 , ( 39
)
The maximal value 4/π 2 is reached when k is smallest. The same value is obtained when Bx, By or a Fourier component of q are considered.

Figure 1 :

 1 Figure 1: Plane Couette flow configuration

Figure 2 :

 2 Figure 2: Bound on dissipation obtained from equations[START_REF] Joseph | Stability of fluid motions I[END_REF] for the velocity spectral condition (blue line), and (28) for the magnetic spectral condition (red lines), for different wall conductances e = 0, 10 -3 , 10 -5 and 10 -7 , plotted here for P m = 10 -5

Figure 3 :

 3 Figure 3: Numerical magnetic upper bound for energy dissipation for different wall thicknesses

Figure 4 :

 4 Figure 4: Numerical upper bound for energy dissipation from the velocity spectral condition and associated spanwise wave-number k y of marginal perturbation

Table 1 :

 1 Extra dissipation due to dynamo action slightly above threshold in experimental dynamosextra dissipation compared to the other two. This experiment is also one where differential velocity is forced deep inside the setup. For the Riga and Cadarache dynamos, a turbine or propellers are set in rotation near the outer boundary of the setup. There can only be a gross comparison to the pure Couette flow configuration presented here, but the case of thick electrically conducting boundaries is also putting shear velocities well inside the electrically conducting domain and results in a huge increase of the dissipation bound. Conversely, it cannot be ruled out that Riga and Cadarache dynamos are still obeying the same classical 'hydrodynamical' law of dissipation, similarly to the Couette flow with electrically insulating walls.

	16]

Acknowledgements

Thanks are due to Alexia Gorecki and Maylis Landeau who have realized projects on closely related topics and provided stimulating discussions. I am also grateful to the members of the Geodynamo team (LGIT lab in Grenoble) for their useful and benevolent comments on this work.

A Bound on magnetic energy for a given amount of Joule dissipation

We show here that the average magnetic energy in the conducting domain is bounded by Joule dissipation. In our dimensionless terms, the following relationship hods:

Let us look for the maximum value of the ratio B 2 / j 2 , which we denote λ. For that particular magnetic field, variations of this ratio are zero for all arbitrary magnetic perturbations δB and associated electrical current perturbations δj (such that curlδB = δj). The first variations around this maximum must vanish:

Integrating by parts on the conducting domain shows that j must vanish on the boundary and that:

in the electrically conducting domain. We then decompose the magnetic field in terms of plane polo ïdal-toro ïdal variables:

Equation (34) leads to q + λ∇ 2 q = 0, p + λ∇ 2 p = 0,

with boundary conditions on z = ±1: q = 0, p = 0, B ′ x = 0 and B ′ y = 0, where a prime subscript denotes differentiation with respect to the z direction. The variables p and q can be expanded in