Statistical and probabilistic modeling of a cloud of particles coupled with a turbulent fluid

This paper exposes a novel exploratory formalism, which end goal is the numerical simulation of the dynamics of a cloud of particles weakly or strongly coupled with a turbulent fluid. Giventhe large panel of expertise of the list of authors, the content of this paper scans a wide range of connexnotions, from the physics of turbulence to the rigorous definition of stochastic processes. Our approachis to develop reduced-order models for the dynamics of both carrying and carried phases which remainconsistant within this formalism, and to set up a numerical process to validate these models. Thenovelties of this paper lie in the gathering of a large panel of mathematical and physical definitionsand results within a common framework and an agreed vocabulary (sections 1 and 2), and in somepreliminary results and achievements within this context, section 3. While the first three sections havebeen simplified to the context of a gas field providing that the disperse phase only retrieves energythrough drag, the fourth section opens this study to the more complex situation when the dispersephase interacts with the continuous phase as well, in an energy conservative manner. This will allowus to expose the perspectives of the project and to conclude.


Introduction
Many applications involve the transport of a disperse phase (particles, droplets, bubbles) coupled with a fluid: spray combustion, fluidized beds, soot dynamics... In the standard case, the evolution of the carrier phase can be described by a deterministic system of equations such as the Navier-Stokes equations. However, in the strongly coupled case the evolution equations are unclosed due to the exchange term with the particles [10,12]. Often, models proposed in the literature only consider the influence of the carrier fluid on the disperse phase and neglect its retroactive consequences, or, at best, limit it to a global balance between the two phases [21]. In particular, these simplifying hypotheses allow to decouple the inaccuracies coming from the approximate resolution of each phase.
But, one of the main difficulties in the derivation of a consistent model for the strongly coupled evolution of a cloud of particles within a turbulent flow, is that inaccuracies arise both from the chaotic behavior of the fluid [9,24,35], and from the initial properties of the particles, such as their starting positions and velocities. Therefore, the proper level of "modeling" consists in making consistant assumptions about the properties of the stochastic processes involved in the global dynamics of both phases. Even if some advances have been made in the field [14], the problematic is far from being closed.
In order to better understand the coupling of the inaccuracies coming from both phases, we split the construction of the fluid dynamic model into four main steps, corresponding to four spatial levels of modeling.
Step-by-step, we then express some links between these levels, in order to better understand the influence of the small scales on the highest level of modeling. Here, one has to understand that this hierarchy of points of view is worth both for the carrier fluid and for the disperse phase. Simply, each of the passing to the limit between each levels does not occur at the same scale for the two phases. Although the carrier fluid is made of nanometric particles, while the dispersed particles seldom reach a micrometrical size, the description of each phase starts at the microscopic level (or molecular level). From there, one can reach reduced order large-scale models rather continuously, by first looking at an intermediate mesoscopic scale, dealing with the law on the presence of the microscopic phase (e.g. the Boltzmann equation), and then consider close-to-equilibrium regimes that we will call the macroscopic scale (e.g. Euler or Navier-Stokes equations). These four different levels of modeling are sketched level-by-level in the following items list: • Microscopic: at the scale of atomes, molecules or particules. Generally speaking, one may say "at the scale of the indivisible". The medium is here modeled by a very large number of ODEs. • Macroscopic: at the scale of the continuum. Fluids (liquid, gaz, spray,. . . ) are now seen as a continuous medium. It is modeled by a system of PDEs. • Mesoscopic: the transition from the micro to the macro scale necessitates an intermediate scale, called "mesoscopic", at which the medium is modeled a statistical manner. At this level, the fluid is modeled by the transport equation of a probability density function (PDF) of particles. • Reduced-Order: despite all the complexity reduction already performed, the simulation of all the macroscopic scales (Direct Numerical Simulation, DNS) is far from being reachable. An additional order reduction is then performed by splitting the solution into a significant part φ and a residual φ : φ = φ + φ . In general, the residual is removed and its action on the resolved part is modeled by a chosen underlying random process.
Throughout this paper the term significant part is kept general on purpose: it could denote one of the numerous choice of decomposition of the macroscopic sought solution into a numerically resolved and an unresolved part, see subsection 2.2.2 for more details. To give an insight of historical context, the usual method is traditionally referred to as a Large Eddy Simulation (LES) of the particulate flow, which means that only the features of the flow at a scale greater than a characteristic cut-off size are computed. The smallest scales, called subscales, need to be modeled from the computed variables in both carrying and disperse phases. As proposed by Pope [25], we chose to place ourselves in a probabilistic formalism where the closure in performed by the definition of a probabilistic process for the residuals. This closure can be seen as a probabilistic mapping between the reduction of the non-linear terms of the solved macroscopic PDEs and the resolved variables, see subsection 2.2.2. As a consequence, defining a subscale model is equivalent to making a choice for this mapping. This is what we are looking for in this project.
An ideal model for the numerical simulation of a turbulent flow loaded with dispersed particles would be a global reduced-order model for both phases, where the residual part would have to be able to take into account the strong coupling between both phases (mass, momentum and energy are exchanged in a bidirectional manner and globally conserved). We think that the formalism introduced in [25], and rapidly sketched in the previous paragraph, is a good starting point. We also believe that the stochastic model of the unresolved fluctuations has a root at the microscopic level in both phases. This is the reason why we then start our exploratory study by considering an idealistic micro/micro modeling with additional stochastic processes on both phases, and then try to derive a global large-scale reduced-order model for the dynamics of the strongly coupled system, which remains reliable, accurate and consistant with the underlying micro/micro description of the physical system. This paper is divided into four sections. In a first section, we give a statistical description at micro and mesoscale which are the beginning of all macroscopic descriptions, with a theorem in the infinite population limit. It explains the link between a system of a large number of ODEs at the microscopic level and a PDE on a Probability Density Function (PDF) of existence of the particles. In section two we describe in a very condensed manner the other levels of continuous description, while staying as consistant as possible. This leads us to a very general definition of turbulence and to the probabilistic framework for the modeling of the subscales in the context described by Pope [25]. In particular we explain the derivation of a reduced-order model for the disperse phase only, when the underlying carrying continuous gas field is supposed to be perfectly known and is not perturbed by the presence of the particles. Section three presents a numerical process intended to validate the reduced-order models possibly created within this micro/micro to reduced-order context, by looking at the statistics missed by the disperse field when the underlying gas velocity field has been reduced (for example filtered). In particular, we show that it seems hard to build a reduced-order turbulent model for the dynamics of a 1D spray, but that the situation improves with higher dimensionality. Finally, section four opens the discussion on the construction of a consistent model for two-way coupled systems. This section being preliminary, this will allow us to expose the perspectives of the current project and to conclude the paper.

Statistical description of the dynamics of a population: from micro-to meso-scale
This section describes the dynamics of a population at micro and mesoscale. This is the beginning of all work implying complex dynamics of turbulent particules-laden flows. This gathering represents a real team effort, especially in finding a common vocabulary between those of us more physics-oriented and those more used to the theory of probability and of stochastic processes. As already said in the previous paragraph, what is written here is valid for both carrying and disperse phases, only the passing to the limit do not occur at the same scales.

Microscopic scale
The studied domain X ⊂ R 3 is filled with a cloud of N identical spherical particles, moving into void or supported by a carrying gas. Assuming that the three degrees of freedom in rotation of each particle can be ignored, the dynamics of the system is described by the 6N parameters (velocity are in R 3 C C C := R 3 ): or equally by the empirical measure or normalized counting measure: . If the set of particles is immersed within an external field G(t, X X X, C C C), interacts with itself following a collision kernel F (X X X, C C C) and each particle is possibly subject to an independent Brownian random process of intensity σ, the phase space (1) evolves with the following system of 6N ODEs: Then, given an initial condition Z Z Z 0 = X X X 0 1 , C C C 0 1 , . . . , X X X 0 N , C C C 0 N , which may be deterministic or stochastic, the empirical measure can be indexed by is the number of particles from the configuration Z Z Z 0 at time t = 0, situated within V X X X and with a velocity belonging to V C C C at time t.

Mesoscopic scale
From now on, the configuration of each particle is denoted by z z z i = (X X X i , C C C i ) ∈ X×R 3 C C C , for all i = 1, . . . , N . The collision kernel F simulates the interaction between the particles and it thus seems fair to have F (−z z z) = F (z z z), which implies F (0) = 0.
Let us consider that the particles are changeable at initial time, which means that their initial distribution µ N 0 ∈ R 2d is invariant by permutation of the N variables. This invariance therefore remains satisfied at any time t > 0 and in particular, the N particles must follow the same one-particle law in R 2d , denoted µ (1) t , which is what we are looking for in this subsection. First, if A is a borelian in R 2d , Then, we recall Z Z Z(t) = (z z z 1 (t), ..., z z z N (t)), later simply noted Z Z Z t , and we introduce the following function and the 2d × N diagonal matrix, denoted Σ, with zero (d times) and σ (d times), repeated N times along the diagonal. Thus, equation (2) can be rewritten For any function Φ : the Itô's formula gives us: denote projection operators on the respective lines of x x x i and c c c i . Taking the expectancy we get: Next, we introduce the following linear form on the measures of R 2d , defined for any Φ such as in (3): Here µ (N ) t is the N -joint law followed by the N particules: it is the law followed by Z Z Z t . Using this dual formulation, we can now extend the definition of the partial derivatives to the measures of R 2d , and we have: Since Φ does not have a compact support in time, integration by part requires to keep the boundary terms and the time partial derivative of µ N t defines as: To sum up, thanks to the Itô's formula, we have obtained a weak form of the equation followed by the law µ So now, we have generalized the results given by Bolley in [5] to a time dependent transport term G. However, equation (4) is a weak formulation of a PDE on the N -particles joint law, when what we are looking for is the equation ruling the one-particle law µ (1) t , which is the marginal of µ (N ) t for particle 1. By integrating Eq. (4) over all the particles but the first one, we get that, in the weak sense, µ (1) t follows: In this expression, µ (2) t is the 2-particles joint probability. In order to close equation (5), we would like to express it as a function of µ (1) t .
To do so, we suppose that the initial data Z Z Z 0 are indistinguishable and all follow the same law f 0 on R 2d . Then, we introduce an intermediate law, as the solution for t > 0 and (x x x, c c c) ∈ R 2d , of the following equation with initial data f 0 : where F and G are now supposed to be Lipschitz functions with respect to the variable x x x ∈ R 2d and G is continuous in the time variable. Next, for i ∈ 1, . . . , N, letz z z i (t) be the solution of the following system with initial dataz z z i (0) = z z z i (0): Then, the fictive particlesz z z i evolve in the field F * f t generated by the distribution f t , while the z z z i particles evolve in the F * µ N t field, generated by the empirical measure µ N t . Itô's formula gives once more the PDE followed byz z z i in the weak sense, and we now wish to show that this measure converges to f t when the number N of particles tends to infinity.
We denote |(x x x, c c c)| = |x x x| 2 + |c c c| 2 and for p > 1 we define The Wasserstein distance of order p between two measures µ andμ of P p (R 2d ) is defined by where Z Z Z andZ Z Z are stochastic variables of law µ andμ respectively. Then, following the lines of [5], Theorem 1.1. we obtain the explicit convergence rates:

3) Let Φ be a Lipschitz function in the second variable, then
In other words, this means that: 1) The one-particle law µ converges to f t in the Wasserstein distance when N → ∞, 2) At the limit of an infinite number of particles, the chaos propagates in time; the particles remain uncorrelated during the whole dynamics: µ . In particular, one recovers the famous molecular chaos assumption of Boltzmann: 3) The weak convergence of the empirical measure µ N t to f t . Finally, equation (5) is now closed rigorously thanks to the molecular chaos propagation in the context of Lipschitz-regular interactions (external G or between particles F ), [5,33]. However, when the interactions are less regular, which is the case for the Boltzmann equation (9) below, an increasing number of positive results let us think that equation (8) remains correct, [20,34]. Nonetheless, no rigorous demonstration is nowadays available.

A population of particles in a turbulent fluid
In the previous section, a general kinetic equation has been derived for a population of "particles" (molecules, droplets, solid particles). As this point, one can be interested in deriving a two-way coupled system of kinetic equations for the carrying fluid and the particles. However, in [7], it has been shown in the context of nanoparticles that such a derivation cannot be performed. Instead, we use the classical strategy of first deriving macroscopic equations for the fluid, and then coupling them to the particle equations, either microscopic or mesoscopic. In the following, we first present the Euler and Navier-Stokes equations that can describe a carrying fluid, with an emphasis on the underlying assumptions at the kinetic level. In a context where dealing with the whole range of scales of the fluid is not accessible, we detail a general strategy for generating large-scale reduced-order models, and we show how it can be taken into account for the description of the particle dynamics at the microscopic level.

Classical theories for macroscopic equations for the fluid
In the context of gaz dynamics, in the limit of an infinite number of particles and when ignoring the stochastic subscale Brownian perturbations for the moment, equation (5) becomes the Boltzmann equation: where Kn = λ L is the Knudsen number, ratio between the mean free path λ and a characteristic size of observation L, and where the quadratic collision operator Q writes: Remark 2.1. When considering a non self-interacting population of particles, its repartition function also follows an equation of the (9) type, where the Knudsen number is infinite: Kn = +∞.

Euler equations
For any PDF f , one can define its microscopic entropy by h = f log f . It can be understood as a local uncertainty rate. Then, the macroscopic entropy reads:

c c c) dc c c, and one can
show that when f is a solution of the Boltzmann equation (9), its macroscopic entropy decreases: dH dt ≤ 0. When the minimum H min is reached, log f must be a collision invariant and this implies that the velocity distribution f is a Maxwellian: The Maxwellian distribution being perfectly defined by its three first moments , the evolution of the Boltzmann equation (9) at isentropic thermodynamic equilibrium H = H min is given by the system of its three first moments, which closes into the Euler equations:

Navier-Stokes-Fourier equations
In the previous paragraph, we have clearly stated that collisions occur everywhere at all time, or, to reformulate, that the Knudsen number remains null: Kn = 0. In reality, it is often very small but strictly positive. Then, we look at near equilibrium regimes by stating ε = Kn and looking for an expansion of f in ε: the Chapman-Enskog expansion. At first order, f = f 0 + εf 1 + •(ε), which, at orders 1/ε and 1, gives in (9): The latest equation is an integral equation in f 1 which might be completely solved. For a monoatomic gas of atoms of mass m and radius r, the three first moments of f verify the following Navier-Stokes equations [6]: and µ and λ are respectively the viscosity and the thermal conductivity. . However, a similar system of PDEs can be obtained by considering the conservative principles of mass, momentum and total energy, added with constitutive equations of the considered fluid, which provide heuristic laws of the viscosity µ and the thermal conductivity λ.

Properties of turbulence
Turbulence is a particular type of flows which can not be rigorously defined. The easiest way to define it is by using the metric of the Reynolds number: Re = , where u f,0 is a characteristic speed of the fluid, L 0 is a characteristic length scale of the system and ν = µ ρ is the kinematic viscosity of the fluid. We will say that a fluid exhibits a turbulent behavior, when its Reynolds number is high. The limit Reynolds number depends on the considered experiment and on the operating condition. However, the flow is generally turbulent when Re >> 10 3 .
Turbulent flows share in common their chaotic behavior. For deterministic systems, there are multiple definitions of chaos, but in this context we choose to say that turbulent flows all are : • highly sensitive to the initial conditions of the system. The present determines the future, but the approximate present does not approximately describe the future. For instance, we say that x 0 is a highly sensitive initial conditions, if for all L 0 > M > 0 and for all δ > 0, there exists another close initial data y 0 and an arbitrary time t > 0 such that • topologically transitive, in the sense that for every pair of non-empty open sets U ⊂ X and V ⊂ X, there is an arbitrary time t > 0 such that From an experimental point of view, some observations have been made on turbulent fluid flows. The main ones are expressed by Kolmogorov [24, p.190].
• At sufficiently high Reynolds number, the small-scale turbulent motions are statistically isotropic. They follow a universal form that is uniquely determined by the viscosity ν and the energy dissipation ε. • The viscosity also defines a cut-off size η K , called the Kolmogorov scale, below which all the inertia of the flow is dissipated. • Between the characteristic length L 0 and η K , there is an intermediate range of scales, called the inertial range, where the statistics of motion have a universal form that is uniquely determined by the dissipation ε and is independent of the viscosity ν. Through dimensional analysis, we get that within this range, the turbulent kinetic energy decreases as: E (|k|) ∝ |k| − 5 3 , with k the wavenumber.

Reduced description of turbulence
It is commonly admitted that the macroscopic Navier-Stokes equations contain the turbulence defined above, in the sense that these equations present solutions which have all the properties listed in paragraph 2.2.1. Nonetheless, in practice the domain size, denoted by |X|, and the dissipative cut-off scale η K , may be separated by many orders of magnitude. In this context, the Direct Numerical Simulation of the Navier-Stokes equations is rapidly unreachable, since the number of needed computational cells will be at least of the order of (η K /|X|) 3 , not speaking about the generally necessary high number of degrees of freedom per cell.
Therefore, while staying very generic, we consider a decomposition of the solution into a significant part and a residual: if φ is a quantity of interest, we consider its reduction φ on the space of significant data and thus write φ = φ + φ . This significant part could be an ensemble average, a filtering, a spatial or a temporal average or even a modal decomposition. The goal is always to reduce the size of the information needed to entirely represent the chosen significant part, hence the name reduced-order model. Now, the reduction operator · is applied directly on the macroscopic equations Eq. (12)- (13). For example, when considering the incompressible version of the Navier-Stokes equation, assuming commutativity between all implied linear operators, one gets: with u f the fluid velocity, ρ f its density (constant for incompressible fluids), ν its kinematic viscosity and p the pressure field.
The main difficulty now lies in the reduction of the non-linear terms. Indeed, nothing indicates that there exists an application giving (u f · ∇) u f as a function of u f . Thus, Eq. (14) is not meaningful in term of the significant unknown u f . To overcome this difficulty, the main idea is to define a more complex application which gives multiple possibilities to the relation between (u f · ∇) u f and u f . This is done by adding a hidden variable ω, which encodes all the complexity of (u f · ∇) u f inside an application F and a space of possibilities Ω in the following way Of course, the definition of F strongly depends on the choice of the reduction operator · . Next, an elegant way to move forward is now to define Ω as a probability space, see [25]. Then, two main techniques emerge : • by drawing many particular ω, thus giving a random modeling of the unknown term (u f · ∇) u f through F, compute many trajectories of the process u f , • considering the statistics or moments of the random variable F, and solve for the evolution of the moments of the random variable u f .
The advantage of the first approach is to preserve the properties of a trajectory of the process u f , which is still the solution of a PDE. Thereby, the random variable u f lies in a large dimensional probability space, which requires a very large number of such succession of draws to hope for some meaningful statistics. On the contrary, solving for the evolution of the means of u f does not preserve the trajectories of the process, but it gives correct estimators and statistics on the general behavior of the gaseous velocity field.

Closures
The obtained reduced-order system as in Eq. (14) is closed by making a calculable choice on F. Three strategies can be found in the literature for this choice, as depicted in [24,29]: • the functional approach: starting from the fact that the regularized version of the flow field will dissipate less energy than the real turbulent flow field does, the unresolved scales can be modeled in a first approximation by an additional diffusion process, consistently with the theory of turbulence described in paragraph 2.2.1: Here, µ turb is an additional turbulent viscosity. In the case of filtering procedures, this viscosity depends on the filter size such that it vanishes for full-resolution [23,32]. As such models can depend on empirical constants, dynamic procedures were also proposed to get the better estimate of theses constants (see [13]). • The structural approach: instead of simply recovering a global property of the unresolved information, structural methods aim at capturing the SGS tensor structure (see [2]). • the "pragmatic" approach: starting from the idea that it is hard to distinguish unresolved scales effects from numerical dissipation, some authors propose to integrate effects of unresolved scales through the numerical schemes (see [15]).

Reduced LES models
The fluid velocity at the location of the particle appears in the expression of the particle acceleration modeled by Stokes drag law: i ∈ [[1, N ]], and τ p being a characteristic relaxation time of the particule toward the underlying velocity field. However, in every LES model existing up to now, only a regularized version of the fluid velocity is computed. Thus, a closure on the fluid velocity seen by the particle is required in order to provide a consistant LES model for the disperse phase. Ideally, this model has to be in agreement with the probability space of the random variable F seen by the inertial particles on the fluid flow.
Up to now, very similarly to the models developed for the fluid flow, the main strategies have been to compensate second order moments of the the particle density distributions by the adjunction of energy in the form of Wiener processes (see [4,11,22,27,30,31]). In its general from, this can be represented by the stochastic differential equation (17): with W W W t a Wiener process, Z Z Z t the state vector of the particle, µ t the drift and σ t the diffusion coefficient. It is to be noted that in most models, the Wiener process only acts on one variable of the particle : either its position, or its velocity, or an other intermediate variable like the velocity seen by the particle. The next section shows that in the context of equation (17), where the closure has been chosen in the form of a Wiener process, the derivation of a mesoscopic equation for the disperse phase is not a major difficulty.

Consistency of modeling approaches with numerical cases
Sections 1 and 2 were mainly focused on providing a meaningful formalism for reduced multiphase flow simulations in agreement with mathematical consistency and physical literature. In this context, we conclude that an appropriate formalism to describe a fluid in a Large-scale reduced order in section 2.2 is the selfconditioned structure proposed by [25] and formalized Eq. (15). In a nutshell, the evolution of the large scale of the flow must be obtained as the expectation of all possible unresolved scales of the flow compatible with the resolved large scales.
Applying this formalism with the full resolution of Navier-Stokes is not easy because it is not straightforward to control large scales and unresolved scales separately. An interesting alternative that has been widely used in the literature is to rely on synthetic turbulence: by means of a summation of analytic modes, and under the constraint of specific spectral distribution and representation, one can expect to reproduce the main characteristics of the turbulence, even without verifying Navier-Stokes equations. In this section, we investigate the use of such analytic representation from 1D to 3D, and we show what is the minimal representation that can be envisaged.

Synthetic turbulence
The synthetic flow field has been designed in order to reproduce somehow the dynamics that could be expected from a self-conditioned LES flow field simulation ( [17,18]). It is represented by a sparse matrix of spectral modes (Eq. (18)) chosen according to the energy density given by Pope's spectrum in Eq. (19) (see [24, p.232]) with Eq. (20). u f (t, z) = N n=0 a n cos (ω n t + k n · z + ϕ n ) The amplitude of the modes is chosen according to the distribution |a n | ∼ N 0, Following [17], the spectral components of the energy spectrum are chosen in order to respect the numerical simulations performed in [16], which show that it seems sensible to approximate E (|k| , ω) by : with a ∈ [0.4, 0.51] depending on the wavenumber and the integral length scale (see [16]). For the numerical simulations, the random number generator chosen is ran2 presented in [26]. The numerical values are chosen such that a = 0.5, u 0 = 1 m.s −1 and k 0 = 1 m −1 . The particle evolution is computed using Runge-Kutta scheme of order four. The evolution of the particles on the fluid is computed by the linearised Stokes drag law in Eq. (16), with the expression of u f given in Eq. (18).
For numerical simplicity, we first start by performing one-dimensional simulations. In one dimension, a realization of the evolution of the particles submitted to a random fluid is given in Fig. 1. Although the initial positons of the particles are random and uniformly distributed on a segment, their trajectories seem very limited. They look more like oscillations around a mean drift rather than dispersion. Furthermore, when observing the evolution of the variance in a one-dimensional space for 10 4 particles, see Fig. 2(a), we see that it seems bounded for this case and that it is highly dependent on the underlying fluid fluctuations.
This kind of behavior is not consistent with the properties of turbulence and the expected behavior of particles in a turbulent flow: we would rather expect a dispersion behavior similar to diffusion (see for instance [28]). Since the stochastic models of the literature have a first order effect on the second order moments of the measure of the disperse phase, it is essential to work on a numerical setup which preserves the basic properties of turbulent flows for realizations of the second order moments of the measure of the disperse phase. Hence, it is of prime importance to understand why such a behavior is observed on the simple fluid model we have chosen if we want to use it for reproducing and understanding the dynamic of inertial particles on fluids described by Navier-Stokes kind of equations.

Simplified one-dimensional case
As explained above, Fig. 1 enlightens an unexpected behavior in one dimensional case. Let us start by looking if it is possible to understand this behavior on a simplified case where the fluid is only represented by one sine. We have the particle evolution in Eq. (22) and the reduced evolution in Eq. (23).
We will prove the following result : Proof. In order to prove this result, we can first study the system (23). This system is autonomous in dimension 2, so by the Poincaré-Bendixon theorem, only three cases are possible : • The trajectories are unbounded, • The trajectories converge to a point, • The trajectories converge to a limit cycle. Let us now try to characterize these behaviors more precisely.
It proves that in finite time, the solution falls under C max . Then we have proved that for all trajectories, there exists a time t max where C i (t max ) < C max . Denote C min = −C max and use again the time t * with symmetric definition, we obtain and δ ≤ (C min − C i (t * ))τ p . It proves that in finite time, the solution rises above v min . Then we have proved that for all trajectories, there exists a time t min where C i (t min ) > C min . Finally we can suppose that for all trajectories, the speed C i ∈ [−C max , C max ] after some transitory time. In fact -with the same procedure-we can prove that C i ∈ [ω − |a |, ω + |a |].
Thus if ω > |a | then it proves that the speed stays strictly greater than ω − |a | > ε > 0, and thus the trajectories cannot be bounded. In order to prove that the particles will follow an increasing signal, we have to study the difference with this linear growing. Denote We can see that V i cannot converge to a constantV , because there is no solution toV = a sin(2πt(V + ω)) (except ω =V = 0). Since V i cannot converge to a constant while staying in a compact, it is non-monotonous. Denote T + a moment where dV i dt changes its sign (without loss of generality, suppose it changes from > 0 to < 0), i.e. a sin(2π(Y i (T + ) + ωT + )) = V i (T + ) =: V + . Thus and in particular Suppose a > 0 to simplify. The quantity dV i dt changes from > 0 to < 0, thus the second derivative is negative, so at a given time, there is a local maximum, and during a period [T + , T + + T ], V i is decreasing and we have also 2π(Y i (T + ) + ωT + ) ∈ [π/2, 3π/2] mod 2π. Or simply 2π(Y i (T + ) + ωT + ) = 2kπ + π/2 + επ with ε ∈ [0, 1].
If V + > 0 we have ε ∈ [0, 1/2]. And since V i is decreasing, and Y i increasing unbounded, there is a moment where 2πY i (t) + 2πωt = 2kπ + π = 2πY i (T + ) + 2πωT + + π/2 − επ. At this moment, V i becomes negative and Y i becomes decreasing. Since V i is bounded, it will reach a minimum (since it cannot converges). Denote this time T − and we are in the symmetric case than previously.
We have proved that there exists two sequences (T n + ) n∈N and (T n − ) n∈N such that T n + < T n − < T n+1 + for all n ∈ N. We can bounded the time (T n − − T n + ) above and below independently of n ∈ N roughly proving that the solution is close to a periodic one. Finally the solution X i is close to a increasing signal having periodic oscillation around its drift, which is incompatible with an expected diffusive behavior.
In this particular case of only one sine, we have performed a transformation which leads to an autonomous system, and hard conclusion with only a discrete set of final positions. With more exciting sines the behavior could be different. But -as it is represented in Fig. 2-even with more exciting sines we do not obtain in 1D a dispersive behavior as expected. It makes a 1D model very dubious.
But, dispersion of particles is greatly influenced by the dimensionality of the underlying space chosen. Although the dynamic in the one dimensional case is very different from the physic we aim at modeling, we expect that when dimensionality is increased, this behavior will change and be most likely similar to diffusion (see Fig.  2), as envisioned by the physic, and as described by the models currently in use in the literature. Let us check this assumption in the following section.

Higher dimensionality
It is possible to observe numerically that by increasing the dimensionality to more than one physical dimension (Figs. 2(b) and 3(a)), the second order moment of f t has a better behavior, i.e. it increases quite monotonously with time, and the particles do not seem to be overly constrained by the underlying fluid flow. The higher the dimensionality, the better the dispersion of the particles. Indeed, one observes in Fig. 2(b) that the dispersion of the particles appears to be much less influenced by the characteristics of the underlying fluid flow than in the 1D case (see Fig. 2(a)), and that the third dimensionality brings even more smoothness (see Fig. 3(a)). The change of behavior between 2D and 3D can also be partly understood by the addition of new topologies for the three-dimensional stationary points as described in [3].
Given these results, it seems relevant to keep on pursuing the simulation effort focusing on the three dimensional configuration.

Towards two-way coupled systems
The next step towards the modeling of particulate flows is to account for the impact of the disperse phase on the turbulent carrier phase, which has strong implications. Let us consider the empirical measure µ N t (Z Z Z) = 1 N N i=1 δ X X Xi(t) δ C C Ci(t) and the following evolution equation In a one-way coupled context the gas phase velocity at the particle location does only depend on the particle position itself and is independent of the others particles as they share the same gas phase. In this context, we satisfy the conditions of Theorem 1.1, i.e. G (t, X X X i , C C C i ) = u f (t,X X Xi)−Ci τp . We can thus state a theorem of convergence towards the law of the process.
In a two-way coupled system, all particles affect the gas phase evolution such that the gas velocity is condi- where κ = mpn l ρ f . The equilibrium solution is then: As a consequence, if we want to study the impact of inhomogeneity of the particulate phase by changing the number of particles but keeping the same physical problem, we need to modify the particle mass m p accordingly, to keep κ m = κdx/L x = N p m p constant.

Particle-laden case with Lagrangian particles
Knowing the sought continuum limit of the particle system, we now investigate the impact of the number of particles, i.e. the impact of the statistical convergence of the randomly-drawn initial condition. We thus simulate the two-way coupled burgers problem by changing the number of particles from 1 to a large number or particles. In Fig. 3a, we compare the time evolution of the gas velocity averaged over a large number of realizations of the initial conditions for different numbers of particles at fixed mass loading. We clearly see the convergence of the Lagrangian simulations towards the homogeneous solution, with a convergence rate of order 1 (see Fig. 3b). This convergence rate is not affected by the number of cells for numerical discretization and by the addition of physical diffusion in the Burgers equation. So even if we do not have a formal proof in the spirit of Theorem 1.1, we still have confidence in the existence of a convergence result, and thus of an Eulerian limit description. Figure 5. 1D burgers problem. Ensemble-averaged gas velocity for different number of particles from 1 to 256, with the Lagrangian tracking (dashed lines) and the "empirical" Eulerian moment method (full lines to compare with the homogeneous limit (black). systems of interest. Finally, we have shown the main limitations in two-way coupled system, proposing some possible solutions to overcome them.