Richard Saurel
email: richard.saurel@univ-amu.fr

Ashwin Chinnayya
email: ashwin.chinnayya@ensma.fr

Quentin Carmouze
email: quentin.carmouze@rs2n.eu

Modelling compressible dense and dilute two-phase flows

Keywords: hyperbolic, two-phase, compressible, shocks

published or not. The documents may come

Modelling compressible dense and dilute two-phase flows

I. Introduction

It is well accepted that hyperbolic models are mandatory to deal with phenomena involving wave propagation. This is the case for multiphase flow mixtures in many situations such as in particular shocks and detonations propagation in granular explosives and in fuel suspensions, as well as liquid-gas mixtures with bubbles, cavitation and flashing, as soon as motion is intense and governed by pressure gradients. This is thus the case of most unsteady two-phase flows situations. Wave propagation is important as it carries pressure, density and velocity disturbances. Sound propagation is also very important as it determines critical (choked) flow conditions and associated mass flow rates. It has also fundamental importance on sonic conditions of detonation waves when the twophase mixture is exothermically reacting [START_REF] Petitpas | Modelling detonation waves in condensed energetic materials: Multiphase CJ conditions and multidimensional computations[END_REF]. Hyperbolicity is also related to the causality principle, meaning that initial and boundary conditions are responsible of time evolution of the solution. When dealing with first-order partial differential equations it means that the Riemann problem must have a solution, and the Riemann problem is correctly posed only if the equations are hyperbolic. However, only a few two-phase flow models are hyperbolic in the whole range of parameters. The [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] model seems to be the only formulation able to deal with such requirement. Its symmetric extension [START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF] facilitates the Riemann problem resolution as the corresponding model involves 7 wave's speeds (instead of 6 in the original version). See also [START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF] for similar conclusions. However, in the dilute limit at least, the acoustic properties of this model seem inconsistent [START_REF] Lhuillier | On the quest for a hyperbolic effective-field model of disperse flows[END_REF]. Indeed with this model, the dispersed phase sound speed corresponds to the one of the pure phase, while this phase is not continuous and unable to propagate sound in reality, at least at a scale larger than particle's one. When the phase is not continuous (dispersed drops in a gas, dispersed bubbles in a liquid), the associated sound speed should vanish, such effect being absent in the formulation. In the low particles concentration limit, the [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF] model is preferred. This model corresponds to the Euler equations with source terms for the gas phase and pressureless gas dynamic equations for the particle phase (see also [START_REF] Zeldovich | Gravitational instability: An approximate theory for large density perturbations[END_REF]. This model is thermodynamically consistent and hyperbolic as well, except that the particle phase equations are hyperbolic degenerate. In this model, contrarily to the BN model, sound doesn't propagate in the particles phase, this behaviour being more physical in this limit. However, the Marble model has a limited range of validity as the volume of the dispersed phase is neglected, this assumption having sense only for low (less than per cent) condensed phase volume fraction. There are thus fundamental differences between these two models:

-The volume occupied by the condensed phase is considered in BN while it is neglected in the dilute model, restricting its validity to low dispersed phase volume fractions. -Condensed phase compressibility is considered in BN while incompressible particles are assumed in the dilute formulation. -Acoustic properties of the BN model are well accepted in the dense domain but seem inappropriate in the dilute limit. Even if these two models can be used in the entire space of two-phase flow variables without yielding computational failure (this is characteristic of thermodynamically consistent hyperbolic models) validity of their results is questionable when they are used out of their range of physical validity. This issue has been clearly understood in [START_REF] Lhuillier | On the quest for a hyperbolic effective-field model of disperse flows[END_REF], [START_REF] Mcgrath | A compressible two-phase model for dispersed particle flows with application from dense to dilute regimes[END_REF] and Houim and Oran (2016) where various attempts to build new formulations are reported. In [START_REF] Lhuillier | On the quest for a hyperbolic effective-field model of disperse flows[END_REF] discussion on the volume fraction equation is done, but no explicit flow model is given. In [START_REF] Mcgrath | A compressible two-phase model for dispersed particle flows with application from dense to dilute regimes[END_REF] a model is given with conditional hyperbolicity. Same issue is present with different cause in Houim and Oran (2016). The aim of the present paper is to build an alternative to the BN model with improved acoustic properties, while remaining unconditionally hyperbolic and thermodynamically consistent. The new model is derived from number density and particle radius (or bubble radius) evolution equations resulting in a volume fraction evolution equation expressed in conservation form with pressure relaxation. Replacing the transport equation in the conventional balance equations of mass, momentum and energy of the phases has dramatic influence on wave's propagation and structure of the equations. The paper is organised as follows. The well-known BN and Marble models are recalled in Section II to present the main alternatives of existing two-phase hyperbolic models. Derivation of the volume fraction equation of the new model is addressed in Section III. The new model is then derived in Section IV and its hyperbolicity demonstrated. Its compatibility with the model of [START_REF] Kapila | Two-phase modeling of deflagrationto-detonation transition in granular materials: Reduced equations[END_REF] is demonstrated as asymptotic limit of the new model, in the limit of stiff mechanical relaxation. Computed results are then examined in Section V, compared to exact and experimental solutions when available. A symmetric variant of the new flow model, aimed to model bubbly liquids, is derived in Section VI and typical solutions are examined. A general model, aimed to address the full spectrum of volume fractions is then derived in Section VII. Conclusions are given in Section VIII.

II. Well-known limit models of two-phase flows

Two hyperbolic models are widely used in the two-phase flow literature and their main characteristics are recalled hereafter.

a) BN type model (1986)

The [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] model is recalled hereafter, in the absence of granular effects ('configuration' pressure and energy) as well as heat and mass transfers. Mechanical relaxation effects only are considered in addition to waves' dynamics. A variant of this model is available as well in [START_REF] Romenski | Compressible two-phase flows: two-pressure models and numerical methods[END_REF], where a conservative formulation is obtained. The symmetric variant of [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shock-bubble interaction[END_REF] is presented hereafter rather than the original BN. The evolution equations for phase 1 read,

1 1 I 1 2 u (p p ) t x ∂α ∂α + = µ - ∂ ∂ 1 1 ( ) ( u) 0 t x ∂ αρ ∂ αρ + = ∂ ∂ (II.1) 2 1 1 1 I 2 1 ( u) ( u p) p (u u ) t x x ∂ αρ ∂ αρ + α ∂α + = + λ - ∂ ∂ ∂ ' 1 1 1 I I I 2 1 I 1 2 ( E) ( ( E p)u) p u u (u u ) p (p p ) t x x ∂ αρ ∂ α ρ + ∂α ′ + = + λ --µ - ∂ ∂ ∂ The evolution equations of the second phase are, 2 2 ( ) ( u) 0 t x ∂ αρ ∂ αρ + = ∂ ∂ (II.2) 2 2 2 2 I 2 1 ( u) ( u p) p (u u ) t x x ∂ αρ ∂ αρ + α ∂α + = -λ - ∂ ∂ ∂ ' 2 2 2 I I I 2 1 I 1 2 ( E) ( ( E p)u) p u u (u u ) p (p p ) t x x ∂ αρ ∂ α ρ + ∂α ′ + = -λ -+ µ - ∂ ∂ ∂ With the following definitions and notations: -k α , k ρ , k u , k E , k
p denote respectively the volume fraction, material density, velocity, total energy and pressure of the phase k (k=1,2).

-The total energy of the phases reads,

2 k k k 1 E e u 2 = + .
-The pressures are given by convex equations of state of the form

k k k k p p ( , e ) = ρ .
-The velocities relax each other to a common equilibrium one at a rate controlled by λ, modelled by conventional drag force correlations and specific interfacial area.

-The pressure relax each other to a common equilibrium one at a rate controlled by µ. Estimates for this relaxation parameter are given in the references above:

I 1 2 A Z Z µ = +
, where I A represents the interfacial exchange area. The specific interfacial area is given by

1 I 1 3 A R α =
if phase 1 represents the dispersed phase made of particles or bubbles of constant radius 1 R . Obviously, more sophisticated models of interfacial area are possible.

-The interfacial variables are estimated by,

' 1 2 1 I I 1 2 p p u u sgn x Z Z ∂α -   = +   ∂ +   with 2 1 2 2 1 1 ' I Z Z u Z u Z u + + = , 2 1 1 2 1 1 2 I 2 1 1 2 1 2 Z p Z p Z Z p sgn (u u ) Z Z x Z Z + ∂α   = + -   + ∂ +   with ' 2 1 1 2 I 1 2 Z p Z p p Z Z + = + ,
where

k k k Z c = ρ represents the acoustic impedance of phase k.
This symmetric formulation of the BN model has some advantages:

-Its extension to more than two phases is quite easy.

-It is able to deal with contact and permeable interfaces [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shock-bubble interaction[END_REF][START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF].

-It involves an extra wave, not aligned with the condensed phase velocity, this property having benefits at least for numerical resolution (Ambroso et al., 2012, Furfaro and[START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF]. This system admits the following mixture entropy equation:

( ) ( ) 1 1 2 2 1 1 1 2 2 2 2 2 2 1 1 1 2 1 2 1 2 2 1 2 1 2 1 2 1 1 2 1 2 1 2 2 2 1 2 1 2 1 2 2 1 2 2 1 2 ( ) s ( ) s ( ) u s ( ) u s t x Z Z Z 1 (p p ) sgn Z (u u ) (u u ) (p p ) T x x Z Z Z Z Z Z Z Z 1 (p p ) sgn Z (u u ) T x x Z Z ∂ αρ + αρ ∂ αρ + αρ + = ∂ ∂    ∂α  ∂α     - + - + λ - + µ -       ∂ ∂ + +   +        ∂α  ∂α   + - + - + λ     ∂ ∂   +   2 2 2 2 1 2 1 1 2 1 2 Z (u u ) (p p ) Z Z Z Z     - + µ -   + +    
Its 7 associated wave speeds are:

I I u λ = , 1 1 u λ = , 2 1 1 u c λ = + , 3 1 1 u c λ = -, 4 2 u λ = , 5 2 2 u c λ = + , 6 2 2 u c λ = -.
This model is consequently hyperbolic, thermodynamically consistent and symmetric. However, the wave speeds are independent of the volume fraction, meaning that in the dilute limit, the sound speed in the condensed phase is unchanged, this behaviour being questionable as this phase is no longer continuous.

b) Dilute two-phase flow model (Marble, 1963)

As the model that follows is no longer symmetric it is necessary to precise the phases. Phase 1 is considered to be the condensed one and the gas phase is denoted by the subscript 2. The 'apparent density' of the dispersed phase is introduced as, 1 1 ( ) ρ = αρ .

In this approach, 1 0.01 α < and volume fraction effects are neglected in the gas phase equations. Phase 1 (dispersed)

1 1 1 u 0 t x ∂ρ ∂ρ + = ∂ ∂ 2 1 1 1 1 2 1 u u (u u ) t x ∂ρ ∂ρ + = λ - ∂ ∂ (II.3) 1 1 1 1 1 e e u 0 t x ∂ρ ∂ρ + = ∂ ∂ or alternatively 1 1 1 1 1 1 2 1 E E u u (u u ) t x ∂ρ ∂ρ + = λ - ∂ ∂ .
Phase 2 (gas)

2 2 2 u 0 t x ∂ρ ∂ρ + = ∂ ∂ (II.4) 2 2 2 2 2 2 2 1 u u p (u u ) t x ∂ρ ∂ρ + + = -λ - ∂ ∂ 2 2 2 2 2 2 1 2 1 E ( E p )u u (u u ) t x ∂ρ ∂ ρ + + = -λ - ∂
∂ This system admits the following mixture entropy equation:

1 1 2 2 1 1 1 2 2 2 1 2 2 s s s u s u (u u )² t x T ∂ρ + ρ ∂ρ + ρ λ - + = ∂ ∂
Its associated wave speeds are:

1 1 u λ = , 2 2 u λ = , 3 2 2 u c λ = + , 4 2 2 u c λ = -. As 1 1 u λ =
is fold three times, the equations of phase 1 are hyperbolic and linearly degenerate, while the ones of the gas phase are strictly hyperbolic. These two models are thus well posed in the sense that they are thermodynamically consistent, frame invariant and hyperbolic. Both models can be solved by Godunov type methods as the Riemann problem has been addressed for both [START_REF] Saurel | Two-phase flows-Second-order schemes and boundary conditions[END_REF][START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF], Deledicque and Papalexandris, 2010[START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF]. However, well posedness is a necessary condition but not a sufficient one for physical validity. In particular, considering again the BN model, the condensed phase sound speed 1 c is well defined as a thermodynamic variable and sound disturbances propagate at the particle or grain level. But sound cannot propagate in the mixture at speed 1 c as the continuum approximation is no longer valid for the condensed phase as soon as the mixture becomes dilute enough. See also [START_REF] Lhuillier | On the quest for a hyperbolic effective-field model of disperse flows[END_REF] and [START_REF] Mcgrath | A compressible two-phase model for dispersed particle flows with application from dense to dilute regimes[END_REF] for further arguments.

III. Alternative volume fraction equations

The volume fraction equation of the BN model is the first equation of System (II.1) and can be derived from averaging method considering the transport of a characteristic function, equal to 1 in a given phase and 0 in the other phase. See for example [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF], [START_REF] Drew | Theory of multicomponent fluids[END_REF]. Let us now consider another point of view as done by [START_REF] Lhuillier | On the quest for a hyperbolic effective-field model of disperse flows[END_REF] and consider liquid drops (or condensed phase particles) suspended in a gas. The radius R 1 of a single spherical compressible liquid drop surrounded by a gas evolves, under acoustic approximation ( u p / c ∆ ≈ ∆ ρ ) with the following transport equation [START_REF] Chinnayya | Modelling detonation waves in heterogeneous energetic materials[END_REF],

1 1 1 2 1 1 d R p p dt c - ≈ ρ , (III.1)
where

1 1 d u dt t x ∂ ∂ = + ∂ ∂ denotes the Lagrangian derivative of phase 1.
Estimate of the pressure relaxation time 1 τ is given by the time needed for an acoustic wave to travel the particle radius,

1 1 1 R c τ ≈ .
For a liquid drop of 1 mm radius suspended in air the pressure relaxation time is therefore of the order of 1 microsecond. This is very small in most practical situations compared to the other characteristic times related to drag, heat exchange, most situations of wave propagation and fluid motion. With this definition (III.1) becomes,

1 1 1 1 2 2 1 1 d R R p p dt c - ≈ τ ρ
Trivial transformation of the former equation implies,

2 1 1 1 1 2 1 2 1 1 1 d V R p p 4 R dt c - ≈ π τ ρ (III.2)
where

3 1 1 4 V R 3 = π denotes the volume of the drop.
In absence of fragmentation and coalescence, the specific number of drops per unit volume obeys the following balance equation:

1 1 1 N N u 0 t x ∂ ∂ + = ∂ ∂ , (III.3)
where 1 N represents the specific number of drops. Multiplying (III.2) by 1 N yields,

1 1 1 1 1 2 2 1 1 1 u 3 p p t x c ∂α ∂α α - + = ∂ ∂ τ ρ (III.4) as 1 1 1 N V α =
. The volume fraction equation is now in conservative form with a pressure relaxation term. It is interesting to consider the symmetric situation of liquid containing spherical bubbles. In this situation the bubble radius evolves according to,

2 2 2 1 1 1 d R p p dt c - ≈ ρ ,
as the acoustic impedance of the less compressible phase ( 1 1 c ρ ) controls the interface velocity.

The specific number of bubbles per unit volume obeys the balance law,

2 2 2 N N u 0 t x ∂ ∂ + = ∂ ∂ ,
and the corresponding volume fraction equation now reads,

2 2 2 2 2 1 2 1 1 2 u 3 p p t x c c ∂α ∂α α - + = ∂ ∂ τ ρ (III.5) with 2 2 2 R c τ ≈ .
We now examine the implications of such volume fraction equations (III.4 and III.5) on the flow model. The analysis begins with a model based on (III.4) to start with a concrete example.

IV. The new model

For the sake of simplicity in the notations and compatibility with (II.1), Equation (III.4) is expressed as

1 1 1 1 2 u μ(p p ) t x ∂α ∂α + = - ∂ ∂ , (IV.1) with 1 2 1 1 1 3 μ c α = τ ρ
, where the estimate for the pressure relaxation time 1 τ has been inserted.

The same mass, momentum and energy equations of Systems (II.1)-(II.2) are reconsidered as,

1 1 ( ) ( u) 0 t x ∂ αρ ∂ αρ + = ∂ ∂ 2 1 1 1 I 2 1 ( u) ( u p) p (u u ) t x x ∂ αρ ∂ αρ + α ∂α + = + λ - ∂ ∂ ∂ 1 1 1 I I 2 1 2 1 ( E) ( ( E p)u) p u (u u ) H(T T ) t x t ∂ αρ ∂ α ρ + ∂α + = - + λ -+ - ∂ ∂ ∂ 2 2 ( ) ( u) 0 t x ∂ αρ ∂ αρ + = ∂ ∂ (IV.2) 2 2 2 2 I 2 1 ( u) ( u p) p (u u ) t x x ∂ αρ ∂ αρ + α ∂α + = -λ - ∂ ∂ ∂ 2 2 1 I I 2 1 2 1 ( E) ( ( E p)u) p u (u u ) H(T T ) t x t ∂ αρ ∂ α ρ + ∂α + = -λ -- - ∂ ∂ ∂
The right hand side of the phase energy equations has been modified with the presence of

1 I p t ∂α ∂ , the interstitial pressure work, present in (II.1-2) differently.
System (IV.2) obviously satisfies mixture mass, mixture momentum and mixture energy conservation, for any model of interfacial pressure I p and interfacial velocity I u . Possible estimates are for example [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shock-bubble interaction[END_REF]:

2 1 1 2 I 1 2 Z p Z p p Z Z + = + , 1 1 2 2 I 1 2 Z u Z u u Z Z + = + , (IV.3) Convective heat exchange ( 2 1
H(T T ) -) has been inserted for the sake of generality where H denotes the product of the specific interfacial area and heat exchange coefficient, related to the Nusselt number. Balance equations (IV.2) are considered not only in the BN formulation, but in any two-phase Eulerian model when the effects of volume fraction are considered. The only point to underline is that the various pressures are distinct at this level. Two questions arise immediately, regarding the fulfilment of the second law of thermodynamics and the hyperbolicity of (IV.1-2). In this aim, the equations are expressed in a set of appropriate variables.

a) Physical variables formulation System (IV.1-2) is expressed with 'physical variables': volume fraction, density, velocity, internal energy and entropy for each phase:

1 1 1 1 1 2 d u μ(p p ) dt x α ∂ + α = - ∂ (IV.4) 1 1 1 1 2 1 d μ (p p ) dt ρ ρ = - - α 1 1 1 I 1 1 2 1 1 1 1 d u p (p p ) (u u ) 1 dt x ( ) x ( ) ∂ - ∂α λ - + = + ρ ∂ αρ ∂ αρ 1 1 1 I 1 I 1 2 1 2 1 I 1 2 1 1 1 1 1 d e p p u (u u )(u u ) H(T T ) μ p (p p ) dt x ( ) ( ) -∂ λ - - - + = - - + + ρ ∂ α ρ αρ αρ 1 1 1 I 1 I 1 2 1 2 1 I 1 1 2 1 1 1 1 1 1 1 1 1 d s p p u (u u )(u u ) H(T T ) μ (p p )(p p ) dt T x T ( ) T ( ) T -∂ λ - - - + = - - - + + ρ ∂ α ρ αρ αρ ( ) 2 2 2 1 1 1 2 2 1 2 2 2 1 2 2 2 2 d u u u u μ(p p ) dt x x x ρ ρ ∂α α ∂ ∂ ρ + - + ρ + ρ = - α ∂ α ∂ ∂ α 2 2 2 I 2 2 2 1 2 2 2 d u p p p (u u ) 1 dt x ( ) x ( ) ∂ -∂α λ - + = - ρ ∂ αρ ∂ αρ ( ) I 1 2 2 2 2 2 1 I 1 1 I 1 2 I 2 2 1 2 1 2 2 2 2 2 2 2 2 2 p u u d e p u p u μp (p p ) (u u )(u u ) H(T T ) dt x x x ( ) ( ) - ∂ α ∂ ∂α - λ - - - + + + = - - ρ ∂ α ρ ∂ α ρ ∂ α ρ αρ αρ ( ) I 2 1 2 2 2 1 I 2 1 1 1 2 I 2 I 2 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 (p p ) u u d s (p p ) u μ(p p )(p p ) (u u )(u u ) H(T T ) dt T x T x T ( ) T ( ) T - - α - ∂ ∂α - - λ - - - + + = - - α ρ ∂ α ρ ∂ α ρ αρ αρ
The second equation of this system is particularly interesting. It means that phase 1 density is independent of velocity divergence. As phase 1 is dispersed, there is no reason that droplet cloud contraction or expansion make the density of that phase vary. In the present formulation, it varies only as a consequence of drop contraction or expansion, due to pressure differential. b) Stiff pressure relaxation limit Former system dramatically simplifies in the stiff pressure relaxation limit, as shown hereafter. Let us consider first-order expansions for the pressures, ( )

2 0 1 k 1 1 k k p c p p ... = ρ + ε + (IV.5)
Where, -ε is of the order of the pressure relaxation time (

1 ε ≈ τ ), tending to zero in most situations ( 0 + ε → ), - 0 k
p and 1 k p are respectively the dimensionless pressures at leading and first order of the Taylor expansion. Inserting these definitions in (III.4) it becomes,

0 1 0 1 1 1 1 1 1 1 2 2 u 3 (p p p p ) t x ∂α ∂α α + = + ε --ε ∂ ∂ ε This equation implies two relations, as ε is arbitrarily small: 0 0 1 2 p p = (IV.6) and, 1 1 1 1 1 1 2 u p p t x ∂α ∂α + = - ∂ ∂ (IV.7)
The equilibrium condition (IV.6) is valid at leading order only and is different of the strict pressure equilibrium condition,

1 2 p p = ,
(IV.8) widely used in two phase flow literature. Such strict equality results in non-hyperbolic models (see for example Guidaglia et al., 2001). Equation (IV.7) means that pressure fluctuations are still present in the flow model, each time the relaxation coefficient μ appears in factor of the pressure differential 1 2 (p p ) -. Let us for example examine the entropy equation of the first phase that becomes, after inserting (IV.5):

( ) ( ) ( )( ) 2 0 0 2 1 1 1 1 I 1 1 I 1 1 1 1 1 1 2 0 0 1 1 0 0 1 1 1 I 1 2 1 2 1 1 I 1 I 1 2 1 2 1 1 1 1 1 c p p c p p d s u u dt T x T x 3c (u u )(u u ) H(T T ) p p (p p ) p p (p p ) T ( ) T ( ) T - - ∂ ∂ + + ε = ∂ ∂ λ - - - - -+ ε - -+ ε - + + ε αρ αρ
With the help of (IV.6) simplifications appear, ( )

2 1 1 2 1 1 I 1 1 1 1 1 1 I 1 2 1 2 1 1 I 1 1 1 1 1 1 c p p d s u 3c (u u )(u u ) H(T T ) (p p )² dt T x T ( ) T ( ) T - ∂ λ - - - + ε = - ε - + + ∂ αρ αρ .
Under the assumption of smooth solutions and as 0 + ε → , it reduces to:

1 1 I 1 2 1 2 1 1 1 1 1 d s (u u )(u u ) H(T T ) dt ( ) T ( ) T λ - - - = + αρ αρ
Therefore all terms involving pressure differential vanish, except those related to first-order pressure relaxation effects (quadratic pressure differential terms vanish). The resulting limit system reads,

1 1 1 1 1 2 d u μ(p p ) dt x α ∂ + α = - ∂
(IV.9)

1 1 1 1 2 1 d μ (p p ) dt ρ ρ = - - α 1 1 1 2 1 1 1 d u p (u u ) 1 dt x ( ) ∂ λ - + = ρ ∂ αρ 1 1 I 1 2 1 2 1 I 1 2 1 1 1 1 d e (u u )(u u ) H(T T ) μ p (p p ) dt ( ) ( ) λ - - - = - - + + α ρ αρ αρ 1 1 I 1 2 1 2 1 1 1 1 1 d s (u u )(u u ) H(T T ) dt ( ) T ( ) T λ - - - = + αρ αρ ( ) 2 2 2 1 1 1 2 2 1 2 2 2 1 2 2 2 2 d u u u u μ(p p ) dt x x x ρ ρ ∂α α ∂ ∂ ρ + - + ρ + ρ = - α ∂ α ∂ ∂ α 2 2 2 2 1 2 2 d u p (u u ) 1 dt x ( ) ∂ λ - + = - ρ ∂ αρ ( ) I 1 2 2 2 2 2 1 I 1 1 I 1 2 I 2 2 1 2 1 2 2 2 2 2 2 2 2 2 p u u d e p u p u μp (p p ) (u u )(u u ) H(T T ) dt x x x ( ) ( ) - ∂ α ∂ ∂α - λ - - - + + + = - - ρ ∂ α ρ ∂ α ρ ∂ α ρ αρ αρ 2 2 I 2 2 1 2 1 2 2 2 2 d s (u u )(u u ) H(T T ) dt ( ) T ( ) T λ - - - = - - αρ αρ
With the help of interfacial variables estimates (IV.3) the entropy equations become,

1 1 2 2 1 2 1 1 1 1 2 1 1 d s Z H(T T ) (u u )² dt ( ) T Z Z ( ) T - λ = - + αρ + αρ 2 2 1 2 1 2 1 2 2 1 2 2 2 d s Z H(T T ) (u u )² dt ( ) T Z Z ( ) T - λ = - - αρ + αρ
Combination of these equations with the mass equations results in the following mixture entropy equation, that guarantees non-negative evolutions,

1 1 1 2 2 2 1 1 1 1 2 2 2 2 2 1 2 1 2 1 1 2 1 2 1 2 s s s u s u Z Z H(T T )² (u u )² t x T T Z Z T T   ∂α ρ + α ρ ∂α ρ + α ρ - λ + = + - +   ∂ ∂ +   (IV.10)
System (IV.9) is consequently entropy preserving.

It is interesting to note that, in the present limit, the internal energy equation of the first phase expresses in conservation form:

1 1 1 1 1 1 1 I 1 2 I 1 2 1 2 1 e e u μp (p p ) (u u )(u u ) H(T T ) t x ∂α ρ ∂α ρ + = - - + λ - -+ - ∂ ∂ (IV.11)
We now check hyperbolicity of the same equations.

c) Hyperbolicity System (IV.9) in absence of relaxation effects, is expressed as,

W W A(W) 0 t x ∂ ∂ + = ∂ ∂ , with, ( ) T 1 1 2 1 1 2 2 W , s , s , , u , , u = ρ α ρ . The Jacobian matrix reads, 1 2 1 1 2 1 1 2 1 1 1 1 1 1 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 u 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 u 0 0 c p 1 0 0 u 0 0 A(W) s 0 0 0 (u u ) u p c 1 0 0 0 0 u s ρ ρ           α      ∂    = ρ ρ ∂      ρ α   - ρ ρ   α α      ∂      ρ ∂ ρ    .
The wave speeds, solution of A I 0 -λ = are,

1 4 1 u - λ = , 5 2 u λ = , 6 2 2 u c λ = -and 7 2 2 u c λ = + .
(IV.12) All roots being real the system is unconditionally hyperbolic. The wave speeds correspond to the one of the dilute model of [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF] (Systems II.2 -II.3) and not those of [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF], as expected.

d) Model summary

The flow model thus consists in System (IV.1-2) with the condition:

1 0 + τ →
(IV.13) Alternatively it can be expressed as,

1 1 1 1 2 u (p p ) t x ∂α ∂α + = µ - ∂ ∂ , with µ → +∞ 1 1 ( ) ( u) 0 t x ∂ αρ ∂ αρ + = ∂ ∂ , 2 1 1 1 1 1 1 1 1 1 I 2 1 u u p p (u u ) t x x ∂α ρ ∂α ρ + α ∂α + = + λ - ∂ ∂ ∂ , 1 1 1 1 1 1 1 2 I 1 2 2 1 2 1 1 2 e u e Z p (p p ) (u u )² H(T T ) t x Z Z ∂α ρ ∂α ρ + = -µ - + λ - + - ∂ ∂ + , 2 2 ( ) ( u) 0 t x ∂ αρ ∂ αρ + = ∂ ∂ , (IV.14) 1 1 1 2 2 2 1 1 1 1 1 2 2 2 2 2 ( u u ) ( u ² p ) ( u ² p ) 0 t x ∂ α ρ + α ρ ∂ α ρ + α + α ρ + α + = ∂ ∂ , 1 1 1 2 2 2 1 1 1 1 1 2 2 2 2 2 ( E E ) u ( E p ) u ( E p ) 0 t x ∂ α ρ + α ρ ∂α ρ + + α ρ + + = ∂ ∂ .
In this formulation there is a single non-conservative equation (the momentum one of the liquid phase) as the momentum of the second phase is deduced from the mixture momentum equation. The conservative internal energy equation for the liquid phase is a consequence of the second equation of (IV.4) and stiff pressure relaxation limit. This is a nice property that may simplify shock conditions determination. Obviously, System (IV.14) can be complemented by mass transfer.

e) Stiff mechanical relaxation limit

We now address both stiff pressure and velocity relaxation limit to check compatibility of the model with the [START_REF] Kapila | Two-phase modeling of deflagrationto-detonation transition in granular materials: Reduced equations[END_REF] one. This is important for the computation of material interfaces with capturing methods. Here, only pressure and velocity relaxation processes are considered. They are considered to relax at infinite rate. The pressure evolution equations read,

2 1 1 1 1 1 1 2 1 p p c u (p p ) t x ∂ ∂ ρ + = - µ - ∂ ∂ α 2 2 2 2 2 2 1 1 2 2 2 2 2 1 2 2 2 p p c u u c u (p p ) t x x ∂ ∂ ρ ∂α + α ρ + + = µ - ∂ ∂ α ∂ α
Taking the difference,

2 2 2 1 2 2 2 1 1 2 2 1 1 2 2 1 2 1 2 2 1 2 p p c u u c c u u (p p ) x x x   ∂ ∂ ρ ∂α + α ρ ρ - - = - + µ -   ∂ ∂ α ∂ α α   .
In the stiff pressure ( 0 0 0 1 2 p p p = = ) and velocity relaxation limits (

0 0 0 1 2 u u u = = ), 2 2 2 2 1 2 2 2 1 1 2 2 1 2 c u (p p ) c c x ρ α ∂ µ - → ρ ρ ∂ + α α
Inserting this result in the volume fraction equation,

2 2 1 1 1 1 2 2 2 2 1 1 2 2 1 2 c c u u c c t x x ∂α ∂α ρ -ρ ∂ + = ρ ρ ∂ ∂ ∂ + α α
, the volume fraction equation of the [START_REF] Kapila | Two-phase modeling of deflagrationto-detonation transition in granular materials: Reduced equations[END_REF] model is recovered.

The mixture sound speed at mechanical equilibrium is thus that of [START_REF] Wood | A Textbook of Sound[END_REF],

1 2 2 2 2 1 1 2 2 1 c c c α α = + ρ ρ ρ
while for the Marble model with stiff velocity relaxation, the mechanical equilibrium sound speed is

2 2 2 2 c c   ρ =   ρ  
, which is very different.

As the Kapila model is recovered in the stiff mechanical relaxation limit, it means that the present flow model is able to compute interfacial flows with the help of stiff velocity and pressure relaxation solvers. This feature is particularly important for the sake of generality of the formulation. The present model is hyperbolic but not symmetric, as sound propagates only with the second phase. It is therefore interesting to compute relevant test problems to examine typical solutions. To do this, an appropriate flow solver is derived in Appendix A.

V. Computed results

Several test problems are addressed, some giving relevant illustrations of model's capabilities, other serving for validation as compared to exact solutions as well as experimental data.

Shock tube tests

The first test corresponds to the simple transport of a volume fraction discontinuity in a flow field in uniform pressure and velocity equilibrium. The method of Appendix A is extended to higher order thanks to the MUSCL algorithm (see for example Toro, 1997). Present computations use the Minmod limiter. The ideal gas equation of state is used to model thermodynamics of the gas phase, while the liquid is modeled by the stiffened gas EOS. These two EOS can be summarized as follows, ( ) Another test is now addressed and corresponds to a two phase shock tube, as shown in Figure 2. In this test, stiff pressure relaxation is not used, while the flow model is hyperbolic and entropy preserving only when stiff pressure relaxation is done. Consequently, the present computations should fail. They however produce results with quite significant pressure disequilibrium, showing robustness of the formulation. (compared to the velocity of the gas phase) even in absence of drag, the pressure term in the momentum equation being responsible for that.

k k k k k k k k p ( ,
These typical profiles are very different to those expected with BN and Marble models (without relaxation terms). For example, the pressure profiles of the phases are very different and even unphysical, but the model and algorithm do not fail.

To recover acceptable pressure evolutions in the condensed phase, stiff pressure relaxation is used. Pressure and velocity relaxation solvers are recalled in Appendix B. The same run as the one defined in Figure 2 is reconsidered hereafter with pressure relaxation and the results are shown in Figure 3. With the help of both velocity and pressure relaxation solvers it is possible to address interfaces separating (nearly) pure liquid and (nearly) pure gas. The aim is to analyze the behavior of the flow model in another limit case, having in mind it has been derived for clouds of droplets, not interfaces. The initial conditions are given in Figure 4 and correspond to a liquid at right set to motion by a pressurized gas at left. The exact solution is available for this test case and used to check accuracy of computations. These various results illustrate method's capabilities in several configurations far from the original design of the model, based on clouds of droplets suspended in air. It is also interesting to note that even if the liquid acoustic wave is absent of the formulation, its dynamics is correctly computed thanks to relaxation effects.

Shock interaction with a fluidized bed -Rogue test

We now consider a test more appropriate to the model. It consists in the fluidization of a particle cloud under shock wave interaction. Such a configuration has been studied experimentally by Rogue et al.

(1998) and is summarized in Figure 6. ) is set in a cross section of the tube, with 2 cm width. The initial solid volume fraction in the particle bed is 0.65. The initial pressure is uniform initially and set at 10 5 Pa. A shock at Mach number 1.3 is created by the expansion of the high pressure gas, equivalent to a shock created by a piston moving at 151 m/s.

In this experiment pressure signals are recorded before and after the particles cloud, to examine reflected and transmitted waves through the granular media as well as its dilution and dispersion. To account for drag effects the following correlation is used, combination of [START_REF] Ergun | Fluid flow through packed columns[END_REF] and [START_REF] Bernecker | Studies in the transition from deflagration to detonation in granular explosives-II. Transitional characteristics and mechanisms observed in 91/9 RDX/Wax[END_REF],

F = ρ 2 d 1 C d u 2 -u 1 u 2 -u 1 ( ) with, C d = 150α 1 Re + 1.75 if α 1 ≥ α cr 150α 1 Re + 1.75 (1-α cr )α 1 α 1 α 2       0.45 if (1-α s ) ≤ α 1 ≤ α cr 150α 1 Re + 0.3 if α 1 ≤ (1-α s )            , 2 2 2 1 1 e 2 u u d R α ρ - = µ
the particulate Reynolds number, cr 0.63 α = and

1 cr s cr 1 1 0.02 -   -α α = +   α   .
The particle diameter appearing in these relations is constant ( 1 d 1.5mm = ) and the gas viscosity is At this stage, more potential than expected appeared with the various shock tube tests in limit configurations, but weakness appeared for particle cloud dynamics in the Rogue test.

A symmetric variant of the new model, with Equation (III.5) instead of (III.4) is thus considered and examined.

VI. Symmetric variant

In the model examined and tested formerly pressure waves propagate with the gas sound speed. This behavior seems incorrect if the flow is mainly liquid, which is the continuous phase with bubbly flows. We thus address the symmetric variant of the previous model on the basis of the volume fraction equation (III.5). Equation (III.5) is thus plugged to the system of balance equations (IV.2). Analysis of the resulting system yields the following wave's speeds: We examine typical solutions on some test problems, as done previously. The same shock tube test case as the one of Figure 5 is considered with liquid at left and gas at right. Computed results are compared against the exact solution in Figure 9 at the same time 250µs. This test is however interesting as the wave dynamics in the gas phase is correctly computed, while it is absent of the non-equilibrium flow model.

It is also interesting to address the Rogue test problem with the symmetric variant model.

Corresponding results are shown in Figure 10, showing improvement of the pressure evolution in the cloud dynamics, but poor wave propagation in the gas phase. This test shows the limits of the pressure relaxation method, to rebuilt an acoustic wave absent of the original hyperbolic system. Accuracy has been lost in the incident and reflected waves compared to the former model while an interesting feature appears regarding the beginning of the pressure rise during particle cloud motion that seems more accurate.

The lack of accuracy in the incident shock that propagates in a single phase gas media is due to the absence of gas sound speed in the eigenvalues. Also, the present formulation is derived for bubbly flows while it is used for a nearly pure gas in the first part of the shock tube. Therefore an attempt for a general formulation is done in the section that follows.

VII. Towards a general formulation

The new model and its symmetric variant are embedded in a general formulation. Parameters 'a' and 'b' are defined as,

fluidization 1 1 if a 0 otherwize  α < α =   (VII.1) b a 1 = -
Parameter fluidization α corresponds to some fluidization limit, for example fluidization 0.5 α ≈ . This parameter has been used in the various computations that will be examined later. The various tests done haven't shown clear dependence to this parameter. Therefore, in this formulation 'a' and 'b' are local constants, but as they vary in space as The general flow model reads,

1 1 1 2 2 1 2 u u a b (p p ) t x x ∂α ∂α ∂α + + = µ - ∂ ∂ ∂ , (VII.2) 1 1 ( ) ( u) 0 t x ∂ αρ ∂ αρ + = ∂ ∂ 2 1 1 1 1 1 1 1 1 1 I 2 1 u u p p (u u ) t x x ∂α ρ ∂α ρ + α ∂α + = + λ - ∂ ∂ ∂ 1 1 1 1 1 1 1 1 1 1 1 1 2 2 I I 1 2 I 2 1 2 1 E u E u p u u p a b p (p p ) u (u u ) H(T T ) t x x x x ∂α ρ ∂α ρ ∂α ∂α ∂α   + + = + -µ - + λ - + -   ∂ ∂ ∂ ∂ ∂   2 2 ( ) ( u) 0 t x ∂ αρ ∂ αρ + = ∂ ∂ 2 2 2 2 2 2 2 2 2 2 I 2 1 u u p p (u u ) t x x ∂α ρ ∂α ρ + α ∂α + = -λ - ∂ ∂ ∂ 2 2 2 2 2 2 2 2 2 2 1 1 2 2 I I 1 2 I 2 1 2 1 E u E u p u u p a b p (p p ) u (u u ) H(T T ) t x x x x ∂α ρ ∂α ρ ∂α ∂α ∂α   + + = - + + µ - -λ -- -   ∂ ∂ ∂ ∂ ∂ 
 It admits the following additional mixture entropy equation,

2 2 2 1 1 1 2 2 2 1 1 1 1 2 2 2 2 2 1 2 1 2 1 1 2 1 2 1 2 1 2 1 2 s s s u s u (u u ) Z Z H(T T ) 1 1 (p p ) t x T T Z Z T T T T     ∂α ρ + α ρ ∂α ρ + α ρ - - + = µ - + + λ + +     ∂ ∂ +    
guaranteeing its thermodynamic consistency. System (VIII.2) can also be written as,

W W A(W) 0 t x ∂ ∂ + = ∂ ∂ , with, ( ) T 1 2 1 1 1 2 2 W s , s , , , u , , u = α ρ ρ , and, (
)

1 2 1 2 2 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 2 2 1 2 2 1 2 2 2 2 2 2 2 2 2 2 2 u 0 0 0 0 0 0 0 u 0 0 0 0 0 c 1 p 0 u 0 0 0 s 0 0 b u b u u 0 b A(W) 0 0 a 0 ( au b u ) 0 b 0 0 a 0 a ( u u ) u a c 1 p 0 0 0 0 u s ρ ρ          ∂    ρ ∂ ρ      ρ α   -ρ - -ρ   α α =   α - α     α ρ ρ - ρ   α α      ∂      ρ ∂ ρ  
 This matrix has a nice structure. Eigenvalues are given by det(A I) 0 -λ = , which results in the following polynomial,

( ) ( ) ( ) 2 2 2 1 1 2 2 1 2 1 2 1 (u )²(u )² ( au b u ) (u )²bc ( a b) u (u )²ac (b a)u 0   -λ -λ - -λ + -λ - -λ + -λ λ + - =   When a=1, b=0 it reduces to, 3 2 1 2 2 (u ) (u )² c 0   -λ -λ - = 
 with the wave speeds of the first model. When a=0, b=-1 it reduces to, ( )

3 2 2 1 1 u (u )² c 0   -λ -λ - =  
with the wave speeds of the symmetric model. Therefore the flow model (VII.2) has the following wave speeds: For numerical computations, the first equation of System (VII.2) is expressed as, ( ) ( )

1 1 u λ = , 2 1 1 u c λ = + , 3 1 1 u c λ = -, 4 2 u λ = , 5 2
1 1 2 2 1 1 1 2 2 a u b u a u u 0 t x x ∂ α + α ∂α ∂ + -α + α = ∂ ∂ ∂ ,
It is solved with the hyperbolic solver of Appendix A based on the Rusanov flux. However, the wave speed estimate, ( )

k k k i 1 i S Max , + = λ λ ,
now involves the six eigenvalues (VII.3).

The volume fraction equation in (VII.2) being non-conservative, appropriate scheme is needed. Similar analysis as the one described in Section V is reused. Details are given in Appendix C. Let us examine typical solutions of the general model on some test problems, as those considered previously. A shock tube test case with liquid at left and gas at right is considered, in the same conditions as the tests in Figures 5 and9. Computed results are compared against exact solution in Figure 11. Results are shown at time 250µs. The numerical solution tends to the exact one and converges quite well. Compared to the original model (Figure 5) and its symmetric variant (Figure 9) improvements are visible on both velocity and density profiles. This last test shows improvements compared to Marble and BN models results:

-reflected and transmitted waves have better accuracy, -pressure rise in the particle cloud shows better accuracy than existing models. Possibly better agreement with experimental data could be obtained by using sophisticated drag force correlation, but this is not the scope of the present work. Even if the BN model and the new one have the same limit model (Kapila's model) in the case of infinite pressure and velocity relaxation, the transient wave dynamics between the pure gaseous shock towards a fan of compression waves has to be well captured and seems to depend on the topology of the two-phase flow. Moreover, it can be expected that the dispersive nature of a compression wave from the BN and the new models is different. Intergranular stress [START_REF] Bdzil | Two-phase modeling of deflagrationto-detonation transition in granular materials: A critical examination of modeling issues[END_REF][START_REF] Saurel | Modelling dynamic and irreversible powder compaction[END_REF] has been considered as a possible effect to improve the computations of Figure 12. These effects have been added to the present formulation and coded, but no noticeable improvement appeared. We have also investigated the effects of the fluidization limit switch fluidization α , taken previously to 0.5. Various estimates have been tested, from 0.1 to 0.9 without noticeable changes to the results.

VIII. Conclusion

A new two-phase hyperbolic and thermodynamically consistent model has been built and typical solutions have been computed. It is able to compute the same flow configurations as the BN model, i.e. interfaces separating pure fluids and non-equilibrium multiphase mixtures. Its acoustic properties sound physical. Moreover, the evolution of the two-phase topology directly influences both number and speed of waves present at a given point of space. The flow dynamics expressed by the model is not only reflected through the change of drag coefficient between phases, through interfacial area evolution. It is expected that two-phase shock waves structure be easier to analyze in the present frame. It is also expected that multidimensional solutions exhibit more differences than present one-dimensional computations, in particular regarding interface instabilities.

   = - -     ∂ ∂ ∂ ∂ ∂     .
The difficulty with this hyperbolic system relies in the non-conservative term U H U, x

∂     ∂   .
For the sake of simplicity the [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves and obstacles[END_REF] approximate Riemann solver is considered. It uses the following estimate for the right facing wave, at a given cell boundary separating cells i and i+1:

( )

k k k i 1 i S Max , + = λ λ .
At a given cell boundary separating left (L) and right (R) states, the approximate flux reads, ( )

* R L R L 1 F F F S U U 2 = + - -     (A.
2) The Godunov scheme for System (A.1) necessarily reads, (  )

n 1 n * * i i i 1/2 i 1/2 i t U U F F tH x + + - ∆ = - - + ∆ ∆ (A.3)
where i H is the numerical approximation of U H U, x

∂     ∂  
, to be determined.

To determine i H , we follow [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] where a flow in uniform mechanical equilibrium is considered: In uniform velocity flows conditions, using the Rusanov flux (A.2) in the Godunov method (A.3), the discrete mass equation results in (A.4) Let us now consider the particular case of uniform density field:

n n n 1,i 1,i 1 1,i 1 + - ρ = ρ = ρ .
Both velocity and density being uniform, the density at the next time step must be invariant:

n 1 n 1,i 1,i + ρ = ρ .
In this context, the mass equation becomes,

[ ] [ ] n 1 n 1 1 i 1 1 i 1 i 1 1 i 1 1 i 1 1 i 1 i 1 u t S t ( ) ( ) ( ) ( ) ( ) 2( ) ( ) 2 x 2 x + + - + - ∆ ∆   α ρ = ρ α - α -α + α -α + α   ∆ ∆   (C.

  state and associated parameters are used in all test problems of the present section. A volume fraction discontinuity separating two mixtures is transported at 100 m/s in a uniform pressure flow field of 0.1 MPa. The initial discontinuity is located at 0.5 m initially and computed results are compared to the exact solution at time 1 ms in the Figure1.

Figure 1 .

 1 Figure 1. Volume fraction transport in uniform pressure and velocity fields. The mesh involves 500 cells and the time step is computed with CFL=0.5. Initial velocities are set to 100m/s and pressures are constant and equal to 10 5 Pa. The volume fraction discontinuity is initially set at 0.5 m. The numerical solution is plotted at 1ms. The exact solution for the volume fraction is presented in dot symbols showing perfect agreement. The numerical solutions are oscillations free.

Figure 2 .

 2 Figure 2. "Smooth shock tube test case". Computations are made in the absence of relaxation terms, with 500 cells and CFL=0.5. Results are shown at time 350µs. Four waves are visible, in spite of the simplified Riemann solver that considers two only. It is interesting to note the discontinuous profile of pressure in the phase 1: no pressure wave is present in this phase. Phase 1 shows slight velocity creation

Figure 3 .

 3 Figure 3. Smooth shock tube computations in the absence of velocity relaxation but with stiff pressure relaxation. Computations made with 500 cells and CFL=0.5. Computed results are shown at time 350µs.All pressures are now equal, modifying significantly the phase 1 velocity profile.

Figure 4 .Figure 5 .

 45 Figure 4. Shock tube with gas-liquid interface: High pressure gas at left and low pressure liquid at right. Computations done with 500 cells and CFL=0.5. Computed results shown at time 150µs. Both velocities and pressure are relaxed, making the interface condition of equal pressures and velocities fulfilled.The same test is considered but with fluids in reverse order: High pressure liquid at left and low pressure gas at right. This test is more severe as maintaining pressure positivity during liquid expansion is hard to manage. Corresponding results are shown in Figure5and compared to the exact solution.

Figure 6 .

 6 Figure 6. Rogue et al. (1998) fluidization shock tube test. A shock tube is filled with gas at density 1.2 kg/m 3 . A dense cloud of nylon particles ( 3 0 1050kg / m ρ =

  is made of nylon particles, treated as compressible material, governed by the SG EOS (V.1) with the same EOS parameters as liquid water. Predictions of the BN, Marble and new model are compared with the same modeling of drag effects given above. The BN model, or more precisely its symmetric version with 7 waves is solved with the Furfaro and Saurel (2015) method. The Marble model is solved with the Saurel et al. (1994) method and the new model is solved with the method presented above. Let us first comment Rogue et al. (1998) experimental data, typical pressure signals being shown in Figure 7.

Figure 7 .Figure 8 .

 78 Figure 7. Experimental pressures signals of Rogue (1998): 1 denotes the incident shock wave, 2 denotes the transmitted shock wave/fan of compression waves, 3 denotes the reflected shock wave on the particles cloud, 4 corresponds to the arrival of the cloud upper front at the pressure gauge location. Computed results with the various flow models are shown in Figure 8. (Marble) (BN)

Figure 9 .

 9 Figure 9. Liquid -gas shock tube test problem solved with the symmetric variant model with stiff pressure and velocity relaxation. Computations are done with 500 cells and CFL=0.5. Results are shown at time 250µs. The numerical solution tends to the exact one but converges more difficultly, as visible on the density graph.

Figure 10 .

 10 Figure 10. Computed pressure profiles of the Rogue test with the symmetric variant model (lines) are compared against pressure records (symbols). Computations are done with 1000 cells and CFL=0.5. Accuracy has been lost in the incident and reflected waves compared to the former model while an interesting feature appears regarding the beginning of the pressure rise during particle cloud motion that seems more accurate.

  wave speeds are not present at any point of space. They change when the volume fraction crosses the fluidization limit ( fluidization α ) somewhere in the domain. In nearly all computational examples considered previously, such instance happens.Let us mention that other guesses have been considered for parameters 'a' and 'b'. For example 1

Figure 11 .

 11 Figure 11. Liquid -gas shock tube test solved with the general model with stiff pressure and velocity relaxation. Computations are done with 500 cells, CFL=0.5 and van Leer limiter in the MUSCL method.Results are shown at time 250µs. The numerical solution tends to the exact one and converges quite well. Compared to the original model (Figure5) and its symmetric variant (Figure9) improvements are visible on both velocity and density profiles.

  It is also interesting to address the Rogue test problem with the general model. Corresponding results are shown in Figure12.

Figure 12 .

 12 Figure 12. Computed pressure signal of the Rogue test with the general model are compared with experimental pressure records. Computations are done with 2000 cells and CFL=0.5. The various transmitted and reflected waves are correctly computed. The pressure evolution in the particle cloud has been improved but is still perfectible.

  energy equation of phase 1 with the corresponding volume fraction equation results in, the specific volume of the considered phase. Considering the relaxed pressure as a constant during time integration, it becomes (this assumption has been analyzed in[START_REF] Saurel | A relaxation-projection method for compressible flows. Part I: The numerical equation of state for the Euler equations[END_REF], equation of state (alternatively given by Equation V.1) is inserted (any other convex EOS can be considered as well), discretization compatible with the mass equation of the same system,

  . The same is done for the volume fraction equation,

	[ is the numerical approximation of n 1 i 1 2 i 1 1 u t ( ) (a b ) (a 2 x + ∆ = α -α + α -α + α b ) 2 i 1 -∆ a x ∆ ∆ a 0 x ∂ = ∂ Rearranging (C.1) with b a 1 n 1 1 i ( ) + α where = -and 2 1 1 α = -α the discrete volume fraction equation becomes, ] [ ] 1 i 1 1 i 1 i 1 S t a ( ) 2( ) ( ) u t 2 x x + -∆ ∆ + α -α + α + ∆ ∆ ∆ , to be determined.	(C.1)
	n 1 1 i ( ) + α	n 1 i ( ) = α -	u t 2 x ∆ ∆	[	1i 1 + α -α 1i 1 -	]	-	u t 2 x ∆ ∆	[	a	i 1 +	-	a	i 1 -	]	+	S t 2 x ∆ ∆	[	1 i 1 ( ) + α	1 i 2( ) ( ) 1 i 1 --α + α	]	u t + ∆	a x ∆ ∆	(C.2)

1,i 1 1,i 1,i 1 2,i 1 2,i 2,i 1 u

Acknowledgements. The authors are grateful to anonymous referees that helped to improve significantly the quality of the manuscript.

Appendix A. Hyperbolic solver

We address derivation of a Godunov type method for System (IV.1-2). In the absence of source terms, it expresses in compact form as,

where, ( )

)

Inserting the Rusanov flux (A.2) in the Godunov method (A.3) for the mass equation of the first phase, the following result is obtained:

4) The same procedure is done for the momentum equation of the same phase:

H must be approximated as,

Considering the balance energy equation of the same phase the following discrete approximation is obtained:

The same analysis as before to maintain mechanical equilibrium provides the following guess:

The flow solver thus consists in (A.3) with (A.2), (A.5) and (A.6).

Appendix B. Stiff pressure and velocity relaxation solvers a) Stiff pressure relaxation

The system to consider for phase 1 during pressure relaxation is the following:

2) Same result is obtained for the second phase,

The saturation constraint

∑ is then considered resulting in the following root for the equilibrium pressure:

where, ( )

and ( )

Once the relaxed pressure is determined with (B.3) the volume fractions at equilibrium are determined with (B.2). b) Stiff velocities relaxation During stiff velocity relaxation, the subsystem to consider reads,

u is defined by (IV.3). As λ → +∞ both velocities relax to the equilibrium one given by the mixture momentum equation,

The velocities of the phases are then reset to this value. The internal energy of the phases are corrected as,


Dr Niki Baccile

Dr Markus Antonietti

Maria-Magdalena Titirici
email: magdalena.titirici@mpikg.mpg.de

Dr Maria-Magdalena Titirici-Corresponding Author

One-step Hydrothermal Synthesis of Nitrogen-Doped Nanocarbons: Albumine Directing the Carbonization of Glucose

Keywords: carbon nanoparticles, hydrothermal carbonization, green chemistry, nanostructured carbon ))

In this work, we present a simple and green one-step pathway towards nitrogen-doped carbon nanostructures with controlled mesoporosity via hydrothermal treatment of glucose in the presence of model proteins. Thus, performing the reaction with different amounts of egg white ovalbumin protein (OvA), carbonaceous nanoparticles or continuous nanosponges with high specific surface areas can be efficiently produced. The nitrogen content of the structures is rather high (up to 8 wt%) and can be kept constant up to 950°C, while oxygen elimination and graphitization of the carbon material occurs. We demonstrate here that sustainable natural resources can be efficiently used in the synthesis of pure highpotential nanomaterials.

Introduction

Carbon-based materials with large surface area or small particle size are made along a number of different synthesis pathways. [1][2][3][4] Activated carbons are generally synthesized from raw biomass but the conditions are usually very harsh, involving high temperatures (600-900°C) and pressures under oxidizing atmospheres (O2, CO2, steam) or chemical treatment with acids (e.g., phosphoric acid), bases (e.g., KOH, NaOH) or salts (ZnCl2). Recently, much softer condensation routes to obtain periodic porous carbon networks [5][6][7][8] were presented. From a cost, processing and sustainability perspective, there is however still space for improvements, as indirect impregnation techniques with mesoporous silica as a template, hydrogen fluoride etching, [5,8] and finally high carbonization temperatures [5] are required. Direct templating of resins (phenol/formaldehyde) [6] is presumably most promising, but restricted to a distinct monomer polarity and functionality and excludes the use of natural products as carbon source.

Energy and sustainability issues are even less well covered in case of nanosized forms of carbon, ranging from nanoparticulate soots and fullerenes to nanotubes. Their synthesis is performed using energy-consuming methods like arc discharge, laser ablation or Chemical Vapor Deposition (CVD), [3] whereas printing soots are still made by incomplete combustion of natural gas, being an element and energy inefficient process. [9,10] Recently, an electrochemical method was developed to produce carbon nano colloids, [11] and carbon nanoparticles with interesting fluorescence [12] and anti-microbial [13] properties from candle soot were presented, underlining the increasing interest towards more sustainable processing approaches to carbon nanostructures.

According to our opinion and excluding activated carbons, only little research has been done to synthesize and recognize the structure of carbon materials based on natural resources. This is hard to understand, as carbon structure synthesis was done from the beginning of civilization on the base of biomass, with the petrochemical age being a late deviation, only. A modern approach towards carbon synthesis based on renewable resources is a significant operation, as the final products represent an important share of materials of the next generation. They already have found applications in different domains ranging from environmental science [14][15][16] to drug delivery [17][18][19] and energy storage [20] depending on their structural, morphological and chemical properties. [21] Synthesizing materials using sustainable approaches is definitely not new but currently revamped for the field of nanomaterials, as shown by recent review papers. [22][23][24] Another issue which has to be addressed in a novel carbon synthesis is functionalization, which is required in most specific applications. [25][26][27][START_REF] Jurewicz | Ammoxidation of brown coals for supercapacitors[END_REF][START_REF] Hsieh | [END_REF] For instance, fuel cell development has turned its attention towards nitrogen-doped carbon nanostructures, as these materials seem to enhance the catalytic activity towards oxygen reduction reactions in fuel cell electrodes. [START_REF] Shao | [END_REF] In addition, N-doped carbons are also interesting as CO2 binders [START_REF] Plaza | Pis in CO2 capture by adsorption with nitrogen enriched carbons[END_REF] or pHresponsive adsorbents. [START_REF] Jia | [END_REF] The problem of the sustainable synthesis in functional carbons was recently revisited by several research groups, [START_REF] Wang | Novel spherical microporous carbon as anode material for Li-ion batteries[END_REF][START_REF] Wang | [END_REF][START_REF] Titirici | [END_REF][START_REF] Titirici | [END_REF][37][38][39] where hydrothermal treatment of biomass in water under relatively mild conditions and, in some cases, in presence of additives provided mesoporous or nanostructured carbon materials. [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55] This technique is known for many years, [56] but the forthcoming raw material change made a re-discovery of this technique timely. In addition, the implementation of a low-cost pathway to recycle byproducts of farmed biomass would represent a way to sequester significant amounts of CO2 [37] in carbon negative products, creating a materials benefit at the same time.

The use of hydrothermal synthesis between 180°C and 220°C allowed to obtain carbon-based powders, nanofibers [57,58] , or sponge-like mesoporous carbons being potentially useful as soil conditioner, ion exchange resins or sorption coals. [START_REF] Titirici | [END_REF] The synthesis proved to be feasible both for simple systems as glucose [38] and for raw biomass side products, such as oak leafs and orange peals. [START_REF] Titirici | [END_REF] In the first case, surface chemistry could also be modified by mean of hydrophilic or hydrophobic postfunctionalization. [59] Along similar lines, the group of Clark showed that a slightly different approach using expanded starch treated with sulfuric acid could provide functional materials with disordered mesoporosity which proved to be satisfactory catalysts in the esterification of succinic acid. [43,60] In these experiments, it was found that some crude biomass mixtures created carbons with exciting nanostructures, [START_REF] Titirici | [END_REF] however bound to the specific biomass with its potentially regional and seasonal fluctuations. It was speculated that this development of nanostructure was due to natural surface stabilizers contained in the biomass, such as tannins and lipids, but also the presence of small amounts of amines and proteins seemed to have a crucial influence. This is why we study in this paper the carbon nanostructures obtained upon hydrothermal carbonization of a mixture of a simple and accessible model protein, ovalbumin (OvA, chicken egg white) and pure glucose as a model carbohydrate. At the same time, we will answer the question whether we are able to introduce the highly attractive nitrogen functionality into the condensed carbon black just by using higher protein loads and co-condensation (formation of "nitrogen blacks"). It will be shown that it is possible to create carbon nanoparticles and inter grown nanostructures with structural sizes between 10 and 50 nm rich in nitrogen just by reacting glucose with ovalbumin as co-reactant under hydrothermal conditions.

Results and Discussion

Characterization of the Hydrothermal Carbons (HC)

Figure 1 shows SEM micrographs of the final carbon materials obtained after hydrothermal treatment of glucose in presence of different concentrations of OvA. It can be observed that all materials exhibit a continuous nanometer-scale sponge-like network structure made up from primary particles ranging between 10 -40 nm. Not all the formed droplets are spherical, but cylindrical tectons are also present, indicating a rather low interface tension in the liquid phases of the reaction. TEM images shown in Figure 2 at lower (a, b) and higher magnification (c, d) nicely confirm SEM data. In sample HC-Alb2, the primary network constituting particles are found to be in the nano-range, between 20 and 50 nm, while HC-Alb1 (the sample with the lowest OvA concentration) shows the existence of more interconnected filaments or particles whose size is bigger and in the order of 70-100 nm. It is indeed remarkable that with increasing OvA content, the carbon structure changes from a micron-sized carbon dispersion to a much finer, cohesive spongelike structure. We assume that OvA (through its reaction and condensation products) indeed stabilizes the surface of much smaller nanoparticles, which throughout the course of reaction turned out to be unstable and to aggregate towards linked networks, however without pronounced coalescence. Aggregation is a common phenomenon in carbon colloids, [61] as the mutual particle-particle interactions are quite strong, and is observed in all our samples as well. Nitrogen adsorption-desorption isotherms allow to quantify the character of the interstitial porous system and were performed for all HC-AlbX bulk materials (Figure 3 a). The presence of both the surface areas and hysteresis loop at P/P0 values lower than 1 for most samples indicate the clear existence of a mesoporous network, whose origin can only come from the inter-particle architecture. Surface areas ranging between 30 and 110 m 2 /g (Table 1) as well as the mesopore sizes between 10 and 40 nm are in good agreement with the morphologies depicted in the electron micrographs presented before. We can exclude intraparticle porosity for the as-synthesized hydrothermal samples, as this is not observed by TEM and nitrogen adsorption experiments. Such a value based on sole mesoporosity is actually among the highest reported in literature for as-synthesized hydrothermal carbons (HC). The BJH pore size distribution for selected samples in the middle of the composition range (Figure 3 b, open circles) is relatively sharp and centered at about 20 nm. The calculated BET constant values (Table 1) for low-content OvA samples are roughly twice the value of the samples obtained using a high OvA amount, suggesting that the surface properties (in terms of type and distribution of functional groups) do significantly change with increasing the OvA content. It is generally accepted that the BET constant reflects the adsorption strength of the adsorbent, and that high values generally mean higher adsorption equilibrium. The reproducibility regarding the synthesis of this material seems to be quite good, as several synthetic repeats of the sample HC-Alb2 always provided, within experimental accuracy, the same characteristics (preservation of the N2 adsorption-desorption behavior given in Figure 1 of Electronic Supplementary Information, ESI).

Surface functionalities present in the OvA-modified carbons were qualitatively examined using FT-IR and XPS (Figure S2a in ESI).

The FT-IR spectra of HC prepared in presence of OvA show adsorption bands corresponding to carbonyl functions (peak at 1699 cm -1-) together with bands corresponding to N-containing groups (1580 cm -1 ) such as amines. [START_REF] Jia | [END_REF] Resonances in the 750-690 cm -1 region can be related to extended C=C bonds together with C-O-C bands. (1020 cm -1 ). [39] However, even if ambiguity in interpretating peaks in the 1580-1600 cm -1 region exists, clear insights on carbonyl functions and C=C groups will be shown and discussed later by mean of 13 C solid state NMR experiments (Figure 6) while indications of nitrogen-containing functionalities will be confirmed by N1s peaks found by XPS (ESI, Figure S2b). All in all, the resulting OvA-modified materials contain structural amino groups in addition to the typical functional groups found in hydrothermal carbons. [38] Besides generating carbon materials with pure mesoporosity at lower temperatures and under sustainable conditions, it is of course also possible to apply temperature treatments for further strengthening of the carbonaceous scaffold. To identify optimal conditions, thermogravimetric (TGA) measurements were performed. The TGA behavior under nitrogen gas of HC-AlbX samples as compared to pure OvA and hydrothermal glucose carbon reveals a number of features (Figure S3 of ESI)). As expected, the pure protein can carbonize without residue and looses 70% of its mass already at about 400°C, while the carbonaceous materials containing OvA show a much higher thermal resistance. At 400°C, they loose only about 25% of their mass, while at 1000°C they still keep up to 40%. Interestingly, hydrothermal carbons synthesized in the presence of OvA are more temperature resistant compared to pure hydrothermal carbon (HC), showing the network strengthening effect of nitrogen on the carbon scaffold. Elemental analysis (EA) (Table 1) shows that the nitrogen content increases from 3 w% to 8 w% with increasing OvA content. TGA and EA therefore clearly suggest that fractions of the amino acid backbone are actually part of the carbonaceous scaffold, as otherwise one would expect an intermediary decomposition behavior between OVA and hydrothermal carbon, but not improved stability. The point of nitrogen incorporation will be discussed in more detail later.

Characterization of the Carbonized Materials

In order to obtain more insights concerning the thermal behavior and increase in carbonization by eliminating oxygen functions, series of samples were further condensed at different temperatures, from 350 to 950°C under an inert gas (N2) atmosphere. Some of the SEM images for selected samples HC-AlbX-Y (X=1, 2, 5; Y= 350, 550, 750°C) can be found in the ESI (Figure S4). The microstructure of all albumine-modified materials seems to be stable against the heat treatment; in particular, no pronounced sintering or coalescence occurs (as it might be expected from tars without structural cross-linking) ,and the overall morphology on larger length scales appears very similar to the primary products.

TEM images in Figure 4 which can provide more insights to changes in the lower nanometer range, indeed show that the mass loss and the coupled creation of additional porosity is taken up on the scale of the carbon scaffold: Samples tempered at 350°C create an exciting and previously unknown "carbon foam" morphology, while further mass loss opens up transport pores between the blobs while all essential constituting elements are still kept. Nitrogen adsorption-desorption experiments for calcined materials were performed on all materials. An exemplary set of data is shown in Figure 5 a-b for sample HC-Alb2 treated at different calcination temperatures. The isotherms show a clear increase in gas adsorption at low relative pressures indicating an increase in surface area, while the microporosity (characterized by the apparent intercept with the Y-axis) is staying remarkably low. The surface areas show a maximum (359 m 2 /g) for the samples calcined at around 550°C, while values for samples treated at higher temperatures are lower. We attribute this feature to the onset of tectonic carbon rearrangements (promoted by nitrogen moieties) which lower the exposed surface of the scaffold structure. This behavior is seen for practical all the compositions, as quantified in Table 1. In good agreement with the TGA behavior and in spite of the "ripening" of the carbon foam structure seen in TEM, treatment at 350°C does not affect the adsorption properties and the overall C/N ratio too much: all N2 adsorption-desorption isotherms appear very similar and are generally within a 10% variation, exception made for the sample with very low OvA content, HC-Alb1 and HC-Alb1-350, where variation reaches 25% (Table 1).

FT-IR experiments (Figure 5, ESI) indicate that treatment at 350°C results in a partial loss of carbonyl groups, as indicated by the decrease of the correspondent peak bound to unsaturated carbons at 1650 cm -1 . [62][63][64] Carbonyl functions bound to saturated carbons (peak at 1699 cm -1 ) disappear only after treatment at 550°C (Figure 5, ESI) while N-containing groups (1580 cm -1 ) [START_REF] Jia | [END_REF] are present up to the highest temperature.

Structural changes with temperature were also investigated using 13 C CP-MAS solid state NMR (Figure 6). Using solid state 13 C CP-MAS NMR, the local structure of hydrothermal carbons obtained from glucose was recently revealed, [START_REF] Titirici | [END_REF] and it was found out that hydroxymethylfurfural (the main decomposition product of glucose) is the main monomer of carbon formation, involving polymerization towards polyfuranes and formation of a conjugated system, while still maintaining a high number of oxygenated functionalities bound to the furan rings. When HC is prepared in presence of OvA, its corresponding 13 C spectrum (HC-Alb2 in Figure 6) shows very similar features in the carbonyl and aliphatic regions as observed for pure HC from glucose, but, interestingly, the amount of delocalized aromatic sp 2 -carbons is much higher (peak at 130 ppm), while the amount of furane moieties is significantly decreased (150 and 110 ppm). This clearly suggests that the existence of reactive nitrogen promotes the formation of a graphitic like structure where N is presumably incorporated in form of C=C-N bonds. This information is strongly supported by XPS analysis shown discussed below. Upon further heat treatment at 550°C, the peak at 130 ppm dominates the spectrum, which is typical for graphitic structures. This clearly indicates that nanoparticles with an aromatic core are formed, while most oxygenated functions have been removed. Besides the structural control of the carbons towards a nanosponge architecture, thermally stable functionalization of the carbonaceous scaffold with nitrogen using OvA as a nitrogen source was a second task being successfully accomplished. Elemental analysis (Table 1) shows that the carbon content increases with temperature from about 67% for bare hydrothermal carbons to 84% while the nitrogen content stays almost constant at relatively high values up to T= 750°C, while the T= 950°C data show only a slight decline. This is a proof that the nitrogen is included within the carbon scaffold in the form of (highly beneficial) aromatic nitrogen units thus contributing to both the electronic and chemical behavior of the material. XPS analysis (ESI, Figure 2 b,c) confirms such an interpretation. The HC-Alb2 material shows the co-existence of amines, pyrrolic and quaternary nitrogen groups while, after treatment at 550°C, a large component at 398.7 eV appears while a shift from 400.1 to 400.7 eV [65] occurs. These data clearly indicate the existence of aromatic amine groups like pyridine, which is in complete accordance with the 13 C solid state NMR data indicating the dominance of aromatic graphitic structures.

Zeta-potential experiments as a function of pH on the HC-Alb2-Y sample series (Figure 7) can further characterize the nitrogen functions. The as-synthesized hydrothermal powders (Figure 7a) clearly show positive zeta-potentials below pH~ 4 indicating the coexistence of positively charged amino groups at the particle's surface with more negatively charged oxygen moieties. In addition, the quantitative zeta-potential values at pH= 2.3 increase with increasing content of nitrogen within the materials, the highest value (+26.6 mV) belonging to HC-Alb10 (8% N) and the lowest (+9.7 mV) belonging to HC-Alb1 (3% N). Pure hydrothermal carbon from glucose shows no positive zetapotential value in the whole examined range, as expected.

Analyzing the effects of calcination on the zeta-potential (Figure 7 b) it is found that the positive character at low pH is not affected by temperature treatment. The material is overall even more basic, as the isoelectric point shifts towards higher pH values. This can be connected with the FT-IR and 13 C solid state NMR observations which show that the oxidic surface functions (COOH, OH) are degraded with temperature. The absolute zeta-potential value in the low pH range does not increase with temperature but stays constant at +18 mV. Again, we attribute that to the conversion of primary and secondary amine groups to aromatic nitrogen moieties generated and being stable at higher temperatures. These results show that OvA takes an active role in the carbonformation reactions: 1) it stabilizes the precipitating carbon in the hydrothermal reaction as nanosized particles; 2) by crossreactions with the condensing monomers it brings up to 8 w% nitrogen to the final carbonaceous material. It is known from previous studies [66][67][68][69] that glucose dehydrates three times to form hydroxymethyl furfural (HMF) in a first step, which is the main intermediate reacting species to form the carbonaceous scaffold. [START_REF] Titirici | [END_REF]38] In presence of proteins and amino acids, the aldehydic group of HMF will react to form Schiff bases and then to Maillard products, as it was previously documented in the chemistry of HMF [70] and furans in general. [71] with amino-group containing molecules such as amino acids. [72] Similar experiments using other common proteins (hemoglobin, gelatin, casein; results not shown) revealed a systematic increase in nitrogen content of the carbon as a function of protein but within this series a sensible reduction of the primary particle size was not found. We have observed during preliminary experiments that use of molecules with accentuated hydrophobic character or amphiphilic molecules act as emulsion stabilizers. As reported in the experimental part, half of the OvA amino acid residues are hydrophobic, and it displays surface activity. Under the reaction conditions, the protein certainly reacts via the Maillard reaction with the polysaccharide, obviously creating an effective stabilizer for the intermediate droplets later transferring to carbon.

Conclusion

We have demonstrated that it is possible to synthesize well defined carbonaceous nanostructures with a developed spongelike mesopore system by using hydrothermal treatment at 180°C of glucose in the presence of ovalbumin. This simple process provides a unique and sustainable pathway to control the primary particle stability in hydrothermal carbonization, resulting in the formation of nanometer-sized sponge scaffolds composed of spherical particles (20-50 nm). Another advantage beside the morphological control is that the protein efficiently introduces nitrogen functions (as much as 8 wt%) into the carbonaceous scaffold, which are known to be very favourable for catalysis, adsorption, energy storage devices, and also to increase the oxidation stability of engineering carbons. Classical characterization methods (nitrogen adsorption-desorption, TEM, SEM, FTIR, XPS and elemental analysis) were used to describe the porous carbon structures. Temperature treatments have shown that all materials keep the nitrogen content almost unchanged up to 950°C, while providing a higher thermal and oxidation stability as compared to pure carbon based hydrothermal scaffolds. This was attributed to a remarkably increased aromaticity.

Last but not least, the whole synthetic process occurs in an energy and atom-saving fashion from cheap and sustainable resources, as employed temperatures in the first step remain below 200°C and neither metals nor surfactants have been used to catalyze and control the reaction. Even if this work has employed pure model compounds, we believe that more complex biomass like starch or chicken egg white might be actually employed in a similar fashion, resulting potentially even in a wide variety of structures.

Experimental Section

Synthesis. For all samples, 13.5 ml of a deionized water solution containing 1.5 g of D(+)-glucose and variable amount of ovalbumin (OvA) was used. D(+)-glucose (CAS: 50-99-7) and ovalbumin (OvA) from chicken egg white (grade V, CAS: 9006-59-1) were used as received (Sigma Aldrich). OvA is the major component of egg white (54%) and is a phosphoglycoprotein with a molecular mass of about 45 kDa [73] and a sequence of 385 amino acids [74] including one disulfidic bond, four sulfhydryl groups, zero to two phosphoryl groups and a carbohydrate chain. [75] Half of OvA amino acids residues are hydrophobic and one third are charged, for the most part acidic, giving to the molecule an isoelectric point (IEP) of 4.5. [76] The mixture was sealed into a glass vial inside a typical PTFE-lined autoclave system and hydrothermally reacted at 180°C for 24h. After reaction, the autoclave is cooled down in a water bath at room temperature. The obtained black solid flake-like powder is then separated from the remaining aqueous solution by filtration and put into an oven at 80°C under vacuum overnight for drying. Further calcination was performed under a N2 flow using a 1°C/min ramp from room up to the desired temperature, at which the sample is left for 4 h. Typical sample notation follows: HC-AlbX, where X stands for the amount of OvA introduced (X=1, 2, 3, 4, 5, 10, respectively for 0.15 g, 0.30 g, 0.45 g, 0.50 g, 0.75 g and 1.50 g). Calcined samples are referred to HC-AlbX-Y, where Y stands for the calcination temperature (Y= 350, 550, 750 and 950) given in degrees celtius. (e) Thermogravimetric analysis. TGA was performed on a Netzsch Iris TG 209 FI instrument under nitrogen gas from room temperature to 1000°C using a 5°C/min ramp.

(f) Infrared spectroscopy. FT-IR spectroscopy was done in ATR geometry on a Varian 1000 FT-IR spectrometer, Scimitar Series (FTS 1000) between 4000-600 cm -1 .

(g) Surface charge analysis. Z-potential measurements were realized on a Malvern Zetasizer Nano ZS instrument. Carbon powders were dispersed and stirred (400 rpm) in water (CNaCl= 0.03 M) at different pH values for three days at 25°C to reach equilibrium. Disposable clear zeta cells (DTS1060c) were used for determing the zeta potential, which was calculated by Smoluchowski's equation via the measurement of the electrophoretic mobility. [START_REF] Hunter | Zeta potential in colloid science : principles and applications[END_REF] Due to the heterogenous nature of our samples, in order to avoid immediate pollution of the electrodes, all solutions were filtered using PTFE 5µm disposable filters and fast monomodal acquisition setting (suitable for high conductivity solution and no distribution plot is required, voltage is set automatically according to the conductivity of the sample) was used. Although the zeta potentials of the anisotropic particles obtained by Smoluchowski's equation may be underestimated, [START_REF] Johnson | [END_REF] they are still useful to qualitatively understand the effects of the surface charge of the HC particles.

(h) Solid state NMR. 1 H and 13 C solid-state Magic Angle Spinning (MAS) NMR experiments have been acquired on Bruker Avance 300 MHz (7 T) spectrometer using the 4 mm zirconia rotors as sample holders spinning at MAS rate νMAS= 14 kHz. The chemical shift reference was tetramethylsilane (TMS; δ= 0 ppm). Proton-to-carbon CP MAS was used to enhance carbon sensitivity: recycle delays for HC-Alb2 and HC-Alb2-550 samples are, respectively, 3 and 5s and TPPM decoupling is applied during signal acquisition. Crosspolarization transfers were performed under adiabatic tangential ramps, [79,80] to enhance the signal with respect to other known methods [81] and CP time tCP= 3 ms was found to be a good compromise in order to have a good overview on all carbon species. Number of transients is 1200. Peak attribution was done after references [82][83][84][85][START_REF] Breitmaier | Carbon-13 NMR spectroscopy : high-resolution methods and applications in organic chemistry and biochemistry[END_REF] and after supporting information.

(i) X-ray Photoelectron spectroscopy. XPS spectra were collected at the facilities of the Laboratoire de Réactivité de Surface (University Pierre et Marie Curie, Paris, France) on a SPECS (Phoibos MCD 150) X-ray photoelectron spectrometer, using a Mg Kα X-ray source (hν = 1253.6 eV) having a 150 W (12 mA, 12.5 kV) electron beam power and a 7X2 mm spot size. The emission of photoelectrons from the sample was analyzed under ultrahigh vacuum conditions (10-8 Pa). High-resolution spectra were collected at pass energy of 10 eV for C1s, O1s, and N1s core XPS levels. After collection, the binding energies were calibrated with respect to either the C-C/C-H components of the C1s peak at a binding energy of 284.8 eV. The spectral decomposition was performed by using Gaussian functions after background subtraction.

Entry for the Table of Contents

Figure 1 .

 1 Figure 1. SEM images for samples HC-AlbX, X=1, 2, 4.

Figure 2 .

 2 Figure 2. TEM images of HC-Alb1 and HC-Alb2. Scale bars is 1000 nm in (a, b) and 200 nm in (c, d).

Figure 3 .

 3 Figure 3. -a) Nitrogen adsorption-desorption isotherms and b) BJH pore size distribution curves for samples HC-AlbX (X=1-10).

Figure 4 .

 4 Figure 4.TEM images for higher temperature samples HC-AlbX-Y, X= 2,5. Y (=350, 750) indicates the calcination temperature. Scale bar is 200 nm.

Figure 5 .

 5 Figure 5. Nitrogen adsorption-desorption isotherms for samples HC-Alb2-Y ( Y= 350, 550, 750, 950); b) BJH pore size distribution curves for samples HC-Alb2-Y.

Figure 6 .

 6 Figure 6. 13 C CPMAS solid state NMR for the HC-Alb-2 before and after heat treatment at T= 550°C. Peak attribution: 1 -ketones and aldehydes; 2 -

Figure 7 .

 7 Figure 7. Zeta-potential experiments on a) HC-AlbX (X= 1, 2, 3, 4, 5, 10) and b) HC-Alb2-Y (Y= 350, 550, 750) materials. HC indicates the hydrothermal carbon synthesized without OvA as reference.

  Elemental analysis. Chemical analysis was performed on a (C, N, O, S, H) Elementar Vario Micro Cube. (b) Scanning Electron Microscopy. SEM images were acquired on a secondary electron and In-lens detectors. (c) Transmission Electron Microscopy. TEM was perfermed on a Zeiss EM 912 instrument equipped with CCD camera and a filament of LaB6 under a 120 kV tension. All samples are microtomed.(d) Porosimetry. N2 adsorption and desorption isotherms were performed at 77 K with a Quadrachrome Adsorption Instrument; BET and BJH methods was used for specific surface area determination and pore size distribution calculations.

  (Please choose one layout) Layout 1: FULL PAPER In this work, we present a simple and green one-step pathway towards nitrogen-doped carbon nanostructures with controlled mesoporosity via hydrothermal treatment of glucose in the presence of model proteins. ((Author(s), Corresponding Author(s)*) this work, we present a simple and green one-step pathway towards nitrogendoped carbon nanostructures with controlled mesoporosity via hydrothermal treatment of glucose in the presence of model proteins. ((Author(s), Corresponding Author(s)*) Page No. -Page No. ((Title))

Table 1 .

 1 Physico-Chemical Characterization of the Carbon Materials

	Material	BET	S(m 2 /g)	N [%w]	C[%w]	N/C
	HC-Alb1	83	74	3.1	67.0	21.7
	HC-Alb1-350	100	100.4	3.2	67.8	20.9
	HC-Alb1-550	397	397.6	3.6	86.4	23.8
	HC-Alb1-750	381	381.4	3.1	84.9	27.5
	HC-Alb1-950	255	358.4	2.5	89.1	36.1
	HC-Alb2	99	112	4.9	67.6	13.7
	HC-Alb2-350	106	106.6	4.8	71.3	15.0
	HC-Alb2-550	359	359.3	4.8	82.2	17.0
	HC-Alb2-750	246	246.1	4.1	83.4	20.3
	HC-Alb2-950	152	152.1	3.2	87.3	27.2
	HC-Alb3	41	81	5.4	65.8	12.3
	HC-Alb3-350	50	89.9	5.9	68.5	11.7
	HC-Alb3-550	-321	217.2	6.0	81.0	13.6
	HC-Alb3-750	-152	115.2	5.4	83.3	15.4
	HC-Alb3-950	-	-	4.5	88.4	19.7
	HC-Alb4	40	34	6.1	66.8	10.9
	HC-Alb4-350	36	36.4	6.3	68.4	10.9
	HC-Alb4-550	106	106.5	6.4	79.4	12.4
	HC-Alb4-750	-1661	60.6	6.0	84.3	14.2
	HC-Alb4-950	64	64.8	4.6	89.0	19.4
	HC-Alb5	20	30	6.6	66.9	10.2
	HC-Alb5-350	50	24.8	7.0	69.5	10.0
	HC-Alb5-550	20	78.3	6.8	81.1	11.9
	HC-Alb5-750	90	19.2	6.3	82.3	13.0
	HC-Alb5-950	-	-	5.1	90.3	17.7
	HC-Alb10	34	56	8.0	65.6	8.2
	HC-Alb10-350	57	48.9	8.2	68.8	8.4
	HC-Alb10-550	33	78.7	8.4	80.0	9.6
	HC-Alb10-750	-203	61.1	7.0	83.0	11.9
	HC-Alb10-950	30	46.2	4.7	84.1	17.9

Acknowledgements

Dr. Lorenzo Stievano and Christophe Méthivier (Laboratoire de Réactivité de Surface, University Pierre et Marie Curie, Paris, France) are gratefully acknowledged for the XPS mesurements.

Supporting Information