Nico Schwan
email: nico.schwan@alcatel-lucent.com

Klaus Satzke
email: klaus.satzke@alcatel-lucent.com

David Griffin
email: dgriffin@ee.ucl.ac.uk

Bertrand Mathieu
email: bertrand2.mathieu@orange-ftgroup.com

Toufik Ahmed

Fostering Cooperation Between Overlay Applications and Underlying Networks

Considering the demanding nature of future high definition (HD) and three dimensional (3D) video communications, application overlays will need increased awareness of underlying network capabilities. We propose a new cross-layer solution to enhance overlay/ISP interactions by providing a comprehensive, media aware, open, and standardized interface between ISPs and application overlays. We further discuss how cooperation between the layers will be improved and describe an example use case to demonstrate the benefit of such interworking. © 2011 Alcatel-Lucent.

(head 1) Introduction

Future networked media environments will be high quality, multi-sensory, multi-viewpoint, and multistreamed and will rely on high definition (HD) and three dimensional (3D) video [1]. For Internet service providers (ISPs), however, this will place unprecedented demands [START_REF]Future Media Internet Research Challenges and the Road Ahead[END_REF] on the networks for high capacity, low latency, low loss and reliable communication paths. Because the traditional solution of pre-provisioning/overprovisioning of sufficient network resources everywhere is not economically viable [1,[START_REF] Velocix | Enabling Digital Media Content Delivery: Emerging Opportunities for Network Service Providers[END_REF] given the huge capacities required, we propose a fundamentally different approach, via a specific interface, to improve cooperation between ISPs and application overlays so that:  Intelligent overlay applications are optimized for end-to-end performance at a global scale, being aware of the actual capabilities of multiple underlying ISPs.

 Specialized network services such as multicast distribution or in-network caching are dynamically invoked in the underlying ISPs by the applications.

 Content, and the way it is generated, accessed, and distributed is adapted on-the-fly to what the networks are able to deliver and according to terminal capabilities.

Due to the lack of a vertical interface that allows cooperative information exchange and service invocation, today's Internet overlay applications, as well as the underlying networks, are using independent and thus often counterproductive content routing strategies. Cooperation will benefit both sides through a mutual optimization of the content delivery. The Internet Engineering Task Force (IETF) has specified the Application-Layer Traffic Optimization (ALTO) protocol [START_REF]Application-Layer Traffic Optimization (ALTO)[END_REF] to begin to bridge this gap. In contrast to ALTO, we propose to establish a much richer interface to exchange additional information and also to invoke services. Projects like NAPA WINE [2] and SmoothIT [START_REF]Seventh Framework Programme[END_REF] are also working to bridge the same gap, but our approach supports more sophisticated applications that are multi-participatory, multi-sourced, and interactive.

(head 1) ENVISION Interface

Current overlay applications are, for the most part, unaware of the underlying network infrastructure since there is no established means for the exchange of cross-layer information. To overcome this limitation, overlay networks try to optimize their topology based on end-to-end measurements [START_REF] Gurbani | A Survey of Research on the Application-Layer Traffic Optimization Problem and the Need for Layer Cooperation[END_REF] or by reverse engineering the network, however this optimization often undermines ISP routing policies [START_REF]Application-Layer Traffic Optimization (ALTO)[END_REF][START_REF] Keralapura | Can ISPs Take the Heat from Overlay Networks?[END_REF]. Further, ISPs today are not part of the value chain of over-the-top applications, partly because there is no standardized way for them to offer value-added services to applications without significant offline pre-provisioning efforts. Our current work focuses on overcoming these limitations by facilitating cooperation between overlay applications and ISP layers via the ENVISION interface. 3. Bandwidth on demand. Edge devices may be connected to multiple access networks. Methods for provisioning capacity-on-demand over aggregated access networks are being investigated. 4. Dynamic QoS mapping. Application QoS requirements may be invoked and mapped dynamically to network capabilities, end user devices, and access networks.

5.

Advertisement/text insertion. Network operators may provide added-value services that also present an opportunity for monetization.

6. Content adaptation. Heterogeneous end user devices and access network capabilities will require dynamic and distributed multimedia content adaptation techniques.

In comparison to the ALTO approach, we propose to extend the scope of the information exchanged with a rich set of data, and with the introduction of a two-way information exchange. Beyond providing applications with network topology information, network operators also may want to provide information about the access network itself (e.g., type, link capabilities, coverage), current network services status (e.g., number of multicast groups, content population), as well as additional metrics such as router load or bandwidth delay. Furthermore, overlays can keep network operators apprised about the applications running within their networks and enable the optimization of data delivery by providing information on users (e.g., estimated total number of users and locations) and on content (e.g., number of different sources, bit rate, and adaptive coding).

The information exchange as well as the invocation of network services will lead to cost reductions for the network operator through more optimized content distribution. Operators will also benefit from access to application-specific information that was formerly unavailable to them, and through new revenue models by offering services to overlay applications. In this way, future network investments can be shared between the ISP and the applications that are using the network. An application, in turn, will be able to offer improved quality of service (QoS) to its users, leading to a win-win situation for all parties.

(head 2) Cooperation Between Layers

The success of emerging networked media environments is driven by services and applications that can be personalized and location-aware, and that can therefore adapt their content according to the user context and network information. The comprehensiveness and content awareness of the ENVISION interface will enable the layers to cooperate for the mutual optimization of both the overlay application and the network resources. One of the primary research challenges for this optimization is to construct an end-to-end knowledge by aggregating and harmonizing information coming from both the application and network layers, and from different network domains. At the same time, the associated management overhead introduced through this cooperation has to be kept to a minimum. In today's network, multimedia content delivery services are centrally managed by provider-side logic with limited or no knowledge of the underlying consumer's network topology and conditions, user requirements, and the usage environment. Offline pre-generation strategies have been used to cope with dynamics and heterogeneous delivery by pre-adaptation of the most popular content. However, this solution cannot meet continuously changing user interests and the huge variation in the capabilities of both heterogeneous terminals and access networks. One goal is thus the improvement of content adaptation services by enabling dynamic, distributed, and on-the-fly content adaptation that can be coordinated through the ENVISION interface. This includes real time encoding and adaptation of content for transmission to a large number of users, processing and mixing of media sources to produce derived content, content customization and tailoring, and distribution of the derived content to the consumer.

Another key challenge is the utilization of network services that are currently hidden from overlay applications. Multicast, for example, is not accessible to today's applications and service providers, mainly because it is not globally deployed. To overcome these limitations, overlay nodes can be classified through interactions via the ENVISION interface according to the underlying network capabilities. This allows the overlay application logic to identify nodes belonging to the same domain where a multicast address can be established for local content distribution, and thereby significantly reduce the load on the network. As an example, the Abilene [START_REF]Internet2 Network[END_REF] core network topology of 11 nodes and 14 links was assumed for the ISP. Assuming equal link metrics, a stream bandwidth of S Mb/s, a single source located at one of the Abilene nodes, and stream consumers at each of the other network nodes, in this particular example we found that a 58.9 percent reduction in total traffic could be achieved through the invocation of native multicast compared to traditional solutions, as detailed in Table I.

(head 2) Bicycle Race Application Use Case

Application use cases were defined to highlight the advantages of the ENVISION solution for crosslayer cooperation. One example is the live coverage of a bicycle race involving a group of athletes in a potentially remote racing environment. Local spectators, professional news and media organizations, and even the athletes themselves generate video footage of the event, e.g., from mobile phones, dedicated video cameras or helmet-cams, from their current location. End users are able to follow the race on their individual terminals, selecting viewpoints and camera angles, following individual cyclists or specific events according to their individual interest. This use case offers a good illustration of our solution for the following reasons:

1. It may span a large geographical region and different network operators. For instance, the dashed cloud in the underlay shown in Figure 2 represents the physical location of the event, where three people are filming the event (the content producers, two connected to ISP A and one to ISP B) and two content consumers are watching it.

2. The competitors and professional camera teams are constantly moving, while spectators' video sources, spread along the race track, are rather fixed.

3. End users (content consumers) can enjoy the service from various devices and access networks and may specify their restrictions/preferences as to the format of the received video, dynamically navigate across the racetrack, or define criteria for feeds of interest. 4. It enables the insertion of additional information, e.g., advertisements or details about the event's location or environment, and can support video chat between viewers.

All of these features impact content delivery and require the exchange of information between actors in order to deliver content with the required QoS, especially since in a live event, the availability of sources can change quickly in the swarm. Typically, trusted and authenticated entities in the overlay (grey boxes) can push or pull information to or from the network via the ENVISION interface (step A-1, D-1). Non-authorized nodes may retrieve non-critical information (C-1). Nodes can request the instantiation of network services (service nodes), e.g., multicast in the ISP domain D, because many peers are watching the same source (step D-2), or the activation of caches in ISP domain A, because viewers want an action replay of a certain scene (step A-2). Some overlay nodes can themselves offer services, such as a content adaptation service (C-1), to fit end user requirements.

(head 1) Conclusion

In this paper, we have introduced the concept of a new interface between overlays and the underlying network infrastructure with the goal of fostering cooperation between those layers to optimize the operation and delivery of future media applications, and thereby improve their quality. We have highlighted potential information that can be exchanged as well as network services that can be invoked through the interface.

Future work will concentrate on refining the interface specifications, developing specific protocols and optimization and adaptation algorithms, and further evaluating the approach through a combination of simulations and testbed experiments.

(head 2) Acknowledgements

This work was supported by the ENVISION project (http://www.envision-project.org), a research project partially funded by the European Commission under its Seventh Framework Programme (contract no. 248565). The authors want to thank all partners for their valuable comments and proposals in various discussions during the ENVISION project.

(head 2) References

Figure 1

 1 highlights this interface (blue arrows) between different networks and the overlay layer. The interface enables the exchange of information and the invocation of network services offered by ISPs. Numerous candidate network services could be offered via the interface. Some examples of those under development include the following. 1. Multicasting. Since applications may be spread over several ISPs, native intra-domain multicast mechanisms can be enhanced with techniques such as hybrid application-network layer multicast or operator-instantiated high capacity nodes in the network. 2. Caching. Employing specialized nodes provided by the ISP or third parties can save bandwidth in the network. Caching can also aid in the provisioning of catch-up services.

Table I .

 I are related to dynamic overlay networks, P2P networks, and information-centric networking.Dr. Mathieu has contributed to several national and European projects and has published more than 30 papers in international conferences, journals, or books. He also serves on the Technical Program Committee for several conferences, and is a senior member of the IEEE and SEE. TOUFIK AHMED is a full professor at the Institut Polytechnique de Bordeaux (IPB) in the ENSEIRB-MATMECA school of engineering. He received his B.Sc. in computer engineering with high honors from the National Institute of Computer Science (I.N.I.) Algiers, Algeria, and an M.Sc. and Ph.D. degree in computer science from the University of Versailles, France. He also obtained the HDR (Habilitation à Diriger des Recherches) degree at University of Bordeaux-1 on adaptive streaming and control of video quality of service over wired/wireless IP networks and P2P architectures. Dr. Ahmed served previously as a visiting scientist at the School of Computer Science at the University of Waterloo, and as a research fellow at the PRiSM laboratory at the University of Versailles. His primary research activities concern quality of service for multimedia wired and wireless networks, end-to-end signaling protocols, P2P networks, and wireless sensor networks. Network service example.

	Distribution mode Client -server*	Maximum link load 4S-8S	Source upload capacity 10S	Total load (sum of link load) 21S-30S
	P2P swarming**	>1S	1S	26.8S
	Native multicast	1S	1S	11S
	* Results depend on source node location
	** Assuming equal load distribution
	P2P-Point-to-point
	[1]	Cisco, "Cisco Visual Networking Index: Forecast and Methodology, 2009-2014," White Paper, June 2,
		2010.
	[2]	European Commission, Seventh Framework Programme, "NAPA WINE," FP7-ICT Project,

<http://www.napa-wine.eu/>. activities

NICO SCHWAN is a research engineer with the Alcatel-Lucent Bell Labs' Service Infrastructure

Research Team in Stuttgart, Germany. He is also an academic at the University


Christophe Chalons

Marc Massot

M Massot † ‡ §

A Vi É

ON THE EULERIAN LARGE EDDY SIMULATION OF DISPERSE PHASE FLOWS: AN ASYMPTOTIC PRESERVING SCHEME FOR SMALL STOKES NUMBER FLOWS

Keywords: disperse phase flows, large-eddy simulation, realizability, asymptotic preserving, gaussian closure AMS subject classifications

come

On the Eulerian Large Eddy Simulation of disperse phase flows: an asymptotic preserving scheme for small

Stokes number flows

1. Introduction. The simulation of disperse phase flows is nowadays of great importance in several applications, such as automotive engines, aeronautical combustors or fluidized beds. Actually, the modelling of such flows relies on the accurate description of both the continuous carrier phase, gaseous or liquid, and the discrete particulate phase, composed of particles or droplets.

In the context of small particles with respect to all carrier phase flow scales, the modelling of the carrier phase could be envisioned at a Mesoscopic level [START_REF] Fox | Large-Eddy-Simulation tools for multiphase flows[END_REF], i.e. the flow around each particle is not resolved and the coupling effects between particles and the carrier phase are modeled using Mesoscopic closures, such as the Stokes law for drag force. At this level, to solve the statistics of the disperse phase, a Population Balance Equation (PBE) on the Number Density Function (NDF) can be used. The NDF represents the probability of having a particle at a certain position of the phase space, the phase space dimensions being the relevant properties of the particles, like their position, velocity, size, temperature... To solve this equation, three approaches are possible:

• Full resolution: the PBE is directly solved by discretizing the entire phase space. This method is the most precise for compactly supported distributions, but is too expensive for unsteady configurations when the phase space of the particles is too large1 ; • Direct Monte Carlo Simulation (DMSC): the NDF is sampled by an ensemble of individual stochastic Lagrangian realizations, which are solved by means of ODEs. This approach is less expensive than the full resolution, but may also be limited by statistical convergence issues when the dimensionality of the phase space is high and many realizations are needed; • Moment methods: instead of solving for the NDF directly, moments of the NDF are solved, which are integrals over the phase space. By reducing the phase space to the physical space only, this method is computationally efficient. However it requires an additional effort in terms of modelling, as the integration step results in a loss of information. In the present work, we are interested in moment methods, because of their computational efficiency with regards to other approaches. One of the main issues with moment methods is the accurate description of the velocity distribution of the particulate phase. Actually, in turbulent flows, the velocity distribution can drastically change with the inertia of the particles, which can be quantified by the Stokes number based on the Kolmogorov time scale. For Stokes number smaller than one, the NDF is monokinetic, i.e. all particles at the same position have the same velocity, and such a distribution can be uniquely determined using zero and first order moments, i.e. density and momentum. For higher Stokes number, particles trajectories may cross, and the velocity distribution is no longer a unique Dirac δ-function, and higher order moments are needed. To handle these higher order moments, several methods can be found in the literature, and can be split into two categories. On one hand, Algebraic-Closure-Based Moment Methods (ACBMM) [START_REF] Alipchenkov | Differential and algebraic models for the second moments of particle velocity and temperature fluctuations in turbulent flows[END_REF][START_REF] Kaufmann | Comparison between Lagrangian and mesoscopic Eulerian modelling approaches for inertial particles suspended in decaying isotropic turbulence[END_REF][START_REF] Masi | Development of an algebraic-closure-based moment method for unsteady eulerian simulations of particle-laden turbulent flows in very dilute regime[END_REF][START_REF] Masi | Algebraic-closure-based moment method for unsteady eulerian simulations of non-isothermal particle-laden turbulent flows at moderate stokes numbers in dilute regime[END_REF] derive closures for the second order moments using physical and/or mathematical assumptions. On the other hand, Kinetic-Based Moment Methods (KBMM) close the system by using a presumed shape for the NDF [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of sectional and sampling methods[END_REF][START_REF] Massot | Lagrangian/Eulerian analysis of the dispersion of evaporating droplets in a non-homogeneous turbulent flow[END_REF][START_REF] Chalons | A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework[END_REF][START_REF] Kah | Eulerian quadrature-based moment models for polydisperse evapo-rating sprays[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF][START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF][START_REF] Vié | On the direct numerical simulation of moderate-stokes-number turbulent particulate flows using Algebraic-Closure-Based and Kinetic-Based Moment Methods[END_REF], which has as many parameters as the number of moments required to be controlled to describe the NDF accurately. The choice between each type of closure is motivated by the structure and the complexity of the encountered PTC, and is directly related to the number of moments.

The moment methods are a powerful tool to simulate academic configurations, but when it comes to complex configurations with a large spectrum of time and space scales, the mesh size may become too large and the computation too expensive to be achieved. To circumvent this issue, the Large Eddy Simulation is a powerful strategy: by filtering the equations in space or frequency domain, the mesh requirements can be significantly reduced. For the gas phase, the problem has been intensively studied, see [START_REF] Smagorinsky | General circulation experiments with the primitive equations. I : The basic experiment[END_REF][START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient[END_REF][START_REF] Nicoud | Using singular values to bluid a subgrid-scale model for large eddy simulations[END_REF] for example. In this work this topic is not addressed, assuming that the gas phase closures affecting the disperse phase are given, and that we are in a sufficiently dilute regime to neglect the impact of the disperse phase on the gas phase, i.e a one-way coupling regime. For the disperse phase, two types of approach can be found in the literature, which differs by the filtering procedure:

• The first approach consists in deriving moment equations of the NDF, and then filtering this system to obtain the filtered moment equations. This method was investigated in [START_REF] Moreau | Development of gas-particle euler-euler LES approach: a priori analysis of particle sub-grid models in homogeneous isotropic turbulence[END_REF], which proposes to use a Smagorinsky-like approach. It had already been applied in academic and industrial applications [START_REF] Boileau | Investigation of two-fluid methods for Large Eddy Simulation of spray combustion in Gas Turbines[END_REF][START_REF] Martinez | Large eddy simulation of fuel sprays using the eulerian mesoscopic approach. validations in realistic engine conditions[END_REF][START_REF] Sanjosé | Fuel injection model for euler-euler and euler-lagrange large-eddy simulations of an evaporating spray inside an aeronautical combustor[END_REF]. The main issue is that the connection between the NDF and the filtered moment equations is lost. Thus, the realizability of the moment equations, that is the moment are always those of a positive NDF, is hard to achieved; • The second approach consists in filtering the kinetic equation, and then deriving moment equations on the filtered NDF. This approach was first envisioned in [START_REF] Pandya | Two-fluid large-eddy simulation approach for particle-laden turbulent flows[END_REF], and has been also studied in the context of the Mesoscopic Eulerian Formalism (MEF) [START_REF] Février | Partitioning of particle velocities in gas-solid turbulent flow into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF] in [START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF]. The main interest of such an approach is that it keeps a clear link between the kinetic level, which is in fact the filtered kinetic equation, and the moment level, which is really helpful to devise realizable methods. Whatever the method used to obtain the LES equations, one of the main issues of the LES modelling is the asymptotic behavior at small Stokes number: when the Stokes number of the particles tends to zero, the resulting system admits a natural asymptotic behavior governed by an advection-diffusion on the number density, as the particles become tracers diffusing in the carrier gas phase. This asymptotic limit is important for real applications, for example in liquid-fuel combustion systems: when the droplets evaporates, they tend to this zero Stokes number, and the right dynamics has to be captured to reproduce the right fuel distribution in the system. For now, LES models have various way of dealing with this limit:

• In [START_REF] Moreau | Development of gas-particle euler-euler LES approach: a priori analysis of particle sub-grid models in homogeneous isotropic turbulence[END_REF], the authors do not treat this limit, because they consider moderate Stokes number flows, for which the limit is not predominant; • In [START_REF] Shotorban | A eulerian model for large-eddy simulation of concentration of particles with small stokes numbers[END_REF], they obtain an advection-diffusion limit, but with an isotropic diffusion coefficient, whereas the real diffusion coefficient, which is the subgrid scale tensor of the gas phase, can be anisotropic; • In [START_REF] Pandya | Two-fluid large-eddy simulation approach for particle-laden turbulent flows[END_REF] and [START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF], the authors get the right limit, that is an asymptotic advectiondiffusion asymptotic limit on the density, with a potentially anisotropic diffusion coefficient. Following these statements, in the present work we will investigate the close link between modelling and numerical methods for the LES approach of Zaichik et al. [START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF]. For this model, in the zero Stokes number limit, the flux and source terms that arises from the closures at the kinetic level become infinite, leading to strong constraints on the time step. Thus, using a global unified scheme for the whole Stokes range and recovering the proper asymptotic limit for small Stokes numbers is not straightforward, as these extreme constraint on the numerics will degrade the quality of the numerical approximation; standard numerical methods for hyperbolic system of conservation laws will not preserve the asymptotic limit or eventually lead to robustness issues.

Consequently, the goal of the present study is to propose a numerical scheme which (1) recovers the advection-diffusion asymptotic limit for the number density, referred as the asymptotic-preserving property [START_REF] Jin | Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | Numerical Schemes for Hyperbolic Conservation Laws with Stiff Relaxation Terms[END_REF] and (2) can get rid of the small time step imposed by the source term and the acoustic waves in a standard time-explicit Godunov-type method. To do so, we rely on the approach proposed by Chalons, Girardin and Kokh in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF], which is based on three ingredients:

• a Lagrange-Projection decomposition [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], that separates the terms responsible for the acoustic waves and the transport waves. This decomposition overcomes the strong restriction on the time step CFL condition coming from the large value of the sound speed when the Stokes number goes to zero by means of a time-implicit treatment of the Lagrangian system, see [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] where such an approach was first proposed; • a relaxation strategy first introduced by Suliciu [START_REF] Suliciu | On the thermodynamics of fluids with relaxation and phase transitions[END_REF] and Jin and Xin [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF], and then further studied by many authors for instance in [START_REF] Coquel | Some new godunov and relaxation methods for two-phase flow problems[END_REF][START_REF] Chalons | Navier-stokes equations with several independant pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Chalons | Relaxation approximation of the euler equations[END_REF][START_REF] Bouchut | reduced stability condition for nonlinear relaxation to conservative laws[END_REF][START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]. The motivation for using a relaxation approach is to circumvent the nonlinearities involved in the pressure law and then to make the time-implicit treatment of the Lagrangian step very low cost; • a USI (Upwinding Sources at Interfaces) approach first initiated by Cargo and Le Roux [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF], Greenberg and Le Roux [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], Gosse and Le Roux [START_REF] Gosse | A well-balanced scheme designed for inhomogeneous scalar conservation laws[END_REF], (see also Gosse [START_REF] Gosse | A priori error estimate for a well-balanced scheme designed for inhomogeneous scalar conservation laws[END_REF], Perthame and Simeoni [START_REF] Perthame | A kinetic scheme for the saint-venant system with a source term[END_REF], Jin [START_REF] Jin | A steady-state capturing method for hyperbolic systems with geometrical source terms[END_REF], R. Botchorishvili, B.Perthame and A.Vasseur [START_REF] Botchorishvili | Equilibrium schemes for scalar conservation laws with stiff sources[END_REF]...) the principle of which is to upwind the sources at interfaces in order to get the expected asymptotic-preserving property (or related well-balanced property). The form of the USI approach used here is the same as the one recently introduced in [START_REF] Chalons | Godunovtype schemes for hyperbolic systems with parameter dependent source. the case of euler system with friction[END_REF] and [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF]. The novelty of the present contribution is related to the particular context of turbulent two-phase flows, which involve singular fluxes and infinite sound speeds in the limit as well as relaxation source terms in the energy equation, a different context compared to the original studies in the field, for which we provide a specific treatment.

The paper is organized as follow. First the modelling approach is presented, based on the kinetic equation of [START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF], and on a Gaussian closure for the moment system [START_REF] Massot | Lagrangian/Eulerian analysis of the dispersion of evaporating droplets in a non-homogeneous turbulent flow[END_REF], arguing on the physical aspect of such a system, and presenting the asymptotic limit. Second, the main features of the numerical scheme are described: (1) a Lagrange-Projection [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] to decouple slow material and fast sound waves, which allows to use an explicit numerical method for material wave for the sake of precision and an implicit method on sound waves to avoid the CFL limitation, (2) a relaxation strategy [START_REF] Chalons | Relaxation approximation of the euler equations[END_REF], to avoid the non-linearity induced by the pressure law, and (3) an HLLC scheme which includes the source terms [START_REF] Gallice | Positive and entropy stable godunov-type schemes for gas dynamics and mhd equations in lagrangian and eulerian coordinates[END_REF], to recover the asymptotic limit of the moment system.

2. Moment methods for the Large Eddy Simulation of particle-laden flows. In the present work, the following assumptions are considered:

• point particles: no effect of the finite size of particles. Particles are smaller than the Kolmogorow length scale of the carrier flow, so that every flow modification induced by the particles is immediately diffused by the micromixing of the turbulence; • high Knuden number: the mean free path of particles is sufficiently large to neglect collisions; • very dilute regime: the impact of the disperse phase on the carrier phase is neglected; • fixed size solid particle: no size change nor breakup are accounted for. Consequently, in the hierarchy of modelling approaches proposed in [START_REF] Fox | Large-Eddy-Simulation tools for multiphase flows[END_REF], here the description of the flow is done at the mesoscopic level: we consider the Number Density Function (NDF) f (t, x, v) where t is the time, x the position and v the velocity of a particle. The state of the particle is then solely described by its position and velocity.

2.1. The Population Balance Equation and the moment methods. The NDF f (t, x, v) satisfies a Population Balance Equation:

∂ t f + v • ∂ x f + ∂ v • (F d f ) = 0, (2.1)
where F d is the drag force applied by the carrier phase on the disperse phase. The drag force is modeled using Stokes law:

F d (v) = u g -v τ p , (2.2)
where τ p is the relaxation time of particles.

In real applications, solving Eq. (2.1) directly is intractable for two reasons: the high dimensionality of the phase space in 3D and the wide range of scales of either carrier or disperse phases. An alternative is the macroscopic Eulerian approach: Eq. (2.1) is integrated over the velocity phase space, and conservation equations on moments are obtained:

M ijk = v i 1 v j 2 v k 3 f (t, x, v)dv, (2.3)
∂ t M ijk + ∂ x   M i+1,j,k M i,j+1,k M i,j,k+1   = - 1 τ p   (i + j + k)M i -u g •   iM i-1,j,k jM i,j-1,k kM i,j,k-1     . (2.4)
This conservation equation lives now in a 3D space. However, for every moment set of order N = i + j + k that is solved, additional higher order moments of order N + 1 are necessary. This additional moment requires a closure law, that determines (1) the accuracy of the method (2) the size N of the moment set.

In turbulent flows, the choice of the closure law is directly related to the occurrence of Particle Trajectory Crossings (PTC). When the Stokes number of the particles with respect to the Kolmogorov time scale St K = τ p /τ K is smaller than one, the velocity distribution is monokinetic i.e. the particulate phase has only one velocity per position in physical space. Several strategies exist in the literature for this range of Stokes number:

• the dusty gas approach [START_REF] Carrier | Shock waves in a dusty gas[END_REF][START_REF] Marble | Dynamics of dusty gases[END_REF]: the particles have the same velocity as the gas phase; • the Equilibrium Eulerian approach [3]: particle velocity is chosen has an expansion around the gas velocity; • the Monokinetic approach [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of sectional and sampling methods[END_REF]: the particle velocity is solved using a conservation equation for the momentum of the particulate phase. This approach had also been envisioned in a volume-averaged sense in [START_REF] Druzhinin | Direct numerical simulations of bubble-laden turbulent flows using the two-fluid formulation[END_REF]. When St K > 1, the velocity distribution becomes multivalued because of PTC, referred as Random Uncorrelated Motion (RUM) in [START_REF] Simonin | On the spatial distribution of heavy particle velocities in turbulent flow: from continuous field to particulate chaos[END_REF], and the closure laws for higher order moments has to be chosen to reproduce the NDF as accurately as possible. Moreover the choice of the NDF is not uniquely determined, and, as stated by the Hamburger Moment problem, for a finite set of moments, an infinite number of distributions is possible. So, two choices have to be made: (1) the number of moments to solve and (2) the assumptions on the unclosed moments or the shape of the NDF. In fact, this second choice splits the literature into two type of methods:

• Algebraic-Closure-Based moment methods (ACBMM) [START_REF] Masi | Development of an algebraic-closure-based moment method for unsteady eulerian simulations of particle-laden turbulent flows in very dilute regime[END_REF][START_REF] Masi | Algebraic-closure-based moment method for unsteady eulerian simulations of non-isothermal particle-laden turbulent flows at moderate stokes numbers in dilute regime[END_REF] 2 : The closure is devised based on the use of a limited information on the moments, like the total energy of the particulate phase, and a series of assumptions on the moments, for example equilibrium assumption on the deviatoric part of the stress tensor, see [START_REF] Kaufmann | Comparison between Lagrangian and mesoscopic Eulerian modelling approaches for inertial particles suspended in decaying isotropic turbulence[END_REF]; • Kinetic-Based moment methods (KBMM) [START_REF] Laurent | A hierarchy of eulerian models for trajectory crossing in particle-laden turbulent flows over a wide range of stokes numbers[END_REF]: the closure is made at the kinetic level, by choosing a presumed equilibrium shape for the NDF. The former are interesting, because it provides closed systems using a limited number of moments (density, momentum, and central energy). However, as the link between the moment and kinetic levels is not straightforward, it faces reliability issues i.e. the presumed moments do not correspond to a positive NDF. The KBMM overcome this issue by relying on a positive presumed shape for the NDF, thus realizability is intrinsically preserved. Moreover, as KBMM generates well-defined systems of equations of hyperbolic or weakly-hyperbolic type with an entropic structure [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Levermore | The Gaussian moment closure for gas dynamics[END_REF][START_REF] Chalons | Beyond pressureless gas dynamics: quadrature-based velocity moment models[END_REF][START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF], the design of realizable numerical methods is easier than for ACBMM, for which the mathematical structure is hard to determine. Thus in this work we will focus on KBMM, and our choice will be driven by the needs of the LES modelling.

The Filtered Population Balance Equation.

In real turbulent applications, the range of scale encountered in the carrier phase is so wide that solving the whole spectrum is totally unreachable with the available computational resources. To avoid this problem, LES methods filter equations in physical or frequency space. The filtering operation for a quantity φ is:

φ = G(x -x 1 , ∆)φ(x 1 , t)dx 1 , (2.5)
where • denotes a filtered quantity. The filtering of the carrier phase will decompose the velocity seen by the particle u g = u g + u ′ g where u g is the filtered carrier phase velocity and u ′ g is the carrier phase velocity fluctuation. For the disperse phase, two methods have been used in the literature:

• Filtering the moment equation 2.4, see [START_REF] Shotorban | A eulerian model for large-eddy simulation of concentration of particles with small stokes numbers[END_REF][START_REF] Moreau | Development of gas-particle euler-euler LES approach: a priori analysis of particle sub-grid models in homogeneous isotropic turbulence[END_REF];

• Filtering the PBE and get moments of the filtered NDF, see [START_REF] Pandya | Two-fluid large-eddy simulation approach for particle-laden turbulent flows[END_REF][START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF]. Theoretically the two methods lead to the same equations, as velocity and physical spaces are independent.However, the latter strategy seems more interesting because:

• the PBE is linear whereas the moment equations are non-linear;

• the filtering at the moment level loses track of the link with the underlying NDF. Consequently, filtering at the kinetic level is our choice for the present work.

Filtering Eq. (2.1), we get:

∂ t f + v • ∂ x f + ∂ v u g -v τ p f = - 1 τ p ∂ v • (u ′ g f ) r , (2.6)
where f is the filtered NDF and (u g f ) r = u g f -u g f is the subgrid correlation between the carrier phase velocity and the NDF of the particles.

To model this subgrid correlation, two contributions have been found in the literature that lead to the same closure:

• In [START_REF] Pandya | Two-fluid large-eddy simulation approach for particle-laden turbulent flows[END_REF], the authors derive the kinetic equation in direct analogy with the work of Reeks [START_REF] Reeks | On a kinetic equation for the transport of particles in turbulent flows[END_REF][START_REF] Reeks | On the continuum equations for dispersed particles in nonuniform flows[END_REF] who use the Lagrangian History Direct Interaction (LHDI); • In [START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF], the authors follow the same strategy but based on the work of Zaichik et al. [START_REF] Zaichik | A statistical model of particle transport and heat transfer in turbulent shear flows[END_REF], which consider the impact of the gas phase on the disperse phase to be a Gaussian process, and thus make use of the Furutsu-Donsker-Novikov Formula [START_REF] Novikov | Functionals and the random-force method in turbulence theory[END_REF]. Following [START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF] , the closure for the subgrid correlation in the kinetic equation is:

- 1 τ p ∂ v • (u ′ g f ) r = ∂ v • µ∂ v f + λ∂ x f , (2.7)
where:

λ ij = g r u τ g,ij , µ ij = τ g,ik f r u τ p + l r u ∂ x k u j , (2.8)
where f r u , g r u and l r u are response coefficients (see [START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF] for details) whose long-time values are:

f r u = 1 1 + St , g r u = 1 St(1 + St) , l r u = 1 St(1 + St) 2 .
(2.9)

Here, the Stokes number is St = τ p /τ r L where τ r L is the Integral time scale of the subgrid eddies.

The resulting kinetic equation in non-dimensional form is:

∂ t * f + v * • ∂ x * f + ∂ v * • u * g -v * St f = ∂ v * • µ * • ∂ v * f + λ * • ∂ x * f , (2.10)
where:

µ * ij = τ * g,ik f r u St δ kj + l r u ∂u * g,j ∂x * k , λ * ij = g r u τ * g,ij . (2.11)
In the following, star exponents are dropped for the sake of clarity.

2.3. Asymptotic limit of the filtered PBE for small Stokes number. To determine the equilibrium state and the asymptotic limit of the filtered kinetic equation Eq. (2.10), The Chapman-Enskog expansion is used [START_REF] Chapman | The mathematical theory of non-uniform gazes[END_REF]. Contrary to the one proposed in [START_REF] Alipchenkov | Differential and algebraic models for the second moments of particle velocity and temperature fluctuations in turbulent flows[END_REF], it is performed in the classical way, assuming a decomposition of the solution into power of St:

f = f 0 + Stf 0 φ 1 + O St 2 .
(2.12)

Rewriting Eq. (2.10)

∂ t f + v.∂ x f = 1 St J (f ), (2.13)
where J (f ) is the relaxation operator:

J (f ) = ∂ v • (v -u g )f + Stµ∂ v f + Stλ∂ x f . (2.14)
Inserting Eq. (2.12) into Eq. (2.13), and grouping terms by powers of St:

∂ t f 0 + v • ∂ x f 0 + St ∂ t f 0 φ 1 + v • ∂ x f 0 φ 1 = 1 St J (f 0 ) + J (f 0 φ 1 ) + O(St). (2.15)
For small Stokes number, the zeroth order of Eq. (2.15) is:

J (f 0 ) = Stλ∂ x f 0 + Stµ∂ v f 0 -(u g -v) f 0 = 0. (2.16)
The solution of such an equation has the form:

f 0 (x, u) = Γ(x, v) ρ (2π) 3/2 |τ g | exp -(v -u g )τ -1 g (v -u g ) T , (2.17)
where ρ = vf dv and Γ(x, v) is an unknown function that is equal to 1 in the case of an homogenous NDF in space. The zeroth order distribution leads to the following equation on density:

∂ t ρ + ∂ x • ρu 0 = 0, (2.18)
where u 0 = vf 0 dv. To find u 0 , the zeroth order moment of Eq. (2. [START_REF] Chalons | Relaxation approximation of the euler equations[END_REF]) is taken and leads to:

u 0 = u g - τ g ρ ∂ x ρ, (2.19)
so that the asymptotic limit at the kinetic level leads to the following asymptotic limit at the moment level:

∂ t ρ + ∂ x • ρu g = ∂ x • (τ g ∂ x ρ) .
(2.20)

One important thing to notice at this point is that this asymptotic limit is not a consequence of the modelling approach we choose, this is the limit that every LES model for the particulate phase has to recover in the low Stokes number limit. It is also worth noticing that the resulting moment system does not depend on any closure for the disperse phase, but only on the closure for the gas phase. Thus, any moment method that treats the number density should have the same asymptotic limit. To achieve the description of the low Stokes number limit, taken the first order moment of Eq. (2.16), the asymptotic limit of the internal energy

σ 0 = 1 2 v ⊗ vf 0 dv -1 2 u ⊗ u is recovered: σ 0 = τ g 1 -∂ x ⊗ u g + ∂ x ⊗ τ g ρ ∂ x ρ , (2.21)
where ⊗ is the symmetric tensor outer product.

2.4. Moment equations for LES. After obtaining the kinetic equation, one has to go up to the moment level, and thus to choose the number of moments to solve and the closure law. Our choice is driven by the asymptotic limit of the PBE. Actually, the equilibrium distribution is a perturbation of a Gaussian distribution. Therefore, the minimal model to reproduce this asymptotic limit has to be a Gaussian distribution. Following [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF], where the Anisotropic Gaussian distribution is used to close the moment system in a DNS context, we will use moments up to second order, that is 10 moments in 3D. The resulting system of equation is the following3 :

∂ t ρ + ∂ x • (ρu) = 0 (2.22) ∂ t (ρu) + ∂ x • (ρuu T + P) = ρ u g -u St (2.23) ∂ t (ρE) + ∂ x • ((ρE + P) ⊗ u) = ρ u g ⊗ u + Stµ -2E St , (2.24)
where

E ij = 1 2 ũi ũj + σij 2
is the total energy tensor, and P ij the pressure tensor:

ρσ ij = [(v i -u i )(v j -u j )] f (t, x, v)dv (2.

25)

P ij = ρ (σ ij + λ ij ) .
(2.26)

2.5. Numerical issues associated to the low Stokes number limit. In the following, the study will be fully performed in 1D for the sake of simplicity, even if all the developments can be straightforward envisaged in 3D. Moreover, we consider an homogeneous gaseous flow field, to focus on the main difficulties of the moment method.

The resulting system of equation is:

       ∂ t ρ + ∂ x ρu = 0, ∂ t ρu + ∂ x ρu 2 + P = -ρ u -u g St , ∂ t ρE + ∂ x (ρuE + uP ) = -ρu u -u g St -ρ 2ǫ -Stµ St (2.27)
where E = u 2 /2 + ǫ is the total energy and P = ρ(2ǫ + λ) is the pressure. Observe from now on that the last equation can be formulated in terms of internal energy and for smooth solutions as follows

∂ t ǫ + u∂ x ǫ + P ρ ∂ x u = - 2ǫ -Stµ St .
Recall that λ and µ are given by

µ = τ g St(1 + St) , λ = τ g St (1 + St) . (2.28)
Therefore, it is clear that the subgrid scale effects represented by τ g = 0 appear in two terms: a relaxation term on the internal energy, which tends to relax the internal energy of the disperse phase towards the subgrid scale internal energy of the gas phase (see paragraph below on the asymptotic analysis), and a pressure-like term via λ, acting like a flux and which propagates the time and space correlation of the gas phase over the disperse phase. Characteristic analysis. Easy manipulations show that (2.27) equivalently writes for smooth solutions

         ∂ t ρ + u∂ x ρ + ρ∂ x u = 0 ∂ t u + (2ǫ + λ) ρ ∂ x ρ + u∂ x u + 2∂ x ǫ = - u -u g St ∂ t ǫ + P ρ ∂ x u + u∂ x ǫ = - 2ǫ -Stµ St . ( 2

.29)

Focusing on the convective part only (without source terms), the characteristic velocities are then given by

λ 1 = u -c < λ 2 = u < λ 3 = u + c,
where the speed of sound c is given by c = √ 6ǫ + 3λ. Observe that the square of the sound speed is made of two contributions, namely the classical one 6ǫ = γ(γ -1)e associated with the perfect gas equation of state with adiabatic coefficient γ = 3, and a new contribution involving the subgrid stress of the gas phase τ g via λ. Importantly, the presence of this new contribution (τ g = 0) makes the sound speed c tend to infinity when the Stokes number St goes to zero. From a physical point of view, it simply means that the smaller the Stokes number St is, the faster the time and space correlations of the gas phase propagate over the dispersed phase. From a numerical point of view, this property is expected to give a strong constraint on the time step for a fully explicit in time Godunov-type method according to the well-known CFL stability condition.

Asymptotic analysis. Here we consider the asymptotic behavior of the model when the Stokes number St tends to 0. Multiplying the last two equations of (2.29) by St (1 + St) and letting St go to 0 gives

u = u g - τ g ρ ∂ x ρ, (2.30) 2ǫ = τ g ∂ x τ g ρ ∂ x ρ + 1 . (2.31)
Combining these two relations with the first equation of (2.29) (written in conservation form) then gives the following asymptotic limit

∂ t ρ + ∂ x ρu g = ∂ x (τ g ∂ x ρ) , (2.32)
2ǫ = τ g ∂ x τ g ρ ∂ x ρ + 1 . (2.33)
The number density then evolves according to an advection-diffusion equation, which is coherent with the asymptotic limit at the kinetic level. If one want to solve small Stokes number flows only, it is easy to directly solve Eq. (2.32). If one want to be able to solve a large range of Stokes numbers with the same method, Eq. (2.27) is needed. However, in the low Stokes number limit, Eq. (2.27) encounter two strong and limiting effects:

• the sound speed c and the source terms tend to infinity: it leads to strong constraints for explicit numerical schemes, as well as high numerical diffusion for material waves u. • the numerical method has to recover the asymptotic limit at low Stokes number, which is not possible with "classical" schemes. In the following, a numerical scheme is proposed to handle these two features of the chosen methodology.

3. Numerical scheme. The aim of this section is to describe a numerical scheme which (1) recovers the advection-diffusion asymptotic limit for the number density, referred as the asymptotic-preserving property [START_REF] Jin | Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | Numerical Schemes for Hyperbolic Conservation Laws with Stiff Relaxation Terms[END_REF] and (2) can get rid of the small time step imposed by the source term and the acoustic waves in a standard time-explicit Godunov-type method. We rely on the approach proposed by Chalons, Girardin and Kokh in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] and based on three ingredients:

• a Lagrange-Projection decomposition [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF];

• a relaxation strategy first introduced by Suliciu [START_REF] Suliciu | On the thermodynamics of fluids with relaxation and phase transitions[END_REF] and Jin and Xin [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF], and then further studied by many authors for instance in [START_REF] Coquel | Some new godunov and relaxation methods for two-phase flow problems[END_REF][START_REF] Chalons | Navier-stokes equations with several independant pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Chalons | Relaxation approximation of the euler equations[END_REF][START_REF] Bouchut | reduced stability condition for nonlinear relaxation to conservative laws[END_REF][START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]; • a USI (Upwinding Sources at Interfaces) approach first initiated by Cargo and Le Roux [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF], Greenberg and Le Roux [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], Gosse and Le Roux [START_REF] Gosse | A well-balanced scheme designed for inhomogeneous scalar conservation laws[END_REF] (see also for related works in a wide range of applications [START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF], [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF], [START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF], [START_REF] Gosse | Asymptotic-preserving and well-balanced schemes for radiative transfer and the rosseland approximation[END_REF], [START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF],

[5], [START_REF] Carillo | Simulation of fluid and particles flows: Asymptotic preserving schemes for bubbling and flowing regimes[END_REF], [START_REF] Degond | An asymptotic preserving scheme for the euler equations in a strong magnetic field[END_REF], [START_REF] Filbet | A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources[END_REF], [START_REF] Després | Design of asymptotic preserving finite volume schemes for the hyperbolic heat equation on unstructured meshes[END_REF], [START_REF] Bouchut | Upwinding of the source term at interfaces for euler equations with high friction[END_REF], [START_REF] Berthon | Late-time relaxation limits of nonlinear hyperbolic systems. A general framework[END_REF][START_REF] Abgrall | Special issue : Numerical approximations of hyperbolic systems with source terms and applications[END_REF], [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic euler and navier-stokes equation[END_REF] and the references therein).

The next three subsections give only a brief description of these three main ingredients as details can be found for instance in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF]. The objective is to provide the reader with the most important update formulas and to focus on the specific treatment of the singular fluxes as well as of the source term associated with the internal energy, which are not included in the upwinding at the interfaces, receives a treatment through a splitting like approach and still satisfied the objective of being asymptotic preserving on energy.

3.1. Lagrange-Projection and source terms decomposition. The objective of the Lagrange-projection is to decompose the full system into two sub-systems using a chain rule argument on the space derivatives. The first one only involves the transport wave, and the second one involves the acoustic waves (and the source terms). This kind of approach can be seen as an operator splitting strategy, and the main interest is to enable the use of different solvers for each subsystem. In this work, we shall consider explicit-explicit solvers but also implicit-explicit solvers in order to get rid of the strong CFL stability restriction imposed by the sound speed. By implicit-explicit, we mean here implicit on the acoustic waves and source terms and explicit on the transport part following the same approach as in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] (see [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]). Note that the source term associated with the internal energy in the last equation of (2.27) being not considered in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF], we shall treat it separetely here using again a splitting strategy. Therefore, we shall end with three sub-systems to be treated numerically. Using the property ∂ x ρuX = u∂ x ρX + ρX∂ x u for X = {1, u, E} in the full system (2.27), we get the following transport system:

   ∂ t ρ + u∂ x ρ = 0 ∂ t ρu + u∂ x ρu = 0 ∂ t ρE + u∂ x ρE = 0, (3.1)
which corresponds to the advection of each field with velocity u. The acoustic and source term part is then given by

       ∂ t ρ + ρ∂ x u = 0 ∂ t ρu + ρu∂ x u + ∂ x P = -ρ u -u g St ∂ t ρE + ρE∂ x u + ∂ x P u = -ρu u -u g St -ρ 2ǫ -Stµ St . (3.2)
Introducing τ = 1/ρ and τ ∂ x = ∂ m this sytem equivalenty writes

       ∂ t τ -∂ m u = 0 ∂ t u + ∂ m P = - u -u g St ∂ t E + ∂ m P u = -u u -u g St - 2ǫ -Stµ St . (3.3)
Now splitting the drag force and energy relaxation terms leads to

       ∂ t τ -∂ m u = 0 ∂ t u + ∂ m P = - u -u g St ∂ t E + ∂ m P u = -u u -u g St , (3.4)
and

     ∂ t τ = 0 ∂ t u = 0 ∂ t ǫ = - 2ǫ -Stµ St . (3.5)
It is now a matter of approximating the three systems (3.4), (3.1) and (3.5).

3.2. Relaxation procedure. As mentioned above, system (3.4) may receive an implicit treatment in order to reduce the CFL stability restriction. In order for this implicit procedure to be cheap we propose to use a relaxation strategy, the objective of which being to overcome the difficulties introduced by the pressure nonlinearities. More precisely, the main idea is to consider the pressure P as a new unknown Π of the system and to solve for an additional equation for this new variable. We propose the following relaxation system

           ∂ t τ -∂ m u = 0 ∂ t u + ∂ m Π = - u -u g St ∂ t E + ∂ m Πu = -u u -u g St ∂ t Π + a 2 ∂ m u = 0. (3.6)
According to the very classical relaxation method, see for instance [START_REF] Suliciu | On the thermodynamics of fluids with relaxation and phase transitions[END_REF], [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF], [START_REF] Chalons | Relaxation approximation of the euler equations[END_REF], [START_REF] Bouchut | reduced stability condition for nonlinear relaxation to conservative laws[END_REF], this system will be solved instead of (3.4) at each time step and with Π = P at initial time (the new relaxation variable Π is said to be at equilibrium). It is worth noticing that the convective part of (3.6) has three constant eigenvalues given by -a, 0 and a so that the associated characteristic fields are linearly degenerate. Note also that a has to be chosen sufficiently large according to the subcharacteristic condition a > max (ρc) in order to avoid instablities in this relaxation procedure. We refer for instance the reader to [START_REF] Chalons | Navier-stokes equations with several independant pressure laws and explicit predictor-corrector schemes[END_REF] and [START_REF] Coquel | Relaxation of fluid systems[END_REF] for a rigorous proof and for more details. To conclude this short section, let us introduce the new variables ←w and -→ w defined by ←w = Π-au and -→ w = Π+au. These quantities are easily shown to be strong Riemann invariants associated with the characteristic speeds ±a and allow to equivalently write (3.6) as follows :

               ∂ t τ -∂ m u = 0 ∂ t - → w + a∂ m - → w = - u -u g St ∂ t ← - w -a∂ m ← - w = u -u g St ∂ t E + ∂ m Πu = -u u -u g St .
(3.7)

3.3. Notations and overview of the numerical scheme. In this section, we introduce some notations and present the proposed numerical method for approximating the solutions of (2.27). This method can be understood as an operator-splitting strategy consisting in solving the three systems (3.4), (3.5) and (3.1) in this order. Recall that (3.4) will be solved using the relaxation system (3.7). Let be given a constant time step ∆t and a constant space step ∆x. Introducing x j+1/2 = (j + 1/2)∆x for j ∈ Z and t n = n∆t for n ∈ N, the approximate solution of (2.27), say u λ (x, t), with u 0 as initial data is classically sought as a piecewise constant function on each slab C n j = [x j-1/2 ; x j+1/2 ) × [t n ; t n+1 ) for (n, j) ∈ N × Z, and we set u λ (x, t) = u n j for all (x, t) ∈ C n j , j ∈ Z, n ∈ N.

Here λ stands for the ratio λ = ∆t/∆x, and for completeness we set

u 0 j = 1 ∆x x j+1/2 x j-1/2 u 0 (x)dx.
Assuming that the approximate solution u λ (x, t n ) is known, we propose to advance it to the next time level t n+1 according to a classical splitting strategy between (3.4), (3.5) and (3.1).

First step (t n → t n+1= ) This step consists in solving the Lagrangian system (3.4), that is to say

∂ t v + ∂ m g(v) = S(v), (3.8
)

with v = (τ, u, E) t , g(v) = (-u, P, P u) t and S(v) = (0, - u -u g St , -u u -u g St ) t . Note
that the change of variables u → v is one-to-one. We consider v λ (x, t n ) (easily calculated from u λ (x, t n )) as initial data. To define u n+1= j , we propose to use a Godunov-type method based on the definition of an approximate Riemann solver for (3.7) as detailed in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF]. The proposed approximate Riemann solver here coincides exactly with the one proposed in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] up to very minor modifications : the gravity is not taken into account in the present paper and the friction terms make the velocity tend to u g here instead of 0. This approximate Riemann solver includes the source terms in its definition in the sense that the definition of the intermediade states actually depends on the source term S(v) and in particular on the Stokes number St. Therefore, the source term S(v) is taken into account at each interface of the mesh, in agreement with the celebrated USI (Upwinding Sources at Interfaces) approach introduced in [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF], [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], [START_REF] Gosse | A well-balanced scheme designed for inhomogeneous scalar conservation laws[END_REF]. Note however that the way the sources are upwinding here is really different and relies on the formalism of consistency in the integral sense introduced by Gallice in [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF], [START_REF] Gallice | Positive and entropy stable godunov-type schemes for gas dynamics and mhd equations in lagrangian and eulerian coordinates[END_REF] for systems of conservation laws with sources. This formalism is nothing but a generalization of the well-known Harten, Lax and van Leer formalism [START_REF] Harten | On upstream differencing and godunovtype schemes for hyperbolic conservation laws[END_REF] for systems of conservation laws. We then refer the reader to [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] for more details.

Time-explicit Godunov-type scheme. Averaging the approximate Riemann solutions defined at each interface leads to the following Godunov-type update formulas

                                       τ n+1= j = τ n j + ∆t ∆m j (u * j+ 1 2 -u * j-1 2 ), - → w n+1= j = - → w n j -a ∆t ∆m j ( - → w n j -- → w n j-1 ) + ∆t a ∆m j-1/2 ∆m j u g -u * j-1 2 St , ← - w n+1= j = ← - w n j + a ∆t ∆m j ( ← - w n j+1 -← - w n j ) -∆t a ∆m j+1/2 ∆m j u g -u * j+ 1 2 St , E n+1= j = E n j - ∆t ∆m j ((up) * j+ 1 2 -(up) * j-1 2 ) + ∆t ∆m j u g ∆m j+1/2 u * j+ 1 2 + ∆m j-1/2 u * j-1 2 2St - ∆t ∆m j ∆m j+1/2 (u * j+ 1 2 ) 2 + ∆m j-1/2 (u * j-1 2 ) 2
2St , (3.9) where we have set ∆m j = ρ n j ∆x, ∆m j+1/2 = ∆m j + ∆m j+1 2 , and

u * j+ 1 2 = St 2aSt + ∆m j+1/2 ( - → w n j -← - w n j+1 +u g ∆m j+1/2 St ), p * j+ 1 2 = - → w n j + ← - w n j+1 2 . (3.10)
Using the relation

u n+1= j = - → w n+1= j -← - w n+1= j 2a ,
we also have

                             τ n+1= j =τ n j + ∆t ∆m j (u * j+ 1 2 -u * j-1 2
),

u n+1= j =u n j - ∆t ∆m j (p * j+ 1 2 -p * j-1 2 ) + ∆t ∆m j ∆m j-1/2 (u g -u * j-1 2 ) + ∆m j+1/2 (u g -u * j+ 1 2 ) 2St , E n+1= j =E n j - ∆t ∆m j ((up) * j+ 1 2 -(up) * j-1 2 ) + ∆t ∆m j u g ∆m j+1/2 u * j+ 1 2 + ∆m j-1/2 u * j-1 2 2St - ∆t ∆m j ∆m j+1/2 (u * j+ 1 2 ) 2 + ∆m j-1/2 (u * j-1 2 ) 2
2St .

(3.11) This scheme is shown to be stable under the following CFL condition

max j∈Z a ∆t ∆m j ≤ 1 2 ,
see again for instance [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF]. Using classical notations, one has thus defined v λ (x, t n+1= ). This piecewise constant solution will be used in the second step as a natural initial condition.

Time-implicit Godunov-type scheme. In order to get rid of the above CFL restriction on the time step ∆t, which becomes stronger and stronger as the Stokes number goes to zero, a time-implicit version can be simply defined as follows,

                                       τ n+1= j = τ n j + ∆t ∆m j (u * j+ 1 2 -u * j-1 2 ), - → w n+1= j = - → w n j -a ∆t ∆m j ( - → w n+1= j -- → w n+1= j-1 ) + ∆t a ∆m j-1/2 ∆m j u g -u * j-1 2 St , ← - w n+1= j = ← - w n j + a ∆t ∆m j ( ← - w n+1= j+1 -← - w n+1= j ) -∆t a ∆m j+1/2 ∆m j u g -u * j+ 1 2 St , E n+1= j = E n j - ∆t ∆m j ((up) * j+ 1 2 -(up) * j-1 2 ) + ∆t ∆m j u g ∆m j+1/2 u * j+ 1 2 + ∆m j-1/2 u * j-1 2 2St - ∆t ∆m j ∆m j+1/2 (u * j+ 1 2 ) 2 + ∆m j-1/2 (u * j-1 2 ) 2 2St , (3.12
) where u *

u * j+ 1 2 = St 2aSt + ∆m j+1/2 ( - → w n+1= j -← - w n+1= j+1 +u g ∆m j+1/2 St ), p * j+ 1 2 = - → w n+1= j + ← - w n+1= j+1 2 . (3.13)
Let us notice that the second and third equations do not depend on τ and E and can thus be solved independenlty. The corresponding system is easily seen to be linear with a pentadiagonal and strictly diagonally dominant matrix. Therefore, ←w n+1= j and -→ w n+1= j are uniquely defined for any ∆t. Then, τ n+1= j and E n+1= j follow explicitly (which makes the overall strategy very low cost) thanks to the first and fourth equations of (3.12).

About the Asymptotic-Preserving property. Let us discuss in this short paragraph the asymptotic behavior of the proposed numerical scheme when the Stokes number St goes to zero. We focus on the explicit-explicit version for simplicity. Let us first observe that λSt goes to τ g when St goes to zero by definition of λ, so that StP goes to τ g ρ by definition of P . We then easily get by definition of ←w and -→ w the following

limit lim St→0 u * j+ 1 2 = 1 ∆m j+1/2 (τ g ρ n j -τ g ρ n j+1 + u g ∆m j+1/2 ), that is to say lim St→0 u * j+ 1 2 = u g - 2τ g ρ n j + ρ n j+1 ρ n j+1 -ρ n j ∆x .
This limit is clearly consistent with the expected mass flux

u = u g - τ g ρ ∂ x ρ, in (2.30).
Regarding the asymptotic behavior of the internal energy (2.33), it will be useful in the next steps of the method to estimate the asymptotic behavior when the Stokes number goes to zero of the following discrete time derivative

St ǫ n+1= j -ǫ n j ∆t . (3.14)
With this in mind, it is first easily shown from the second equation of (3.11) that Stu n+1= j goes to zero with the Stokes number (the calculations are left to the reader). Then letting St go to zero in the third equation of (3.11) shows that (3.14) asymptotically behaves like

St - ∆t ∆m j ((up) * j+ 1 2 -(up) * j-1 2 ) + ∆t ∆m j u g ∆m j+1/2 u * j+ 1 2 + ∆m j-1/2 u * j-1 2 2St - ∆t ∆m j ∆m j+1/2 (u * j+ 1 2 ) 2 + ∆m j-1/2 (u * j-1 2 ) 2 2St .
Using the properties that u *

-(up) * j-1 2 = 1 2 u * j+ 1 2 + u * j-1 2 p * j+ 1 2 -p * j-1 2 + 1 2 p * j+ 1 2 -p * j-1 2 u * j+ 1 2 -u * j-1 2 ,
we easily get that (3.14) is asymptotically consistent with -τ g ∂ x u when St goes to zero.

Second step (t n+1= → t n+1-) The second step consists in solving the transport step (3.1) with v λ (x, t n+1= ) as initial data, which can be easily and equivalently transformed in terms of the u = (ρ, ρu, ρE) t variables. Following [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] (see also [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF]), we consider a very classic upwind and time-explicit numerical scheme given by

X n+1- j = X n+1= j + ∆t ∆x u * ,+ j-1/2 X n+1= j-1 + u * ,- j+1/2 -u * ,+ j-1/2 X n+1= j -u * ,- j+1/2 X n+1= j+1 (3.15)
where α ± = (α ± |α|)/2 for any α and X ∈ {ρ, ρu, ρE}. This scheme is shown to be stable under the following CFL condition

∆t ∆x u * ,+ j-1/2 -u * ,- j+1/2 ≤ 1. (3.16)
About the Asymptotic-Preserving property. Let us go on with the asymptotic-preserving property. The first equation of (3.15) writes

ρ n+1- j = ρ n+1= j + ∆t ∆x u * ,+ j-1/2 ρ n+1= j-1 + u * ,- j+1/2 -u * ,+ j-1/2 ρ n+1= j -u * ,- j+1/2 ρ n+1= j+1
, that we can combine with the first equation of (3.9), namely

τ n+1= j = τ n j 1 + ∆t ∆x (u * j+ 1 2 -u * j-1 2 ) ,
to give the conservative update formula

ρ n+1- j = ρ n j - ∆t ∆x {ρu} j+ 1 2 -{ρu} j-1 2 , (3.17)
with numerical flux given by

{ρu} j+ 1 2 = ρ n+1= j u * ,+ j+1/2 + ρ n+1= j+1 u * ,- j+1/2 .
We have just seen in the previous step that u * Third step (t n+1-→ t n+1 ) The last step consists in solving (3.5) with v λ (x, t n+1-) as initial data. Compared to the work of [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF], this step is the main difference, as no internal energy relaxation exists in their work. Here we choose to decouple this relaxation step because the leading order for the asymptotic limit is on the mean velocity component, such as it is not mandatory to include this step in the HLLC solver. To define v λ (x, t n+1 ), it is a natural idea to use a classical pointwise implicit evaluation of the source term. More precisely, it amounts to set

j+ 1 2 is consistent with u g - τ g ρ ∂ x ρ
       τ n+1 j = τ n+1- j , u n+1 j = u n+1- j , ǫ n+1 j = ǫ n+1- j -∆t 2ǫ n+1 j -Stµ St ,
or equivalently

       τ n+1 j = τ n+1- j , u n+1 j = u n+1- j , ǫ n+1 j = St(µ∆t + ǫ n+1- j ) St + 2∆t . ( 3
= St ǫ n+1= j -ǫ n j ∆t + St ǫ n+1- j -ǫ n+1= j ∆t + Stµ.
Since lim St→0 Stµ = τ g (1 + ∂ x u g ), using the results obtained in the previous two steps then clearly gives that 2ǫ n+1 j is asymptotically consistant with -τ g ∂ x u + τ g + τ g ∂ x u g , that is to say with (2.33) by the relation u = u g -τ g ρ ∂ x ρ.

3.4. Main properties. We gather in this section the main properties of the proposed algorithm. We focus on the implicit-explicit version. Theorem 3.1. Under the CFL condition (3.16) and provided that a is chosen sufficiently large, the implicit-explicit in time numerical scheme is well defined and satisfies the following stability properties:

(i) it is a conservative scheme for the density ρ. It is also a conservative scheme for ρu and ρE when the source terms are omitted, (ii) the density ρ n j is positive for all j and n > 0 provided that ρ 0 j is positive for all j, (iii) it is asymptotic preserving.

Sketch of the proof. Property (i) has been proved in the course of the description of the second step for the density, see (3.17). The proof is similar for ρu and ρE. Property (ii) is obtained from standard manipulations [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]. The asymptoticpreserving property has been proved in the previous subsections.

Remark. The validity of an entropy inequality has been proved in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] for the explicitexplicit version of the scheme. Regarding the implicit-explicit version, we refer the reader to [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF].

Numerical results.

Here we propose a test case which highlights the effects of the subgrid scale of turbulence on the disperse phase. The domain is [-1, 1]. The initial state for the disperse phase is a spatially Gaussian distribution ρ(t = 0, x) = 1+ exp -x 2 /(2σ 2 0 ) , where σ = 0.01, which mimics the dispersion of particles occurring in a turbulent field. At time t = 0, the particles are at rest, i.e. u(t = 0, x) = 0 and ǫ(t = 0, x) = 0. The gas phase is constant in time and space, and is decomposed into a mean velocity u g = 0, which in fact generates no fluxes, and a subgrid scale energy τ g = 1/10.

For the explicit schemes, two different explicit constraints are imposed on the time step through the source terms and the CFL number:

∆t source ≤ St 2 , ∆t CF L ≤ CFL ∆x c = CFL ∆x √ 6ǫ + 3λ . ( 4.1)
At time t = 0s, the time step is then:

∆t = min St 2 , CFL 3τ g ∆x St(1 + St) . (4.2)
In Figs. 4.1-4.2, results at time t = 0.2s for the density, the mean velocity and the internal energy are plotted against the position, for the Explicit non-AP and AP schemes, and for St = 0.1 and St = 0.01, using 100 cells. The Stokes number is not sufficiently small to use the asymptotic solution as a reference. So a 2000-cell solution is use as a reference, for which a sufficient mesh convergence has been verified for the sake of the comparisons. For these Stokes number and this space discretization, the differences between non-AP and AP schemes are not obvious, as the contribution of the subgrid flux and source terms are not predominant.

In Figs. 4.3-4.4, results at time t = 0.2s for the density, the mean velocity and the internal energy are plotted against the position, for the Explicit non-AP and AP schemes, and for St = 0.001 and St = 0.0001, using 100 cells. Here the Stokes numbers are small enough for the time step to be CFL-driven. The Stokes number is now sufficiently small to use the asymptotic solution as a reference, the statement having been verified. Now the differences between non-AP and AP scheme are clearly highlighted: for both Stokes numbers, the AP scheme matches the asymptotic solution whereas the non-AP scheme comes up with a significant deviation which increases while the Stokes number decreases. The quality of the results are also quantitatively assessed in Fig. 4.5, in which the error on the density against the number of cells is plotted for the non-AP and the AP schemes and St = 0.0001. It shows two orders of magnitude between the two schemes, definitely demonstrating the necessity of AP schemes in such regimes, and the quality of the proposed explicit one.

However, even if the explicit AP scheme obtains good results for every Stokes number, it is constrained by the explicit description of acoustic waves. To avoid such a constraint, the Implicit scheme is a solution. We recall that by Implicit, we mean that the acoustic and source terms are addressed implicitly, where the material wave are still addressed explicitly.

In Figs. 4.6-4.7, the results for the Implicit AP and Implicit non AP schemes are presented, for St = 0.0001 and for a time step which is respectively 10 times and 50 times larger than the explicit time step. It is shown, even using big time steps comparing to the explicit one, increasing the error, the Implicit AP scheme still matches the asymptotic solution, where the non AP one obtains less precise results. The qualitative evaluation of the error on the density plotted in Fig. 4.8 demonstrates it again, the AP scheme being impacted by larger time steps but keeping a low error comparing to the explicit one.

Conclusions.

In this work, a comparison of existing fully Eulerian strategies for LES has been presented. This comparison highlighted two important aspects of each method: the way the moment equations are closed and the asymptotic behavior at small Stokes number.

Concerning the closure of the moment equations, on the one side ACBMM can not always ensure the realizability of the moments, because of the lack of a clear link between the moment closure and an underlying NDF. On the other side, KBMM are intrinsically realizable, as long as the underlying presumed NDF is realizable, and accurate and stable numerics are used.

Regarding the available method for particle-laden flow LES, a Kinetic-Based LES method based on the work of [START_REF] Zaichik | An eulerian approach for large eddy simulation of particle transport in turbulent flows[END_REF] and [START_REF] Pandya | Two-fluid large-eddy simulation approach for particle-laden turbulent flows[END_REF] has been chosen. The resulting moment method has two main drawbacks: when the Stokes number tends to zero, the source terms become infinite and in the case of Euler equations, the sound speed too. Moreover the system of equations must tend to a clear advection-diffusion asymptotic limit. To handle it, a new numerical scheme has been designed based on the work of [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] and constituted of a Lagrangian-Projection, a relaxation method and a HLLC scheme with source terms, with a specific treatment of the peculiarities of the models we work with, that is the singular behavior of the fluxes and the inclusion of energy relaxation. The resulting scheme is proven to be Asymptotic-Preserving for the small Stokes number limit in 1D cases.

The next step is to use the strategy for the simulation of multi-dimensional flows, especially by considering an Anisotropic Gaussian closure [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] to recover the full details of the asymptotic limit. The adaptation of the scheme to space varying Stokes number (for example in evaporating flows) has also to be investigated.

j+ 1 2 and Stp * j+ 1 2

 21 are respectively consistent with u = u g -τ g ρ ∂ x ρ and τ g ρ, together with the Leibniz relation (up) * j+ 1 2

Fig. 4 . 1 .

 41 Fig. 4.1. Stokes number=0.1: Comparison between refined solution with 2000 cells (black full line), Explicit non-AP scheme (squares) and Explicit AP scheme (crosses) with 100 cells at time t = 0.2s.

Fig. 4 . 2 .Fig. 4 . 3 .

 4243 Fig. 4.2. Stokes number=0.01: Comparison between refined solution with 2000 cells (black full line), Explicit non-AP scheme (squares) and Explicit AP scheme (crosses) with 100 cells at time t = 0.2s.

Fig. 4 . 4 .

 44 Fig. 4.4. Stokes number=0.0001: Comparison between asymptotic analytic solution (black full line), Explicit non-AP scheme (squares) and Explicit AP scheme (crosses) with 100 cells at time t = 0.2s.

Fig. 4 . 5 .

 45 Fig. 4.5. Stokes number=0.0001: Error on the number density comparing to the asymptotic analytic solution using the Explicit non-AP scheme (squares) and the Explicit AP scheme (crosses).

Fig. 4 . 6 .Fig. 4 . 7 .

 4647 Fig. 4.6. Stokes number=0.0001: Comparison between asymptotic analytic solution (black full line), Implicit non-AP scheme (squares) and Implicit AP scheme (crosses) with 100 cells at time t = 0.2s and ∆t = 10∆t explicit .

Fig. 4 . 8 .

 48 Fig. 4.8. Stokes number=0.0001: Error on the number density comparing to the asymptotic analytic solution using the Implicit non-AP scheme (squares) and the Implicit AP scheme (crosses) for ∆t = 10∆t explicit (full lines) and ∆t = 50∆t explicit (dot lines).

6 .

 6 Acknowledgments. The support of the France-Stanford Center for Interdisciplinary Studies through a collaborative project grant (PIs: P. Moin and M. Massot) is also gratefully acknowledged. The post-doctoral stay of A.Vié has also been supported by the ANR Sechelles (PIs S. Descombes and M. Massot) and DIGITEO MUSE Project (PI M. Massot).

  .[START_REF] Chalons | Beyond pressureless gas dynamics: quadrature-based velocity moment models[END_REF] About the Asymptotic-Preserving property. Let us first recall that the first equation (2.32) has been proved at the end of the second step and remains valid since ρ is not expected to be modified in this last step. Regarding now the second equation (2.33), we first write(3.19) under the following equivalent form

	St	ǫ n+1 j ∆t -ǫ n j	+ 2ǫ n+1 j

For example for 3D simulations where the phase space is at least 6D, 3D for the position and 3D for the velocity of the particles.

These groups were also referred as differential and algebraic model respectively in[START_REF] Alipchenkov | Differential and algebraic models for the second moments of particle velocity and temperature fluctuations in turbulent flows[END_REF]

In[START_REF] Pandya | Two-fluid large-eddy simulation approach for particle-laden turbulent flows[END_REF], the authors derive equations for zero-to-second order moments and close the fluxes by using a zero third order central moments assumption, but they do not rely on the link with the kinetic level, see[START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] for details.