domain: phys spi math nlin sde

text:
Thomas Fonty
email: thomas.fonty@edf.fr

Martin Ferrand

Agnès Leroy

Antoine Joly

Damien Violeau

Mixture model for two-phase flows with high density ratios: A conservative and realizable SPH formulation

Keywords: Smoothed Particle Hydrodynamics, two-phase mixture flows, large density ratios

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Introduction

Multi-component flows with high density ratio play a prominent role in many engineering applications and imply complex strong flow dynamics (e.g. air-water turbulent mixing in hydraulic jump or over dam spillways as in [START_REF] Wan | Numerical simulation for the air entrainment of aerated flow with an improved multiphase SPH model[END_REF]). The SPH (Smoothed Particle Hydrodynamics) method therefore appears as a natural way to deal with such cases, due to its ability to model such highly deformed flows (e.g. dam break or wedge entry in water detailed in [START_REF] Gong | Two-phase SPH Simulation of Fluid-Structure Interactions[END_REF]).

Multi-fluid SPH models have been extensively studied to take advantage of the good ability of the method in tracking the interface between different phases. Particular attention was paid to high density ratio cases applied to air-water cases such that dam-break or rising-bubble cases in [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF], or Rayleigh-Taylor instability and gravity currents for moderate density ratios in [START_REF] Grenier | Modélisation numérique par la méthode SPH de la séparation eau-huile dans les séparateurs gravitaires[END_REF], [START_REF] Monaghan | A simple SPH algorithm for multi-fluid flow with high density ratios[END_REF]. However, accurately modelling multi-fluid phenomena with the usual SPH approach in the air-water case requires choosing a particle discretisation of less than the size of an air bubble or water drop, which leads to prohibitive computational cost at the scale of practical interest for engineering applications.

Local instant formulation for multi-phase flows being generally beyond computational capability, multi-phase averaged models have been developed, as in [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow -Second Edition[END_REF]. Due to the averaging process, the interface between phases is no longer explicitly tracked. One gets a continuum in which phases are followed through a volume fraction representation. This kind of model encompasses two main approaches1 :

• Two-fluid models. Each phase is described by a continuity and a momentum equations.

Interfacial interactions introduce additional terms in these equations. They are closed by constitutive laws, depending on the flow regime, but their determination can prove to be tedious;

• Mixture models. The flow is seen as a single-fluid flow with one continuity and one momentum equation that rule the evolution of mixture quantities (mixture density and velocity, combinations of phase properties to be defined) complemented by an additional equation for the mass conservation of one phase. Some additional terms linked to the relative velocity between phases, computed through a closure law depending on the flow regime, appear in 1. INTRODUCTION these equations. In the absence of relative velocity, one has a homogeneous fluid model. These models are of particular interest for modelling flows with small-scale interfaces (e.g. dispersed air phase in a water flow) and are usually implemented using the finite volume approach, as in [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF]. A prominent point consists in solving the volume fraction evolution properly as underlined by [START_REF] Rusche | Computational fluid dynamics of dispersed two-phase flows at high phase fractions[END_REF]. However, free surface tracking is not an easy task with meshbased approaches, so that these kinds of multi-phase models were also developed in SPH to deal with sediment-laden flows with a two-fluid model in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] or separated flows with mixture models in [START_REF] Grenier | Modélisation numérique par la méthode SPH de la séparation eau-huile dans les séparateurs gravitaires[END_REF]. It introduces an original way of dealing with phases in SPH: instead of having different sets of particles for each phase, each individual SPH particle is now able to carry the different phases with their respective volume fractions. First studies employing volume fraction representation in SPH were done in an astrophysical framework to study gas/dust mixtures in [START_REF] Monaghan | SPH simulation of multi-phase flow[END_REF] and [START_REF] Monaghan | Implicit SPH drag and dusty gas dynamics[END_REF], following a two-fluid model, further applied to sedimentation in [START_REF] Kwon | Sedimentation in homogeneous and inhomogeneous fluids using SPH[END_REF]. The unique velocity field for the motion of SPH particles makes the mixture model enticing. It is close to the usual single-phase set of equations and can therefore be implemented with limited effort. It is the approach taken here. The reduced number of equations to solve is expected to limit the computations compared to two-fluid approaches with a different set of particles for each phase, as described in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] or [START_REF] Kwon | Sedimentation in homogeneous and inhomogeneous fluids using SPH[END_REF]. Mixture models have already been tested in SPH:

• [START_REF] Price | Two-phase mixtures in SPH -A new approach[END_REF] proposed a mixture model applied to dust settling in protoplanetary discs by rewriting the two-phase equations through the introduction of new variables, combinations of phase variables, to overcome some issues inherent to the two-fluid approach (overdamping at high drag, artificial trapping of dust particles);

• In the graphics community, a mixture model with a volume fraction representation was implemented to capture a wide range of physical phenomena in [START_REF] Ren | Multiple-fluid SPH Simulation Using a Mixture Model[END_REF] but without quantitative validation. This work includes a closure for diffusion velocities (velocity of the phases relative to the mass centre, equal to the difference of the phase velocities and massweighted mixture velocity);

• A mixture model with a volume fraction formulation and a diffusion of phases between particles following a Fick law has been implemented and successfully tested on gravity currents 2. MATHEMATICAL MODEL in [START_REF] Cueille | Modélisation par Smoothed Particle Hydrodynamic des phénomènes de diffusion présents dans un écoulement[END_REF], extending to the Lagrangian framework the Eulerian model proposed by [START_REF] Chanteperdrix | Modélisation et simulation numérique d'écoulements diphasiques à interface libre. Application à l'étude des mouvements de liquides dans les réservoirs de véhicules spatiaux[END_REF];

• Following [START_REF] Cueille | Modélisation par Smoothed Particle Hydrodynamic des phénomènes de diffusion présents dans un écoulement[END_REF], a mixture model with a volume fraction formulation without phase exchanges between particles was implemented and compared to a multi-fluid formulation in [START_REF] Grenier | Modélisation numérique par la méthode SPH de la séparation eau-huile dans les séparateurs gravitaires[END_REF], especially on a Rayleigh-Taylor instability case. It highlighted the diffusion of the interfaces triggered by the mixture model. The multi-fluid formulation of SPH proved to perform better on the cases tested.

In this paper, the physical and numerical formulations chosen for a new mixture model are presented. It is designed to be able to deal with high density ratio flows and versatile physics. Special attention is paid to the conservativity, realizability (respect of physical bounds) and accuracy of this resolution, issues that have raised some noticeable interest as in [START_REF] Rusche | Computational fluid dynamics of dispersed two-phase flows at high phase fractions[END_REF]. Then, numerical validation on several cases is detailed, with comparisons to reference solutions in the literature.

The cases comprise an oil-water and an air-water separation cases, an analytical case similar to a two-phase mixture Poiseuille flow, a Rayleigh-Taylor instability case and a sand dumping case.

Mathematical model

Let us consider a two-phase flow (e.g. air bubbles rising or sediments falling within water). This section aims at detailing the notations and the mixture equations adopted in this work to model such a flow. We consider the control volume2 presented in Figure 1. The presence of phases is taken into account through volume fractions, as illustrated schematically in Figure 2. The two phases are denoted by α and β.

Notations

The mixture quantities detailed in Table 2 are computed based on the value of some phase quantities listed in Table 1. Several choices are possible for the phase quantities to be considered: in view of achieving the conservation of phase quantities at high density ratios, a volume formulation has been retained here for the phase description (volume fraction3 ) and mixture velocity. The mixture velocity is indeed defined here with respect to the volume centre, and is also called volumetric flux, due to weighting by volume fractions. A weighting by mass fractions would correspond to a definition with respect to the mass center. A mass formulation using mass fractions would lead to a simpler set of mixture equations but trigger some numerical issues for high density ratio flows. To preserve the symmetry of the system, we choose to work with relative velocities instead of drift/diffusion velocities.

Mixture quantities

As a consequence of the above-mentioned definitions, one has the following useful relations:

α + β = 1
(1)

v α = j + βv r and v β = j -αv r (2)
The Lagrangian derivative for a particle moving at the mixture velocity j is denoted by:

d dt = ∂ ∂t + j • ∇ (3)

MATHEMATICAL MODEL

Volume V k Mass m k Density ρ k = m k /V k Kinematic viscosity ν k Dynamic viscosity µ k Velocity v k Pressure p k

Mixture model

Following [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow -Second Edition[END_REF], the averaged model for two-phase flows gives the following set of continuous phase equations:

∂αρ α ∂t + ∇ • (αρ α v α ) = 0 (4) ∂αρ α v α ∂t + ∇ • (αρ α v α ⊗ v α ) = -∇ (αp α ) + ∇ • (αT α ) + αρ α g (5)
where T α = µ α (∇v α + t ∇v α ) stands for the viscous stress tensor and g the gravity vector. We neglected here the interfacial terms. Combining these relations and using the mixture model variables introduced in 2.1, one gets:

∂ρ ∂t + j • ∇ρ = -ρ∇ • j + ∇ • αβ ρ β -ρ α v r (6) ∂αρ α ∂t + j • ∇ (αρ α ) = -αρ α ∇ • j -∇ • (αρ α βv r ) (7)
Under the hypothesis of constant phase densities (i.e. ρ α = cst and ρ β = cst), one gets the volume fraction equation:

∂α ∂t + j • ∇α = -∇ • (αβv r ) (8)
Summing this relation on both phases leads to a divergence-free mixture velocity field. This conservation equation, if we except the relative velocity term, is somewhat similar to the Volume of Fluid approach proposed by [START_REF] Hirt | Volume of Fluid (VOF) method for the dynamics of free boundaries[END_REF] or the topological equation introduced in 7-equations models for two-phase flows as detailed in [START_REF] Baer | A Two-phase Mixture Theory for the Deflagration-to-Detonation Transition (DDT) in Reactive Granular Materials[END_REF], neglecting compressibility and

MATHEMATICAL MODEL

Volume V = V α + V β Mass m = m α + m β Density ρ = m/V = αρ α + βρ β Volume fractions α = V α /V and β = V β /V Mass fractions Y α = αρ α /
= Y α v α + Y β v β Relative velocity v r = v α -v β Drift velocity v α -j Diffusion velocity v α -u
Mixture pressure p thermodynamical effects. A volume fraction equation has been proposed in [START_REF] Chanteperdrix | Modélisation et simulation numérique d'écoulements diphasiques à interface libre. Application à l'étude des mouvements de liquides dans les réservoirs de véhicules spatiaux[END_REF] and includes additional terms due to compressibility effects that might be considered in further studies. In this work we choose to use the usual momentum equation using mixture velocities and densities for the sake of simplicity, as we focus on volume fraction and density resolution in a first attempt at deriving an SPH mixture model:

dj dt = - 1 ρ ∇p + 1 ρ ∇ • T + g (9)
where T = ρν (∇j + t ∇j) with ν = αν α + βν β the mixture kinematic molecular viscosity. [START_REF] Chanteperdrix | Modélisation et simulation numérique d'écoulements diphasiques à interface libre. Application à l'étude des mouvements de liquides dans les réservoirs de véhicules spatiaux[END_REF] indicates that there is no preferential argument to choose either a mixture kinematic or dynamic viscosity. The choice done here is motivated by a momentum equation derivation in this volume-based framework that makes a mixture kinematic viscosity appear naturally. An equation consistent with the volume-based formulation would require additional careful work (as explained in [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow -Second Edition[END_REF], including an additional viscous terms and a contribution of the form

∇ • (αβv r ⊗ v r )).

State equation

We assume the mechanical equilibrium of partial pressures as pressure relaxation time is generally small compared to other characteristic times of the flow, as illustrated in [START_REF] Labois | Modélisation des déséquilibres mécaniques dans les écoulements diphasiques : approches par relaxation et par modèle réduit[END_REF] in case 3. NUMERICAL MODEL of bubbly flows:

p = p α = p β (10)
In the standard WCSPH (Weakly Compressible SPH) approach, the mixture is considered as barotropic and the pressure is calculated based on the density values, using the Tait equation of state as reported in [START_REF] Cole | Underwater explosions[END_REF]. We adapted the state equation to a two-phase mixture with a unit polytropic index (common assumption, as in [START_REF] Marrone | δ-SPH model for simulating violent impact flows[END_REF] for example) and a sound speed computed following [START_REF] Grenier | Modélisation numérique par la méthode SPH de la séparation eau-huile dans les séparateurs gravitaires[END_REF]. We linked the pressure to a ratio of volumes instead of densities:

p = αρ α (c α ) 2 + βρ β c β 2 V 0 V -1 + p B ( 11
)
where V 0 is the reference volume equal to δr 2 in 2D, δr 3 in 3D, δr being the particle discretization.

This formula reduces to a one-phase formulation for null or unit volume fraction: in single-phase SPH, when particle masses are constant, the ratio of densities ρ/ρ 0 can be identified to the ratio of volumes V 0 /V . In the following, if not specified, we assume that the background pressure p B is null and sound speeds are chosen so as to ensure that c k > 10 max v max , √ gH where v max is the maximum mixture velocity and H the maximum height of the flow under gravity. As in classical WCSPH, this ensures the compressibility effect is purely numerical and the sound speeds c α and c β aim at ensuring relative density fluctuations within 1%, as explained by [START_REF] Monaghan | Simulating free surface flows with SPH[END_REF].

Numerical model

The fluid domain Ω is discretized with a set of SPH fluid particles F. The discrete value A b = A (r b , t) denotes the quantity A at the particle position vector r b and time t.

SPH operators

The continuous SPH interpolation is performed with the Wendland C 2 kernel introduced in with a volume computed with respect to a reference volume through a fraction θ. The reference volume for vertices and segments writes:

V a = b∈F V b w ab b∈F w ab (12)
It is equal to the actual volume for fluid particles. On the other hand, θ is defined as the angle between two connected segments divided by 2π in 2D for vertices, 1/2 for segments and 1 for fluid particles. Segments and vertices contribute to accurately compute fields and their derivatives close to the boundaries through boundary and volume integrals, respectively. Vertices are indeed part of the fluid and allow one to impose field quantities. If not stated otherwise, multi-fluid operators adapted to this framework are employed following [START_REF] Ghaitanellis | Application of the unified semi-analytical wall boundary conditions to multi-phase SPH[END_REF]. The SPH pressure gradient approximation is therefore:

(∇p) a ≈ G a {p b } = 1 γ a V a b∈F ∪V θ b V 2 a p a + V 2 b p b ∇w ab - 1 γ a V a s∈S 1 V s V 2 a p a + V 2 s p s ∇γ as (13)
In this operator, the renormalization factor γ a presented in [START_REF] Ferrand | Unsteady open boundaries for SPH using semi-analytical conditions and Riemann solver in 2D[END_REF] was used to account for potential truncated kernel supports near boundaries:

γ a = Ωa∩Ω w (r a -r) dr (14)
where Ω a denotes the kernel support of the particle a. Moreover, following [START_REF] Ferrand | Unsteady open boundaries for SPH using semi-analytical conditions and Riemann solver in 2D[END_REF], the contributions to the gradient of γ are:

∇γ as = s w(r)dS n s (15)
where n s is the inward unit normal to the boundary segment s. These values can be computed analytically as explained in [START_REF] Ferrand | Unsteady open boundaries for SPH using semi-analytical conditions and Riemann solver in 2D[END_REF].

Except in subsection 4.2, the SPH Laplacian approximation used to compute the viscous stress tensor writes:

(∇ • (µ∇j)) a ≈ L a {µ b , j b } = 1 γ a b∈F ∪V θ b V b 4µ a µ b µ a + µ b j ab r 2 ab r ab • ∇w ab - 1 γ a s∈S 2µ a j as r 2 as r as • ∇γ as (16)
where r ab = |r ab |, j ab = j a -j b and µ = ρν being the dynamic molecular viscosity. Following [START_REF] Violeau | Dissipative forces for Lagrangian models in computational fluid dynamics and application to smoothed-particle hydrodynamics[END_REF] and [START_REF] Español | Smoothed dissipative particle dynamics[END_REF], the alternative formulation used in case of variable viscosities is:

L a {µ b , j b } = 1 γ a b∈F ∪V θ b V b 2µ a µ b µ a + µ b (4 (j ab • e ab ) e ab + j ab ) r ab • ∇w ab r 2 ab - 1 γ a s∈S 2µ a j as r 2 as r as • ∇γ as (17)
where e ab = r ab r ab . This formulation is used in subsection 4.2.

Integral balance for a single-phase flow

The usual SPH discretisation of volume fraction equation ( 7) leads to non conservative and non realizable solutions (i.e. volume fractions may take non-physical values). We intend to derive a realizable phase volume equation, conservative with respect to the relative velocity contribution, through a finite-volume like approach. In this prospect, we will consider a Voronoi tesselation that discretizes the fluid, as illustrated in Figure 4. The cells in this Figure are a representation of the SPH particles, used to derive the subsequent equations. During a simulation, the exact shape of the particles is never actually known (or used). a|b stands for the interface index between the cells a and b, the outwards-oriented surface vector associated to this interface being denoted S a|b . For the sake of clarity, we first apply this approach to a single phase medium to see how one can recover the WCSPH volume conservation equation:

dV dt = V ∇ • j (18)
After this, a similar reasoning shall be applied to the two-phase case. To start with, we only consider particles that do not interact with the boundaries of the domain.

NUMERICAL MODEL

Leibniz' rule

Let us apply the Leibniz' rule for a scalar field A on a particle a seen as a material volume Ω a moved and deformed by the fluid velocity j:

d dt Ωa(t) AdV = Ωa(t) ∂A ∂t dV + ∂Ωa(t) Aj • dS (19)
For A = 1, equation ( 19) gives:

dV a dt = ∂Ωa(t) j • dS (20)
We therefore recognize an integral form of equation ( 18). For a closed surface:

∂Ωa(t) dS = 0 (21)
Subtracting j a • (21), one gets:

dV a dt = ∂Ωa(t) (j -j a ) • dS (22)

Discrete SPH approximation

Let us proceed with a finite volume formulation of the integrals in equation ( 22). For each neighbouring cell, one has: Equation ( 23) is exact while equation ( 24) is a definition of the interface mixture velocity j a|b . By summing over all neighbours (i.e. all the a|b interfaces), one gets:

a|b dS = S a|b (23) a|b j • dS ≡ j a|b • S a|b (
dV a dt = a|b j a|b -j a • S a|b (25)
In the above equation, j is the velocity at the continuous level, while j a is the velocity of the particle a, which is a discrete field, constant over the support of the particle (see Figure 5).

The similarity between SPH-ALE and Finite Volume (FV) approaches have been underlined in [START_REF] Neuhauser | Development of a coupled SPH-ALE/Finite Volume method for the simulation of transient flows in hydraulic machines[END_REF], among others. In the SPH formalism, this analogy is based on the following relation:

S a|b = 2V a V b ∇w ab (26)
to determine the virtual surface vector of the interface between particles. An important difference, as highlighted by [START_REF] Neuhauser | Development of a coupled SPH-ALE/Finite Volume method for the simulation of transient flows in hydraulic machines[END_REF], between FV and SPH, is the definition of a neighbour in each method: only cells that share (part of) a face with the considered cell in FV are considered as neighbours, whereas all SPH particles within the kernel support are considered as neihgbours (their number thus depends on the smoothing length). This difference of neighbour definition requires an adequate handling of the weighting of the neighbouring SPH particle contribution, which is done through the presence of the kernel gradient in equation ( 26). A limit of this approach is that the property b S a|b = 0, that is true in FV, does not hold in SPH (as zeroth consistency is not fulfilled). This notation will be used throughout this document.

Choosing a centred value to approximate j a|b , equation ( 25) becomes:

dV a dt = -V a b∈F V b (j a -j b ) • ∇w ab (27)

NUMERICAL MODEL

We recognize a discrete SPH approximation of ( 18), the right hand side being an SPH discrete divergence operator as described in [START_REF] Violeau | Fluid mechanics and the SPH method: theory and applications[END_REF].

Integral balance for a two-phase flow

One can now make a similar reasoning in the case of a two-phase flow. Due to the mixture model we have adopted, we can consider that we have a single fluid, but whose characteristics will vary depending on the local volume fraction α.

Leibniz rule

For A = α, equation ( 19) becomes:

d dt Ωa(t) αdV = Ωa(t) ∂α ∂t dV + ∂Ωa(t) αj • dS (28)
Using ( 8) to compute the right-hand side volume integral and using Gauss' theorem, we get:

dV α a dt = Ωa(t) [α∇ • j -∇ • (αβv r )] dV (29)
Expanding α at the zeroth order, we make the following approximation:

Ωa(t) α∇ • jdV = α a Ωa(t) ∇ • jdV (30)
Subtracting α a j a • (21) leads to:

dV α a dt = α a ∂Ωa(t) (j -j a ) • dS - ∂Ωa(t)
αβv r • dS (31)

Discrete SPH approximation

The finite volume approximation of the integrals leads to:

dV α a dt = α a a|b j a|b -j a • S a|b - a|b (αβv r ) a|b • S a|b (32)
where we identify two contributions: the first term corresponds to a divergence of the mixture velocity field α∇ • j and the second to a contribution from the relative displacements of phases ∇ • (αβv r ). The interface a|b quantities still need to be defined in equation (32). Four principles guide this definition:

• Conservativity: the fluxes must be symmetrical with respect to the particle labels (a, b);

NUMERICAL MODEL

• Realizability: phase volumes V α and V β should remain positive;

• Reduction to single phase model: for α = 1 and v r = 0, we want to recover single-fluid WCSPH equations;

• Symmetry with respect to phases: we can interchange α ↔ β without modifying the equations.

Guided by these principles, we propose the following relations on both phase volumes:

dV α a dt = α a dV a dt -2V a b∈F V b α a β b v r a|b • ∇w ab + + α b β a v r a|b • ∇w ab - ( 33
)
dV β a dt = β a dV a dt + 2V a b∈F V b α a β b v r a|b • ∇w ab + + α b β a v r a|b • ∇w ab - ( 34
)
where [x] + = max (x, 0) and [x] -= min (x, 0). The relative velocity is computed with a centred approach:

v r a|b = 1 2 (v r a + v r b ) (35)
Equation ( 34) was obtained from (33) by changing α into β, since these relations are symmetric (this switch changes the sign of the relative velocity). There are several important choices for the derivation of equation ( 33) from equation ( 32):

• The first term of the right hand side of equation ( 32) can be identified as the first term of equation ( 33) using the relation (25). The computation of the temporal derivative of the total volume relies on the use of an exact time integration scheme, following what is suggested by [START_REF] Ferrand | Unsteady open boundaries for SPH using semi-analytical conditions and Riemann solver in 2D[END_REF]. It proves to be better at simulating the mixture at rest (the usual continuity equation resolution with antisymmetric divergence led to spurious convection cells which appeared to be linked to the accumulation of numerical errors). Starting from the SPH volume interpolation:

γ a V a = b∈(F ∪V) θ b w ab (36)
one can update the total volume through:

γ n+1 a V n+1 a - γ n a V n a = b∈(F ∪V) θ n+1 b w n+1 ab -θ n b w n ab (37)
• The second term of the right hand side of equation ( 32) is computed through an upwind formulation, symmetric with respect to phases, that ensures the realizability of the scheme

NUMERICAL MODEL

(under a CFL condition, see section 3.7). The flux v r a|b • S a|b at the midpoint of the pair of particles (a, b) is split into its positive and negative parts. The factor αβ then takes a different form depending on the sign of the flux: in α a β b , α is upwinded with respect to v r and β is convected by (-v r ). The choice of alternate indices a and b for the volume fractions ensures the symmetry with respect to the phases and the conservation of the total quantity of each phase.

A similar formulation can be found in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF]) that details a two-phase model for sediment laden flows with an asymmetric treatment of the liquid and solid phases. The volume fraction of sediment is updated with an equation looking like (33) but with β taken to 1 and a usual SPH divergence of the fluid velocity field to compute the total volume variation. Instead of a β equation, an equation on βρ β is solved. Our approach has the advantage to keep a symmetrical treatment of both phases and to be conservative with respect to the relative velocity contribution. The model can therefore handle only water or air phase within a particle when the approach detailed in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] assumes constant water mass in particles, which may trigger some significant volume variations. We are only solving three equations instead of four, which decreases the computational cost of the approach. A mixed upwind/downwind operator in a finite volume framework is proposed in [START_REF] Larreteguy | A bounded upwind-downwind semidiscrete scheme for finite volume methods for phase separation problems[END_REF] to solve a volume fraction equation with drift velocities and a massweighted mixture velocity (α is upwinded whereas the drift velocity encompassing β and v r is downwinded). Such a distinction between upwind and downwind, due to a stencil limited to a pair of particles in SPH, does not seem possible to implement.

As in usual WCSPH, the total volume is not conserved due to weakly-compressibility. However, the antisymmetry of the discrete terms concerning the relative velocity ensures that the relative phase motion does not affect the total volume. Introducing the wall renormalization factor γ of equation ( 14), the phase volume equation (33) becomes:

dV α a dt = α a dV a dt -2 V a γ a b∈F V b α a β b v r a|b • ∇w ab + + α b β a v r a|b • ∇w ab - ( 38
)
the temporal variation of V a being computed from (37), the symmetry of the relative velocity term is somewhat broken by this factor, as always with the USAW technique but it did not prove to introduce significant discrepancies. No boundary term is considered so as to ensure a no-flux condition at the boundaries. Let us finally underline that in the numerical implementation phase volumes have been adimensionalized by the reference volume V 0 to avoid the accumulation of

NUMERICAL MODEL

numerical errors that triggered some spurious pressure profiles.

Volume diffusion

For constant relative velocities, checkerboard effects clearly appear in simulations as a consequence of the colocated nature of SPH. One possible way to circumvent this issue is to introduce volume diffusion when phase volumes are updated. The derivation of this additional term is detailed in Appendix A.

Particle characteristics

Once phase volumes are computed, one can assess the other particle quantities according to their definitions (here particle labels a were dropped for the sake of simplicity):

• m = (αρ α + βρ β )V 0 • ρ = m/V • α = V α /V
Special attention shall be paid to the particle mass computation. In single-fluid SPH, particle masses are usually constant (except in the Arbitrary Lagrangian-Eulerian approaches such as in [START_REF] Vila | On particle weighted methods and Smoothed Particle Hydrodynamics[END_REF]), so that they do not vary due to compressibility effects; only the density does. To recover this behaviour, we introduce the reference volume V 0 to compute particle masses, so that they only vary due the volume fraction variations. With this choice, we indeed have a varying density due to the fluid (weak) compressibility. This way of calculating the particle quantities constitutes a discrete approximation of the original continuous equations ( 6) and (8). To check that, one can deduce from (33) the discrete governing equations of the mixture density and volume fraction:

dρ dt = - ρ a V a dV a dt -ρ α -ρ β V 0 V 2 a b α a β b v r a|b • S a|b + + α b β a v r a|b • S a|b - ( 39
)
dα dt = - 1 V a b α a β b v r a|b • S a|b + + α b β a v r a|b • S a|b - ( 40
)
the sum on the right hand side being an SPH approximation of ∇ (αβv r ), as previously explained.

NUMERICAL MODEL

Closure law for the relative velocity

The closure law for the relative velocity is flow-dependent and one can refer to the extensive literature, for example [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow -Second Edition[END_REF], to find adequate closures depending on the twophase flow regime. The relations used in this article assume constant particle sizes. More complex models solving transport equations for particle diameters exist but we will focus on the simple closures described hereafter. Let us note v 0 the rising velocity of a single, small spherical bubble of diameter d p within an infinite medium at rest under gravity, according to Stokes' law:

v 0 = ρ α -ρ β d 2 p 18µ β g (41)
Under the assumption that a local equilibrium can be reached over a short spatial length scale, one can write an algebraic expression of the relative velocity according to [START_REF] Manninen | On the mixture model for multiphase flow[END_REF]. Several cases are presented in this work:

• A constant relative velocity:

v r = v 0 (42)
• In [START_REF] Brethour | Drift model for two-component flows[END_REF], the relative velocity is variable and linked to the pressure gradient through:

v r = V p K p α ρ α -ρ β ρ ∇p ( 43
)
where the drag coefficient K p is computed as:

K p = 1 2 A p ρ β C d |v r | + 24ν β d p (44)
with A p = 1 4 πd 2 p and V p = 1 6 πd 3 p the cross sections and volumes of the bubbles. To account for the particle-particle interactions, the Richardson-Zaki correlation can be employed to compute the effective relative velocity:

v r eff = α ζ v r ( 45
)
where ζ is computed using the Reynolds number Re = d p |v r |/ν β as:

             Re < 1 ζ = 4.35/Re 0.03 1 < Re < 500 ζ = 4.45/Re 0.1 500 < Re ζ = 2.39 (46)

NUMERICAL MODEL

To close (43), we use here ∇p = ρg in the separation case 4.1.1, assuming hydrostatic equilibrium.

• Finally, we propose a relative velocity including a volume fraction gradient:

v r = v 0 -K ∇α α ( 47
)
where the coefficient K > 0 accounts for turbulent-like diffusion effects. The presence of a gradient requires some numerical precautions: to avoid double summation that can lead to instabilities, a second order SPH operator is used. Substituting ( 47) in ( 32), one gets:

dV α a dt = α a dV a dt - b∈F (αβv 0 ) a|b • S a|b + K b∈F (β∇α) a|b • S a|b (48)
The relative velocity contribution in the phase volume equation is divided in two parts: the v 0 part is treated as in (33) whereas the diffusion term in K is discretized separately as a usual SPH Laplacian introduced by [START_REF] Morris | Modeling low Reynolds number incompressible flows using SPH[END_REF]:

dV α a dt = α a dV a dt - b∈F α a β b v 0,a|b • S a|b -+ α b β a v 0,a|b • S ab + +K b∈F β a + β b 2 α a -α b r ab r ab • S a|b r ab (49)

CFL condition

The resolution of the phase volume equation ( 33) introduces a condition on the time step for positivity of the phase volume (i.e. realizability of the numerical scheme). Let us consider the relative velocity of general form (47). It is associated to the numerical Péclet number

P e σ = v max σ K (50)
where v max = max (|v 0 |) and σ is the kernel standard deviation described in [START_REF] Violeau | On the maximum time step in weakly compressible SPH[END_REF]. Appendix B details the derivation of the sufficient condition that guarantees the positivity

of V α,n+1 a given the positivity of V α,n a , V β,n a : CF L = v max δt σ ≤ ξ η + 1 P e σ -1 (51)
The non-positivity of V α can trigger some instabilities so that this sufficient condition of positivity is interpreted as a condition of stability. For an advection-diffusion of the volume fraction α along the vertical axis, we got numerically ξ ≈ 1.1 and η ≈ 0.8 for P e σ ∈ [0.0012, 0.12] (optimized with respect to the L 1 norm of the relative error). The fitting curve is displayed in Figure 6.

NUMERICAL MODEL

Boundary conditions

In what follows, except in subsection 4.2, the volume fraction at boundaries is computed with a Neumann condition using the following approximation:

α v = b∈F V α b w ab (52)
This value of the volume fraction of vertex particles is used in some of the operators described above.

Time marching scheme

The detailed resolution of our system is as follows (particle labels are dropped here):

1. Relative velocity update using 3.6:

v r,n = g (ρ n , α n , β n ) (53)
2. Operator splitting for the momentum equation: velocity update

• Potential force step:

j * = j n + δt - 1 ρ n ∇p n + g (54)
4. NUMERICAL VALIDATION

• Viscous force step:

j n+1 = j * + δt 1 ρ n ∇ • T * with T * = ρ n ν n ∇j * + t ∇j * (55)
3. Particle position update:

r n+1 = r n + δtj n+1 (56)
4. Total volume update using (37):

V n+1 = f (V n , r n , r n+1 ) (57)
5. Phase volume computations with ( 33) and ( 34):

V α,n+1 = V α,n + α n V n+1 -V n -δtV n ∇ • (α n β n v r,n ) (58) V β,n+1 = V β,n + β n V n+1 -V n + δtV n ∇ • (α n β n v r,n ) (59)
6. Computation of other flow features using 3.5:

m n+1 , ρ n+1 , α n+1 (60)
7. Pressure computation using state equation ( 11):

p n+1 = p α n+1 , β n+1 , ρ n+1 (61)
To comply with the saturation condition (1), we deduce β = 1 -α when needed in the implementation.

Numerical validation

All our application cases are two-dimensional with particles initially placed on a Cartesian lattice.

Separation

The separation of liquid dispersions, the sedimentation in a settling tank, etc. are common industrial separation problems. In this subsection, we investigate two applications of this separation phenomenon illustrated in Figure 7. This test case is of particular interest to check the realizability and conservativity of the numerical scheme. 20 as pressure values are quite small (for the relative error computation, the minimum pressure in the domain has been subtracted, as a constant pressure offset was present throughout the domain). In the bulk of the fluid domain, the error is close to 1%, as expected in WCSPH.

4. NUMERICAL VALIDATION x z α = 1 α = 0 z x V β V α z α = 1 α = 0

NUMERICAL VALIDATION

Two-phase mixture Poiseuille flow

The computation domain is a rectangle of height 2e and length e, periodic along the x direction as illustrated in Figure 12. The fluid is submitted to a longitudinal force F mimicking a pressure gradient as in a Poiseuille flow. However, no gravity is introduced as it proved to trigger some small instabilities in the velocities and pressure profiles (the gravity effect is all included in the definition of the relative velocity). Starting from an homogeneous mixture α = α 0 at rest, the flow should converge towards the analytical dimensionless steady state4 detailed in Appendix C, obtained with a relative velocity given by (47) at constant mixture kinematic viscosity ν α = ν β . The volume fraction at vertices v of normal n v is computed thanks to a second order approximation method for Robin boundary condition as detailed in [START_REF] Mayrhofer | Investigation of wall bounded flows using SPH and the unified semi-analytical wall boundary conditions[END_REF]. The no-flux condition at the boundary: where n is the normal to the wall is approximated using:

v r = 0 ←→ α (v 0 • n) - ∂α ∂n = 0 (
α v = b∈F V α b ∆ v Y v (1 -(v 0 • n v ) (r vb • n v ) /K) -Z v (r vb • n v ) 2 w vb (63) ∆ v = X v Y v -Z 2 v (64) X v = b∈F V b (1 -(v 0 • n v ) (r vb • n v ) /K) 2 w vb (65)
Y v = b∈F V b (r vb • n v ) 4 w vb (66)
Z v = b∈F V b (1 -(v 0 • n v ) (r vb • n v ) /K) (r vb • n v ) 2 w vb ( 67
)
This flow is analogous to a Rouse flow described in [START_REF] Rouse | Modern conceptions of the mechanics of turbulence[END_REF] as a function of the dimensionless particle size δr/e. This relative error for a field A compared to the analytical profile A th is computed as:

||A -A th || 2 ||A th || 2 = b∈F (A b -A th (z b )) 2 b∈F A th (z b ) 2 (68)
The convergence slope is approximately of order 1/2 for the volume fraction and 1 for the longitudinal velocity. The upwind formulation of the phase volume equation might explain this weak order of convergence for the volume fraction.

As a second step, one can consider a case with a varying mixture kinematic viscosity. Only the longitudinal velocity profile is modified: its expression is detailed in Appendix C. All parameters remain unchanged compared to previous configuration, except the viscosity ν α = 0.001 m 2 /s. We get, after convergence, the longitudinal velocity profile displayed in Figure 17 (the volume fraction profile is unchanged). Again, analytical and numerical solutions agree very well, even though, due to slight discrepancies on the volume fraction profiles, one can notice that the velocities are a bit smaller than expected in the upper part of the flow (the kinematic viscosity also varies with the volume fraction in this case, and the largest absolute errors on the volume fraction are made in the upper part of the flow, therefore influencing the shear stress, and consequently the quality of the

NUMERICAL VALIDATION

p = αρ α (c α ) 2 + βρ β c β 2 γ ρ αρ α + βρ β γ -1 + p B ( 69
)
where γ = 7 is the polytropic index. Sound speeds are c α = 14 m/s and c β = 10 m/s to comply with the idea used for example in [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF] that one should choose values such that ρ α (c α ) 2 = ρ β c β 2 . As spurious fragmentation occurs at the interface between fluids, a small repulsive force was added in the pressure gradient, following [START_REF] Grenier | An Hamiltonian interface SPH formulation for multi-fluid and free surface flows[END_REF]. This force consists of an additional term in the pressure gradient when the neighbouring particle belongs to another phase. Its expression was extended to the mixture model:

G a {p b } = 1 γ a V a b∈F V 2 a p a + V 2 b p b ∇w ab + 1 γ a V a b∈F (α a β b + α b β a ) V 2 a |p a | + V 2 b |p b | ∇w ab ( 70
)
where the factor (α a β b + α b β a ) is equal to 1 when a and b belong to distinct phases, and 0 otherwise as long as phases remain separated, as in this Rayleigh-Taylor case without relative velocity. The relevance of this term for intermediate volume fractions is questionable. As suggested by [START_REF] Grenier | An Hamiltonian interface SPH formulation for multi-fluid and free surface flows[END_REF], is taken equal to 0.01. This test case aims at checking the numerical behaviour for separated phases: the multi-fluid behaviour should be recovered. Figure 19 shows the temporal evolution of the phases with the finest resolution (300 particles per unit length). Figure 20(a)

displays a convergence study of the shape of the interface at the adimensional time t * = t g/e = 5 in which one can see the convergence of the interface shape. in the curling up region of the mushroom-shaped heads.

Introducing a constant relative velocity (42) chosen arbitrarily with d p = 0.07 m, interfaces first diffuse and rapidly converge towards a separated state as in section 4.1. Simulations without and with relative velocity are displayed in Figure 21. One can see that without relative velocity (i.e.

configuration equivalent to the multi-fluid model), phases are not fully separated as some drops remain in both phases and are under-resolved because of the size of SPH particles (a steady state is not reached even at t * = 500 even though phases are then more clearly separated). An inaccurate pressure gradient computation for those small size structures might explain this behaviour, making our mixture approach useful in such circumstances. In this section, the performances of the presented mixture model are assessed on a case of sand dumping described in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF]. This case implies enriched physics compared to the previous sections, with the presence of a free surface and a physically-based relative velocity closure.

NUMERICAL VALIDATION

α denotes the sediment phase, and β the water phase. The configuration is a square water tank

(1 m×1 m) with a free surface as described in Figure 22. In the middle of the domain, just under the free surface, a sand block is released initially with a volume fraction α 0 equal to the maximum sediment concentration. Many different configurations were considered by [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF]) but we will focus on the tests for which experimental results are available. The sand cloud is initially W 0 = 2 or 4 cm wide and H 0 = 2.5 cm high. The three sizes of sediment d p considered are 0.8 mm, 1.3 mm and 5 mm that correspond roughly to settling velocities of 12.6 cm/s, 19.61 cm/s and 49.52 cm/s, respectively. The particle discretization is δr = 0.005 m.

In order to have a fair comparison with the two-fluid approach developed in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF], several changes were introduced to deal with the complete physics of this case: • Computation of the sediment viscosity following Alihan and Sleath's formula [START_REF] Alihan | Sediment Transport in Oscillatory Flow Over Flat Beds[END_REF]:

ν α = 1.2 ρ β ρ α α sm α 1 3 -1 -2 ν β (71)
where α sm is the maximum sediment volumetric concentration, that was set to 0.607 in our computations (instead of 0.606 in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF]) to avoid singularity of the viscosity at the initialization (indeed the sediment cloud starts at this concentration initially).

• Computation of a turbulent viscosity, added to the physical one, following a Smagorinsky model5 :

ν t = (Cδr) 2 ||S|| 1 - α α sm n (72)
where C = 0.1 and n = 5. Please note that we are using here the mixture rate-of-strain tensor defined as:

S = 1 2 ∇j + t ∇j and ||S|| = √ 2S : S (73)
that differs from the phase rate-of-strain tensors used for distinct phase turbulent viscosities computations in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF]. Moreover, a sensitivity analysis has been conducted in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] to optimize the results with adequate choice of the Smagorinsky coefficients and power law, choices that might not be optimal in the framework of our mixture model.

CONCLUSION

Conclusion

In this paper, a two-phase SPH mixture model suited to high density ratios was introduced.

Thanks to a finite volume-like reasoning, a realizable numerical scheme, conservative with respect 5. CONCLUSION et al., 1990) experimental results for dp = 1.3 mm and W 0 = 4 cm. If we assume V α,n > 0 and V β,n > 0 for every particle (therefore V n > 0), then we have V α,n+1 > 0 (and symmetrically V β,n+1 > 0) under the condition that:

CONCLUSION

Figure 1 :Figure 2 :

 12 Figure1: Control volume (green area) in two-phase flow (here a gas phase in red, liquid phase in blue).

Wendland ( 1995 )

 1995 denoted w. For a pair of particles (a, b), w ab = w (r ab ) where r ab = r a -r b and the kernel gradient writes ∇w ab = ∇ a w (r a -r b ). Unified Semi Analytical Wall (USAW) boundary conditions described in Ferrand et al. (2017) are used. In this framework, boundaries ∂Ω of the domain are meshed by segments s ∈ S of size δr (δr being the initial interparticle spacing) connected at vertices v ∈ V as illustrated in Figure 3. Those vertices are truncated fluid particles

Figure 3 :

 3 Figure 3: SPH kernel and intersecting boundary.

Figure 4 :

 4 Figure 4: Voronoi diagram.

Figure 5 :

 5 Figure5: Mixture velocity field j. ja is the mixture velocity of the particle a, constant within the volume of the particle.

Figure 6 :

 6 Figure 6: Advection-diffusion: limit of stability of the volume fraction profile (stability region below the curve).

Figure 7 :Figure 8 :

 78 Figure 7: Separation case: initial and final states.

Figure 9 :

 9 Figure 9: Separation of an oil-water dispersion: evolution of the upper and lower interface positions. Comparison of experimental data of Jeelani & Hartland (1998) and numerical results (Flow-3D R results are taken from Brethour & Hirt (2009)).

Figure 12 :

 12 Figure 12: Two-phase mixture Poiseuille flow: initial and final states.

  Figures 15 and 16 respectively, plotting the instantaneous L 2 relative error for the steady-state field

Figure 15 :

 15 Figure 15: Two-phase mixture Poiseuille flow: convergence study on the volume fraction profile.

Figure 16 :

 16 Figure 16: Two-phase mixture Poiseuille flow: convergence study on the longitudinal velocity profile.

Figure 17 :

 17 Figure 17: Two-phase mixture Poiseuille flow with varying kinematic viscosity: longitudinal velocity profile at steady-state.

Figure 19 :

 19 Figure 18: Rayleigh-Taylor instability: initial and final states.

  Figure20: Comparison of the interfaces for the Rayleigh-Taylor instability case at t * = 5. On the left, a convergence study is displayed for 75 (red line), 150 (blue line) and 300 (black line) particles per unit length. On the right, a comparison is made between Level-Set (red line) and SPH (blue line) results of[START_REF] Grenier | An Hamiltonian interface SPH formulation for multi-fluid and free surface flows[END_REF] and the present SPH model (black line) for the same resolution.

  Figure 21: Rayleigh-Taylor instability: evolution of the volume fractions from the light (blue) to heavy (red) phases without (upper series) and with (lower series) relative velocity for a discretization of 300 particles per unit length.

Figure 22 :

 22 Figure 22: Sand dumping: sand cloud features.

Figure 23 :

 23 Figure23: Comparison of the cloud frontal velocity for the mixture model,[START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] numerical results and[START_REF] Nakasuji | Dynamic behaviour of sand clouds in water[END_REF] experimental results for dp = 0.8 mm and W 0 = 2 cm.

Figure 24 :

 24 Figure24: Comparison of the cloud frontal velocity for the mixture model,[START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] numerical results and[START_REF] Nakasuji | Dynamic behaviour of sand clouds in water[END_REF] experimental results for dp = 5 mm and W 0 = 2 cm.

Figure 25 :

 25 Figure 25: Comparison of the cloud width for the mixture model, (Shi et al., 2017) numerical results and (Nakasuji

Figure 26 :

 26 Figure 26: Comparison of the cloud width for the mixture model, (Shi et al., 2017) numerical results and (Nakasujiet al., 1990) experimental results for dp = 5 mm and W 0 = 4 cm.

Figure 27 :Figure 28 :

 2728 Figure 27: Comparison of the fluctuations of the free surface for the mixture model and (Shi et al., 2017) numerical results for dp = 0.8 mm and W 0 = 2 cm.

Figure 29 :

 29 Figure 29: Sediment cloud for the finest resolution δr = 0.0025 m, dp = 0.8 mm and W 0 = 2 cm.

  term of the left hand side of equation (B.6) is always positive.We now search for an upper bound of the left hand side term of the unequality (B.7interpolation framework, we have according to[START_REF] Violeau | On the maximum time step in weakly compressible SPH[END_REF]: the C 2 Wendland kernel in 2 dimensions. Noting v max = max (|v 0 |), an upper bound can be written: ξ are unknown. Volumes do not vary a lot so that we may consider that the ratio dilogarithm function 6 . C 1 and C 2 are deduced from the no-slip condition at walls:C 1 = -Li 2,r(1)+ Li 2,r (-1) -P e (ln r (1) + ln r (-1)) 2P e 2 + P e (ln r (1) -ln r (-1)) (C.8)C 2 = [P e ln r (1) + Li 2,r(1)] [ln r (-1) -P e] 2P e + ln r (1) -ln r (-1) + [P e ln r (-1) -Li 2,r (-1)] [ln r (1) + P e] 2P e + ln r (1) -ln r (separation of phases, (C.6) gives a condition on α 1 : 0 ≤ α 1 ≤ exp (-P e) (C.11) And therefore a condition on the initial uniform volume fraction α 0 using (C.10): 0 ≤ α 0 ≤ 1 -exp (-2P e) 2P e (C.12) As a second step, one can consider a case with a varying mixture kinematic viscosity. Noting R = ν α -ν β /ν β , only the longitudinal velocity profile is modified: r) P e 2 [P e (R ln R (z * ) -r ln r (z * )) z * +RLi 2,R (z * ) -rLi 2,r (z * ) -r ln R (z * ) + r ln r (z * ) +C 1 (R -r) P e 2 z * -P e (R ln R (z * ) + r ln r (z * )) + C 2 (C.13) where the additional notations ln R (z * ) = ln (1 + α 1 R exp (P e z * )) and Li 2,R (z * ) = Li 2 (-α 1 R exp (P e z * ))are used. C 1 and C 2 are deduced from the no-slip condition at walls:C 1 = -(r + P eR) ln R (-1) + r (P e + 1) ln r (-1) 2 (R -r) P e 2 -P e (R (ln R (1) -ln R (-1)) + r (ln r (1) -ln r (-1))) + -R (Li 2,R (1) -Li 2,R (-1)) + r (Li 2,r (1) -Li 2,r (-1)) 2 (R -r) P e 2 -P e (R (ln R (1) -ln R (-1)) + r (ln r (1) -ln r (-1))) +(r -P eR) ln R (1) + r (P e -1) ln r (1) 2 (R -r) P e 2 -P e (R (ln R (1) -ln R (-1)) + r (ln r (1) -ln r (-1If |x| ≤ 1, the series expression writes Li 2 (x) = n x n n 2

Table 1 :

 1 Phase quantities (k = α or β)

Table 2 :

 2 Mixture quantities

  ρ and Y β = βρ β /ρ Mixture velocity w.r.t. the volume center j = αv α + βv β

	Mixture velocity w.r.t. the mass center	u

In view of the hydraulic applications, energy considerations are neglected in what follows.

The control volume models an SPH particle that will be introduced in the next section.

We denote α and β the volume fractions as well as the phase names, without any risk of confusion.

Obtained with an incompressibility hypothesis.

Some prefer using the smoothing length h in place of δr. It is a matter of convenience about the definition of C.However, we emphasize that, following[START_REF] Dehnen | Improving convergence in smoothed particle hydrodynamics simulations without pairing instability[END_REF], h is meaningless and should be replaced by the kernel standard deviation.

Acknowledgment

This work was funded by EDF R&D and the French Research Agency (CIFRE grant agreement #2016-0362).

C 2 = -C 1 (R -r) P e 2 -P e (R ln R (1) + r ln r (1))

-P e [R ln R (1) -r ln r (1)] -RLi 2,R (1) + rLi 2,r (1) + r ln R (1) -r ln r (1) (C.15)

• A relative velocity closure is deduced from [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] continuous equations. It includes the drag effect, an hindering factor computed with Richardson and Zaki formula [START_REF] Richardson | Sedimentation and fluidisation: Part I[END_REF] and a turbulent diffusion of volume fraction based on a gradient transport law:

where the drag coefficient is computed according to the Reynolds number Re = ||v r ||d p /ν β following:

Re 1 + 0.15Re Iterations are made to converge to a well-defined value of the relative velocity, due to the implicit relation implied by the drag. The Schmidt number Sc is taken equal to 1.

• Combining the momentum equations of both phases under the constant phase density hypothesis, one can obtain a momentum equation for the mixture, distinct from (9). The pressure gradient term is modified and a convective transfer term appears. The additional viscous terms are neglected, as they trigger some unphysical behaviours at the beginning of the simulation and do not change the results described in this paper. Hence equation ( 9) was changed into:

A different state equation is used in [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] but the tests made did not highlight any significant modification of the following results. One can follow temporally the evolution of the cloud, defined as the particles with a sediment volume fraction α equal to at least 5% of the maximum value in the domain at the considered instant, through its vertical position Z, its width W and its frontal velocity w c . Moreoever, one can follow the free surface fluctuations at the release point z w : this must be handled with care however, as the flucutations after the initial jump are of the order of the particle size. The adimensionalization coefficients used in the graphs are

Numerical results of the mixture model are compared with SPH results of [START_REF] Shi | Development of a two-phase SPH model for sediment laden flows[END_REF] and experimental results of [START_REF] Nakasuji | Dynamic behaviour of sand clouds in water[END_REF] for the frontal velocity of the cloud in Figures two-velocity approach. The pressure gradient computation for high density ratios needs more investigation: indeed, the multi-fluid operator did not appear sufficient to prevent numerical instabilities close to the interface in the separation cases. A momentum equation consistent with the volumebased formulation shall also be further investigated. Second order schemes using slope limiters can be studied to improve the convergence order.

Appendix A. Volume diffusion

We follow the derivation of [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF], but apply it to the volume equation ( 18) instead of the density equation. Using the momentum equation ( 9), the discrete volume equation writes:

where D a and G a stands for the discrete divergence and gradient operators and δt the time step.

Using the Laplacian operator L a , which is slightly different from the divergence of a gradient in SPH, the volume diffusion term F a writes:

It is written discretely as

The easiest approach is now to introduce the usual volume diffusion in both V α and V β equations, but weighted respectively by α and β. The new phase volumes computed by ( 33) and ( 34) being denoted V α, * , V β, * , using a volume diffusion coefficient Λ (generally Λ = 0.1), we finally write at time n + 1:

Derivation of a condition for positiveness of phase volumes

In the following, the particle a considered is far from boundaries (no boundary term) and we neglect the volume diffusion term. We search for a sufficient condition of realizability of the phase volume using equation ( 49):

Let us note

Replacing α, β and S a|b by their definitions:

Appendix C. Analytical solution of the two-phase mixture Poiseuille flow At steady state, the system of equations to solve ( 8) and ( 9) becomes:

Under the longitudinal periodicity condition and projecting the momentum equation along the longitudinal axis, this system becomes:


Martin Virte
email: mvirte@b-phot.org

Marc Sciamanna
email: marc.sciamanna@supelec.fr

Emeric Mercier

Krassimir Panajotov
email: kpanajot@b-phot.org

Experimental and theoretical analysis of limit cycle bistability in a free-running VCSEL

Keywords: Vertical-Cavity Surface-Emitting Laser, Nonlinear Dynamics, Semiconductor Laser

Laser diodes typically behave like damped oscillators: they are generally expected to only show damped relaxation oscillations toward a stable fixed point. In vertical-cavity surface-emitting lasers (VCSELs), the picture appears to be quite different as polarization dynamics can be experimentally observed including bifurcations to selfpulsation and even chaos. Physically, the circular geometry of VCSELs makes the polarization selection very weak and, thus, the additional degree of freedom can enable complex dynamical behavior in the laser diode.

Here we report on a new dynamical behavior in a free-running VCSEL: we observe a bistability between two limit cycles oscillating around two distinct elliptical polarization states whose main axes are symmetrical with respect to the polarization at threshold. Although the existence of two symmetric elliptical polarizations and the associated limit cycles are predicted by the San Miguel, Feng and Moloney (SFM) model, the hysteresis cycle observed experimentally highlights the importance of asymmetry in the dynamics from the elliptically polarized states. We demonstrate that this behavior can be accurately reproduced in theory within the SFM framework when taking into account a small misalignment between the phase and amplitude anisotropies of the laser cavity. Our results bring new light into VCSEL polarization dynamics and provide a very good qualitative agreement with the bifurcation scenario predicted by the SFM model.

INTRODUCTION

Compared to regular edge-emitting lasers, Vertical-Cavity Surface-Emitting Lasers (VCSELs) exhibit key advantages such as a circular beam, a significant threshold reduction and on-chip testing capability. Nevertheless many experimental studies during the past 20 years reported severe polarization instabilities including polarization switching, [START_REF] Choquette | Temperature dependence of gainguided verticalcavity surface emitting laser polarization[END_REF][START_REF] Choquette | Gain-dependent polarization properties of vertical-cavity lasers[END_REF][START_REF] Van Exter | Polarization fluctuations in vertical-cavity semiconductor lasers[END_REF][START_REF] Panajotov | Polarization switching in VCSEL's due to thermal lensing[END_REF] self-pulsating dynamics, [START_REF] Ackemann | Characteristics of polarization switching from the low to the high frequency mode in vertical-cavity surface-emitting lasers[END_REF][START_REF] Sondermann | Experimental and theoretical investigations on elliptically polarized dynamical transition states in the polarization switching of verticalcavity surface-emitting lasers[END_REF][START_REF] Olejniczak | Polarization switching in quantum-dot vertical-cavity surface-emitting lasers[END_REF][START_REF] Olejniczak | Polarization switching and polarization mode hopping in quantum dot vertical-cavity surface-emitting lasers[END_REF][START_REF] Virte | Bifurcation to nonlinear polarization dynamics and chaos in vertical-cavity surface-emitting lasers[END_REF] and polarization chaos; 10 see Ref. [START_REF] Panajotov | Polarization dynamics of vcsels[END_REF] for a recent review. Since those instabilities are major drawbacks for VCSEL applications, they motivated numerous studies and several mechanisms have been suggested to explain the polarization dynamics, including spin-relaxation processes, in the so-called spin-flip model or SFM. [START_REF] San Miguel | Light-polarization dynamics in surface-emitting semiconductor lasers[END_REF][START_REF] Martin-Regalado | Polarization properties of vertical-cavity surface-emitting lasers[END_REF] In the frame of these studies, elliptically polarized (EP) states have been reported experimentally, [START_REF] Sondermann | Experimental and theoretical investigations on elliptically polarized dynamical transition states in the polarization switching of verticalcavity surface-emitting lasers[END_REF] but there have been only few investigations of their underlying dynamics. [START_REF] Prati | Analysis of elliptically polarized states in vertical-cavitysurface-emitting lasers[END_REF] Yet, it was recently demonstrated that these EP states play an essential role in the bifurcation scenario creating polarization chaos in a solitary VCSEL. [START_REF] Olejniczak | Polarization switching in quantum-dot vertical-cavity surface-emitting lasers[END_REF][START_REF] Olejniczak | Polarization switching and polarization mode hopping in quantum dot vertical-cavity surface-emitting lasers[END_REF][START_REF] Virte | Deterministic polarization chaos from a laser diode[END_REF] According to the SFM model as described in Ref., [START_REF] Martin-Regalado | Polarization properties of vertical-cavity surface-emitting lasers[END_REF] the pitchfork bifurcation destabilizing the linear polarization stable at threshold creates two symmetrical EP states. In theory, the bifurcation scenario will then be strictly identical for both orientations: the ellipticity evolution will be the same -although with different orientations of the main axis, the two Hopf bifurcations destabilizing the EP states will appear at the same current and will create two identical limit cycles, i.e. self-pulsating solutions, oscillating at the same frequency. Because of this equivalence, the laser should settle at random on one or the other EP state when the current is increased. Experimentally, however, we observe that the laser always selects the same orientation. Moreover, reproducible switchings between the two orientations have been reported, [START_REF] Olejniczak | Polarization switching and polarization mode hopping in quantum dot vertical-cavity surface-emitting lasers[END_REF] which confirm that the strict equivalence as seen in theory is not a viable option.

Here we unveil an even more surprising phenomenon: a bistability between two limit cycles oscillating around two EPs with main axes symmetrical with respect to the linear polarization at threshold. The two cycles exhibit different frequencies and amplitudes, and we unveil a clear hysteresis cycle when changing the injection current of the laser. Theoretically, we show that this behavior can be accurately reproduced by breaking the symmetry between the two EP states. In particular, when including a small misalignement between the phase and amplitude anisotropies considered in the SFM framework, we obtain a convincing agreement between experimental results and theoretical predictions. This contribution is organized as follows: the experiment and experimental results are described in section 2. In section 3, we present simulation results based on an extended SFM model, including a misalignment between the cavity anisotropies, which shows good qualitative agreement with the experimental findings. Finally, we discuss the results obtained and the perspectives for future work in section 4.

EXPERIMENTAL RESULTS: OBSERVATION OF A BISTABILITY OF LIMIT CYCLES

The VCSEL analyzed in this experiment is a single-mode InGaAs submonolayer quantum dot VCSEL emitting at 990 nm stabilized in temperature at 22 • C. The device has been described in details in Ref. [START_REF] Hopfer | Single-mode submonolayer quantum-dot vertical-cavity surface-emitting lasers with high modulation bandwidth[END_REF] The laser is collimated with an aspheric lens and coupled into a multimode fiber for measurement. Reflections from the fiber front facet are prevented by using an optical isolator and polarization measurement are performed by rotating a half-waveplate. In this section, we define the polarization at 0 • as the linear polarization at threshold. The output power is measured using a Newport Power Meter, and the time-series are recorded using a fastphotodetector (NewFocus 1554-B, 12GHz), an electronic amplifier (NewFocus 1422LF, 20GHz) and a large bandwidth oscilloscope (Agilent DS0X92504A, 25GHz, 80 GS/s). The Radio-frequency spectra are obtained through a Fast Fourier Transform of the recorded time-series. At threshold the laser starts emitting on the vertical polarization P0, see Fig. 1-(b), but an increase of the current destabilizes this steady-state and creates two symmetrical EP states, shown as two ellipses. To discriminate these two solutions, we need to observe the polarizations at ±45 • where the contrast between the two EPs is maximized. This is what we do in fig. 1-(a), where we plot the polarization resolved output power evolution for increasing (solid line) and decreasing injection current (dashed line). Thus we observe a clear switching between the two solutions -for current around 2.2 mA, similar to what as been reported in Ref. [START_REF] Olejniczak | Polarization switching and polarization mode hopping in quantum dot vertical-cavity surface-emitting lasers[END_REF] When decreasing the current, we unveil an hysteresis cycle between the two EP states for currents between 1.8 mA and 2.2 mA. However these two solutions are not stationary as can be seen from the time-series measurement given in fig. 2. Based on the bifurcation diagrams generated from experimental time-traces, we observe two distinct limit cycles. The first one, the one stable at lower current, shows a larger amplitude and a frequency around 6.25 GHz. When the current is increased and reaches 1.95 mA, we observe a sudden dynamical change, see Fig. 2-(a, b), the amplitude of the oscillations decreases and the frequency goes up to 6.55 GHz. Then, when the current is decreased, the oscillations are strongly reduced and the main peak of the FFT spectrum vanishes, see Fig. 2-(c, d). Under 1.77 mA, the system finally settles back to the first limit cycle hence closing the hysteresis loop. Because the device is very sensitive to external conditions such as stress from the probes, the range of current at which the bistability appears can vary between different measurements. In particular, the hysteresis is observed at two different ranges of current between Fig. 1-(a) and Fig. 2.

SFM MODEL: ANISOTROPY MISALIGNMENT BREAKS EP STATE SYMMETRY

As stated earlier, in the SFM model described in Ref., [START_REF] San Miguel | Light-polarization dynamics in surface-emitting semiconductor lasers[END_REF][START_REF] Martin-Regalado | Polarization properties of vertical-cavity surface-emitting lasers[END_REF] no distinction can be made between the two EP states, and the system is therefore expected to choose randomly between the two. On the other hand, we experimentally report on a biased situation where the two elliptical polarizations are not equivalent, as demonstrated by the hysteresis between the two limit cycles. Hence to model this new behavior we need to identify a mechanism that could break the symmetry between the EP states. We demonstrate here that a misalignment between the phase and amplitude anisotropy has such an effect. Based on this modified SFM model, we then show that a good agreement between theoretical prediction and experimental results can be achieved.

The SFM model considers the right and left circular polarizations as two different emission processes with separated carrier reservoirs; the coupling between them appears through the spin-flip relaxation processes that equilibrate carrier population. The cavity anisotropies then select the main orthogonal axes of polarization and introduce an asymmetry between the two, but the symmetry of the EP states remain. Yet, the introduction of a misalignment between the phase and amplitude anisotropies can break this symmetry. The impact of such misalignment on the steady-states of the SFM model has already been investigated previously, [START_REF] Travagnin | Role of optical anisotropies in the polarization properties of surface-emitting semiconductor lasers[END_REF][START_REF] Travagnin | Linear anisotropies and polarization properties of vertical-cavity surface-emitting semiconductor lasers[END_REF] but no studies were conducted to observe the impact on the dynamical states of the laser. This modified SFM model then writes as follow:

dE + dt = κ.(1 + iα).(N + n -1).E + -(i.γ p + (cos(2.θ) -i.sin(2.θ)).γ a ).E - (1)
dE - dt = κ.(1 + iα).(N -n -1).E --(i.γ p + (cos(2.θ) + i.sin(2.θ)).γ a ).E + (2)
dN dt = µ -N -(N + n).|E + | 2 -(N -n).|E -| 2 (3)
dn dt = -γ s .n -(N + n).|E + | 2 + (N -n).|E -| 2 (4)
With E ± the field for the left (+) and right (-) circular polarization, N the normalized total carrier population and n the normalized carrier population difference between the two reservoirs. The field decay rate is κ, the spin-flip relaxation rate is γ s , the injection current is µ and α is the linewidth enhancement factor. The cavity anisotropies are modeled by the phase anisotropy or birefringence γ p , the amplitude anisotropy γ a and θ, the angle between the phase and amplitude anisotropy axes -i.e. the axis of maximum frequency and the axis of maximum losses. All parameters and variables are dimensionless and the time is normalized by the carrier lifetime. For numerical convenience, we reduce eqs. 1-4 it to a set of 5 real equations considering only the phase difference between the two electrical field for the left and right circular polarization similarly to what was done in Ref. [START_REF] Erneux | Two-variable reduction of the San Miguel-Feng-Moloney model for vertical-cavity surface-emitting lasers[END_REF] To reproduce the behavior reported experimentally we use the following parameter values: κ = 600, α = 3, γ s = 100, γ p = 35, γ a = -10 and θ = 0.023. As seen in fig. 3-(a) and (c), at threshold the laser first emits linearly polarized light, but no pitchfork bifurcation appears to destabilize the state stable at threshold. Instead, we observe a smooth increase of the state ellipticity until a Hopf bifurcation occurs around J ≈ 3.3 hence creating a self-pulsating solution. The newly created limit cycle oscillates around the now unstable EP state: we call this solution LC+. For increasing injection current, a sequence of bifurcation occurs at J = 4.8 that destabilize LC+ and a new limit cycle (LC-) is unveiled. LC-exhibits almost perfectly reversed polarization properties compared to LC+, but, as can be seen from Fig. 3-(e), its frequency is slightly different. We estimate the gap to be grossly around 0.1, i.e. 100 M Hz for a carrier lifetime of 1 ns. From this point, if we decrease the current we unveil a relatively large bistability region between LC-and LC+. We observe that the amplitude of LC-oscillations decreases along with the current; the laser even reaches a steady-state for a tiny range of current. Then at J = 3.5, the laser settles back to LC+, hence closing the hysteresis cycle.

CONCLUSION

To sum up, we report a peculiar limit-cycle bistability in a free-running VCSEL where the two cycles oscillate around two different elliptically polarized states with main axes symmetrical with respect to the linear polarization at threshold. We also demonstrate that this behavior can be qualitatively well reproduced when considering a small misalignment between the phase and amplitude anisotropies in the laser cavity. Our work therefore provides further evidence of the validity of the SFM model as a framework to describe the nonlinear dynamics and polarization properties of VCSELs. We also highlight a new dynamical behavior in a free-running VCSEL and bring new light on the role of the anisotropies in the selection and dynamics of VCSEL polarization. In addition, we motivate additional experimental and theoretical investigations for a better understanding and a better evaluation of the impact of these anisotropies.

Figure 1 .

 1 Figure 1. Experimental observations of the bistability. (a) Polarization LI curve at +45• (thick, red) and -45 • (thin, black) for an increasing (solid) and a decreasing injection current (dashed). (b) Scheme of the different polarization states invoveld: linear polarization at threshold is the vertical gray line P0, and the two EP states around which the limit cycles oscillate are the two ellipses symmetrical with respect to P0. The dashed lines show the polarizations used for the measurement of the LI curve: +45 • (thick red) and -45 • (thin black).

Figure 2 .

 2 Figure 2. Identification of the limit cycles. (a, c) Plot of the time-series extrema versus increasing (top, a) and decreasing current (bottom, c). (b, d) FFT of the time-series versus increasing (top, b) and decreasing current (bottom, d). Time series considered for all panels are recorded at 0 • , i.e. polarization at threshold, with DC removed.

Figure 3 .

 3 Figure 3. Simulated bifurcation diagram for increasing and decreasing injection current. (a, b) Mean output power for polarization at ±45 for increasing (a) and decreasing injection current (b). (c, d) Extrema of the polarization resolved time-series at ±45 • in red (+) and black (-) for increasing (c) and decreasing current (d). (e) Normalized frequency evolution for increasing (black) and decreasing (red) injection current.

ACKNOWLEDGMENTS

The authors acknowledge support from the Conseil Régional de Lorraine, Fondation Supélec, FWO-Vlaanderen, the METHUSALEM program of the Flemish government, and the interuniversity attraction poles program of the Belgian Science Policy Office (grant no. IAP P7-35 "Photonics@be").