Romain Nguyen Van Yen

Eric Sonnendrücker

Kai Schneider

Marie Farge

Éric Sonnendrücker

PARTICLE-IN-WAVELETS SCHEME FOR THE 1D VLASOV-POISSON EQUATIONS

whether they are published or not. The documents may come

Introduction

The state of a hot plasma out of local thermodynamic equilibrium is characterized by a particle distribution function depending on position, velocity, and time. When the effects of collisions are neglected, this distribution satisfies the partial differential equation (PDE) known as the Vlasov equation. In the presence of simplifying hypotheses and/or special symmetries, the complexity of the problem can be reduced to permit a resolution of the Vlasov equation using the Eulerian PDE discretization approach which requires a mesh in position and velocity. Recent progresses have pushed back the limit of the feasible to include for example the 5D gyrokinetic approximation [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Idomura | New conservative gyrokinetic full-f Vlasov code and its comparison to gyrokinetic δ-f particle-in-cell code[END_REF][START_REF] Wang | Gyro-kinetic simulation of global turbulent transport properties in tokamak experiments[END_REF]. However the full equation with six variables plus time remains a major computational challenge, which could take as much advantage of better numerical methods as of further machine development and parallelization.

In the past, even more stringent limitations in memory size and computing power have driven physicists to develop specific simulation methods, which approximate the plasma by a collection of macro-particles [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Hockney | Computer Simulation using Particles[END_REF] and can be linked to Monte-Carlo methods [START_REF] Aydemir | A unified monte carlo interpretation of particle simulations and applications to non-neutral plasmas[END_REF]. Although the precision of Monte-Carlo methods is limited by sampling error, they are known to be competitive for problems involving integrals over many dimensions. For the specific problem of solving the Vlasov equations, it is well known that particle methods offer a more straightforward implementation and an easier parallelization than their Eulerian counterparts. However, the debate over the best compromise in terms of overall development and simulation cost versus precision remains unsettled. Current outstanding applications of particle methods include magnetic confinement fusion [START_REF] Heikkinen | Full f gyrokinetic method for particle simulation of tokamak transport[END_REF][START_REF] Jolliet | A global collisionless PIC code in magnetic coordinates[END_REF], and laser-plasma acceleration [START_REF] Cormier-Michel | Unphysical kinetic effects in particle-in-cell modeling of laser wakefield accelerators[END_REF]. But in the last decade, there has been a tendency to challenge some results obtained using legacy particle methods on the grounds that the simulations were underresolved [START_REF] Nevins | Discrete particle noise in particle-in-cell simulations of plasma microturbulence[END_REF].

Plasma simulation using particles is dominated by derivatives of the particle-in-cell (PIC) method [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Hockney | Computer Simulation using Particles[END_REF], to whom a vast amount of research has afforded the status of a reference. PIC schemes are obtained by convolving the empirical density corresponding to the distribution of marker particles by a fixed kernel, also known as charge assignment function [START_REF] Hockney | Computer Simulation using Particles[END_REF]. Their numerical properties are well understood empirically and they also benefit from an elegant physical interpretation [START_REF] Langdon | Theory of plasma simulation using finite-size particles[END_REF]. The most widely used PIC method uses a piecewise linear kernel. Smoother kernels have been studied in the past [START_REF] Okuda | Splines and high order interpolations in plasma simulations[END_REF] but computational limitations have made their use relatively rare, except for specific applications where precision and reduction of noise are of high concern.

In [START_REF] Nguyen Van Yen | Wavelet based density estimation for noise reduction in plasma simulation using particles[END_REF], we have studied another way of estimating the density from the particle positions, called waveletbased density estimation (WBDE). WBDE was first introduced in the more general context of nonparameteric density estimation [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF], as an alternative to the classical kernel density estimation (KDE) method. Instead of being convolved with a kernel, the empirical density is projected onto a finite dimensional linear space spanned by a family of orthogonal wavelets. The projection space is determined from the data itself, which allows for a refined representation around sharp features, and could make the method more precise than PIC for a given computational cost.

The aim of this paper is to present a new scheme for the 1D Vlasov-Poisson equations, that we call particlein-wavelets (PIW) because it relies crucially on wavelet expansions of the Dirac delta functions corresponding to every particle. WBDE is used for the density estimation step in the PIW scheme. The information pertaining to each particle thus gets spread to neighboring wavelets, in the same manner as it gets spread among neighboring grid points in the PIC case. The general philosophy of this work has been to build the scheme from scratch, taking as much advantage of wavelets as possible, so as to have a proof of principle of various properties specific to wavelets. We are focusing mostly on understanding the precision of the scheme and the conservation properties, as opposed to optimizing numerical efficiency.

Other methods have been proposed, and sometimes successfully implemented, to improve PIC simulations, but they all rely on the application of a linear operator to the particle distribution function in order to reduce the noise. Linear denoising operators are optimal when the regularity of the signal to be denoised is homogeneous in space. However, the distribution function within collisionless plasmas is known to develop sharp features. Therefore a nonlinear thresholding approach could be more efficient, and this work is a first step in this direction.

We first recall some background, on the Vlasov-Poisson equations in one space dimension on the one hand, and on wavelets on the other hand. Reference semi-Lagrangian Eulerian and PIC schemes are also described. After that we proceed to define the various steps of the PIW (particle-in-wavelets) wavelet-based scheme. Finally, numerical results from the three schemes are presented and the properties of the new PIW scheme are assessed in detail.

Background

Vlasov-Poisson equations

The Vlasov equation governs the evolution of the particle distribution function f (x, v, t) of a given species in a collisionless plasma. In the 1D case which we focus on, v varies in R, and for simplicity, we assume periodicity in the x direction, i.e. x ∈ T := R lZ , where l is the period. Γ = T × R is the phase space for one particle in the plasma. We restrict ourselves to the case where the ion density is assumed to remain constant, and the magnetic field is neglected. The system is then fully described by the Vlasov-Poisson equations, which can be written in non-dimensional form as follows:

∂ t f + v∂ x f + ∂ x φ∂ v f = 0 (1)
∂ xx φ + 1 -l R f (x, v, t)dv = 0 (2)
where φ(x, t) is the electric potential. In the following, we denote by E(x, t) = -∂ x φ(x, t) the electric field, and by ρ(x, t) = R f (x, v, t)dv the charge density. Note that we have adopted here the following slightly uncommon normalizations for ρ and f :

Γ f (x, v, t)dxdv = 1, (3)
T ρ(x, t)dx = 1. ( 4
)
We shall also need the definitions of two physical quantities associated to the Vlasov-Poisson problem, namely the electrostatic energy

V(t) = 1 2 T |E(x, t)| 2 (5)
and the total energy

E(t) = 1 2 Γ v 2 f (x, v, t)dxdv + V(t), (6)
the latter being constant for solutions of the Vlasov-Poisson equations. The goal of this paper is to describe a new numerical scheme for solving (1-2) together with an initial condition for f . We will use the fact that any solution of (1) is constant along the characteristic trajectories defined by:

x ′ (t) = v(t) v ′ (t) = -E (x(t), v(t), t) (7)
Existing approaches to the same problem are divided into two main categories, that were already mentionned in the introduction:

• particle solvers, which integrate the characteristic equations [START_REF] Daubechies | Two-scale difference equations. I. Existence and global regularity of solutions[END_REF] and discretize the potential on a grid in the x variable, • Eulerian solvers where the full distribution function f is discretized on a phase space mesh or grid. The first category is especially appealing when the phase space has dimension 4 or more, since it is then very costly to use an Eulerian approach. However, the discretization of f using particles is subject to statistical sampling noise which normally makes it of rather low precision [START_REF] Nevins | Discrete particle noise in particle-in-cell simulations of plasma microturbulence[END_REF]. In particular, for the 1D problem that we are considering here, we expect Eulerian solvers to be much more precise for a given computational cost. This allows us to take the solution obtained with a high precision Eulerian solver as a reference for measuring the convergence of particle methods. We have chosen a well validated semi-Lagrangian (SL) solver using cubic spline interpolation [START_REF] Sonnendrücker | The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF].

It is a reasonable requirement that any new particle method should offer some improvement over the legacy PIC approach. As a contestant to compare with our newly proposed PIW method, we therefore take a PIC solver, which uses a triangular, piece-wise linear charge assignment function to estimate the charge density, and Fourier transforms to solve the Poisson equation.

When assessing the convergence of particle solvers, the loading scheme is especially important since it controls the initial amplitude of the noise. Provided that the initial data f 0 is factorized

f 0 (x, v) = f 0,x (x)f 0,v (v)
we ensure that the initial noise is small using the following classical approach. First, compute a 2D Hammersley sequence [START_REF] Hammersley | Monte Carlo Methods[END_REF], (a i , b i ) i∈N , Then set as initial particle positions

x i (0) = F -1 0,x (a i ) v i (0) = F -1 0,v (b i )
where F -1 0,x and F -1 0,v are the inverse mappings corresponding to the cumulative distribution functions F 0,x and F 0,v associated respectively to f 0,x and f 0,v , and are in practice approximated using a pre-generated lookup table and polynomial interpolation. With this initialization scheme, called "quiet start", the initial error on the particle distribution function is expected to scale like N -1 p , as opposed to N -1/2 p for a random start.

Wavelets

In this section we recall some notions on wavelets. We keep the same notations as in [START_REF] Nguyen Van Yen | Wavelet based density estimation for noise reduction in plasma simulation using particles[END_REF]. More background may be found in the literature (see e.g. [START_REF] Mallat | A wavelet tour of signal processing[END_REF]). The construction is done in R/lZ with l = 1, and the adjustements for l = 1 will be made further down the moment they are required. Let ψ be a 1-periodic wavelet generating an orthogonal basis of L 2 (R/Z), and ϕ be the associated scaling function.

The dilated and translated scaling functions and wavelets are denoted

ϕ λ (x) = 2 j λ 2 ϕ(2 j λ x -i λ ) (8)
ψ λ (x) = 2 j λ 2 ψ(2 j λ x -i λ ) (9)
where λ = (j λ , i λ ) is a multi-index characterizing the scale j λ and position i λ of a wavelet. We denote by V J the 2 J -dimensional linear space spanned by ϕ (J,i) 0≤i≤2 J -1 . An important property of the function ϕ is that it satisfies the two-scale recurrence equation [START_REF] Mallat | A wavelet tour of signal processing[END_REF]:

ϕ(t) = √ 2 S-1 i=0 h i ϕ(2t -i) (10)
where (h i ) 0≤i≤S-1 is called the scaling function filter.

The scaling function and wavelet coefficients of a function u ∈ L 2 (T) are defined by Only the values for positive k are given. The values for negative k are obtained from the conditions that r (1) is odd and r (2) is even.

u λ = u | ϕ λ (11)
u λ = u | ψ λ (12)
where • | • denotes the usual scalar product in L 2 (T). For any positive integer L, u can be decomposed as follows:

u = 2 L -1 i=0 u (L,i) ϕ (L,i) + ∞ j=L 2 j -1 i=0 u (J,i) ψ (J,i) (13)
where the first sum is the projection of u on V L , and the second sum contains all the additional details needed for the full reconstruction of u.

In the following we shall also need to represent first and second order differential operators in the space V J . The results necessary to do so were derived in [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF], and we summarize them here for completeness. For any u ∈ V J we have, using [START_REF] Heikkinen | Full f gyrokinetic method for particle simulation of tokamak transport[END_REF]:

u ′ = 2 J -1 i=0 u (J,i) ϕ ′ (J,i)
The Galerkin representation of the first order derivation operator in V J is thus defined by the matrix G such that:

G i1,i2 = ϕ (J,i1) | ϕ ′ (J,i2) = 2 J T dxϕ 2 J x -i 1 ϕ ′ 2 J x -i 2 = T dyϕ (y + i 2 -i 1 ) ϕ ′ (y) := r (1) i2-i1
In [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF], it is shown that the matrix elements r

i2-i1 can be determined by solving a linear system with coefficients depending only on the h i , and that they satisfy the antisymmetry relation r

-i = -r (1)
i , so that applying the matrix G is numerically equivalent to applying a centered finite difference operator. G is thus an antisymmetric matrix. The same procedure can be followed to obtain the Galerkin representation of the second derivative in V J , defined by the matrix L such that:

L i1,i2 = ϕ (J,i1) | ϕ ′′ (J,i2)
and which corresponds to a symmetric filter r

(2) i .

In the following we make use of only one orthogonal wavelet basis, namely the C 1 Coiflets with S = 6 introduced at the end of [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets II. Variations on a theme[END_REF], which we dub "R-Coiflet 6". This wavelet is more regular than the normal Coiflet 6 wavelet, hence the R which stands for "regular". The scaling function filter coefficients (h i ) are recalled in Table 1. Except otherwised noticed, the results presented below also hold for any other orthogonal wavelet basis. The filters r (1) and r (2) corresponding to the R-Coiflet 6 wavelet basis are provided in Table 2 for completeness. Note that 6 is the minimal filter length for which both r (1) and r (2) are well defined, which is why we have chosen this wavelet family. Another candidate having the same filter length was the Daubechies minimal phase family, but we have discarded it because of its lack of symmetry.

Description of the PIW scheme

Wavelet-based density estimation

Let us briefly recall the wavelet-base density estimation (WBDE) algorithm for estimating the charge density from the position of a finite number N p of particles. A more detailed account, illustrated by several applications, may be found in [START_REF] Nguyen Van Yen | Wavelet based density estimation for noise reduction in plasma simulation using particles[END_REF]. WBDE, first introduced in [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF] belongs to the class of density estimators constructed from orthogonal projections [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. Its main originality is that the members of the orthogonal family on which the density is projected are chosen in an adaptive manner, depending on the observed particle positions.

In this work, we have focused on the Vlasov-Poisson problem, for which only the spatial distribution of particles matters when estimating the electric field, while the velocity distribution is irrelevant. Therefore, we apply the WBDE algorithm to the spatial distribution only, which is one dimensional. In the units we are working with, this spatial distribution is identical to the charge density ρ. Defining the empirical density associated to the particle positions (x n ) 1≤n≤Np by

ρ δ (x) = 1 N p Np n=1 δ(x -x n ) ( 14
)
where δ is the Dirac measure, the idea of the WBDE algorithm is to project f δ on the finite dimensional subspace spanned by all the wavelets at scales coarser than L and part of the wavelets at scales j such that L ≤ j ≤ J, where the scales L and J are defined as follows:

L = log 2 N 1 1+2r 0 p (15) J = log 2 N p log 2 N p (16)
In these formulas, ⌈•⌉ and ⌊•⌋ denote respectively the lower and upper integer parts, and r 0 is the order of regularity of the wavelet, which is 1 in our case. The wavelets that are kept for L ≤ j ≤ J are determined according to the amplitude of the corresponding coefficients in ρ δ , in the spirit of the classical wavelet denoising algorithms, see [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF][START_REF] Donoho | Density estimation by wavelet thresholding[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The threshold at scale j is defined by:

Θ j = K j N p
where K is a constant which in principle depends on the regularity of the solution. In the following, we have fixed K = 1. The threshold was derived rigorously in [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF], and as detailed in [START_REF] Nguyen Van Yen | Wavelet based density estimation for noise reduction in plasma simulation using particles[END_REF], the N

-1 2 p
dependence can be justified heuristically from the asymptotic distribution of the wavelet coefficients when estimating a probability density from N p independent realizations.

To treat the case l = 1, we can rescale ρ δ to obtain a normalized distribution on R Z :

x → l -1 x (17)

ρ δ → lρ δ (18)
and notice that this change of variables leaves the scaling function and wavelet coefficients [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Hammersley | Monte Carlo Methods[END_REF] unchanged. Therefore at this stage the algorithm can be applied harmlessly as if l = 1.

In summary, the following steps are required to actually implement the estimation:

(1) compute some approximation ρ S (J,i) of the scaling function coefficients of ρ δ at scale J, by which the projection of ρ δ on the space V J can be expressed:

ρ S = 2 J -1 i=0 ρ S (J,i) ϕ (J,i) (19)
(2) apply the fast wavelet transform algorithm to compute the wavelet coefficients of ρ S at scales j such that L ≤ j ≤ J:

ρ S = 2 L -1 i=0 ρ S (L,i) ϕ (L,i) + J-1 j=L 2 j -1 i=0 ρ S (j,i) ψ (j,i) (20)
(3) apply the nonlinear threshold operator :

ρ W = 2 L -1 i=0 ρ S (L,i) ϕ (L,i) + J-1 j=L 2 j -1 i=0 ½ [0,+∞[ | ρ S (j,i) | -Θ j ρ S (j,i) ψ (j,i) (21)
where ½ [0,+∞[ is the characteristic function of the interval [0, +∞[. The remaining difficulty in the above procedure is its first step, to which we dedicate the next section entirely. If steps (2-3) are omited, we obtain a scheme which we call Linear PIW (L-PIW), and that we will use below for comparison purposes. In that case, the value of L has no influence, and we allow J to be chosen as a free parameter.

Approximation of empirical scaling function coefficients

To perform step (1) of the above WBDE algorithm, some approximate values ρ S λ of the coefficients ρ δ λ must be defined. According to [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF] and [START_REF] Hockney | Computer Simulation using Particles[END_REF], their exact values are:

ρ δ λ = 1 N p Np n=1 ϕ λ (x n ), (22)
which is difficult to compute because there is no analytical expression available for ϕ λ . In [START_REF] Nguyen Van Yen | Wavelet based density estimation for noise reduction in plasma simulation using particles[END_REF] we proposed to use the approximation ρ S (J,i) = 2 -J 2 ρ H 2 -J i , where ρ H is a histogram constructed over a grid of size 2 J in each direction. To increase the accuracy in the computation, we adopt here the following alternative. We known that the scaling function ϕ is supported on the interval [0, S -1], which can be discretized using 2 JLT regularly spaced points s k = 2 -JLT (S -1)k, at which we can compute approximations of the values of ϕ φk ≃ ϕ(s k ) using for example the iterative algorithm described in [START_REF] Daubechies | Two-scale difference equations. I. Existence and global regularity of solutions[END_REF]. The φk consitute a look-up table of approximations to the ϕ(s k ), which needs to be computed only once for any chosen wavelet basis. Now assume that a particle is located at position x and that we want to approximate

ϕ (J,i) (x) = 2 j 2 ϕ(2 j x -i).
This can be nonzero only if 0 ≤ 2 j xi ≤ S -1, that is, for S values of i. For each of these values of i we find the unique integer k such that s k ≤ 2 j xi < s k+1 , we let α = 2 j xix k x k+1x k and we compute a linear interpolation from the values stored in the lookup table, φλ (x) = (1α) φk + α φk+1 .

There remains to sum over all the particles to compute the scaling function coefficient of the charge density which we use in the following:

ρ S λ = 1 N Np n=1 φλ (x n ) ( 23
)
For each particle, we need to make S lookup operations in the table, so that the overall cost is proportional to SN p . For the size of the table, we have used J LT = 18. Note that when generalizing this algorithm to d dimensions, the required lookup table will remain one dimensional because the d-dimensional scaling functions are factorized into products of 1D scaling functions. Therefore there will be no increase in memory requirements for this stage of the algorithm.

Wavelet-Galerkin Poisson solver

Once the wavelet coefficients of the electron charge density ρ have been obtained by WBDE, the electric potential needs to be computed by solving the Poisson equation [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF]. A wavelet preconditioning technique has been previously developed to solve (2) in 3D using finite differences in the context of a PIC code [START_REF] Terzić | Particle-in-cell beam dynamics simulations with a wavelet-based Poisson solver[END_REF]. Here, we prefer to use a full Galerkin discretization of (2) in the wavelet basis, as introduced in Section 1.2:

l -2 Lφ S = 1 -ρ W ( 24
)
where φ S and ρ W stand respectively for the column vectors (φ S λ ) and (ρ W µ ), φ S stands for the unknown electric potential, and the l -2 prefactor comes from the rescaling of the domain (Eq. 17).

In order to better exhibit the properties of the scheme as it is likely to be implemented in the future for two or three spatial dimensions, or in the case of adaptive grids, we have decided to use an iterative method to solve [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF], even though a direct method would be faster in the present one-dimensional case. As such, this system is ill-conditionned and solving it directly using an iterative method would be inefficient. Fortunately, there is a well known technique to improve the condition number. The main idea is to solve for the wavelet coefficients instead of the scaling function coefficients, and to use a diagonal preconditioner in wavelet space. First ( 24) is rewritten

l -2 F LF † φ S = F (1 -ρ W ) ( 25
)
where F is the matrix form of the wavelet transform operation, and • † stands for transposition. Now ( 25) can be efficiently solved using the conjugate gradient method with the diagonal preconditioner defined by:

D λ,µ = δ λµ 2 j λ
The wavelet coefficients of the potential at the previous time step are used as initial guess to speed-up the convergence. Once φ S has been obtained, an inverse wavelet transform is applied to yield φ S .

Note that an alternative and non equivalent approach would be to write down a Galerkin discretization of (2) in the basis made of the wavelets whose coefficients are non-zero in the expansion (21) for ρ W . We would thus obtain an approximation φ W to the potential which has the same support in wavelet space, allowing for an adaptive strategy. In the present approach, φ S has a priori a full support in wavelet space. As we will see below, this has an important consequence on the self-force.

Interpolation and particle push

Like in the standard PIC method, it is essential that the interpolation method be compatible with the charge assignment scheme in order to avoid self-forces [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF]. The scaling function coefficients of the electric field can be obtained by Galerkin projection of the gradient operator in the scaling function basis (see Sec. 1.2):

E S λ = ϕ λ , ϕ ′ µ φ S µ = lGφ S
A natural first attempt is to directly interpolate E S at the particle positions to obtain the forces exerted on the particles. For this purpose we could use the approximation:

E S (x) = λ∈Λ E S λ φλ (x) (26)
but as we now show this does not guarantee the vanishing of self-forces. Indeed, let us assume for the rest of this section that there is a single particle located at position x. From ( 23), we see that, prior to the nonlinear thresholding step, the approximation by the scheme of the charge density is simply given by ρ S λ = φλ (x). Hence using the matrix notations defined above we have

E S (x) = lρ S † Gφ S = l -1 ρ S † GL -1 ρ W . (27)
Now since G and L are both circulating matrices, they commute, and using the fact that G is antisymmetric and L is symmetric, we obtain that GL -1 is also antisymmetric. If we would replace ρ W by ρ S above, we would immediately get that E S (x) = 0, but in general ρ W = ρ S because of the nonlinear thresholding step. Fortunately a simple remedy exists to enforce the vanishing of self-forces in all cases. Denote by M the matrix of the linear operator defined diagonally in wavelet space by the mask corresponding to the thresholding operation. By definition we have that ρ W = M ρ S so that ( 27) is equivalent to

E S (x) = l -1 ρ S † GL -1 M ρ S
which leads us to define a new candidate for the electric field reconstruction:

E W = lM Gφ S (28)
that is, E W is obtained from E S by discarding the same wavelet coefficients that were discarded when going from ρ S to ρ W . The self-force corresponding to E W is obtained as before:

E W (x) = l -1 ρ S † M GL -1 M ρ S
but now, using the fact that M is symmetric and GL -1 is antisymmetric, we obtain as desired

E W (x) = 0.
Once we know how to interpolate the electric field at the particle positions, the advancement in time of ( 7) is done using a classical Verlet integrator.

Numerical results

Landau damping

To begin with we consider a Landau damping test case, for which the initial condition is given by

f 0 (x, v) = l -1 (1 + ε cos(kx)) 1 √ 2π e -v 2 2 ,
where k = 0.5, l = 4π, ε = 10 -3 . The small value of ε implies that the evolution is dominated by the linear transport effect. The reference solution used for comparison was obtained with the SL solver, using 4097 grid points in the x and v directions. The timestep was δt = 0.025 for the reference solution, and δt = 0.05 for the others.

The time evolution of the electrostatic energy V(t) is shown in Fig. 1 (left). Thanks to the high order interpolation properties of the wavelets, the PIW solution matches the decay of the electrostatic energy better than the PIC one. The time evolution of total energy E(t) for the PIW and reference scheme is shown in Fig. 1 (right). Note that the kinetic energy appearing in E(t) is approximated by a discrete sum over the particles. Variations remain small for the two schemes, and are especially low for the PIW scheme with 2 19 particles.

This classical test case serves to validate the PIW method, but since the solution remains very smooth it is not very well suited to provide more insight into the effect of nonlinear thresholding. We therefore move to the more challenging two-stream instability.

Two-streams instability

The initial condition for this test case is

f 0 (x, v) = l -1 (1 + ε cos(kx)) 1 √ 2π e -(v-v 0 ) 2 2 + e -(v+v 0 ) 2 2 ,
where k = 0.2, l = 10π, ε = 10 -2 , and v 0 = 2.4. The reference distribution function is shown at t = 10 and t = 30 in Fig. 2. It was obtained using the SL solver with 4097 grid points in the x and v directions. The timestep was δt = 0.0125 for the reference solution, and δt = 0.05 for the others. For this test case the electrostatic energy (Fig. 3, left) first oscillates for a few periods, roughly up to t = 10, and then undergoes an exponential growth which leads to a strongly nonlinear behavior. The PIC and PIW schemes with 2 19 particles allow for a good approximation of the electrostatic energy up to t = 30. To assess the quality of the approximations of the electric field in the PIC and PIW schemes, we focus on the relative L 2 error:

δ E (t) = E(x, t) -E ref (x, t) E ref (x, t) ,
where E ref (x, t) is the electric field computed using the reference SL solver. To allow for the computation of δ E (t), E(x, t) is interpolated on the grid corresponding to E ref (x, t). In the PIW case, E(x, t) is first reconstructed on a grid from its scaling function coefficients using the methods presented in [START_REF] Genovese | Daubechies wavelets as a basis set for density functional pseudopotential calculations[END_REF]. δ E (t) is shown in Fig. 4 (left) at t = 30, as a function of the number of particles. For all considered methods, the scaling lies between N -1 p and N -1 2 p . PIW offers an improvement of roughly a factor 3 over PIC, over a wide range of particle numbers. Linear PIW with J = 7 is a good candidate but starts to saturate for N p = 2 19 probably because the spatial discretization error is reached. In contrast, nonlinear PIW remains uniformly better for any number of particles, thanks to its adaptive properties. In the right panel of Fig. 4, δ E (t) is plotted versus the total CPU time (in seconds) needed for the integration of the equations on the time interval [0, 30]. The same computer was used in all cases, but the comparison between the absolute values of the CPU time should be made with caution since the languages that were used are different (Fortran for the PIC code, and C++ for the PIW code). Anyway, an interesting observation which can be made is that the decrease of the error with computing time for the PIW code is uniformly N -1 p over the whole range of considered parameters. For the PIC and L-PIW methods, in contrast, a saturation is observed when the number of particles is too low for a given spatial resolution.

To get a feeling of how the error we make on the electric field, we plot it as a function of x for t = 30 (Fig. 5). It seems that the error can be qualitatively decomposed into the following three components:

(1) a large scale component, which is probably due to a phase shift between the PIW and reference solutions,

(2) localized peaks, where some sharp features in the electric field are not properly resolved, (3) noise spread out over the whole domain.

Our goals with nonlinear thresholding was to preserve the sharp features that develop due to nonlinear interactions.

Discussion

The PIW approach was proposed as an alternative to PIC for plasma simulation using particles, and successfully implemented to solve the simple one dimensional Vlasov-Poisson toy model. The main conceptual difference between PIC and PIW is that the empirical distribution of particles is convolved with a kernel for the former, whereas it is projected onto a finite dimensional linear space for the latter. Thanks to the use of wavelets to generate the projection space, the denoising step which can be included in PIW and which we have studied is genuinely different from those that are developed for PIC schemes. Indeed, the wavelet representation can adapt locally to the regularity of the density, whereas linear approaches such as Krook operators [START_REF] Mcmillan | Long global gyrokinetic simulations: Source terms and particle noise control[END_REF], Fourier filtering [START_REF] Jolliet | A global collisionless PIC code in magnetic coordinates[END_REF] or coarse graining [START_REF] Chen | Coarse-graining phase space in delta f particle-in-cell simulations[END_REF], cannot.

The adaptive denoising strategy is especially useful when the density develops sharp features, which is known to occur generically as soon as nonlinear effects become important. To better understand this effect, we have compared a linear Landau damping test case, where the density remains homogeneously smooth, to a two stream instability test case. We have shown that for the first, PIW behaves as a higher order PIC approach, offering a constant improvement in terms of L 2 error, but no improvement in rate of convergence. In contrast, for the second test case, fine resolution was required to resolve sharp features developing in the electric field.as a result of nonlinear effects. The PIW scheme is able to automatically increase the resolution locally when the number of particles is sufficient to guarantee that the estimation of the density will not get corrupted by noise. This adaptivity guarantees that the error decays like the inverse square root of computing time over a wide range of particle numbers. The topic of discretization error control using wavelets has been studied elsewhere, see e.g. [START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF] for a review in the context of computational fluid dynamics.

For the simple 1D case considered here, the improvement observed was roughly a factor 3 of precision. The increase in computational cost of the PIW scheme with respect to the PIC scheme is proportional to the length of the wavelet filter, which we have fixed to 6 in this study. Extending the full PIW scheme to more realistic 2D or 3D geometries is a much more ambitious challenge. A first undertaking, complementary to the present study, could be to interface wavelet denoising using WBDE in an existing production code in order to further support the possible benefits of the approach.

h0 = 0 .

 0 03858077774789 h2 = -0.07716155549577 h4 = 0.74568755893443 h1 = -0.12696912539621 h3 = 0.60749164138568 h5 = 0.22658426519707

Figure 1 .

 1 Figure 1. Comparison of PIW, PIC and reference solvers for the Landau damping test case. Left: electrostatic energy as a function of time. Right: kinetic energy as a function of time.

Figure 2 .

 2 Figure 2. Two streams in stability test case. Particle distribution function obtained using the reference solver at t = 10 (left) and t = 30 (right).

Figure 3 .

 3 Figure 3. Comparison of PIW, PIC and reference solvers for the two-stream in stability test case. Left: electrostatic energy as a function of time (Eq. 5). Right: total energy as a function of time (Eq. 6).

  ,

Figure 4 .

 4 Figure 4. Two streams instability test case. L 2 error on the electric field at t = 30, as a function of number of particles (left) and computing time (right).

Figure 5 .

 5 Figure 5. Two streams in stability test case. Error on the electric field as a function of space at t = 30.

Table 1 .

 1 Scaling function filter coefficients for the C 1 Coiflet family with 2 vanishing moments introduced in[START_REF] Daubechies | Orthonormal bases of compactly supported wavelets II. Variations on a theme[END_REF] and used throughout this paper.

		i	0	1	2	3	4
	r	(1) i	0	-0.619368727403235 -1.403784906417689	2.427684222146337 -0.839028531550099
	r	(2) i	0.622642595976898 -0.389477243531951	0.179082676222473 -0.326853461357942	0.225926730678971

Table 2 .

 2 Stencils for first and second derivative operators in the R-Coiflet 6 scaling functions basis.

thank Diego del-Castillo-Negrete for his inspiration in the pursuit of PIC denoising, and Matthieu Haefele for developing and supporting the VLASY platform. * * The authors acknowledge financial support from the Euratom-CEA association and the French Federation for Fusion Studies. This work, supported by the European Communities under the contract of Association between EURATOM, CEA and the French Research Federation for fusion studies, was carried out within the framework of the European Fusion Development Agreement. The views and opinions expressed


Sergey Gavrilyuk

Multiphase Flow Modeling via Hamilton's Principle

We present here a variational approach to derivation of multiphase flow models. Two basic ingredients of this method are as follows. First, a conservative part of the model is derived based on the Hamilton principle of stationary action. Second, phenomenological dissipative terms are added which are compatible with the entropy inequality. The variational technique is shown up, and mathematical models (classical and non-classical) describing fluid-fluid and fluid-solid mixtures and interfaces are derived.

Introduction

For any vectors a, b we shall use the notations a • b for their scalar product (the line vector (covector) is multiplied by the column vector) and a ⊗ b for their tensor product (the column vector is "multiplied" by the line vector). Linear transformations of vectors into vectors (second order tensors) are denoted by uppercase letters. The divergence of a linear transformation A is the covector divA such that, for any constant vector a div(A) • a = div(Aa).

In particular, for any vector field v

div(Av) = divA • v + tr A ∂v ∂x div(u ⊗ v) = vdivu + ∂v ∂x u
In the last formula we have been identified the line vectors and the column vectors. The identical transformation is denoted by I, and the gradient operator by ∇. The elements of a linear transformation A in a given base e j , j = 1, 2, 3 are denoted be A i j : Ae j = A i j e i where i means lines and j columns. The repeated Latin indices imply summation. In particular, (a ⊗ b) i j = a i b j .

If f (A) is a scalar function of A, the matrix of ∂f ∂A is defined by the formula ∂f ∂A

j i = ∂f ∂A i j .
The definition of the divergence of a linear transformation used in France is opposite to the one used in USA where the divergence of A is defined as follows div(A) • a = div(A T a)

for any constant vector a. Once an appropriate definition has been selected, it should remain unchanged during calculations.

Variations

The results of this Section are well known and can be partially found, for example, in the book by [START_REF] Berdichevsky | Variational Principles of Continuum Mechanics[END_REF]. Let D(t) be a volume of the physical space occupied by a fluid at time t. A particle is labeled by its initial position X in the reference space D(t 0 ), the motion of the continuum is defined as a diffeomorphism from D(t 0 ) into D(t) :

x = ϕ t (X).

The velocity u and the deformation gradient F are defined by v = ∂ϕ t (X) ∂t , F = ∂x ∂X = ∂ϕ t (X) ∂X ,

Let

x = Φ t (X, ε), be a one-parameter family of virtual motions of the medium such that

Φ t (X,ε)| ε=0 = ϕ t (X).
Here ε is a scalar defined in the vicinity of zero. We define the virtual displacement δx (t, X) associated with the virtual motion :

δx (t, X) = ∂Φ t (X,ε) ∂ε ε=0 . (1)
Let us consider any quantity f defined in the Eulerian coordinates (or in the Lagrangian coordinates)

(t, x) -→ f (t, x) or (t, X) → f (t, X)
We abuse the notation by using the same letter f for completely different functions, because the last one is obtained from the first one by replacing x = ϕ t (X). To avoid a possible confusion, we will often indicate the arguments of these functions. The corresponding virtual families of f in the Eulerian coordinates will be denoted by f (t, x, ε) and f (t, X, ε). The relation between them is as follows :

f (t, X, ε) = f (t, Φ t (X, ε), ε) (2)
Let us define the Eulerian variations of f

δf (t, x) = ∂ f ∂ε ε=0
and the Lagrangian variation :

δf (t, X) = ∂ f ∂ε ε=0
Differentiating [START_REF] Alvarez-Samaniego | Large time existence for 3D water waves and asymptotics[END_REF] with respect to ε we obtain the following relation between two types of variations :

δf = δf + ∂f ∂x δx = δf + ∇f • δx (3)

Variation of density

The density ρ (t, X) is determined by

ρ det F = ρ 0 (X),
where ρ 0 (X) is the reference density. Lemma

δF = ∂δx ∂X , δ (det F ) = det F tr F -1 δF = det F div (δx)
Here the variation δx defined by ( 1) is considered as a function of the Eulerian coordinates. In the following, all operations from the vector analysis (div, ∇, ... ) will be applied only to functions of the Eulerian coordinates. Proof By definition, F (t, X, ε) = ∂Φ t (X, ε) ∂X .

Differentiating it ones with respect to ε we obtain the first formula. The proof of the second formula we present below is due to P. Casal. Let w 1 , w 2 , ..., w n be n constant vectors. We denote by det (w 1 , w 2 , ..., w n ) the determinant of the matrix having columns w 1 , w 2 , ..., w n . For any linear transformation A its determinant and its trace satisfy the following relations: det (Aw = tr(dAA -1 ) det (A) det(w 1 , w 2 , ..., w n ) = det (A) tr(A -1 dA) det(w 1 , w 2 , ..., w n ).

Hence, d (detA) = det A tr A -1 dA .

It is now sufficient to take d = δ and to apply this result to det (F ). The Lemma has been proved. In particular, the Lemma and (3) imply : δρ = -ρdiv (δx) , δρ = -div (ρδx) .

Variation of entropy

The entropy per unit mass η is conserved along trajectories, so it is a function only of the Lagrangian coordinates : η = η 0 (X). Hence, its Lagrangian and Eulerian variations are : δη = 0, δη = -∇η • δx.

Variation of velocity

The variation of the velocity is direct :

δv = ∂δx ∂t , δv = Dδx Dt - ∂v ∂x δx,
where the material derivative D Dt is defined as

D Dt = ∂ ∂t + v • ▽.

Hamilton's principle

Let us consider the Hamilton action :

a = t1 t0 Ldt,
where [t 0 , t 1 ] is a time interval and L is the corresponding Lagrangian. The Hamilton principle is : for any virtual motion satisfying

Φ t0 (X,ε) = ϕ t0 (X), Φ t1 (X,ε) = ϕ t1 (X), (4)
Φ t (X,ε) = ϕ t (X), X ∈ ∂D(t 0 ), (5)
and imposed constraints, the variation of the Hamilton action is stationary : δa = 0.

Let us remark that (4), [START_REF] Drew | Theory of multicomponent fluids[END_REF] imply

δx| t=t0 = δx| t=t1 = 0, δx| ∂D(t) = 0.
We will recall first the Lagrangians for classical models of continuum mechanics and the corresponding constraints.

Lagrangian of an incompressible fluid

L = D(t) ρ |v| 2 2 dD
The density ρ and the velocity field v verify the constraints :

Dρ Dt = 0, div (v) = 0.
Lagrangian of a compressible fluid

L = D(t) ρ |v| 2 2 -e (ρ, η) dD
The density ρ, the velocity field v and the entropy η verify the constraints :

Dρ Dt + ρdiv (v) = 0, Dη Dt = 0.
The specific internal energy e (ρ, η) is a given function satisfying the Gibbs identity :

θdη = de + pd 1 ρ ,
where θ is the temperature and p is the pressure.

Lagrangian of an elastic body

L = D(t) ρ |v| 2 2 -e (G, η) dD.
Here

G = F F T -1 = F T -1 F -1
is the Finger tensor. In the case of isotropic solids the energy is a function of the invariants of G. This Lagrangian is suitable for the Eulerian description of solids. Another possibility commonly used in the literature is to consider the specific energy in the form e (C, η) where C = F T F is the right Cauchy-Green deformation tensor.

4 Governing equations of classical continua Hence δρ = -∇ρ • δx, div (δx) = 0.

To satisfy the last constraint we take

δx = rot (A)
where A is a smooth vector field vanishing at the boundary D(t) :

A| ∂D(t) = 0.
The variation of the Hamilton action is

δa = d dε ε=0 t1 t0 D(t) ρ |v| 2 2 dD = t1 t0 D(t) δ ρ |v| 2 2 dD = t1 t0 D(t) δρ |v| 2 2 + ρv• δv dD = t1 t0 D(t) - |v| 2 2 ∇ρ • δx + ρv• Dδx Dt - ∂v ∂x δx dD = t1 t0 D(t) -▽ ρ |v| 2 2 • δx + ρv• Dδx Dt dD = t1 t0 D(t) -▽ ρ |v| 2 2 • δx + ρv• ∂δx ∂t + ∂δx ∂x v dD = t1 t0 D(t) ∂ (ρv•δx) ∂t + div (ρv ⊗ v) •δx- ρ |v| 2 2 δx dD - t1 t0 D(t) ∂ρv ∂t + div (ρv ⊗ v) •δxdD
Using the Gauss-Ostrogradski formula and taking into account the fact that δx vanishes at ∂D(t) and δx| t=t0 = δx| t=t1 = 0, we finally get

δa = - t1 t0 D(t) ∂ρv ∂t + div (ρv ⊗ v) •δxdD
Let us remark that for any function f

t1 t0 D(t) ▽f • rot (A) dD = t1 t0 D(t) div (f rot (A)) dD = t1 t0 ∂D(t) f rot (A) •ndD = 0 if rot (A) •n = 0.
Here n is the external unit normal vector to ∂D(t). But the last condition is verified if the field A is orthogonal to ∂D(t) :

A • s = 0 for any vector s tangent to ∂D(t). But A| ∂D(t) = 0. Hence δa = - t1 t0 D(t) ∂ρv ∂t + div (ρv ⊗ v) •δxdD = - t1 t0 D(t) ∂ρv ∂t + div (ρv ⊗ v) •rot (A) dD.
Hence, there exists a scalar function p called pressure such that

∂ρv ∂t + div (ρv ⊗ v) = -▽ p.
So, we have obtained the Euler equations of incompressible non-homogeneous fluids :

∂ρv ∂t + div (ρv ⊗ v+pI) = 0, Dρ Dt = 0, div (v) = 0.
When the mass forces f are added, the equations become

∂ρv ∂t + div (ρv ⊗ v+pI) = ρf , (6)
Dρ Dt = 0, div (v) = 0.
In the case of potential forces (f = -∇Ψ (x)) the equations admit the energy conservation law :

∂ ∂t ρ |v| 2 2 + Φ + div ρv |v| 2 2 + Φ + pv = 0.

Compressible fluids

Let us consider the Hamilton action for compressible fluids :

a = t1 t0 Ldt = t1 t0 D(t) ρ |v| 2 2 -ρe(ρ, η) dD The constraints Dρ Dt + ρdiv (v) = 0, Dη Dt = 0,
can by rewritten in the Lagrangian coordinates in the form

ρ det F = ρ 0 (X) , η = η 0 (X). Hence δρ = -div(ρδx), δη = -∇η • δx.
The variation of the Hamilton action is :

δa = d dε ε=0 t1 t0 D(t) ρ |v| 2 2 -ρe(ρ, η) dD = t1 t0 D(t) δ ρ |v| 2 2 -ρe(ρ, η) dD = t1 t0 D(t
-ρδx• ▽p ρ + ρv• ∂δx ∂t + ∂δx ∂x v dD = t1 t0 D(t) - ∂ρv ∂t + div (ρv ⊗ v) + ▽p •δxdD = 0.
In this derivation, we have used the Gauss-Ostrogradski formula and have taken into account the fact that δx vanishes at ∂D(t) and δx| t=t0 = δx| t=t1 = 0. So, we have obtained the Euler equations of compressible fluids :

∂ρv ∂t + div (ρv ⊗ v+pI) = 0, Dρ Dt + ρdiv (v) = 0, Dη Dt = 0.
These equations admit the energy conservation law :

∂ ∂t ρ |v| 2 2 + e(ρ, η) + div ρv |v| 2 2
+ e(ρ, η) + pv = 0.

Isotropic elastic bodies

Here G = F F T -1 = F -T F -1
is the Finger tensor. The tensor G is more convenient for the Eulerian formulation of governing equations of isotropic elastic materials. Let us recall that the material is isotropic if for any orthogonal transformation O e OGO T , η = e (G, η) .

Hence, the internal energy depends only on the invariants of G. In particular, we choose such invariants in the form :

J i = tr G i , i = 1, 2, 3.
Let X = (X α ) be the Lagrangian coordinates, α = 1, 2, 3, x = x i be the Eulerian coordinates, i = 1, 2, 3. We introduce the curvilinear cobasis

e α = ∇X α , α = 1, 2, 3
which is dual to the natural curvilinear basis

e α = ∂x ∂X α , α = 1, 2, 3.
In particular, the vector e α is the α -th column of F -T :

F -T =    ∂X 1 ∂x 1 ∂X 2 ∂x 1 ∂X 3 ∂x 1 ∂X 1 ∂x 2 ∂X 2 ∂x 2 ∂X 3 ∂x 2 ∂X 1 ∂x 3 ∂X 2 ∂x 3 ∂X 3 ∂x3    = e 1 ,
Since DX α Dt = 0

we have by taking the gradient of that equation :

∂e α ∂t + ∇ (v • e α ) = 0, rote α = 0. (8)
In particular, equations [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF] and [START_REF] Gavrilyuk | A new form of governing equations of fluid arising from Hamilton's principle[END_REF] imply that

DG Dt + G ∂v ∂x + ∂v ∂x T G = 0.
Remark 1. It is interesting to note that the left hand side of this equation can be considered as an objective derivative. Indeed, let O(t) be a time dependent orthogonal transformation : O(t)O(t) T = I. Consider new variables denoted with primes :

x ′ = O(t)x, v ′ = dO(t) dt x+O(t)v.
The aim is to show that the equation for

G ′ = O(t)GO T (t)
will be the same as for G. We get :

DG Dt = ∂G ∂t + ∂G ∂x v = ∂ O T (t)G ′ O(t) ∂t ′ + ∂ O T (t)G ′ O(t) ∂x ′ ∂x ′ ∂t + ∂ O T (t)G ′ O(t) ∂x ′ ∂x ′ ∂x O T (t)v ′ -O T (t) dO(t) dt O T (t)x ′ = ∂ O T (t)G ′ O(t) ∂t ′ + ∂ O T (t)G ′ O(t) ∂x ′ dO(t) dt O T (t)x ′ + ∂ O T (t)G ′ O(t) ∂x ′ v ′ - dO(t) dt O T (t)x ′ = ∂ O T (t)G ′ O(t) ∂t ′ + ∂ O T (t)G ′ O(t) ∂x ′ v ′ = dO T (t) dt G ′ O(t) + O T (t)G ′ dO(t) dt + O T (t) D ′ G ′ Dt ′ O(t). Hence DG Dt + G ∂v ∂x + ∂v ∂x T G = = dO T (t) dt G ′ O(t) + O T (t)G ′ dO(t) dt + O T (t) D ′ G ′ Dt ′ O(t) +O T (t)G ′ O(t) O T (t) ∂v ′ ∂x ′ O(t) -O T (t) dO(t) dt + O T (t) ∂v ′ ∂x ′ O(t) -O T (t) dO(t) dt O T (t)G ′ O(t) = O T (t) D ′ G ′ Dt ′ + G ′ ∂v ′ ∂x ′ + ∂v ′ ∂x ′ T G ′ O(t).
Finally,

D ′ G ′ Dt ′ + G ′ ∂v ′ ∂x ′ + ∂v ′ ∂x ′ T G ′ = O(t) DG Dt + G ∂v ∂x + ∂v ∂x T G O T (t).
Remark 2. The differential constraint rote α = 0 in ( 8) is compatible with the evolution equation for e α : if rote α vanishes at time t = 0, it vanishes at all t > 0. This constraint is automatically satisfied in the onedimensional case. In the multi-dimensional case equation ( 8) can also be replaced by

∂e α ∂t + ∇ (v • e α ) = -rote α ∧ v. (9)
We will obtain now the expression for the stress tensor. Let e be the specific internal energy : e = e (G, η)

where η is the specific entropy. Let us consider the variation of the internal energy E i at fixed value of η :

δE i = δ D(t) ρe (G, η) dD
We have

δE i = δ D(t) ρe (G, η) dD = D(0) ρ 0 δe (G, η) dD 0 = D(0)
ρ 0 tr ∂e ∂G δG dD 0 .

The matrix ∂e/∂G is symmetric. Since

δG = δF -1 T F -1 + F -1 T δF -1 = -F -1 δF F -1 T F -1 -F -1 T F -1 δF F -1 = - ∂δx ∂x T G -G ∂δx ∂x
we obtain

δE i = - D(0) 2ρ 0 tr ∂e ∂G G ∂δx ∂x dD 0 = D(t)
tr σ ∂δx ∂x dD.

The tensor

σ = -2ρ ∂e ∂G G
is called the Murnaghan stress tensor [START_REF] Godunov | Elements of Continuum Mechanics[END_REF][START_REF] Godunov | Elements of Continuum Mechanics and Conservation Laws[END_REF]. It is symmetric, if e is isotropic. Indeed, in the isotropic case, e depends only on the invariants of G :

J i = tr G i , i = 1, 2, 3
In particular, the determinant of G (denoted by |G|) can be expressed in terms of J i in the form

|G| = 2J 3 -3J 1 J 2 + J 3 1 6
It can be proved that

∂J i ∂G = iG i-1 , i = 1, 2 , 3.
In particular (see also the Lemma)

∂ |G| ∂G = |G| G -1 . Hence σ = -2ρ ∂e ∂G G = -2ρ ∂e ∂J 1 I + 2 ∂e ∂J 2 G + 3 ∂e ∂J 3 G 2 G = -2ρG ∂e ∂G = σ T
Here the density ρ is expressed by

ρ = ρ 0 |G| 1/2
EOS formulation in separable form It is natural to present the energy in a separable form : the energy is the sum of a "hydrodynamic" part and an "elastic" part : e = e h (ρ, η) + e e (g)

where

g = G |G| 1/3 .
The hydrodynamic part of the energy e h (ρ, η) can be taken in the form of stiffened gas equation of state :

e h = p + γp ∞ ρ (γ -1) , p + p ∞ = A exp η -η 0 c v ρ γ , η 0 = const, A = const (10) Here γ is the polytropic exponent, p ∞ is a constant, c v is the specific heat at constant volume.
The elastic part of the internal energy e e depends only on g. The tensor g has a unit determinant, so it is unaffected by the volume change. This idea to take the arguments of the internal energy in this form was first proposed by Gouin and Debieve (1986) (see also [START_REF] Plohr | Linearized analysis of Richtmyer-Meshkov flow for elastic materials[END_REF], but for the dependence of the energy on the right Cauchy-Green tensor. The simplest example of the elastic energy is

e e (g, η) = µ s 4ρ 0 tr (g -I) 2 = µ s 4ρ 0 J 2 |G| 2/3 - 2J 1 |G| 1/3 + 3 ( 11
)
where µ s is the shear modulus. The stress tensor will be then

σ = -2ρ ∂e ∂G G = -pI + S, tr(S) = 0, (12)
where the deviatoric part is

S = -µ s ρ ρ 0 1 |G| 2/3 G 2 - J 2 3 I - 1 |G| 1/3 G - J 1 3 I
and the thermodynamic pressure is

p = ρ 2 ∂e h ∂ρ .
The elastic part of the energy has no influence on the pressure, it is determined only by the hydrodynamic part. In the case of small displacements, we obtain the classical Hooke law.

Governing equations

The governing equations can be written in the form :

∂e α ∂t + ∇ (v • e α ) = 0, rote α = 0, ∂ρ ∂t + div (ρv) = 0, ∂ρv ∂t + div (ρv ⊗ v -σ) = 0, ∂ ρ e + v•v 2 ∂t + div ρv e + v • v 2 -σv = 0.
For the energy in the form (10), 11), the equations are hyperbolic at least in the one-dimensional case. The proof will be given in a more general situation involving solid-fluid mixtures.

Multiphase flow modeling: general definitions

Consider a multiphase mixture of two compressible components. Each ath component (a = 1, 2) has its own phase average characteristics: the velocity v a ; the density ρ a ; entropy η a , the internal energy per unit mass e a (ρ a , η a ), the temperature θ a etc. The apparent densities ρa are defined by : ρa = α a ρ a . The mixture density is :

ρ = α 1 ρ 1 + α 2 ρ 2 .
Here α a is the volume fraction of the a-th component,

α 1 + α 2 = 1.
The definition of the volume fraction is as follows. At each point (t, x) consider a representative mixture volume Ω : Ω = Ω 1 ∪ Ω 2 where Ω a is the volume occupied by the α -th phase. Consider the characteristic function of the a -th phase :

χ a (t, x ′ ) = 1, if (t, x ′ ) ∈ Ω a 0, if (t, x ′ ) / ∈ Ω a Then α a (t, x) = χ a = Ω χ a (t, x ′ ) dω ′ V ol(Ω) = V ol(Ω a ) V ol(Ω)
The phase average characteristics are defined in standard way :

ρ a (t, x) = Ωa ρ ′ a (t, x ′ ) dω ′ V ol(Ω a ) v a (t, x) = Ωa ρ ′ a (t, x ′ ) v ′ a (t, x ′ ) dω ′ Ωa ρ ′ a (t, x ′ ) dω ′
etc. Other types of average can be introduced: time average, space-time average etc. (see, for example, Drew and Passman (1998), [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow[END_REF]). The usual technique in obtaining the governing equations of multiphase flows is to multiply the system of balance laws by the characteristic function and average it. In particular, the mass conservation laws become :

∂ (α a ρ a ) ∂t + div (α a ρ a v a ) = 0.
However, applying the same averaging procedure to the conservation laws of momentum and energy, we immediately arrive to the closure problem : the balance equations contain more unknowns than equations. A possible remedy to that is to formulate the Lagrangian of the system directly in terms of average variables. And then to proceed in usual way by writing the corresponding Euler-Lagrange equations. The principal difficulty in applying such an approach is : how to determine this Lagrangian? At least two approaches can be used. The first is a phenomenological formulation of such a Lagrangian. The validation of such an approach should be established by comparing to experiments. The second possibility is an asymptotic modeling where the development with respect to a small parameter allows us to obtain an approximate Lagrangian from the exact one. Also, it is not sufficient just to obtain the Euler-Lagrange equations : dissipation terms consistent with the second law of thermodynamics should be introduced to describe real physics.

Equilibrium one-velocity model

Quite often one can suppose that the average velocities of phases are equal : v 1 = v 2 = v. We introduce the volume kinetic energy of the system :

T = 2 a=1 α a ρ a |v a | 2 2 = 2 a=1 α a ρ a |v| 2 2 = ρ |v| 2 2 ,
and the potential energy

W = 2 a=1 α a ρ a e a = ρ 2 a=1 Y a e a , Y a = α a ρ a ρ , Y 1 + Y 2 = 1.
The variable Y a is called mass fraction. The energies e a (ρ a , η a ) verify the Gibbs identity:

θ a dη a = de a + p a d( 1 ρ a )
where p a are the pressures. The Lagrangian of the mixture is :

L = D(t) ρ |v| 2 2 -e dD
where we assume the additivity of the mixture energy defined as :

e = Y 1 e 1 Y 1 ρ α 1 , η 1 + Y 2 e 2 Y 2 ρ α 2 , η 2 .
The variation of the Hamilton action can be found under the usual constraints :

∂ρ ∂t + div (ρv) = 0, DY 1 Dt = 0, Dη a Dt = 0, a = 1, 2.
The novelty is that a new independent variable appears, the volume fraction, which is not related to the virtual displacements. Hence, an additional variation should be taken. The variations of the density, mass fraction and entropies are given by : δρ = -div(ρδx), δη a = -∇η a • δx,

δY a = -∇Y a • δx.
The corresponding variation of the Lagrangian can be found in the same way as for classical compressible fluids :

δa = t1 t0 D(t) - ∂ρv ∂t + div (ρv ⊗ v) + ▽p •δxdD = 0 where p = α 1 p 1 + α 2 p 2
is the mixture pressure. It gives the momentum equation :

∂ρv ∂t + div (ρv ⊗ v) + ▽p = 0.
Doing so, we have supposed that the volume fraction is fixed. Now, it is necessary to find the variation of the Lagrangian with respect to α 1 :

δL = D(t) ρ -Y 1 ∂e 1 ∂ρ 1 - Y 1 ρ α 2 1 -Y 2 ∂e 2 ∂ρ 2 Y 2 ρ α 2 2 δα 1 dD = D(t) (p 1 -p 2 ) δα 1 dD = 0.
Since it is valid for any δα 1 , it follows from here that

p 1 -p 2 = 0.
The equilibrium of velocities implies also the pressure equilibrium : p = p 1 = p 2 . This is an algebraic equation for the volume fraction. The energy equation is a consequence of the mass conservation laws, the momentum equation and the entropy equations :

ρ e + |v| 2 2 t + div ρv e + |v| 2 2 + pv = 0. ( 13
)
6.1 What is the sound speed in multiphase flow models?

The sound speeds of pure phases are determined as :

c 2 a = ∂p a (ρ a , η a ) ∂ρ a ηa=const , a = 1, 2.
They correspond to the phase velocity of propagation of linear perturbations of the corresponding Euler equations linearized at rest. The corresponding sound speed c w in a two-fluid mixture is called Wood sound speed and is defined as

1 ρc 2 w = α 1 ρ 1 c 2 1 + α 2 ρ 2 c 2 2 . ( 14
)
Consider first an euristic method of derivation of this formula. The definition of the specific mixture volume is :

τ = Y 1 τ 1 + Y 2 τ 2 , where τ = 1 ρ , τ a = 1 ρ a , a = 1, 2.
Let us differentiate it with respect to the mixture pressure p = p 1 = p 2 :

dτ dp = - 1 ρ 2 c 2 w = Y 1 dτ 1 dp 1 + Y 2 dτ 2 dp 2 = - Y 1 ρ 2 1 c 2 1 - Y 2 ρ 2 2 c 2 2 .
It implies immediately [START_REF] Godunov | Elements of Continuum Mechanics and Conservation Laws[END_REF]. Another way, more rigorous one, is to calculate the characteristic speeds of the governing system. For simplicity, we consider only one-dimensional case (v = (u), x = (x)). The equations are :

Dτ Dt -u x = 0, Du Dt + p x ρ = 0, DY 1 Dt = 0, Dη 1 Dt = 0, Dη 2 Dt = 0.
The mean specific volume is given by :

τ = τ (p, η 1 , η 2 , Y 1 ) = Y 1 τ 1 (p, η 1 ) + Y 2 τ 2 (p, η 2 ) .
Obviously,

Dτ Dt = Y 1 Dτ 1 (p, η 1 ) Dt + Y 2 Dτ 1 (p, η 1 ) Dt = - Y 1 ρ 2 1 c 2 1 + Y 2 ρ 2 2 c 2 2 Dp Dt = - α 1 τ ρ 1 c 2 1 + α 2 τ ρ 2 c 2 2 Dp Dt .
Hence, the system is

Dp Dt + ρc 2 w u x = 0, (15)
Du Dt + p x ρ = 0, DY 1 Dt = 0, Dη 1 Dt = 0, Dη 2 Dt = 0.
20 System (15) can be rewritten in terms of unknowns U = (p, u, Y 1 , η 1 , η 2 ) T :

U t + A (U) U x = 0,
where

A (U) =       u ρc 2 w 0 0 0 1 ρ u 0 0 0 0 0 u 0 0 0 0 0 u 0 0 0 0 0 u      
The characteristic values of the matrix A are

λ 1,2 = u ± c w , λ 3,4,5 = u.
On can see that the system is hyperbolic (obviously, the eigenvalue λ = u of multiplicity three has exactly three eigenvectors).

The Wood sound speed c w has very interesting physical properties. Indeed, let us fix the phase densities and entropies and plot the Wood sound speed as a function of the volume fraction.

c w = 1 (α 1 ρ 1 + α 2 ρ 2 ) α1 ρ1c 2 1 + α2 ρ2c 2 2 .
Its minimal value is much lower than the sound speeds in either medium. For example, for the "air-water" mixture with c 1 = 330m/s, ρ 1 = 1, 24 kg/m 3 , c 2 = 1500m/s, ρ 2 = 1000 kg/m 3 the minimal sound speed is about 23m/s. This very surprising fact has also been validated experimentally [START_REF] Micaelli | Propagation d'ondes dans les écoulements diphasiques à bulles à deux constituants[END_REF].

The model [START_REF] Green | On the Theory of Water Waves[END_REF] we have derived is usually called "5 equations model" (see, for example, [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials : reduced equations[END_REF] for the derivation of this model from a more general two-velocity model). As we have seen, it is hyperbolic. However, it is not in conservative form. The last renders questionable a formal determination of the Rankine-Hugoniot relations for this model. A possible solution to this problem can be found in [START_REF] Saurel | Shock jump relations for multiphase mixtures with stiff mechanical relaxation[END_REF].

Nonequilibrium one-velocity model

The introduction of dissipation is always a phenomenological procedure which however should satisfy the entropy inequality. A natural relaxation equation can be written to replace the algebraic relation p 1 -p 2 = 0. Introducing the dissipation function

D = 1 2µ Dα 1 Dt Figure 1.
The minimal value of the Wood sound speed is much lower than the sound speeds in either components.

one can replace the equilibrium condition by the following one :

p 1 -p 2 = ∂D ∂ Dα1 Dt = 1 µ Dα 1 Dt .
Here µ is some constant. The physical sense of this equation is clear : the volume fraction changes due to the pressure difference. This equation is the simplest one. In more general case it may also contain inertia terms (see [START_REF] Gavrilyuk | Mathematical and Numerical Modeling of Two-Phase Compressible Flows with Microinertia[END_REF] and Saurel, Gavrilyuk and Renaud (2003) for details). The equations for entropies should satisfy the entropy inequality in the form :

∂ (ρη) ∂t + div (ρηv) ≥ 0 ( 16
)
where η is the mixture entropy defined by

ρη = α 1 ρ 1 η 1 + α 2 ρ 2 η 2 .
We take the equations for the entropies in the form :

α a ρ a θ a Dη a Dt = f a , a = 1, 2
where f a should be chosen to satisfy the entropy inequality [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF]. The production terms f a should also be compatible with the energy equation [START_REF] Godunov | Elements of Continuum Mechanics[END_REF].

Developing the energy equation one can obtain

De Dt + p D Dt 1 ρ = 0.
Or

α 1 ρ 1 θ 1 Dη 1 Dt + α 2 ρ 2 θ 2 Dη 2 Dt = µ (p 1 -p 2 ) 2 ≥ 0. ( 17
)
To satisfy [START_REF] Iordansky | On the equations of motion of the liquid containing gas bubbles[END_REF] one can choose

α a ρ a θ a Dη a Dt = f a = (p a -p I ) Dα a Dt , a = 1, 2. ( 18
)
where the interface pressure p I can be taken as :

p I = β 2 p 1 + β 1 p 2 , β 1 + β 2 = 1, β a ≥ 0. ( 19
)
Here β i are some constants. In particular, the interface pressure p I can be taken as the pressure at the contact discontinuity obtained as the solution of the linearized Riemann problem with initial pressures p 1 , p 2 (Saurel, [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF]. For that case

β 2 = Z 2 Z 1 + Z 2 , β 1 = Z 1 Z 1 + Z 2 ,
where Z a are acoustical impedances of pure phases. With such a choice, the entropies of each phase increase :

α 1 ρ 1 Dη 1 Dt = 1 θ 1 (p 1 -p I ) Dα 1 Dt = µβ 1 θ 1 (p 1 -p 2 ) 2 ≥ 0, α 2 ρ 2 Dη 2 Dt = 1 θ 2 (p 2 -p I ) Dα 2 Dt = µβ 2 θ 2 (p 1 -p 2 ) 2 ≥ 0.
Summing these equations we obtain the entropy inequality [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF]. Finally, the corresponding 1D non-equilibrium model is :

∂ρ ∂t + ∂ (ρu) ∂x = 0, (20)
Du Dt + 1 ρ ∂p ∂x = 0, p = α 1 p 1 Y 1 ρ α 1 , η 1 + α 2 p 2 Y 2 ρ α 2 , η 2 , Dα 1 Dt = µ (p 1 -p 2 ) , DY 1 Dt = 0, Dη 1 Dt = µβ 1 α 1 ρ 1 θ 1 (p 1 -p 2 ) 2 , Dη 2 Dt = µβ 2 α 2 ρ 2 θ 2 (p 1 -p 2 ) 2 .
The equations are also hyperbolic, however the characteristic speeds are now :

λ 1,2 = u ± c f , λ k = u, k = 3, 4, 5, 6.
Here c f is the "frozen" sound speed :

c 2 f = Y 1 c 2 1 + Y 2 c 2 2 . One can prove that c 2 f > c 2 w . (21)
When µ → ∞, we recover in this asymptotic limit the equilibrium model [START_REF] Green | On the Theory of Water Waves[END_REF]. The condition ( 21) is known as Whitham's stability condition (or Whitham's subcharacteristic condition) [START_REF] Whitham | Linear and Nonlinear Waves[END_REF]. The Wood sound speed and the "frozen" sound speed are the velocities of propagation of long and short waves, respectively, for [START_REF] Karni | Multi-component flow calculations by a consistent primitive algorithm[END_REF].

Bubbly fluids

In applications (for example, shock wave propagation in bubbly fluids, high velocity impacts etc.) it is necessary to take into account inertia effects. In particular, in such media, the internal energy of a mixture is a function not only of the flow parameters but also their derivatives. Consider a particular case of fluids containing gas bubbles ("bubbly fluids"). For simplicity, we will suppose that the liquid phase is incompressible : only gas bubbles are compressible. We will use the index "1" for a liquid and the index "2" for a gas. The kinetic energy of an incompressible fluid of density ρ 10 = const due to single bubble oscillations is [START_REF] Iordansky | On the equations of motion of the liquid containing gas bubbles[END_REF], [START_REF] Kogarko | On the model of cavitating liquid[END_REF], Wijngaarden (1968)) :

2πR 3 ρ 10 dR dt 2
where R is the bubble radius. If N 2 is the number of bubbles per unit volume of the mixture satisfying the conservation law

∂N 2 ∂t + div (N 2 v) = 0,
the corresponding energy per unit volume of N 2 moving bubbles (having the same mean velocity as a surrounding fluid) will be :

2πR 3 N 2 ρ 10 DR Dt 2 , D Dt = ∂ ∂t + v • ∇,
The bubble volume fraction is :

α 2 = 4 3 πR 3 N 2 .
The Lagrangian of the bubbly fluid can be written in the form :

L = D(t) ρ |v| 2 2 + 3 2 α 2 ρ 10 DR Dt 2 -ρY 2 e 2 Y 2 ρ α 2 , η 2 dD (22)
We consider the case where the mass concentration of bubbles Y 2 = α2ρ2 ρ , the number of bubbles per unit mass n 2 = N2 ρ and the gas entropy η 2 are constant. The bubble radius and the gas density can be expressed through the average density ρ as :

4 3 πR 3 = 1 ρ - Y 1 ρ 10 /n 2 , ρ 2 = Y 2 1 ρ - Y 1 ρ 10 -1 . ( 23
)
We rewrite [START_REF] Li | Linear stability of solitary waves of the Green-Naghdi equations[END_REF] in the following generic form :

L = D(t) ρ |v| 2 2 -W ρ, Dρ Dt dD (24)
In our particular case

W ρ, Dρ Dt = ρ Y 2 e 2 (ρ 2 , η 2 ) -2πR 3 n 2 ρ 10 DR Dt 2
The usual constraint should be respected :

∂ρ ∂t + div (ρv) = 0.
Let us define the variation of such a generic Lagrangian [START_REF] Makarenko | A second long-wave approximation in the Cauchy-Poisson problem[END_REF]. Obviously,

δa = t1 t0 D(t)  δ ρ |v| 2 2 - δW δρ + ρv• δv - ∂W ∂ Dρ Dt δ Dρ Dt   dDdt = t1 t0 D(t)     ρv- ∂W ∂ Dρ Dt ∇ρ   • δv+ δρ |v| 2 2 - δW δρ   dDdt + t1 t0 D(t)   - ∂ ∂t   ∂W ∂ Dρ Dt δρ   -div   ∂W ∂ Dρ Dt v δρ     dDdt = t1 t0 D(t)   |v| 2 2 - δW δρ δρ +   ρv- ∂W ∂ Dρ Dt ∇ρ   • δv   dDdt + t1 t0 D(t)   - ∂ ∂t   ∂W ∂ Dρ Dt δρ   -div   ∂W ∂ Dρ Dt v δρ     dDdt = t1 t0 D(t)   |v| 2 2 - δW δρ δρ +   ρv- ∂W ∂ Dρ Dt ∇ρ   • δv   dDdt.
since the conservative terms are vanishing at the boundary ∂ ([t 0 , t 1 ] × D(t))

. Here we have introduced the classical notation for the variational derivative of W : Finally, the momentum equation is :

δW δρ = ∂W ∂ρ - ∂ ∂t   ∂W ∂ ∂ρ ∂t   -div ∂W ∂ (▽ρ) v = ∂W ∂ρ - ∂ ∂t   ∂W ∂ Dρ Dt   -div   ∂W ∂ Dρ Dt v   = ∂W ∂ρ -ρ D Dt   1 ρ ∂W ∂ Dρ
ρ Dv Dt + ∇ ρ δW δρ -W = 0.
The governing equations can also be rewritten in conservative form :

∂ρ ∂t + div (ρv) = 0, (25)
∂ρv ∂t + div (ρv ⊗ v) + ∇p = 0. ( 26
)
where the pressure is defined by

p = ρ δW δρ -W. ( 27
)
It admits also the energy conservation law :

∂ ∂t ρ |v| 2 2 + E + div v ρ |v| 2 2 + E + p = 0 (28)
where

E = W - Dρ Dt ∂W ∂ Dρ Dt .
If W depends only on ρ, then W = E and equations ( 25), ( 26) coincide with the classical Euler equations of barotropic fluids.

Let us calculate the pressure equation ( 27) in explicit form for the bubbly fluids. We introduce the specific potential w :

W = ρw, w = Y 2 e 2 (ρ 2 , η 2 ) -2πR 3 n 2 ρ 10 DR Dt 2 .
Obviously,

p = ρ 2   ∂w ∂ρ - D Dt   ∂w ∂ Dρ Dt     .
Consider now w ρ, Dρ Dt as a function of R and DR Dt . Then

p = - 1 4πR 2 n 2 ∂w ∂R - D Dt ∂w ∂ DR Dt = - 1 4πR 2 n 2 -p 2 3 4π 3 R 3 n 2 2 4π 3 R 4 n 2 -6πR 2 n 2 ρ 10 DR Dt 2 - 1 4πR 2 n 2 D Dt 4πR 3 n 2 ρ 10 DR Dt = - 1 4πR 2 n 2 -p 2 3 4π 3 R 3 n 2 2 4π 3 R 4 n 2 -6πR 2 n 2 ρ 10 DR Dt 2 +4πR 3 n 2 ρ 10 D 2 R Dt 2 + 12πR 2 n 2 ρ 10 DR Dt 2 = p 2 + ρ 10 3 2 DR Dt 2 + R D 2 R Dt 2 .
The equation

p = p 2 + ρ 10 3 2 DR Dt 2 + R D 2 R Dt 2 (29)
is known as the Rayleigh-Lamb equation describing oscillations of a single bubble in an infinite fluid. Model ( 25), ( 26) and ( 29) is called in the literature Iordansky-Kogarko-Wijngaarden model [START_REF] Iordansky | On the equations of motion of the liquid containing gas bubbles[END_REF], [START_REF] Kogarko | On the model of cavitating liquid[END_REF], van Wijngaarden (1968)).

Dispersive shallow water flows as an analogue of bubbly flows

Much work has been done in the derivation of relatively simple mathematical models of long non-linear water waves. One of popular models is the Green-Naghdi model obtained in the one-dimensional case by [START_REF] Su | Korteweg-de Vries Equation and Generalizations. III. Derivation of the Korteweg-de Vries Equation and Burgers Equation[END_REF] and in the multi-dimensional case by [START_REF] Green | On the Theory of Water Waves[END_REF]Naghdi (1974, 1976) within the context of a homogeneous one-layer fluid. In the literature, this model is usually called Green-Naghdi model (GN model or GN system). A derivation of the GN model based on the variational formulation of the Euler equations was done by [START_REF] Miles | Weakly dispersive nonlinear gravity waves[END_REF] (see also [START_REF] Salmon | Hamiltonian Fluid Mechanics[END_REF][START_REF] Salmon | Lectures on Geophysical Fluid Dynamics[END_REF]). A mathematical justification of the GN model was done by [START_REF] Makarenko | A second long-wave approximation in the Cauchy-Poisson problem[END_REF] and Alvarez-Samaniego and Lannes (2008). Solitary wave solutions of the GN model were obtained by [START_REF] Su | Korteweg-de Vries Equation and Generalizations. III. Derivation of the Korteweg-de Vries Equation and Burgers Equation[END_REF]. The linear stability of solitary waves has been proved by [START_REF] Li | Linear stability of solitary waves of the Green-Naghdi equations[END_REF]. A criterium of stability of shear flows for the GN model was proposed by [START_REF] Gavrilyuk | Linear stability of parallel inviscid flows of shallow water and bubbly fluid[END_REF]. A wide class of multi-dimensional solutions and approximate solutions of nonlinear multi-dimensional GN model has been found by [START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water equations[END_REF] and [START_REF] Teshukov | Three-Dimensional Nonlinear Dispersive waves on Shear Flows[END_REF]. Unsteady undular bores were described by [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF].

It is interesting to note that such a model can be viewed as an analogue of "bubbly fluids" flow model. The Green-Naghdi equations describing shallow water flows over a plane bottom are :

∂h ∂t + div (hv) = 0, ∂hv ∂t + div hv ⊗ v+ gh 2 2 + h 2 3 D 2 h Dt 2 I 2 = 0.
Here h is the position of free surface, v is the average with respect the vertical direction 2D velocity field, g is the gravity acceleration, and I 2 is the two-dimensional identity transformation. The Green-Naghdi model is formally described with the same type of Lagrangian :

L = D(t) h |v| 2 2 -W h, Dh Dt dD
where the function W is :

W (h, Dh Dt ) = gh 2 2 - h 6
Dh Dt 2 .
Obviously,

p = h δW δh -W = gh 2 2 + h 2 3 D 2 h Dt 2 .
Such an analogy between "bubbly fluids" and dispersive shallow water flows is very useful. Indeed, the parameters of waves propagating in bubbly fluids can not be observed, they can only be measured, while the position of the free surface (playing the role of average density) can be at least easily observed.

9 Two-velocity one-pressure model Consider the case where the averaging velocities of phases are not equal. The Lagrangian of the mixture is :

L = D 2 a=1 α a ρ a |v a | 2 2 - 2 a=1 α a ρ a e a dD
Here D is not a material volume, this notion has no sense in the case of two velocity continua. We will use the Eulerian variations to define the governing equations. The imposed constraints are

(α a ρ a ) t + div (α a ρ a v a ) = 0, (α a ρ a η a ) t + div (α a ρ a v a η a ) = 0.
Or, in terms of apparent densities ρa :

(ρ a ) t + div (ρ a v a ) = 0, (ρ a η a ) t + div (ρ a v a η a ) = 0.
At least three types of independent variations should be considered. The first one is the variation with respect to the volume fraction giving us the equilibrium condition :

p 1 = p 2 = p, p = α 1 p 1 + α 2 p 2 .
The second and the third variations in terms of virtual displacements δx a (which can be introduced in the same way as in the case of one-velocity continua) give us usual expressions of variations of the density, entropy and velocity :

δ ρa = -div(ρ a δx a ), δη a = -∇η a • δx a , δv a = D a δx a Dt - ∂v a ∂x δx a , D a Dt = ∂ ∂t + v a • ∇
The Euler-Lagrange equations correspond to the classical momentum equations for each phase :

ρa D a v a Dt + α a ∇p a = 0.
Taking into account the equilibrium condition they can also be rewritten in non-conservative form :

∂α 1 ρ 1 v 1 ∂t + div (α 1 ρ 1 v 1 ⊗ v 1 + α 1 p 1 I) = p∇α 1 , ∂α 2 ρ 2 v 2 ∂t + div (α 2 ρ 2 v 2 ⊗ v 2 + α 2 p 2 I) = p∇α 2 .
Summing them we obtain the equation of the total momentum :

2 a=1 α a ρ a v a t + div 2 a=1 α a ρ a v a ⊗ v a + pI = 0.
The total energy is :

2 a=1 α a ρ a e a + |v a | 2 2 t +div 2 a=1 α a ρ a v a e a + |v a | 2 2 + α a p a v a = 0.

Analysis of the two-velocity one-pressure model

For simplicity, consider the one-dimensional isentropic case ( η a = const ). The equations are :

(α 1 ρ 1 ) t + (α 1 ρ 1 u 1 ) x = 0, (30)
(α 2 ρ 2 ) t + (α 2 ρ 2 u 2 ) x = 0, D 1 u 1 Dt + p x ρ 1 = 0, D 2 u 2 Dt + p x ρ 2 = 0.
The pressures are then functions of densities:

p 1 ρ1 α 1 = p 2 ρ2 α 2 = p. (31)
Using the expression for the differential of the volume fraction obtained from the equilibrium condition ( 31)

dα 1 = c 2 1 α1 dρ 1 - c 2 2 α2 dρ 2 c 2 1 ρ1 α 2 1 + c 2 2 ρ2 α 2 2 = c 2 1 α1 dρ 1 - c 2 2 α2 dρ 2 c 2 1 ρ1 α1 + c 2 2 ρ2 α2
, and rewriting the equations in variables (ρ 1 , ρ2 , u 1 , u 2 ), one can obtain the following characteristic polynomial for ( 30) :

(λ -u 1 ) 2 (λ -u 2 ) 2 - α 2 ρc 2 w ρ 2 (λ -u 1 ) 2 - α 1 ρc 2 w ρ 1 (λ -u 2 ) 2 = 0.
Transforming it to the form

1 = α 2 ρc 2 w ρ 2 1 (λ -u 2 ) 2 + α 1 ρc 2 w ρ 1 1 (λ -u 1 ) 2
and looking for the minimal value of the right-hand side, one obtain the following criterion of hyperbolicity (existence of four real eigenvalues):

(u 2 -u 1 ) 2 > ρc 2 w α 1 ρ 1 1/3 + α 2 ρ 2 1/3 3
.

However, this model is not hyperbolic for small relative velocity. One can prove that the corresponding volume energy

E = α 1 ρ 1 e 1 + α 2 ρ 2 e 2 = ρ1 e 1 + ρ2 e 2
has the following differential :

dE = h 1 dρ 1 + h 2 dρ 2 .
where h i , i = 1, 2 are the phase enthalpies :

dh i = dp i ρ i = dp ρ i . Also, ∂ 2 E ∂ ρ2 1 = ρc 2 w ρ 2 1 , ∂ 2 E ∂ ρ2 2 = ∂h 2 ∂ ρ2 = c 2 2 ρ 2 1 α 2 + ρ2 α 2 2 ∂α 1 ∂ ρ2 = ρc 2 w ρ 2 2 , ∂ 2 E ∂ ρ1 ∂ ρ2 = ρc 2 w ρ 1 ρ 2 .
In particular this implies that the corresponding Hessian matrix is degenerate:

det

∂ 2 E ∂ ρ2 1 ∂ 2 E ∂ ρ1∂ ρ2 ∂ 2 E ∂ ρ1∂ ρ2 ∂ 2 E ∂ ρ2 2 = 0.
In mechanics, the non-convexity of the internal energy shows up nonstability of the mechanical system. In our case, the degeneracy of the internal energy is responsible for the non-hyperbolicity of the governing equations and, as a consequence, the ill-posedness of the corresponding Cauchy problem.

Two-velocity two-pressure model

The following non-equilibrium model obtained by adding relaxation terms in the pressure equilibrium model is called Baer-Nunziato model (1986) (or BN model) :

∂α 1 ρ 1 ∂t + div (α 1 ρ 1 v 1 ) = 0, ( 32
)
∂α 2 ρ 2 ∂t + div (α 2 ρ 2 v 2 ) = 0, ∂α 1 ρ 1 v 1 ∂t + div (α 1 ρ 1 v 1 ⊗ v 1 + α 1 p 1 I) = p I ∇α 1 + λ (v 2 -v 1 ) , ∂α 2 ρ 2 v 2 ∂t + div (α 2 ρ 2 v 2 ⊗ v 2 + α 2 p 2 I) = p I ∇α 2 -λ (v 2 -v 1 ) , ∂ ∂t α 1 ρ 1 |v 1 | 2 2 + e 1 + div α 1 ρ 1 v 1 |v 1 | 2 2 + e 1 + α 1 p 1 v 1 = = -p I ∂α 1 ∂t + λv 1 • (v 2 -v 1 ) , ∂ ∂t α 2 ρ 2 |v 2 | 2 2 + e 2 + div α 2 ρ 2 v 2 |v 2 | 2 2 + e 2 + α 2 p 2 v 2 = = -p I ∂α 2 ∂t -λv 1 • (v 2 -v 1 ) , D I α 1 Dt = µ(p 1 -p 2 ), D I Dt = ∂ ∂t + v I • ∇.
Here p I is the "interface" pressure, v I is the "interface" velocity, and λ is the friction coefficient. In the classical BN model the following duality is supposed p I = p 2 and v I = v 1 , or p I = p 1 and v I = v 2 . The right-hand side in the momentum equations is the sum of "nozzling terms" p I ∇α a ( the word "nozzling" is to underline the analogy between these equations and those describing gas flow in the duct of variable cross section) and Stokes type terms λ (v 2v 1 ). Other choice of "interface" variables is possible (see, for example, Saurel, Gavrilyuk and Renaud, 2003). Such a system of equations satisfies the entropy inequality and has only real eigenvalues.

Equilibrium diffuse interface model of solid-fluid interactions

Pioneering works by [START_REF] Karni | Multi-component flow calculations by a consistent primitive algorithm[END_REF], [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations : A quasi-conservative approach[END_REF] and Saurel and Abgrall (1999) have shown the attractivity the diffuse interface approach for modelling interfaces between ideal fluids having different thermodynamic characteristics. In such an approach the interface between different materials is considered as a diffusion zone : it is a mixture of two components where the volume fraction of the phases varies continuously from zero to one. To assure the non-degeneracy of the governing equations, one also supposes that even in the bulk (where, formally, we have only pure phases ) a negligible quantity of the other phase is present. Hence, the interface can be viewed as an effective fluid mixture. We will generalize this approach to the case of solid-fluid interfaces by using Hamilton's principle [START_REF] Gavrilyuk | Modeling wave dynamics of compressible elastic materials[END_REF], [START_REF] Favrie | Diffuse solid-fluid interface model in cases of extreme deformations[END_REF]).

Consider the following Lagrangian for a solid-fluid mixture :

L = Ωt ρ v • v 2 -e dΩ
where the average density and the average volume energy are defined as ρ = α s ρ s + α g ρ g , ρe = α s ρ s e s + α g ρ g e g .

The indices "s" and "g" mean "solid" and "gas", respectively. The mass conservation laws can be rewritten in terms of the mass fractions Y s and Y g :

DY s Dt = 0, DY g Dt = 0, Y s = α s ρ s ρ , Y g = α g ρ g ρ , Y s + Y g = 1.
Also, the entropies of each phase are conserved :

Dη s Dt = 0, Dη g Dt = 0.
The energy of the solid phase e s is a function of the entropy η s and the mixture deformation tensor

G = 3 β=1 E β ⊗ E β where E β = ∇X β
and X = (X β ) are the Lagrangian coordinates of a mixture particle. The mixture density is related to G by the formula :

ρ(t, X) = ρ 0 (X) |G| 1/2
We use here capital letters E β to distinguish the mixture case from the pure solid one. The gas energy e g is a function of the gas density ρ g and the gas entropy η s .

To apply the Hamilton principle we find the variations of ρ, v, Y s , Y g , η g , η s as functions of virtual displacements δx at fixed Lagrangian coordinates : δρ = -ρdiv (δx) , δv = ∂δx ∂t ,

δη i = 0, δY i = 0. ( 33
) δE β = - ∂δx ∂x T E β .
As usually, the symbol δ means the variation at fixed Lagrangian coordinates. The last formula giving the variation of E β can be found in the following way. There exist potentials X β , β = 1, 2, 3 (the mixture Lagrangian coordinates) such that

DX β Dt = 0.
It follows from the definition of X β that its Lagrangian variation is zero :

δX β = 0. Since δF -1 = -F -1 δF F -1 = -F -1 ∂δx ∂X F -1 = -F -1 ∂δx ∂x
we get from the definition of E β :

δE β = δ∇X β = δ F -T ∇ X X β = -F -1 ∂δx ∂x T ∇ X X β = - ∂δx ∂x T F -T ∇ X X β = - ∂δx ∂x T E β .
In particular, the variation of G =

3 β=1 E β ⊗ E β is given by : δG = δ   3 β=1 E β ⊗ E β   = - ∂δx ∂x T G -G ∂δx ∂x . (34)
The energy of the solid is taken in separable form

e s =
Analogous considerations give a simpler variation for the gas phase :

δe g = - 1 ρ g tr σ g ∂δx ∂x - p g α g ρ g δα s = p g ρ g divδx - p g α g ρ g δα g (36)
Now, we are ready to take the variation of the Hamilton action. We use formulas [START_REF] Su | Korteweg-de Vries Equation and Generalizations. III. Derivation of the Korteweg-de Vries Equation and Burgers Equation[END_REF], [START_REF] Van Wijngaarden | On the equations of motion for mixtures of liquid and gas bubbles[END_REF], [START_REF] Whitham | Linear and Nonlinear Waves[END_REF], the mass conservation law ρdΩ = ρ 0 (X) dΩ 0 to present the variation of a in the form :

0 = δa = δ t2 t1 Ωt ρ v • v 2 -e dΩdt = δ t2 t1 Ω0 ρ 0 (X) v • v 2 -e dΩ 0 dt = t2 t1 Ω0 ρ 0 (X) δv • v -δe dΩ 0 dt = t2 t1 Ω0 ρ 0 (X) ∂δx ∂t •v + Y s 1 ρ s tr σ s ∂δx ∂x + p s α s ρ s δα s dΩ 0 dt - t2 t1 Ω0 ρ 0 (X) Y g p g ρ g divδx - p g α g ρ g δα g dΩ 0 dt = t2 t1 Ω0 ρ 0 (X) ∂δx ∂t •v dΩ 0 dt + t2 t1 Ω0
ρ 0 (X) ρ tr (α s σ s -α g p g ) ∂δx ∂x + p s δα s + p g δα g dΩ 0 dt

Since δα s + δα g = 0, it implies the equilibrium condition :

p s -p g = 0. (37)
Here

p s = - 1 3 tr (σ s ) , p g = - 1 3 tr (σ g ) .
Turning back to the Eulerian coordinates we transform the last variation to the following one :

0 = δa = t2 t1
∂ρv ∂t + div (ρv ⊗ v -(α s σ s + α g σ g )) = 0. ( 38
)
The entropy equations

Dη s Dt = 0, Dη g Dt = 0
complemented by the mass and momentum balance laws imply the mixture energy conservation law:

∂ ∂t ρ e + v • v 2 + div ρv e + v • v 2 -(α s σ s + α g σ g ) v = 0. (39)
Finally, the equilibrium solid-gas governing equations are

∂E β ∂t + ∇ E β • v = -rotE β ∧ v, (40a)
rotE β = 0, (40b)
∂ (α g ρ g ) ∂t + div (α g ρ g v) = 0, (40c)
∂ρv ∂t + div (ρv ⊗ v -(α s σ s + α g σ g )) = 0, (40d)
Dη s Dt = 0, Dη g Dt = 0, (40e)
p s = p g . (40f)
Equations ( 40) admit the energy conservation law (39).

12 Non-equilibrium diffuse interface model of solid-fluid interactions

For numerical reasons it is preferable to use a non-equilibrium model where the algebraic equation for the volume fraction p s -p g = 0 is replaced by a differential equation. This procedure will be done in the same way as for fluid-fluid mixtures. The relaxation equation will be :

p s -p g = 1 µ
Dα s Dt Now, we postulate the mixture energy equation :

∂ ∂t ρ e + v • v 2 + div ρv e + v • v 2 -(α s σ s + α g σ g ) v = 0. (41)
As in the case of fluid-fluid mixtures, we can transform the energy equation (41) to the following one : As in the case of fluid mixtures, the same equations ( 18) for the entropy production can be used :

α s ρ s θ s Dη s Dt + α g ρ g θ g Dη g Dt = (p s -p g ) Dα s Dt . ( 42
) Indeed, 0 = ∂ ∂t ρ e + v • v 2 + div ρv e + v • v 2 -(α s σ s + α g σ g ) v = ρ De Dt -tr (α s σ s + α g σ g ) ∂v ∂x = ρ
α s ρ s θ s Dη s Dt = (p s -p I ) Dα s Dt , α g ρ g θ g Dη g Dt = (p g -p I ) Dα s Dt .
Obviously, they are compatible with (42). Finally, the non-equilibrium model is :

∂ E β ∂t + ∇ E β • v = -rotE β ∧ v, (43a)
∂ (α g ρ g ) ∂t + div (α g ρ g v) = 0, (43b)
∂ρv ∂t + div (ρv ⊗ v -(α s σ s + α g σ g )) = 0, (43c)
α s ρ s θ s Dη s Dt = (p s -p I ) Dα s Dt , α g ρ g θ g Dη g Dt = -(p g -p I ) Dα s Dt , (43d)
Dα s Dt = µ (p s -p g ) , µ > 0. ( 43e
)
The geometric equations (43a) and the entropy equations (43d) can also be rewritten in different forms which are more convenient for numerical computations :

∂E β ∂t + ∂E β ∂x v+ ∂v ∂x T E β = 0, (44a)
∂ ∂t (α s ρ s e s ) + div (α s ρ s e s v) -α s tr σ s ∂v ∂x = -p I µ (p s -p g ) , (44b)
∂ ∂t (α g ρ g e g ) + div (α g ρ g e g v) -α g tr σ g ∂v ∂x = -p I µ (p g -p s ) . (44c)

Hyperbolicity of the non-equilibrium model

Let us show that the non-equilibrium model ( 43) is hyperbolic. We present here the proof of hyperbolicity in the one-dimensional case. We denote v = (u, v, w)

T . Let E β = (A (β) , B (β) , C (β) ). Relation rotE β = 0 implies that only components A (β) vary, the other components are constant. We can always suppose that initially the vectors E β coincide with the Cartesian basis. Hence,

B (2) = C (3) = 1 , B (1) = B (3) = C (1) = C (2) = 0. The matrix G becomes : G =   A (1) 2 + A (2) 2 + A (3) 2 A (2) A (3) A (2) 1 0 A (3) 0 1   .
Obviously,

|G| = A (1) 2 , J 1 = tr(G) = A (1) 2 + A (2) 2 + A (3) 2 + 2, J 2 = tr(G 2 ) = A (1) 2 + A (2) 2 + A (3) 2 2 +2 A (2) 2 + A (3) 2 + 2.
We remark that the invariants J 1 and J 2 depend only on A (1) and Z = A (2) 2 + A (3) 2 . The one-dimensional governing equations become :

∂A (1) ∂t + ∂ A (1) u ∂x = 0, (45a)
∂A (2) ∂t + ∂ A (2) u + v ∂x = 0, (45b)
∂A (3) ∂t + ∂ A (3) u + w ∂x = 0, (45c)
∂ρ ∂t + ∂ (ρu) ∂x = 0, ( 45d
)
∂ρY g ∂t + ∂ (ρY g u) ∂x = 0, (45e)
∂ρu ∂t + ∂ ρu 2 -α s σ s11 -α g σ g11 ∂x = 0, (45f)
∂ρv ∂t + ∂ (ρuv -α s σ s12 -α g σ g12 ) ∂x = 0, (45g)
∂ρw ∂t + ∂ (ρuw -α s σ s13 -α g σ g13 ) ∂x = 0, ( 45h
)
∂ρα g ∂t + ∂ (ρα g u) ∂x = ρµ (p g -p s ) , µ > 0, (45i)
∂ (ρη s ) ∂t + div (ρη s v) = (p s -p I ) Y s θ s Dα s Dt , (45j)
Let us recall a general result proved by Godunov, Friedrichs and Lax : if a system of conservation laws

V t + G(V) x = 0 (48)
admits an additional conservation law ϕ (V) t + ψ (V) x = 0 where ϕ (V) is a convex function of V, the system (48) is hyperbolic (see [START_REF] Godunov | Elements of Continuum Mechanics and Conservation Laws[END_REF] for detail). In our case this function ϕ is the total volume energy defined from ( 46) :

E = ρ e + u 2 + v 2 + w 2 2
and V = A (1) , A (2) , A (3) , ρ, ρY g , ρu, ρv, ρw, ρα g , ρη s , ρη g T . Let us remark that the equations for the volume and mass fractions, and the entropies evolve along contact characteristics. Hence, for hyperbolicity it is sufficient to check the convexity of the energy with respect to a lower number of variables V = A (1) , A (2) , A (3) , ρ, ρu, ρv, ρw T . The total energy is in the form :

E = E k (ρ, ρu, ρv, ρw) + E h (ρ, α g , Y g , η g , η s ) + E e s A (1) , A (2) , A (3) where the kinetic, hydrodynamic and elastic parts of the volume energy are defined as : (1) , A (2) , A (3) = ρY s e e s (g) = Y s ρ 0 |G| 1/2 e e s (g) .

E k (ρ,
The energy E is a convex function of V, if

E k + E h ρ = u 2 + v 2 + w 2 2 + Y g ε h g (τ g , η g ) + Y g ε h s (τ s , η s )
is convex with respect to (τ = 1/ρ, u, v, w), and

E e s = Y s ρ 0 |G| 1/2 e e s G |G| 1/3
is convex with respect to A (β) . The function E k + E h /ρ is convex with respect (τ = 1/ρ, u, v, w) if the hydrodynamic energies of pure phases e h g and e h s are convex with respect to τ g and τ s , respectively. This is the case, for example, of the stiffened gas equation of state we use for applications.

We prove the convexity of E e s in a special case where e e s (g) = µ s 4ρ s0 tr (g -I)

2 = µ s 4ρ s0 J 2 |G| 2/3 - 2J 1 |G| 1/3 + 3 . ( 49
)
Here ρ s0 is the solid reference density. Let us show the convexity of E e s in the case (49). This energy is in the form : If the energy E e s is convex with respect to A (1) and Z , and ∂E e s /∂Z > 0, it will be convex with respect to A (1) , A (2) and A (3) . Obviously, -A (1) 2/3 + 1 > 0 since for any positive x we have x 3 -x + 1 > 0. Finally, to prove the convexity of E e s with respect to variables A (1) , Z we have to prove that the corresponding 2 × 2 Hessian matrix is positive definite in all domain of parameters Z > 0, A (1) The proof can be done by direct calculations. The symmetrization method assures the hyperbolicity but does not give the characteristic eigenvalues in explicit form.

Applications

We made a comparison with experimental results on cylindrical gel samples (cylindrical drops of jelly-like materials) impacting a rigid surface (see Figure 2). Jelly-like materials remain elastic even for extreme deformations (between 100 % and 200 % ). This regime was recently studied experimentally by Luu Lia-Hua and Forterre (2009). A cylindrical gel sample (carbopol) impacting a rigid hydrophobic surface, first spreads until some limit size and then a full elastic recoil is observed which may be followed by a complete rebound. We reproduced this quasi-reversible behavior numerically. The separable form of the energy equation was also used with µ s = 85P a. The comparison during the spreading phase at a given time instant is shown in Figure 3.

Air Drop of jelly

  div (v) = 0 can be written in the form ρ = ρ 0 (X) , det F = 1.

e 2 , e 3 = 3 α=1 e α ⊗ e α , G - 1 = 3 α=1e

 3313 α j = (e α ) j = ∂X α ∂x j . The scalar product of these vectors satisfies e β •e α = δ β α where δ β α are the Kronecker symbols. With these definitions G = α ⊗ e α .

2 + 2Z 2 + 2 A ( 1 ) 1 / 3 - 2 A ( 1 ) 1 / 3 A ( 1 ) 2 + Z 2 + 2 + 3 where Z 2 =

 221132113122232 A (2) 2 + A (3) 2.

Figure 2 .

 2 Figure 2. A cylindrical drop of jelly is impacting a rigid surface.

Figure 3 .

 3 Figure 3. Comparison of the spreading phase of the cylinder impacting a rigid surface. Experiments by L. -H. Luu and Y. Forterre (2009) are shown on the top, computed results are shown on the bottom.

  )

				+	t1 t0	D(t)	ρv•	∂δx ∂t	+	∂δx ∂x	v +	∂e(ρ, η) ∂η	∇η • δx dD
							t1
						=
							t0		D(t)
					δρ	|v| 2	2	-e -	p ρ	+ ρv• δv-	∂e(ρ, η) ∂η	δη dD
		= -	t1 t0	D(t)	div(ρδx)	|v| 2	2	-e +	p ρ	dD
	+	t1 t0	D(t)	ρv•	Dδx Dt	-	∂v ∂x	δx +	∂e(ρ, η) ∂η	∇η • δx dD
			=	t1 t0		D(t)	ρδx • ▽	|v| 2	2	-e +	p ρ	dD
	+	t1 t0	D(t)	ρv•	Dδx Dt	-	∂v ∂x	δx +	∂e(ρ, η) ∂η	∇η • δx dD
				= -	t1 t0	D(t)	ρδx • ▽ e +	p ρ	dD

  -σ 11 u -σ 12 v -σ 13 w) = 0, where σ ij = α s σ sij + α g σ gij ,

											Dα g Dt	.	(45k)
	System (45) admits the energy conservation law :
	∂ ∂t	ρ e +	u 2 + v 2 + w 2 2	+	∂ ∂x	ρu e +	u 2 + v 2 + w 2 2	(46)
	ρe = ρ Y g e g	Y g ρ α g	,	ρη g ρ	+ Y s e h s	Y s ρ 1 -α g	,	ρη s ρ	+ Y s e e s	G |G| 1/3

∂ (ρη g ) ∂t + div (ρη g v) = (p g -p I ) Y g θ g