An Overview of the Lagrangian Dispersion Modeling of Heavy Particles in Homogeneous Isotropic Turbulence and Considerations on Related LES Simulations

: Particle tracking is a competitive technique widely used in two-phase ﬂows and best suited to simulate the dispersion of heavy particles in the atmosphere. Most Lagrangian models in the statistical approach to turbulence are based either on the eddy interaction model (EIM) and the Monte-Carlo method or on random walk models (RWMs) making use of Markov chains and a Langevin equation. In the present work, both discontinuous and continuous random walk techniques are used to model the dispersion of heavy spherical particles in homogeneous isotropic stationary turbulence (HIST). Their efﬁciency to predict particle long time dispersion, mean-square velocity and Lagrangian integral time scales are discussed. Computation results with zero and no-zero mean drift velocity are reported; they are intended to quantify the inertia, gravity, crossing-trajectory and continuity effects controlling the dispersion. The calculations concern dense monodisperse spheres in air, the particle Stokes number ranging from 0.007 to 4. Due to the weaknesses of such models, a more sophisticated matrix method will also be explored, able to simulate the true ﬂuid turbulence experienced by the particle for long time dispersion studies. Computer evolution and performance since allowed to develop, instead of Reynold-Averaged Navier-Stokes (RANS)-based studies, large eddy simulation (LES) and direct numerical simulation (DNS) of turbulence coupled to Generalized Langevin Models. A short review on the progress of the Lagrangian simulations based on large eddy simulation (LES) will therefore be provided too, highlighting preferential concentration. The theoretical framework for the ﬂuid time correlation functions along the heavy particle path is that suggested by Wang and Stock.


Introduction
In order to simulate multiphase flows, both Eulerian and Lagrangian models have been developed and rendered more efficient over the last four decades [1][2][3][4].Lagrangian particle tracking methods are particularly developed for turbulent dispersion problems, at first for dilute multiphase flows [5,6] and more dense gas-solid flows [5,6] in the industrial and engineering field or to simulate atmospheric pollution of tracers or heavy particles [7].However, the dynamic response of a heavy particle to a turbulent flow is different from that of a tracer behaving like a fluid element because of inertia and gravity effects.Insofar as the Lagrangian approach pioneered by Taylor [8] for turbulent diffusion is no longer valid, it must be reformulated.A plethora of models have been proposed to simulate the turbulent dispersion of heavy particles.Most are based on a Monte-Carlo "process."In the present paper, Monte-Carlo process is used as an abbreviation for the well-known Monte-Carlo computational method, using repeated random sampling to numerically analyze the turbulent dispersion process.The average ensemble over many trajectories of discrete particles is representative of their dispersion in a turbulent carrier flow.The trajectories are obtained by solving numerically the equation of motion, taking into account the fluid velocity along the trajectories of the discrete particles.How the fluid velocity "seen" or "experienced" along the trajectory can be simulated is the crucial question in all Lagrangian models.Although experiments can provide Eulerian-at a fixed pointinformation, dispersion needs information in a Lagrangian frame, and the relationships between Eulerian and Lagrangian statistics or an accurate value of the Kolmogorov constant are not yet available [9,10].
The Lagrangian models, well-adapted to the study of dilute gas-solid flows, are interesting for modeling the dispersion of the particles emitted from a point source even in complex geometry configurations.These models are the most appropriate ones for studying particle dispersion in turbulent flows insofar as they allow to take into account relatively easily phenomena such as interactions between particles, wall-particles collisions or the physical laws of each particle (combustion, evaporation, rotation, etc.).Lagrangian stochastic (LS) models used for predicting the turbulent dispersion of the particles can be divided into two important classes: the eddy interaction model (EIM, [11]) and the random flight or walk model (RWM, [7,[12][13][14]).The main difference between these two kinds of Lagrangian models is the method they use to statistically generate the turbulent fluid velocity in the particle surrounding, since this fluid velocity along the particle path is necessary to solve the equation of motion and calculate the trajectory.Other classifications are indeed possible.However, they are all influenced by the way the fluid velocity is simulated; models will have different abbreviations or terminology, whether or not they have been developed early for turbulent diffusion, homogeneous or shear turbulence or heavy particles, for EIM, RWM, NLWs (normalized Langevin models), GLW (generalized Langevin equations), SDEs (stochastic differential equations), DRW (discrete/discontinuous random walk) or CRW (continuous random walk) models.
A Monte-Carlo (MC) process will be associated to the EIM-DRW class; the fluid Lagrangian velocity is generated in this case through independent random numbers, and the fluid velocity is kept constant for a given time of the order of the eddy lifetime.The CRW models are related to Markov chains as well as a Langevin equation, and make use of random numbers with a one-step memory to simulate the fluid velocity "continuously".Multistep Markov chains have also been developed over the last 30 years and can be classified as matrix methods.They are more efficient but need more computer resources.
Numerous studies in turbulent diffusion and dispersion use models based on Monte-Carlo methods or Markov chain methods, and the validation of these models is generally performed on a few available experimental test flow cases in grid turbulence [15,16] or the center region of a pipe flow [17].These flows are well adapted to simulate an ideal turbulent field, homogeneous isotropic stationary turbulence (HIST).The main objectives of the present work are to present an overview of their efficiency to simulate the dispersion of heavy particles in HIST, since Lagrangian modeling is a tool also developed in more complex flow configurations.The reference will be the theoretical contribution of Wang and Stock [18][19][20].It will be shown that most models are far from satisfactory, even in such an ideal theoretical turbulent field as in the case of dilute gas-solid particle-laden flows.
The article is organized as follows.Section 2 gives a summary of the theory on turbulent diffusion, a recall of the turbulent scales, equation of motion and parameters governing turbulent dispersion of heavy particles.Section 3 is dedicated to eddy interaction and classical random walk models.Numerical results of modeling turbulent dispersion in HIST with these models are developed and discussed in Section 4. Finally, Section 5 is dedicated to a short overview of the progress on random walk models based on the generalized Langevin equation and LES modeling and introduce preferential concentration.

Theoretical Background
Turbulent diffusion, also called passive dispersion, refers to the spreading in turbulent flows of scalar quantities such as heat, light seeding particles, chemical species, vectors, e.g., momentum, or tensors, e.g., second-order velocity correlation functions.On the other side, turbulent dispersion, or particle dispersion, consists of the spreading, due to turbulence, of discrete solid particles or droplets exhibiting inertia and a mean relative fluid-particle velocity due to gravity.In a turbulent flow, the particle motion is influenced by two mechanisms that make the particle dispersion more complex than the fluid diffusion, one relevant to the turbulent carrier flow, another related to the particle properties.Sometimes, turbulent dispersion includes diffusion of particles having such small dimensions that inertia vanishes.In any case, effects of molecular diffusion are masked by the turbulent effects, a realistic assumption, except in the near-wall regions, where molecular mixing becomes important.

Taylor's Turbulent Diffusion Theory and Batchelor's Generalization
The seminal theory of turbulent diffusion was first formulated by Taylor [8], in analogy to Brownian motion transposed to continuous movements of fluid particles or passive markers.Taylor presented a first attempt at a theory describing the turbulent diffusivity.In a turbulent field homogeneous in space and stationary in time, he studied the Lagrangian turbulent diffusion of fluid particles from a purely kinematic point of view, and proposed limiting cases for small and large diffusion times.He considered the continuous random motion of a tagged fluid packet (eddy).The fundamental result was a relation between the mean lateral square fluid point displacement y 2 2 (t), the mean-square fluctuating eddy velocity u 2  2 in the direction transverse to the mean flow and the Lagrangian velocity correlation of the fluid R f L22 (τ): the overbar holds for time statistics, e.g., the time average equal to ensemble average under the ergodic assumption in stationary turbulence.For small times, the correlation R f L22 (τ) is close to 1, from which it follows that y 2 2 (t) = u 2 2 t 2 .For long times, the velocity fluctuations become uncorrelated and a Lagrangian fluid integral time scale T f L22 can be defined as and mathematically one can show that By analogy to Einstein's formulation of diffusion by Brownian motion, turbulent diffusivity for the fluid particles is then given by For long time dispersion, a turbulent diffusivity can then be introduced by as well as a Lagrangian integral lengthscale Λ L22 = u 2 2 T L22 , a space length on which a marker substantially follows a single direction (direction 2 here).
Well-known papers by Kampe de Feriet [21,22] and Frenkiel [23,24] put the fundamental work of Taylor on turbulent diffusion and statistical theory of turbulence in highlight.
Batchelor [25,26] extended Taylor's results in a more rigorous study on three-dimensional stationary homogeneous non-isotropic turbulence, the diffusion of fluid particles being statistically described by tensors.The stationary Lagrangian velocity autocorrelation of a single fluid particle is given by where u k and u kRMS = u 2 k are the fluid velocity fluctuation and the root mean square (rms) of the velocity fluctuations in direction k.The mean squared dispersion can be written as a tensor by and the Lagrangian time integral scales are Fluid turbulent diffusivity can now be defined in terms of a generalized diffusion tensor: Again for short times, dispersion coefficients increase linearly with time.For long times, the integral of the Lagrangian correlation coefficients are constant and equal to an Lagrangian integral time scale T f Lij .Therefore, for long times, the particle diffusivity tensor can be written as In terms of turbulent energy spectra, the expression of the mean-square displacement is given by which shows that turbulent diffusion increasingly depends on the lower frequency components of the velocity power spectral density distribution as the time of diffusion increases.

Toward Turbulent Dispersion
Taylor's and Batchelor's formalism can be extended to turbulent dispersion of discrete solid particles in two-phase flows.There is no contradiction with their theory, provided that the quantities are interpreted as particulate characteristics, i.e., the displacement and velocity of fluid particles are replaced by those of the discrete solid particle.Turbulent dispersion is controlled by both the properties of the particles and those of the fluid flow-mean and turbulent flow-carrying the dispersed phase.
Particles can be characterized by their density ρ p and their size or diameter d p .For the carrier fluid phase, apart from its density ρ f and dynamic viscosity µ f , one has to consider both the convection by the mean velocity U and the time and length scales of the turbulent flow.Furthermore, both Eulerian and Lagrangian scales can govern dispersion.While diffusion and dispersion of very small particles are controlled by Lagrangian scales, Eulerian quantities rather dominate the dispersion of larger particles.

Scales of Turbulent Motion
The Eulerian and Lagrangian scales characterizing fluid turbulence are identified in several textbooks [27][28][29][30] or fundamental papers [31].Fluid turbulence can be described in terms of Lagrangian fluid particle properties, rms fluid velocities u , sometimes noted as u 0 in the paper, time correlation functions R f Lij (τ), integral (or macroscales) times scales T f Lij or lengths scales Λ L = u 0 T L , and a microscale τ L , also called Taylor scale.Eulerian scales can be more easily obtained by experiments, the measuring probes are fixed in space, and time or space correlations can be calculated.In a fixed frame, the Eulerian integral time scale, T E , as well as the microscale τ E are obtained from the Eulerian correlation coefficient at a fixed point in the fluid with a time lag τ: and T E is a good measure of the longest correlation time for the turbulent velocity at a fixed point.To obtain the time microscale, one can expand the Eulerian correlation coefficient around zero: This parabola intersects the x-axis at τ = τ E .The microscale τ E is a measure of the rapid time changes in u(t).The Eulerian integral scales characterize the energy-containing eddy length scales, they are obtained in terms of normalized two-point velocity (space) correlations.
L Eij is a good measure of the longest space correlation for two fluctuating velocity components at two points a distance r → e i apart.Of interest are two space correlations in homogeneous turbulence, In isotropic stationary turbulence, dr and if f (r) = exp(−r/L E11 ), one easily shows that L E22 = L E11 /2 and that g(r) has negative loops.
Eulerian length microscales λ f , λ g can be defined, too, as for the microscale τ E .Particularly in low level turbulent flows, (u' << U), the frozen pattern hypothesis of Taylor can be applied: Eulerian correlations can also be measured in a coordinate frame moving at the average velocity of the flow.In that case, and the integral Eulerian space scale will be represented by L E11 = (u 2 ) 1/2 T mE = u 0 T mE = UT E in most cases.Further, turbulence can be described in terms of eddy size and eddy lifetime [30].The mean flow can be scaled by U, a characteristic of the local mean veloc- ity, the turbulence by turbulent velocity u' and a turbulent kinetic energy k ∝ 3u 2 0 /2 in homogeneous isotropic turbulence.The size of the largest eddies is related to the physical boundaries of the flow and to an integral length scale L f , and it is generally assumed that L f = L E = u 0 T mE = UT E .The size of the largest eddies, L f can be defined by the kinetic energy of the flow, k ∝ 3u 2 0 /2, and the energy dissipation rate ε, L f = k 3/2 /ε, and an integral time scale T 0 is commonly referred to as the large eddy turnover time-T 0 = L f /k 1/2 ∝ k/ε.The large energy-containing eddies lose their kinetic energy within the turnover time and the energy transfer rate is ε ∝ k/T 0 ∝ k 3/2 /L f .The smallest eddies can be characterized by the smallest dissipation scales, according to Kolmogorov, and depend on the dissipation rate and viscosity: where ν is the kinematic viscosity of the fluid, and ε is the turbulent kinetic energy dissipation rate.The Eulerian microscale (or Taylor scale) λ marks the transition from the inertial subrange to the dissipation range [27,31], the length scale for which viscous dissipation begins to affect the eddies: Note: More generally, space-time correlations certainly play a role in dispersion, too.Literature data on relationships between Eulerian and Lagrangian scales of turbulence have been developed in some review papers [9,10], but are still contradictory, generally Λ L ≈ L E = L f , T L ≈ T mE , but turbulence structure parameters can be introduced, for instance, m = u 0 T mE /L f , the ratio of the moving scale T mE to the eddy turnover time Although it is often accepted that for simple turbulent flows (grid turbulence, channel flows) Λ L ≈ L E = L f and T L ≈ T mE ,m ≈ 1, significant differences can be observed.Taylor's Lagrangian macroscale can be related to the turbulent energy dissipation rate by the "Universal Kolmogorov constant" C, T L = 2•u 2 0 Cε , however, the constant C can vary between 3 and 10 depending on the Reynolds number of the flow [10].
Since the dispersion of discrete particles depends on the properties of the fluid and particles, one can of course classify heavy particles according to their diameter (characteristic size) compared to the different turbulent scales, the Kolmogorov scale or the smallest eddy size and the integral scale L. In solid-gas flows or particle-laden flows, it is often assumed that the particle size is smaller than the Kolmogorov scale η.
Nevertheless, other fluid and particle properties, such as fluid and particle densities, control inertia and gravity effects in turbulent dispersion.The following section will focus on the equation of motion of a single heavy particle in a fluid.

Equation of Motion
The trajectory of a discrete particle in a flow is governed by its equation of motion.A sphere moving in a fluid is subjected to gravity, viscous and pressure forces, and other non-stationary forces.Several textbooks and articles on multiphase flows or reviews have provided a summary of the different equations of motion of an isolated sphere [32][33][34][35][36].
One of the most recent forms for the equation of motion for an inertial sphere was given by Maxey [35] or Maxey-Riley [36].If a particle has a diameter d p and a mass m p , its motion is given to a good approximation by Newton's equations: where Here For turbulent dispersion in gas-solid flows, it is reasonable to consider uniquely the drag and gravity forces, eventually buoyancy, and for heavy particles buoyancy can be neglected as well as non-stationary forces, Basset history forces, Faxen term, added mass.The equation of motion of a spherical particle can then be simplified to an equation generally used in gas-solid flows still including buoyancy: or with gravity in direction 3: where → v (t) is the instantaneous particle velocity, → u p (t) is the instantaneous fluid velocity at the position of the discrete heavy particle, ρ = ρ p /ρ F is the ratio of particle to fluid specific mass, ) is the Stokesian particle relaxation time, or the particle aerodynamic response time (assuming Stokes drag applies) and the factor f is a correction to small particle Reynolds numbers and is given, for instance, by the drag correlation of Schiller and Naumann [37] and generally used in gas-solid flows [32][33][34]38,39].
the ratio of the effective drag to the Stokesian drag (non-Stokesian correction factor), Re p is the particle Reynolds number, is the mean drift or fall velocity due to gravity g in direction 3 (v ch > 0).
For dilute gas-solid flows, Equation ( 21) can be honestly reduced to the simplified equation of motion without buoyancy, with The instantaneous slip velocity v i (t) − u p i (t) is of importance because it fixes the particle Reynolds number Re p , the nonlinear factor f and depends on the fluid velocity The equations can be nondimensionalized using u 0 for the velocity, T L or T mE for time, and L mE = u 0 T mE for length.Here, for convenience, the moving Eulerian integral time scale, T mE , is chosen.Multiplying Equation ( 22) by T mE u 0 yields These dimensionless equations of motion contain three parameters: • the Stokes number St = St mE , based on the Eulerian moving scale T mE is a measure of the relative importance of the particle inertia; it characterizes the particle's response to the turbulent fluid velocity fluctuations; • the drift parameter γ, which is a dimensionless drift velocity related to the turbulence level, given by • the drag correction factor f , a function that increases with Reynolds number (fluidparticle drift velocity or particle size).
The three parameters are not independent since the relative drift parameter γ and the Stokes number St depend on the particle relaxation time for a given turbulence; in case of no gravity, the particle motion is controlled by inertia (IE e.g., Inertia effect).With gravity, both inertia and drift will control what is called "crossing-trajectory effects"-CTE.Further nonlinear drag can modify the behavior since the parameter f is correlated to the drift velocity.
A fourth parameter has to be considered, for the dependence of γ on the relaxation time: it can be defined as an eddy Froude number [40], the ratio of eddy convective forces to gravitational forces

Note on Stokes Numbers
Stokes numbers are the ratio of the response time of the particle, τ s or τ p = τ s / f to a time scale of the fluid.τ p is a corrected response time due to non-Stokesian drag.Diffusion of tracers or dispersion of small particles is controlled by Lagrangian scales, but for larger and heavy particles dispersion is governed by Eulerian scales.That is why several definitions of the Stokes numbers exist, mainly because particles are carried by the mean flow (even with mean drift) but are dispersed by turbulence, but care has to be made with all publications mentioning a Stokes number: The motion of a single spherical particle in a homogeneous isotropic turbulent field differs from that of a fluid point due to both the characteristic properties of the particle (size, density, inertia and free-fall velocity) and the turbulent flow field (time and length scales).Considering heavy particles, in general, the Stokes number St and the drift term γ = v ch u 0 = gτ s u 0 control the trajectory and three effects control the dispersion.Inertia, the first mechanism, inhibits the particle response to high-frequency turbulent fluctuations of the carrier flow.Gravity can generate a relative mean drift velocity v ch and particles will not be kept trapped by the successive eddies but will cross them; this was first pointed out by Yudine [41] as the effect of crossing trajectories (CTE) or overshooting.Even in HIST, the Lagrangian fluid velocity correlation tensor is not spherical; dispersion in gravity direction (3 in the present text) is different from that in the lateral directions and this constitutes the continuity effect (CE) as introduced by Csanady [42].Spatial correlation in gravity direction is larger than in the lateral directions, in particular, if the drift parameter increases.
Although not unique, a classification is proposed in Table 1 considering the particle size, the Kolmogorov length and the Stokes number St = τ p /τ L where τ p is the particle relaxation time and τ L is the Lagrangian time microscale.This classification is based on physics and is, in the spirit of Reynolds-averaged Navier-Stokes (RANS) computations, certainly appropriate.Cases 1 and 2 refer to the assumption that d p > η K .In these cases, on a computational level, depending on the grid meshing, the filter scale D in case of LES and the particle size, another aspect will directly concern the equation of motion.The point-approximation for the particle and fluid velocities in the equation of motion can be no more valid, if d p > D in LES, and is anyway not valid in direct numerical simulation (DNS).Cases 1 and 2 concern large particles, bubbles or droplets.In these cases, the turbulence of the carrier flow will be modified, a complex problem in multiphase flows.When the particle is larger than the smallest turbulence Kolmogorov scales, the fluid-particle interactions are extremely complex.The influence of the higher wave number components decreases but at the same time the small-scale structures may be altered by the presence of the particles.Only a part of the spectrum-larger-scale structures-will contribute to the dispersion, and the particles will further be perturbed by turbulence of the surrounding flow, exerting non-uniform and unsteady perturbations on the particles.If the drift parameter is non-negligible, inertia (IE), CE and CTE effects will all influence dispersion.Large and heavy particles (case 1) will only respond easily to the largest energy-containing structures, IE, CE and CTE interact strongly.Large but less inertial particles (case 2) will be carried by all structures larger than their size d p .
The cases 3 and 4 hold for particles whose size is smaller than the dissipation scale, typical of classical gas-solid flows, particularly in atmospheric turbulence.In dilute twophase flows, the carrier flow properties are not modified by the particulate phase.
Turbulent diffusion (case 4) concerns tracer particles or small droplets.The entire spectrum of turbulent energy of the carrier fluid then participates in the turbulent diffusion of the particles, which easily recognize all the turbulent fluid fluctuations.Furthermore, the particle will travel in a uniform velocity field during its residence time (the eddy lifetime) within each specific eddy, without inertia.Then the particles will respond to the entire energy spectrum of the fluid, and it is reasonable to expect that the root mean-square turbulent velocity of the particle coincides with that of the fluid, as the Lagrangian time scales and the diffusion coefficients, too: Case 3 is typical of heavy particles in gas-solid flows and can concern liquid-solid flows too.If τ p > τ L and d p < η K the response to rapid velocity fluctuations is damped or cut-off, although the particles move in and through all possible-sized eddies.In the absence of gravity, with an increase of the inertia, e.g., for small dense particles, the rms particle velocity will decrease v 2 < u 2 , whereas it is expected that the particle Lagrangian integral time scale T p L will eventually increase, depending on the relationships between the fluid Eulerian and Lagrangian scales.However, the effect of inertia on particle dispersion is not obvious since according to the statistical theory of turbulent dispersion, the particle diffusivity is the product of the mean-square fluctuating particle velocity v 2 and the Lagrangian time scale T p L , although a decrease of the dispersion coefficient is intuitively expected with inertia effects.In the presence of gravity, CTE and CE effects must be taken into account, but how they simultaneously modify the dispersion process is quite unclear.Nevertheless, the CTE influence is opposite to that of inertia; the particle integral time scale is reduced since the particle will, by overshooting, lose correlation.In that case, it is expected that Gravity can give rise to appreciable relative mean velocity between the particle and the carrier fluid.The particles will not remain within an eddy but migrate from one to another.An overshooting or free-fall time scale can be defined that represents the rate of change of the eddy structures by the particle.This characteristic time for the crossing-trajectories effect could be: where v ch is the relative fluid-particle velocity V-U or the free-fall velocity and l e an appropriate turbulent eddy length scale.Nevertheless, there are still two crucial questions.First, a particle moving in a turbulent field is influenced by the instantaneous fluid velocity at its position X(t), which means that the particle will interact with a fluid spectrum along its trajectory; this spectrum is not necessarily the fluid Lagrangian or Eulerian spectrum.Secondly, the last effect controlling turbulent dispersion is less clear.If in homogeneous isotropic stationary turbulence (HIST), it is evident that the dispersion tensor should be spherical without gravity; intuitively, in such a turbulent flow, how gravity perhaps modifies this dispersion tensor due to a continuity effect induced by gravity has to be well understood.The following section summarizes analytical studies including inertia and gravity effect before the advent of numerical methods.It will focus on the most important contributions to better understand the influence of inertia and gravity, CTE-CE, which control the turbulent dispersion of heavy particles in HIST.

Short Review of Analytical Approaches
A one-dimensional theory of discrete particle dispersion was proposed by Tchen [43].In Tchen's theory, it is supposed that no overshooting takes place: if the particle has a small size, it will remain captured by the same fluid elements along its trajectory.This assumption is hardly satisfied, since external forces, gravity, injection conditions can cause relative mean velocity between a discrete particle and the fluid particle.This is particularly true for heavier particles and when the ratio between the material densities of the particle and fluid is large.Tchen proposed an equation of motion that, once simplified, reduces to linear drag: τ s where τ s is the Stokesian relaxation time.By introducing the fluid and particle Lagrangian energy-spectrum density functions, a Fourier transform solution of the one-dimensional equation of motion without gravity [44,45] can be obtained.The main results are (known as the Tchen-Hinze model [46]) and the Schmidt number = T L T L +τ s , valid for small relaxation times and no overshooting.The particle mean-square velocity and diffusion coefficient should decrease with inertia.Insofar as drag is a damping force, the dynamics of heavy particles is reduced to the small wave numbers of the spectrum, i.e., the larger structures.
Early studies [47,48] on turbulent diffusion in the atmosphere demonstrated that heavy particles do not move surrounded by the same fluid particle.When they fall, they cross many turbulent structures and therefore move from eddy to eddy.The concept of slip fluid-particle velocity or overshooting was thereby introduced in order to describe the deterministic gravitational drift drawing particles out of the strongly correlated fluid environment and thereby reducing the particle diffusivity, as will be shown.
The effect, which might be called an "effect of crossing trajectories," is described as such: when falling, a heavy particle crosses trajectories of air particles so that it interacts consecutively with different air particles.As a result, the succession of velocities of a heavy particle does not coincide with individual changes of the velocity of an air particle [41] (p.186).Yudine stated that the vertical dispersion process depends upon the terminal velocity in three ways: the terminal velocity determines the vertical displacement of the center of dispersion of the particle; this effect is easily accounted for by introducing a convective term into the diffusion equation [49,50]; -since the terminal velocity is a measure of the inertia, the particle does not completely follow the high-frequency fluctuations of the turbulent fluid velocity; thus, Yudine did not separate inertia and gravity effects; -if it has an appreciable settling velocity, a particle will fall from one eddy to another, whereas a fluid point will remain in the same eddy throughout the lifetime of the eddy; this is one of the first papers mentioning "overshooting".
Yudine formulated upper and lower limits for the change in the dispersion coefficient due to the free-fall velocity.In fact, the development of that study was upon space-time correlations including the free-fall velocity.Particularly, for large terminal velocity v ch , the dispersion coefficient should take an asymptotic form inversely proportional to the terminal velocity v ch , D P ∝ (v ch ) −1 .
Csanady [42] accounted for the crossing-trajectory effect caused by gravity and estimated the reduction of dispersion rate of the heavy solid particles attributed to their trajectory across the turbulent eddies.Csanady included parameters controlling the dispersion of heavy particles.
The first one is the ratio v ch /w , the vertical drift parameter γ, which can be considered as a measure of the crossing-trajectories effects, where v ch is the free-fall velocity and w' is the vertical fluid turbulent intensity w 2 1/2 .The second one, β, is a ratio between the Lagrangian and Eulerian integral time scales (Section 2.2.1).A functional form for the vertical velocity correlation w(0)w(τ) consistent with similar shapes for Eulerian spacetime and Lagrangian fluid point correlations was proposed: and L E = L Ez = L E33 is the vertical integral Eulerian length scale.The asymptotic diffusivity in the vertical direction (subscript z or 3) is easily calculated by Batchelor's relations: At large free-falling speeds, both the dispersion coefficient and the variance become proportional to (v ch ) −1 , as already stated by Yudine [41].
Csanady [42] introduced the concept of "Continuity effects," which were found to reduce the lateral dispersion.Continuity effects refer to the relations existing between the space correlation functions via the turbulent part of the continuity equation of the fluid.Therefore, lateral correlations must exhibit negative loops as a consequence of backflow [51] (p.18).In terms of diffusion coefficients, Csanady's study leads to corrected relations for long diffusion times that could be applied to atmospheric dispersion of heavy particles; the asymptotic long time Schmidt number is given by: where D f is the asymptotic fluid diffusivity.β is probably different in the three directions, due to continuity effects.If we consider exponential space correlations in isotropic turbulence, as stated in Section 2.2.1: , the transverse correlation lengths being one-half of the longitudinal one.As a consequence, in lateral directions x and y (or 1 and 2): Later, Lumley [52] (p.305) reconsidered Csanady's results and proposed a few changes.Lumley showed that Csanady's model gives the subscripts x and z corresponding to horizontal and vertical dispersion of the particles in the directions that are perpendicular and parallel to the gravity field.In the limit of a large drift parameter ( γ = v ch /w → ∞ ), the lateral dispersion coefficient is one-half of that in the vertical direction, a characteristic of the continuity effect and loss of correlation.
If one assumes that the particle and fluid velocity variance are equal in terms of Lagrangian time scales, as derived from Csanady's analysis, two-particle Lagrangian integral scales can be then obtained: for the gravity direction, and for the lateral case, with β = w T L /L E = T L /T E indicating a loss of correlation because of gravity.Nevertheless, Lumley stated that the value β x = 4/3 is much less well determined than β z = 2/3 and should be determined via a comparison with adequate experimental data.
Theoretical contributions [53][54][55] by spectral analysis of Tchen's equation followed in terms of transfer functions, confirming that for heavy particles only drag and gravity will control dispersion.Of interest are two results by Meek and Jones [54]: .
Fluids 2021, 6, 145 13 of 48 The inertial effects can be significant, particularly when v ch < u 0 .The inertial effects increase the particle Lagrangian time scale compared to that of the fluid if there is no crossing-trajectories effect.

•
the ratio of the fluctuating velocity variances is .
As expected, free-fall effects reduce dispersion and as drift velocity increases its correlation decreases (CTE effects).
The ratio of the Lagrangian particle and fluid macroscales is written as F , which indicates that if inertia effects are dominant, they tend to increase the particle integral scale in particular if the drift parameter is small.Inertial effects tend to modify the drift effects, acting to give enhanced correlation and, thus, enhanced dispersion.For large settling velocity, the Lagrangian particle macroscale T p L varies again as 1/v ch .Thus inertia and CTE can have opposite contributions to turbulent dispersion, which is determinant in understanding the mechanisms of dispersion of heavy particles.
Reeks [40] brought a significant contribution to the understanding of inertia effects.Reeks discussed the particle dispersion in a stationary homogeneous isotropic turbulence in both the presence and absence of gravity, using the method initially developed by Phythian [56] for fluid dispersion.His results indicate that without settling velocity (without gravity), the asymptotic particle dispersion can be higher than the fluid diffusion if the structure parameter m = 1, a notable contribution that was experimentally verified earlier in several cases by Snyder and Lumley [15].With increasing Stokes number, the turbulence experienced by a heavy particle is much more correlated to a Eulerian frame, considering that Eulerian scales can be greater than Lagrangian scales.
At the same time Pismen and Nir [57], according to Kraichnan's theory [58], explicitly demonstrated that the particle long time diffusion coefficient can become larger than the fluid Eulerian diffusivity (Figure 3 in Ref [57]) if the particle relaxation time is significant in the absence of crossing-trajectory effects.The case of additional free-fall velocity was then introduced in a subsequent paper by Nir and Pismen [59].Free drift velocity is a source of anisotropy and often present in experimental studies where particles are submitted to gravity and attain a steady terminal sedimentation velocity.A second parameter β * = γ, is discussed under gravity effects.They showed that for small values of β * , the particle transverse correlation functions exhibit negative loops as a consequence of the continuity effects.

General Description
The first stochastic models used to simulate the turbulent dispersion of particles were the so-called eddy lifetime or discrete turbulent eddy concept models and were based on the Monte-Carlo process.The impulse was given by the combustion community; even developing Eulerian/Eulerian models, a Lagrangian approach treating the behavior of non-fluid particles seemed essential because it concerns droplets with inertial effects.The impulse was always the result of the growth of computer performance able to solve such problems within short times.The computational methods to solve an equation of motion being available as well as random numbers generation methods, the harder work was to simulate the turbulence field and its interaction with discrete particles along their trajectories.
The emerging concept in eddy interaction models is that the interactions between a particle and a succession of fluid eddies can be characterized by three parameters: an eddy instantaneous velocity, an eddy lifetime and an eddy size, all relevant to the structure of turbulence, the Reynolds stress tensor and the energy spectrum, among others, combined with Lagrangian and Eulerian time and length scales.The velocity and, in the most general case, the lifetime and the length of the eddy are all independent random variables.The turbulent velocity is sampled randomly from a Gaussian probability distribution function with standard deviation Even if the principle of the eddy interaction models is to allocate to the randomly sampled eddy physical properties such as a length or lifetime, the main difficulty is to statistically simulate the turbulent fluid velocity along the particle trajectories, sometimes referred to as the turbulence "seen" or "experienced" by the particle.All the models proposed in the literature try with more or less success to optimize an eddy interaction time, the transit time a particle needs to cross an eddy, or a correlation length (the eddy size), but without modeling space-time correlation functions.
Both the inertial particle (P) and the fluid eddy (F) are tracked simultaneously (Figure 1).At the initial instant (t), the particle and the fluid element are at the same location; then, on account of their different natures, they will separate.Their respective velocities v p (t + ∆t), u f (t + ∆t) and the fluid velocity at the particle position u p f (t + ∆t) will nevertheless be correlated during a process time T EI M depending on two situations: • the particle (P) leaves the fluid eddy (F) for entering into a new one when the eddy (F) lifetime τ * is elapsed (condition 1), or when the distance between the particle and the eddy center does exceed the eddy length λ * a radial dimension (condition 2).
problems within short times.The computational methods to solve an equation of motion being available as well as random numbers generation methods, the harder work was to simulate the turbulence field and its interaction with discrete particles along their trajectories.
The emerging concept in eddy interaction models is that the interactions between a particle and a succession of fluid eddies can be characterized by three parameters: an eddy instantaneous velocity, an eddy lifetime and an eddy size, all relevant to the structure of turbulence, the Reynolds stress tensor and the energy spectrum, among others, combined with Lagrangian and Eulerian time and length scales.The velocity and, in the most general case, the lifetime and the length of the eddy are all independent random variables.The turbulent velocity is sampled randomly from a Gaussian probability distribution function with standard deviation Even if the principle of the eddy interaction models is to allocate to the randomly sampled eddy physical properties such as a length or lifetime, the main difficulty is to statistically simulate the turbulent fluid velocity along the particle trajectories, sometimes referred to as the turbulence "seen" or "experienced" by the particle.All the models proposed in the literature try with more or less success to optimize an eddy interaction time, the transit time a particle needs to cross an eddy, or a correlation length (the eddy size), but without modeling space-time correlation functions.
Both the inertial particle (P) and the fluid eddy (F) are tracked simultaneously (Figure 1).At the initial instant (t), the particle and the fluid element are at the same location; then, on account of their different natures, they will separate.Their respective velocities • or when the distance between the particle and the eddy center does exceed the eddy length λ * a radial dimension (condition 2).The time step T EI M of the Monte-Carlo process will be, if condition 1 is verified, the eddy lifetime τ * or, if condition 2 is respected, the transit time τ R taken by the inertial particle to cover a distance larger than the eddy length λ * , that is, X p (t) − X f (t) > λ * .The Monte-Carlo time is given by T EI M = min(τ * , τ R ).Thus the eddy-particle interaction time is limited by the eddy lifetime or the transit time needed for the particle to cross the eddy.The eddy lifetime τ * is related to a Lagrangian integral time scale T L or to a spectrum of the fluid turbulent flow, the eddy length λ * to an Eulerian length scale L E of the turbulence.This kind of model should at least simulate part of the crossing-trajectory effects (CTEs), and, if possible, the inertia effects by solving the equation of motion with a time step ∆t << T EI M .In the literature, the EIMs differ in the way the different time scales T L , τ * , ∆t or length scales L E , λ * are chosen.In these stochastic trajectory models, the random fluid elements interacting with a particle have a length and time scales , λ * , τ * , and the fluid velocities can be generated with constant ∆t or random time steps δt i .The velocity correlation functions are dependent on the choice of the type of time step probability density function f (δ t i ) as discussed by Wang and Stock [60].In some cases, the Monte-Carlo process uses exponential distributions for the random eddy lifetime and size, with mean values equal to the Lagrangian integral time scale T L and to the Eulerian transverse length scale L E .

Historical Background
The basic ideas of a stochastic process for turbulent dispersion of particles were probably first proposed by Hutchinson et al. [61], and later by Brown and Hutchinson [62].Eddies would be characterized by both a mean decay time and a "contact" time between an eddy and a particle, which could not exceed this decay time; a stochastic model simulated the random displacement of the particle colliding with eddies.Hotchkiss and Hirt [63] proposed a particle tracking scheme to simulate material transport with reference to a diffusion equation.
Later, Yuu et al. [64] assumed the particles would be trapped in successive energycontaining eddies without overshooting as in Tchen's theory [43], each eddy having a constant lifetime and a random velocity fluctuation.Dukowicz [65] coupled a trajectory model to the SOLA code developed at Los Alamos by Hirt et al. [66] for liquid fuel sprays.The Monte-Carlo model was extended by Gosman and Ioannides [67,68].The characteristic size of an eddy was set to the dissipation length scale λ * = λ e = C 3/4 µ k 3/2 /ε and the eddy life time was τ * = τ e = λ e /u where C µ is an empirical constant of the k − ε turbulence model TEACH-T developed at Imperial College.The eddy interaction time is then given by τ int = min(τ e , τ R ) where τ R is the transit time.The eddy transit time can be easily obtained by a linearized equation of motion: where v R is the slip velocity and the particle-eddy interaction time is therefore obtained classically by It appears that the eddy interaction time for inertial particles cannot exceed the eddy lifetime, which is the fluid-particle interaction time or the eddy lifetime with an upper limit T L , the Lagrangian fluid integral scale.As a consequence, the heavy particle dispersion can never exceed the fluid-particle diffusion, a deficiency of most of the EIMs.
Chen and Crowe [74] simulated the experiment in a pipe flow of Arnason [75,76] and for the near isotropic homogeneous grid turbulence in Snyder and Lumley's experiments [15]; however they calibrated the model by changing the constant C µ to fit the last experiment.Durst et al. [77] developed a simple EIM based on the previous work of Gosman and Ioannides [67,68], and the equation of motion was solved analytically for short time steps.Milojevic [78] optimized the model for various test flow configurations [15,16,75,76].It was modified for particulate two-phase systems with a size distribution by Sommerfeld and coworkers [79,80].Mostafa et al. [81][82][83][84] coupled a k − ε model to the EIM to simulate experimental studies in turbulent evaporating sprays by optimizing the initial constants of the k − ε to functions.They used the relationship T L = 0.35k/ε.Govan et al. [85] developed a simplified EIM to simulate their experiments in gas-solid pipe flows.
Ormancey and Martinon [86,87] proposed a 2-D model and introduced a procedure in order to calculate the fluid velocity at the particle position, taking into account a space correlation based on Frenkiel's functions [88].A relevant constant time step ∆t has to be used to numerically solve the equation of motion by a classical time discretization of the ODE's.To account for the Lagrangian time correlations, eddy lifetimes are sampled according to a Poisson process.A variable y is sampled with a uniform probability density on [0, 1] and for each time step, two cases are possible.If y > ∆t/T Eulerian step is added to make a space correction to u f (t) to obtain the fluid velocity u P f (t) at the particle position, where r is the distance between the particle (P) and the eddy (F) at time t.This fluid velocity is sampled according to a Gaussian probability-density function with In the case where ∆t << T f L , it can be demonstrated that the process obeys an exponential fluid correlation function: The procedure can be divided into three steps, • a first estimate of the fluid velocity with the criteria y < ∆t/T f L , • a spatial step that calculates the fluid-particle velocity u p f (X p (t), t) at time t at X p (t) with a correlation function: • the integration of the equation of motion with u p f (X p (t), t) and v (t) to obtain the new discrete particle velocity.
Part of this method was used later by Hajji et al. [89] and compared to a spectral direct simulation model for heavy particles in HIST; nonlinear drag was also investigated for large Stokes numbers with drift velocity.
In some cases, the Monte-Carlo process used exponential distributions of the random eddy lifetime and size, with their mean values equal to the Lagrangian integral time scale T L and to the Eulerian transverse length scale L g .Kallio and Reeks [90] modeled particle deposition in an inhomogeneous boundary layer with random time scales drawn from an exponential probability density function (PDF).Burnage et al. [91][92][93][94] proposed an EIM based on Poisson and exponential distributions for the random eddy lifetime and size.
Most complete papers on the eddy interaction models and characteristics of the associated Monte-Carlo process are those written by Wang and Stock [18][19][20] and Graham [11,[95][96][97] with mathematical arguments to clear the problem of self-consistency as pointed out by Kallio and Reeks [90], but also to make a real distinction between the standard eddy lifetime τ * or τ e also defined as the fluid-particle interaction time T f pit by Graham, and a maximum interaction time T Max , which can be randomly sampled, too [97].T Max is independent of T f pit and its choice determines a relationship between the turbulence structure parameter β = u T L /L E and the ratio of the Lagrangian and Eulerian integral time scale T L /T E .It enables finite-inertia particles to disperse faster than fluid particles.Clearly, the particle-eddy interaction time is given by where the factor 2 holds for self-consistency.Wang and Stock [60] analyzed the relationships between different time step probability density functions f (t e ), their mean value T and the Lagrangian integral fluid correlation time T L .investigated the performance of the EIM in several of his papers and proposed enhancements in order to take into account the three effects, CTE, inertia and continuity, in inhomogeneous flows.He completed the study of Wang and Stock [5] and provided relations between the EIM, Lagrangian autocorrelations and Eulerian (fixedpoint) correlations.These improved EIMs predicted that inertia effects without gravity can in some cases increase the long time dispersion of heavy particles, as was already obtained theoretically by Reeks [40] and Pismen-Nir [57], by the experiment of Wells and Stock [16] or by DNS [98] or LES [99].The best agreement was obtained with the analytical results of Reeks in HIST.Furthermore, Graham [97] extended his work to random length and time scales, with valuable results for non-homogeneous flows.Chen and Pereira [100] tested the Graham model; Chen [101] used the Graham/Gosman-Ioannides EIM for inhomogeneous anisotropic turbulence obtained by a Reynolds stress model and studied the inertia effects.
The advantage of the EIMs is that no time or space correlation has to be introduced in order to obtain the fluid velocity at the particle position.However, it is necessary to provide a time and length scale to characterize the interacting eddies.However, the method used for the time step or the eddy lifetime of the EIM, the time step probability density functions f (t e ), is related to the form of the Lagrangian fluid autocorrelation R f L (τ) and T f L .Thus, with a PDF function All the proposed models generally link the eddy characteristics (lifetime, size) to the turbulent kinetic energy and the dissipation rate and the Lagrangian and Eulerian scales with help of constants.
The constants A and B are adjusted to let the models simulate the classical test cases [102].However, the method fixing the time step of the EIM, the time step probability density functions f (t e ), is related to the form of the Lagrangian fluid autocorrelation R f L (τ) and to the integral time scale T f L .For instance, the way the constants A, B, C T are calibrated to fit a specific experiment is questionable, since they will in fact fix a virtual constant C * T .The corresponding time scale, T L = C * T 2k 3ε , will not be the Lagrangian fluid integral scale of a fluid particle T f L but the effective Lagrangian fluid integral scale of a discrete particle T p L along its trajectory.In that case the model will be optimized anyway to fit the simulation of the experimental test case.
Another advantage is that the EIM can easily be implemented in Eulerian codes of turbulence, and they can give a first qualitative and even quantitative idea of the dispersion.An open question in the EIM still exists in case of time and length distributions: since physically it is expected to have a relation between eddy size and lifetime, a small lifetime is expected for small dissipative eddies.A relationship between eddy size and lifetime should exist; the eddy turnover time t l of a size l eddy can be related to the time it takes for that size eddy to be reduced to the Kolmogorov scale, with an integral length scale L. Such relation should be

Classical Approach
The stochastic differential equation (SDE) modeling the behavior of fluid velocities in Gaussian turbulence, verifying Kolmogorov's similarity theory for locally isotropic turbulence is Langevin's equation contains both a damping or drift term, which accounts for the friction exerted by the surrounding fluid on the heavy particle, and a random diffusion term, which models the fluid impulses, dW, being a Wiener white noise process, a stochastic process of zero mean and variance equal to the time interval dt.Equation ( 31) is the original Langevin equation for homogeneous stationary one-dimensional turbulent diffusion with no mean flow, and a Lagrangian correlation function R L (t) = exp(−t/T L ).The discrete form of the Langevin's equation (or Markov chain equation) correlates the fluid velocity fluctuation u f (t + ∆t) to its value u f (t) at the previous step: where u 0 is the rms value of the fluid velocity, R L (∆t) is the Lagrangian autocorrelation and ξ is a statistically independent, Gaussian, normalized and centered random number.
A classical Lagrangian autocorrelation function of the fluid velocity is simply given by this process: R L (∆t) = exp(−∆t/T L ).
∆t is the time step and T L the Lagrangian integral time scale of the turbulence.Note: For the simulation of the diffusion in a non-homogeneous turbulence, Equation ( 31) is not valid and Langevin's equation including drift correction terms, as suggested by Legg and Raupach [108] must be used.This will be discussed in Section 5.
The trajectory of fluid-discrete particle pairs is the elementary approach to model the fluid velocity seen by a discrete particle and is divided in two steps (Figure 2).Initially the fluid (F) and discrete (P) particles are at the same location → X f p (t) at time t. and their respective velocities u f (t) = u p f (t) and v(t) are known.After a short time step they are located at X f (t + ∆t) and X p (t + ∆t).The first step uses the discrete Langevin equation for the fluid element to estimate the eddy fluid velocity at X f (t + ∆t) and is a Lagrangian process in time: The idea of the second step is based on the distance ∆s = that separates the pair at time t + ∆t and uses a Eulerian space correlation coupled to the Langevin equation to obtain the fluid velocity at → X p (t + ∆t) associated to a new fluid element (F').Then, a new fluid element (F') is centered on the inertial particle (P), and the procedure begins again.
A combined two-step relation was first proposed for simplicity [117,118].This version directly coupled the time and space correlations as In case of a heavy particle with a mean fall velocity of ch v , t v s ch Δ ≈ Δ and a large drift parameter , a Eulerian-Lagrangrian relationship [12,13].This was suggested by Walklate [13] and can be compared to Csanady's formulation: Hunt and Nalpanis [119] suggested a different correlation function optimizing Snyder's measurements [15]: Simulations in HIST [120] with the above equation were performed for heavy particles and compared to other simulations obtained with a Frenkiel zero correlation (m = 0) and the theoretical results of Wang and Stock.
Furthermore, much effort was dedicated to anisotropic turbulence.Berlemont et al. [121][122][123] proposed a sophisticated formulation; it will be briefly discussed in Section 3.3 since it is a multistep random method.Zhou and Leschziner [124] introduced a temporal correlation approach taking into account anisotropic Lagrangian time scales, including This value is used to solve the equation of motion and to calculate the velocity → v (t + 2∆t) and the new position of the fluid element (F') and of the discrete particle → X p (t + 2∆t).A new fluid element is associated to the discrete particle at that new location.At each step a Cartesian frame is defined, with the origin at the fluid element position, and the first axis joins the inertial particle position to take into account the specificities of the directional correlation functions.
Then, a new fluid element (F') is centered on the inertial particle (P), and the procedure begins again.
A combined two-step relation was first proposed for simplicity [117,118].This version directly coupled the time and space correlations as In case of a heavy particle with a mean fall velocity of v ch , ∆s ≈ v ch ∆t and a large drift parameter γ = v ch u 0 = gτ s u 0 << 1: [12,13].This was suggested by Walklate [13] and can be compared to Csanady's formulation: Hunt and Nalpanis [119] suggested a different correlation function optimizing Snyder's measurements [15]: . Simulations in HIST [120] with the above equation were performed for heavy particles and compared to other simulations obtained with a Frenkiel zero correlation (m = 0) and the theoretical results of Wang and Stock.
Furthermore, much effort was dedicated to anisotropic turbulence.Berlemont et al. [121][122][123] proposed a sophisticated formulation; it will be briefly discussed in Section 3.3 since it is a multistep random method.Zhou and Leschziner [124] introduced a temporal correlation approach taking into account anisotropic Lagrangian time scales, including directional and temporal correlation coefficients between u i (t + ∆t) and u j (t) and the covariance matrix of u j (t).Burry and Bergeles [125] and Lu et al. [126-130] generalized the method to anisotropic flows, with temporal and spatial correlations, according to the anisotropy of the Reynolds stresses.Lu used a Eulerian velocity autocorrelation instead of a Lagrangian one, Mashayek [131] compared results obtained with Lu's Eulerian stochastic model with the direct numerical simulation (DNS) data of Mei et al. [132].
The proposed models are related to complex generalized Langevin equations use relationships , relationships often similar to those of the EIMs, or more complex correlation functions in case of tensors.As for the EIMs, the different constants are tuned to best simulate turbulent diffusion and other experiments and test cases, and are compared to DNS results.Other papers have compared the performance of different Lagrangian stochastic models, too.For the Snyder and Lumley experiment [15], an extensive comparison of different EIM and RWM results obtained by Berlemont [122], Gosman and Ioannides [67,68], Ormancey and Martinon [86,87], Shuen et al. [69], Lu et al. [127], Nir and Pismen [57], and Meek and Jones [54] was proposed by Huilier et al. [133].Apart the choice of the constants A, B or C T , some other factors can influence the simulations and justify discrepancies, the number of computed trajectories, the choice of the random number generator or the time step ∆t.Lain and Grillo [134] simulated the Wells and Stock experiment [16] and a jet flow [135].
Within the last two decades, research has been focused on enhancements of RWM and on the formulation of Lagrangian stochastic models and PDF models [131, for the transport of particles in more complex turbulent fields, and substantial work was also addressed on the spurious mean particle drift in inhomogeneous turbulence as introduced by McInnes and Bracco [102] or Legg and Raupach [108].Zhou and Leschziner [136] presented more sophisticated models able to integrate directional/cross correlations (full correlation tensors R ij (∆ t) or the directional correlations, and the anisotropy of turbulence scales).These models should be able to simulate the dispersion of particles in anisotropic turbulence, shear flows or boundary layers.It is not the intention of the present paper to detail these approaches, but some of these methods are summarized and discussed by Pozorski and Minier [139], Bocksell and Loth [142,149], Minier et al. [157] or Tanière and Arcen [158].The models were also modified according to the true fluid velocity along a particle path and the classical fluid integral time scales T Lij where replaced by integral scales integrating inertia and gravity in relation with Csanady's formulas (Section 2.2.4).Most models often used also the theoretical results by Wang and Stock developed in Appendix A.

CRW-Matrix Methods
Several models used multistep Markov chains, which is the fluid velocity along the trajectory of the discrete particle at time t is evaluated with more than one previous calculated velocity, that is Berlemont et al. [121][122][123] were the first to propose a matrix method able to simulate the fluid velocity along a part of the trajectory for any time correlation function.The aim was to generate a sequence u(∆t), u(2∆t), . . .u(i∆t), . . ., u(j∆t) . . .verifying Lagrangian Frenkiel family correlation functions and the stress tensor.At the same time for a given time step, a Eulerian correction was applied; once the distance between the two particles exceeded the lateral Eulerian scale, a new fluid particle was generated.However, the fluidparticle velocity field was that of a tracer, because related to the classical fluid-particle correlation, the Eulerian process should simulate the CTE and inertial effects.This kind of method can be classified as a hybrid RWM-matrix method since they use a multistep Markov chain simulation of the fluid-particle velocity.
In Berlemont's formalism, the fluid-particle fluctuating velocities are given by where the vector Y has uncorrelated variables with a Gaussian distribution and the matrix B is related to the Lagrangian time correlation tensor and a matrix A by The set of equations can be rewritten as a generalized Langevin equation under matrix form where the matrix M includes the effects of previous time steps, vector d n is sampled from Gaussian PDF, and the fluid velocity at the new particle position is then obtained with Frenkiel space correlation functions.A matrix method, different from that earlier proposed by Berlemont et al., has been developed by the research group at Strasbourg University [159][160][161].The idea is to generate a random sequence of fluid velocities for the whole trajectory.The time series algorithm is able to reproduce stochastically any fluctuating velocity field in HIST, at least if the autocorrelation tensor R P F,ii (τ) is known.Starting with the following system, where at each time step, the fluid fluctuating velocity is a linear combination of random numbers, statistically independent, centered and normalized ξ k with a Gaussian probability density function: and by multiplying successively all the equations of the system with the first equation, and averaging the obtained products, we obtain: .
This yields the first column of the low triangular matrix A. Then the initial set is multiplied by the second equation and so on, and the whole matrix A can be built.Matrix A is unique for each chosen Lagrangian autocorrelation.A given set of random numbers ξ k will represent a fluid velocity along a particle trajectory (indeed, a linear system has to be solved for each trajectory) and will be used to solve n times the momentum equation of motion.It is clear that for long time diffusion T, this method requires computer resources since the dimension of the matrix will be n = T/dt.Nevertheless, correlation is typically lost after a time T max = 5T L and the triangular matrix A can be replaced by a band matrix B (Figure 3), the terms of that matrix can be used to calculate the fluid velocity and the number of floating point operations will be reduced.Pétrissans et al. [162] later proposed a more effective model, less computer timeconsuming, based on different time series, an autoregressive-moving-average ARMA (2,1) process, using a second-order time series and a first moving average.The model was able to generate the Wang and Stock's Lagrangian correlation functions and predicted the effect of negative loops.

Numerical Context
Most of the Lagrangian models, among which are either EIM or RWM, use experimental data or test flows to verify their ability and performance to simulate turbulent dispersion.First of all, they must be able to simulate turbulent diffusion in homogeneous isotropic turbulence as predicted by the analysis of Taylor; furthermore, they are generally tested and optimized to predict the dispersion of particles or droplets in various configurations, grid turbulent flows [15,16], horizontal or vertical pipe flows [17,75,76,[163][164][165][166], round jets [84, [167][168][169][170][171][172] or in mixing shear layers [173].Last, they are often used to simulate the dispersion in HIST, and their performance is tested on theoretical or numerical results obtained by other techniques.One of the aims of the present paper is to verify if these methods can take into account the inertia and gravity effects.
Concerning the flow and particle properties in HIST, we focus on the impact of inertia and drift effects governing the dispersion of spherical particles, with a density of ρ p = 1000 kg/m 3 and a diameter d in the range of 10 to 500 microns.The spheres are released from a point source and disperse in a homogeneous isotropic stationary turbulence (HIST), the fluid flow being air under standard conditions (density ρ F = 1.275 kg/m 3 and kinematic viscosity ν F = 1.36 × 10 −5 m/s 2 ) without gravity and a gravity field.The scale of the turbulent fluid is u 0 = 0.131 m/s and the Lagrangian integral time scale is set to T L = 0.091 s.The experimental data by Sato and Yamamoto [174] suggested that for a grid turbulence T L = 0.3 − 0.6 T mE where T mE is the moving Eulerian integral time scale.Burnage and Huilier [175] performed a diffusion experiment with small droplets released isokinetically from a tube in a grid decaying turbulence, using LDA and a light diffusion technique to measure the turbulence and concentration field and found T L = 0.5 − 0.95 T mE [13].In the present study, we suppose that T L /T mE = 0.356 (Appendix A).As a consequence, in the present simulations, the Stokes number, St = τ s /T mE , will range from 0.007 to 4.4, τ s being the Stokesian particle relaxation time (details are summarized in Table 2).All the statistics that will characterize the particulate dispersion, namely, turbulent energies (variances) v 2 ii , Lagrangian particle integral time scales T p,ii and long time particle dispersion coefficients D P,ii (∞) = v 2 ii • T p,ii , are calculated for 50,000 trajectories for each particle size (Stokes number), a number of trajectories sufficient to ensure statistical stability and convergence.The equations of motion as simplified for heavier particles (section equation X) are solved numerically by a fourth order Runge-Kutta algorithm, and the time step is set from 10 −4 s to 10 −3 s with particle relaxation times from 1.2 × 10 −3 s (d = 20 µm) to 0.8 s (d = 500 µm).
The drift parameter γ = ν ch /u 0 due to gravity g = 9.8 m/s 2 varies from 0.02 for the smallest particles to 15 for the largest.The long time dispersion characteristics are calculated for 4 s trajectories, a time that largely exceeds by a factor of 5 to 10 the particle relaxation time of the heavier particles and the Lagrangian scales.Besides other tests on nonstationary forces (added mass, Basset term), lift forces confirmed that statistically there is no influence on long time dispersion of heavy particles [93], in agreement with other results in the literature.

EIM-DRW Results
We consider here an extension of the EIM initiated by Burnage and Moon [91] and furthermore developed by Huilier et al. [93,133,[175][176][177].The eddy lifetime and the eddy length are distributed randomly according to exponential probability distribution functions with average values respectively equal to the fluid Lagrangian integral time scale T L in Section 4.1.1 and to the Lagrangian time proposed by Wang and Stock in the Appendix A in Section 4.1.2;for the length scale related to the eddy lifetime, the transverse Eulerian integral length scale L E was used.The instantaneous fluid eddy's velocity was given by the sum of a mean and a fluctuating velocity.Each component of the fluctuating velocity was evaluated from a Gaussian probability distribution with zero mean and a variance equal to u 2 0 , the rms value u 0 being also called "the eddy velocity scale."In this study, the assumption is made that during the interaction between an eddy and the particle, the fluctuating fluid velocity remains constant.Besides, it is supposed that the particle stays in a same fluid eddy as long as: (1) the distance between the eddy center and the heavy particle position does not exceed the random eddy length (no overshooting condition), (2) the interaction time of the particle does not exceed the random eddy lifetime.
If either of these two conditions becomes invalid, the heavy particle will be centered on a new random fluid eddy.

Classical Dispersion Without Gravity (g = 0)
Without any external force field (no gravity, g = 0), one can reasonably assume that the particle Reynolds number is less than one, since the mean fluid-particle slip velocity In that case, the drag law is linear (Stokesian, f = 1).All classical EIMs consider that the turbulence along a particle path has T L as Lagrangian time scale.Only inertial effects have to be taken into account with no gravity term, although crossingtrajectory effects cannot be excluded since due to inertia particles can cross small eddies rapidly.The general trend is that as expected, the turbulent particle integral scale increases with the Stokes number (Figure 4) and particle turbulence decreases with inertia (Figure 5).But very inertial (heavy) particles disperse much less than fluid elements (Figure 6), the dispersion coefficients being normalized by the fluid diffusion coefficient equal to u 2 0 T L .This is in contradiction with now well recognized results, either theoretical [52] or numerical results (obtained by direct or large eddy simulations).In conclusion, the inertial effect is not well simulated by this type of model based on T L .Simulations were performed with the same model without overshooting [176], they simply verified Tchen's theory, which is obsolete for heavy particles since crossing-trajectory effects are not considered.

Modified EIM with Wang and Stock Correction (No Gravity)
A similar simulation but with Wang and Stock's relation (Equation (A1) in Appendix A) gave much more promising results for the turbulent energies (Figure 7).In particular, Figure 8 indicates that the ratio of the particle dispersion coefficient to the fluid diffusivity is greater than unity for the large Stokes number; this result is in accordance with the theories of Reeks [40] or Pismen and Nir [57], with the experimental observations of Wells and Stock [16], with the direct numerical simulation of Squires and Eaton [98] and with the large eddy simulation of Deutsch [178].Only for large Stokes numbers are the particle

Modified EIM with Wang and Stock Correction (No Gravity)
A similar simulation but with Wang and Stock's relation (Equation (A1) in Appendix A) gave much more promising results for the turbulent energies (Figure 7).In particular, Figure 8 indicates that the ratio of the particle dispersion coefficient to the fluid diffusivity is greater than unity for the large Stokes number; this result is in accordance with the theories of Reeks [40] or Pismen and Nir [57], with the experimental observations of Wells and Stock [16], with the direct numerical simulation of Squires and Eaton [98] and with the large eddy simulation of Deutsch [178].Only for large Stokes numbers are the particle dispersion coefficients underestimated; this is due to the effects of nonlinear drag on the Lagrangian particle integral time scale and the particle velocity scales [176].Comparisons between results obtained with Stokes' drag and nonlinear drag were also proposed by Wang and Stock [6] and also showed an increase of the velocity variance with use of the nonlinear drag mainly due to the difference between the Stokesian relaxation times (and so the related Stokes number) being larger than the nonlinear terms when the particle Reynolds number Re p > 1 (Table 2).dispersion coefficients underestimated; this is due to the effects of nonlinear drag on the Lagrangian particle integral time scale and the particle velocity scales [176].Comparisons between results obtained with Stokes' drag and nonlinear drag were also proposed by Wang and Stock [6] and also showed an increase of the velocity variance with use of the nonlinear drag mainly due to the difference between the Stokesian relaxation times (and so the related Stokes number) being larger than the nonlinear terms when the particle Reynolds number 1 Re > p (Table 2).

Dispersion Under Gravity Effect
Fluid-particle slip velocity under gravity can be non-negligible for heavy particles, so that the particle Reynolds number can be much larger than 1 and a non-Stokesian drag  dispersion coefficients underestimated; this is due to the effects of nonlinear drag on the Lagrangian particle integral time scale and the particle velocity scales [176].Comparisons between results obtained with Stokes' drag and nonlinear drag were also proposed by Wang and Stock [6] and also showed an increase of the velocity variance with use of the nonlinear drag mainly due to the difference between the Stokesian relaxation times (and so the related Stokes number) being larger than the nonlinear terms when the particle Reynolds number 1 Re > p (Table 2).

Dispersion Under Gravity Effect
Fluid-particle slip velocity under gravity can be non-negligible for heavy particles, so that the particle Reynolds number can be much larger than 1 and a non-Stokesian drag law must be used, since the fall velocity increases with particle size.In the present case, a drag law given by Schiller and Naumann [37] was proposed.This law is generally ac-

Dispersion under Gravity Effect
Fluid-particle slip velocity under gravity can be non-negligible for heavy particles, so that the particle Reynolds number can be much larger than 1 and a non-Stokesian drag law must be used, since the fall velocity increases with particle size.In the present case, a drag law given by Schiller and Naumann [37] was proposed.This law is generally accepted, as other laws bring no real modifications.The present results consider again that the heavy particle is controlled by a fluid Lagrangian scale equal to T L .Although the influence of the time step and number of total trajectories on the results was investigated, the results are mixed and scattered for the larger Stokes numbers at least for the time scales (Figure 9) and the particle dispersion coefficients (Figure 10).The normalized variances of the particle velocity fluctuations (Figure 11) and the dispersion coefficients are overestimated for the larger Stokes numbers (heavy particles) although decreasing as predicted by Wang and Stock's theory, even if the decrease is induced by the overshooting (crossing-trajectory) and inertia.But the continuity effect is not restored.

. Conclusion on EIM-DRW Monte-Carlo Methods
It is a matter of fact that Lagrangian numerical simulations based on EIM-Monte-Carlo processes cannot predict the dispersion of heavy particles with or without gravity.The shortcomings of the method are due to the fact that the heavy particle will interact with an eddy during a given time, the maximum of which is normally the eddy life (taken as the Lagrangian scale T L ).Major modifications were brought to Monte-Carlo approaches by Huang et al. [179], Graham [95-97], Launay et al. [176,177], among others, and some eddyparticle interaction models tried to better take into account the space-time characteristics of the effective fluid turbulence controlling the particle motion, to better simulate the so-called "true turbulence experienced by the particle'." The weakness of the EIM-DRW method is that the correlation functions and Lagrangian integral time scales are related to the way the time step or eddy lifetime is chosen [5,60].Furthermore, the eddy length and lifetime are not correlated, so that no turbulence spectrum is really simulated.In addition, the different models proposed are generally tested on well-known flow cases [18,19] and optimized by calibrating a few parameters.Several papers compared different EIMs [120,133] and showed that even for these test cases, some EIM approaches are deficient and the differences in the numerical results are tremendous.It follows that EIM-DRW models are certainly not appropriate to more complex gas-solid flows, shear flows, anisotropic flows, since even in case of HIST deficiencies have been observed.That is why scientists have focused more and more of their research on Markovian methods for two decades.

Markovian Methods
We consider here the models as defined in Section 3.2 and two kinds of simulations will be presented here:

•
The Markovian model developed by a two-step space-time approach (Section 3.2.2) with the following equations and the method proposed by Lu et al. [126][127][128][129][130]: where ξ ti are normalized Gaussian random numbers.

•
A more classical Markovian model with a Langevin equation based on the correlation functions proposed by Wang and Stock (Equation (A2)):

CRW-Lu Model
In the 1990's it seemed that a modified EIM influenced by a Markov chain model could bring some progress.That is why a space-time Langevin equation was explored.The idea was to couple Eulerian space correlations to Lagrangian time correlations; the idea was not so bad and was explored with some hope and used to simulate grid turbulence experiments [15,16].The method as proposed by Lu et al. [126][127][128][129][130] was applied to HIST, but gave poor results.Without gravity and for small particles, which behave like small inertial particles, the particle Lagrangian time was 15% less than the Lagrangian fluid integral time T L .All possible parameters, among them the time increment ∆t, were changed and gave no better results.Presently, simulations with this model are still being undertaken.The particle variance seems to yield results similar to an EIM based on T L .It is clear that the effect on underestimating the particle dispersion is directly related to the integral time scales, and the space correction is certainly not valid, even if the approach gives satisfactory results for grid turbulence experiments, due to calibration of constants related to turbulent kinetic energy and dissipation.As for many models it is now clear that just simulating these experiments and modifying them for anisotropic flows is not the best choice.

•
Dispersion without gravity (g = 0) As shown in Figure 12, the normalized variances of the fluctuating velocity of the particles are generally underestimated, except for the smaller Stokes numbers (less inertial particles).The same holds for the Lagrangian integral time scales (Figure 13) as well as for the normalized dispersion coefficients (Figure 14).Thus, the model based on Markovian chains is not suitable since the inertial effect (IE) of the particles cannot be correctly simulated.

•
RWM-Lu with gravity effects (g = 0) Fluids 2021, 6, x FOR PEER REVIEW 31 of 53 • Dispersion without gravity (g = 0) As shown in Figure 12, the normalized variances of the fluctuating velocity of the particles are generally underestimated, except for the smaller Stokes numbers (less inertial particles).The same holds for the Lagrangian integral time scales (Figure 13) as well as for the normalized dispersion coefficients (Figure 14).Thus, the model based on Markovian chains is not suitable since the inertial effect (IE) of the particles cannot be correctly simulated.As shown in Figure 12, the normalized variances of the fluctuating velocity of the particles are generally underestimated, except for the smaller Stokes numbers (less inertial particles).The same holds for the Lagrangian integral time scales (Figure 13) as well as for the normalized dispersion coefficients (Figure 14).Thus, the model based on Markovian chains is not suitable since the inertial effect (IE) of the particles cannot be correctly simulated.As shown in Figure 12, the normalized variances of the fluctuating velocity of the particles are generally underestimated, except for the smaller Stokes numbers (less inertial particles).The same holds for the Lagrangian integral time scales (Figure 13) as well as for the normalized dispersion coefficients (Figure 14).Thus, the model based on Markovian chains is not suitable since the inertial effect (IE) of the particles cannot be correctly simulated.The turbulent energy of the heavy particles is well predicted by the model except for intermediate Stokes numbers, where it is somewhat underestimated (Figure 15).One can note that overall, this model yields better results for the particle turbulent energy than those obtained by an EIM-Monte-Carlo model (Figure 11).Again, long time dispersion coefficients are systematically much smaller (up to a factor of 2) than Wang and Stock theory for the smaller Stokes numbers (Figure 16).They decrease with larger particle inertia and for larger Stokes numbers, the ratio of transverse to gravity direction dispersion coefficients D p,11 /D p,33 is close to a value of 0.5 in accordance with Csanady's early results [42].Such a Markovian model seems to be able, at least partially, to simulate effects of crossing trajectories and continuity, but the estimates of the dispersion coefficients are significantly small.

•
RWM-Lu with gravity effects (g ≠ 0) The turbulent energy of the heavy particles is well predicted by the model except for intermediate Stokes numbers, where it is somewhat underestimated (Figure 15).One can note that overall, this model yields better results for the particle turbulent energy than those obtained by an EIM-Monte-Carlo model (Figure 11).Again, long time dispersion coefficients are systematically much smaller (up to a factor of 2) than Wang and Stock theory for the smaller Stokes numbers (Figure 16).They decrease with larger particle inertia and for larger Stokes numbers, the ratio of transverse to gravity direction dispersion coefficients is close to a value of 0.5 in accordance with Csanady's early results [42].Such a Markovian model seems to be able, at least partially, to simulate effects of crossing trajectories and continuity, but the estimates of the dispersion coefficients are significantly small.

Langevin-WS Model
The Langevin model based on the Wang and Stock correlation functions seems more promising.Without gravity, this model gave results close to the theory of Wang and Stock, in accordance with results presented for the EIM in Section 4.1.2.We present here only the simulations with the gravity-CTE-CE effects (Figures 17-19).The following equation is used with correlation functions given by Equation (A2): • RWM-Lu with gravity effects (g ≠ 0) The turbulent energy of the heavy particles is well predicted by the model except for intermediate Stokes numbers, where it is somewhat underestimated (Figure 15).One can note that overall, this model yields better results for the particle turbulent energy than those obtained by an EIM-Monte-Carlo model (Figure 11).Again, long time dispersion coefficients are systematically much smaller (up to a factor of 2) than Wang and Stock theory for the smaller Stokes numbers (Figure 16).They decrease with larger particle inertia and for larger Stokes numbers, the ratio of transverse to gravity direction dispersion coefficients is close to a value of 0.5 in accordance with Csanady's early results [42].Such a Markovian model seems to be able, at least partially, to simulate effects of crossing trajectories and continuity, but the estimates of the dispersion coefficients are significantly small.

Langevin-WS Model
The Langevin model based on the Wang and Stock correlation functions seems more promising.Without gravity, this model gave results close to the theory of Wang and Stock, in accordance with results presented for the EIM in Section 4.1.2.We present here only the simulations with the gravity-CTE-CE effects (Figures 17-19).The following equation is used with correlation functions given by Equation (A2):

Langevin-WS Model
The Langevin model based on the Wang and Stock correlation functions seems more promising.Without gravity, this model gave results close to the theory of Wang and Stock, in accordance with results presented for the EIM in Section 4.1.2.We present here only the simulations with the gravity-CTE-CE effects (Figures 17-19).The following equation is used with correlation functions given by Equation (A2): It is clear that this model is efficient, but in other models as well as the last, which will be presented in Section 4.3, the support is conditioned to an appropriate relationship of the correlation functions to the Stokes number and the drift parameter, in this case the Wang and Stock model detailed in Appendix A. The influence of the time step was not investigated but was always less than 1/10 of the typical time scales ∆t < min (T L /10, τ s /10).It is clear that this model is efficient, but in other models as well as the last, which will be presented in Section 4.3, the support is conditioned to an appropriate relationship of the correlation functions to the Stokes number and the drift parameter, in this case the Wang and Stock model detailed in Appendix A. The influence of the time step was not investigated but was always less than 1/10 of the typical time scales     It is clear that this model is efficient, but in other models as well as the last, which will be presented in Section 4.3, the support is conditioned to an appropriate relationship of the correlation functions to the Stokes number and the drift parameter, in this case the Wang and Stock model detailed in Appendix A. The influence of the time step was not investigated but was always less than 1/10 of the typical time scales     It is clear that this model is efficient, but in other models as well as the last, which will be presented in Section 4.3, the support is conditioned to an appropriate relationship of the correlation functions to the Stokes number and the drift parameter, in this case the Wang and Stock model detailed in Appendix A. The influence of the time step was not investigated but was always less than 1/10 of the typical time scales

Conclusions on Markovian Models
The prediction of the dispersion of particles, based on Markov one-chain Lagrangian models, seems encouraging.Two of the three main effects, crossing-trajectory and continuity, are rather well reproduced, except for the inertial effect, which enables particles under zero gravity to disperse more than fluid elements, a result obtained theoretically by Reeks [40], for instance.Furthermore, the long time dispersion coefficients are underpredicted both with and without gravity effects.Nevertheless, Huilier [120], on the basis of Launay's thesis [159] and complementary simulations, clearly showed in a previous paper the important influence of various autocorrelation functions R P F,ii (∆t, ∆s) [20,23,24,120] on the dispersion coefficients of heavy particles and one-step Markovian simulations with adequate autocorrelations can bring enhanced results.

Matrix Method and Related Considerations
The matrix method proposed by Launay has no real obvious gap, it is a statistical method based on any correlation function and therefore simulates the fluid velocity along the path of a discrete particle if the function is known.The following results are of course in agreement with the Wang and Stock theory.The main problem persists with large Stokes numbers because it is not obvious that the Wang and Stock theory is guaranteed by such large Stokes numbers.In order to check that the proposed Lagrangian method also correctly restores the inertia effect (increased diffusivity of high-inertia particles), the dispersion of particles without any drift velocity was first studied.In this case the parameter f was taken as one in the Equation ( 9), id linear drag.The matrix model is self-consistent, in accordance with the model proposed by Wang and Stock, for the particle Lagrangian integral time scales (Figure 20) and the particle turbulence (Figure 21).Again, as in Section 4.1.2with the modified EIM, previous theoretical works (Reeks [40], Pismen and Nir [57]), experiments (Wells and Stock [16]) and numerical studies (Squires and Eaton [98], Deutsch [178], Graham [96,97]) are confirmed for heavy particles with large inertia that disperse faster than the fluid elements.Figure 22 shows that the ratio of the particle to the fluid dispersion (u 2 0 T L ) coefficients is greater than unity for large Stokes numbers: this indicates that the inertia effect is well taken into account by the proposed method for predicting the turbulent dispersion of the heavy particles.The drop-off for Stokes numbers larger than 3 (last point in Figure 22) was never explained up to now; in fact, referring to Launay [159], the particle integral time scale is under-predicted, although the particle turbulence is in accordance with the theory of Wang and Stock.Either a non-linear drag for these very inertial particles should be considered; but calculations based on a Monte-Carlo and a non-linear drag law also underestimated the theoretical values for Stokes numbers larger than 3 as shown by Matt and Huilier [180].Trajectory calculations over a longer path time (typically T * ∼ 20 • τ s ) should then give better results.tion 3) as well as in the two directions orthogonal to the free-fall direction (directions 1 and 2).The agreement between the numerical simulations and the theoretical expressions is very good and the effect of crossing trajectories and continuity effects is clearly seen.

Matrix Method with Gravity (g = 0)
Figures 23-25 yield the numerical results obtained with the proposed Lagrangian model.The whole of the numerical results is again compared with the theory developed by Wang and Stock.In Figure 25, the analytical and calculated long time dispersion particle's coefficients D i p (∞) decrease with the Stokes number St in free-fall direction (direction 3) as well as in the two directions orthogonal to the free-fall direction (directions 1 and 2).The agreement between the numerical simulations and the theoretical expressions is very good and the effect of crossing trajectories and continuity effects is clearly seen.Besides, it is shown in Figure 25 that the particle's dispersion in the free-fall direction is larger than in the two transverse directions by a factor close to 2. Besides, it is shown in Figure 25 that the particle's dispersion in the free-fall direction is larger than in the two transverse directions by a factor close to 2.   Besides, it is shown in Figure 25 that the particle's dispersion in the free-fall direction is larger than in the two transverse directions by a factor close to 2.

Discussion and Present Trends of Modeling Approaches
The present study only gave an overview of the relative performance of Lagrangian particle tracking methods for the turbulent gas-particle transport in HIST.The EIM was historically extensively summarized, the RWM partly.The main advantage of these methods is that they help to better understand the mechanisms and concepts of inertia and

Discussion and Present Trends of Modeling Approaches
The present study only gave an overview of the relative performance of Lagrangian particle tracking methods for the turbulent gas-particle transport in HIST.The EIM was historically extensively summarized, the RWM partly.The main advantage of these methods is that they help to better understand the mechanisms and concepts of inertia and crossingtrajectory effects.Although presenting some other advantages, a crucial disadvantage is that these methods require many trajectory calculations to achieve statistical convergence.
The EIM, sometimes also called ELT (Eddy Life Time), is often coupled to k-eps or RANS turbulence models, is easy to implement and gives a first idea of the influence of inertia and gravity on the behavior of discrete particles in interaction with turbulent structures.But it seems, as demonstrated, that EIM hardly can predict accurate behavior of heavy particles in HIST and certainly not in complex flow configurations, although these methods were extensively applied to polydisperse sprays in combustion.
Markovian methods using time and space correlations are still discontinuous random walk methods, and they have been developed as a possible improvement of the EIM methods.Mainly the experiments by Snyder and Lumley [15] and Wells and Stock [16] served as test cases to optimize these models and calibrate some constants.As shown in Section 4.2, they were not able to simulate correctly turbulent dispersion in HIST, too.
The meteorology community working on air pollution and turbulent diffusion in the atmosphere certainly brought a decisive contribution.A major effort was devoted to optimize DRW and CRW models in inhomogeneous turbulence with the help of the Langevin equation.Initially the Langevin equation was used to model the Brownian velocity fluctuations.The stochastic equation was then applied to homogeneous turbulence by Obukhov [181] in terms of a Fokker-Plank equation related to a Gaussian probability distribution [182,183].The CRW-Langevin equation first served to model turbulent dispersion in the tracer limit, when particles are considered as passive contaminants or fluid elements.In homogeneous stationary one-dimensional Gaussian turbulence the Langevin equation reduces to equation When applied to inhomogeneous flow, such as in boundary or shear layers, the classical models, among which is the original Langevin equation, predict that particles migrate from regions of high turbulence intensity (core or bulk region) to the regions where the intensity is minimum (wall) and those even in the tracer limit.This is due to the sampling of fluctuation velocity, inducing a spurious drift component of the mean flow and introducing an unphysical accumulation and deposition velocity into the wall regions.
Early Wilson et al. [184,185] proposed a correction related to the effect of turbulent acceleration of the particles for the vertical velocity component w (t), in terms of a vertical bias velocity w (z) = σ z T L ∂σ z /∂z added to the RHS of the Langevin equation, σ z being the root-mean square of the vertical velocity.Several papers [186,187] refer to a problem of consistency between the classical Langevin equation and the equations of the mean flow field when particles are so small that they behave like fluid eddies.Classical models do not obey the well-mixed criterion as defined by Thomson [7], that is: "A tracer field which is initially well mixed in an inhomogeneous turbulent field should remain well mixed".
The spurious drift term is related to the fluid acceleration by averaging The second term is non-zero in case of inhomogeneous turbulence, and the drift correction velocity is then given by ∆u i = u j ∂u i ∂x j ∆t in the tracer limit.In case of inertial particles, Bocksell and Loth [149] proposed a spurious drift correction given by so that the classical Langevin equation with spurious drift correction reads: an equation that meets the well-mixed condition and can be reasonably used except in strongly inhomogeneous turbulent flows such as boundary layers.More generally, a generalized Langevin equation model, as proposed by Haworth and Pope [188], was developed and is in active use.Pope [189] initially established PDF methods for reactive flows.He also fixed the consistency conditions [190].Several simplified versions exist [191] and one simplified equation is .
From a scientific point of view, RANS methods have been progressively replaced by more accurate methods (LES, DNS), even if the computational cost is higher.With the progress in computer performance, the Eulerian-Lagrangian models and particle tracking methods have focused progressively on a large eddy simulation (LES), coupled to stochastic particle subgrid scale (SGS) models, proposing Langevin equations for the fluid seen along the heavy particle trajectory in order to solve the momentum equation for particle-laden flows .
One of the more detailed particle SGS models is that proposed by Berrouk et al. [209][210][211], based on the Wang and Stock equations (Appendix A), integrating inertia, CTE, continuity and anisotropy; the filtered fluid field (p * , u * f i ) is calculated by an LES closed by a classical fluid SGS model, and the fluid velocity u pi seen by the particle is obtained by a Langevin equation solved by a first-order Euler scheme; the Langevin equation is given by Continuity, CTE and inertia appear in the drift term.T * i (T L , β, St) is related to the Lagrangian time scale of the subgrid fluctuations with the Wang and Stock correction for inertia and Csanady's expressions for mean drift and transverse directions in case of gravity.More details are given in Berrouk et al. [210] (pp.8-10) for the closure of the drift and diffusion terms, the last being modified for anisotropy and non-stationarity through a diagonal, non-isotropic diffusion matrix C * 0,i .Both parameters (T * i , C * 0,i ) depend on the ratio between the Lagrangian and Eulerian time scales or β = T L /T mE , which was kept constant in all the studies proposed by Wang and Stock [18][19][20] and equal to 0.356, referring to a turbulence structure parameter m equal to 1 (Appendix A).The model was used to simulate the dispersion experiments of Arnason in a pipe flow [75,76].Nevertheless, as mentioned, the normalized time scale β or the Kolmogorov constant C 0 varies considerably as given in the literature [9,10].Last, a main persisting problem is related to the choice of the value of the turbulence structure parameter m = T mE u 0 /L f .Several benchmark tests exist and LES results are sometimes compared to DNS models, which normally provide exact numerical flow fields, to evaluate the predictive capability of particle dispersion in more complex turbulent flows.For instance, an international benchmark [227] was proposed for LES Eulerian-Lagrangian simulations in a particle-laden channel flow, gathering a large database, and was compared to DNS results.Marchioli [228] recently presented a review of LES in turbulent dispersed flows, which was complementary to other reviews and guidelines [157,197].Objectives at least should be able to model physical mechanisms such as erosion or deposition, settling, saltation of aerosols or solid particles in bounded flows, migration in shear layers and boundary layers, preferential concentration and turbophoresis.Particles tend to segregate in decreasing turbulence level region wall areas, and the rate of deposition is enhanced on surfaces.Caporaloni et al. [229] and Reeks [230] independently called this phenomenon "Turbophoresis," analogous to the thermophoresis or Ludwig-Soret effects.Furthermore, dense particles tend to accumulate and cluster in specific convergence regions and are centrifuged away from vortex zones.The influence of coherent vortex structures on particle dispersion was recognized early by Crowe et al. [231][232][233], Chein and Chung [234], and studied by Eaton and coworkers [235][236][237][238][239] as well as the research group of Cartellier or Bourgoin [240][241][242][243][244][245][246][247], among others.Preferential concentration has since been extensively reported in the literature during the last two decades, and the physical clustering mechanisms have been identified, namely, sweep-stick, heavy particles migrating to low acceleration regions [248,249] and centrifugal expulsion effects, inertial particles accumulating in low-vorticity zones of the carrier flow.
To conclude, other review papers on modeling particle-laden flows or multiphase flows are of interest [250][251][252][253][254][255][256][257], including, in some cases, collision and turbulence modification by particles in dense gas-particle flows.Some general and complete reviews should be mentioned, too; the relevant review by Balachandar and Eaton [258] on different aspects of turbulent dispersed multiphase flows focusing on preferential concentration and different coupling mechanisms (e.g., interphase, two-way momentum coupling, turbulence modulation), and a general review of Lagrangian statistical properties of particles by Toshi and Bodenschatz [259].
Furthermore, aside from the interest brought in the literature to single particle Lagrangian structure functions, with comparison between experiments and DNS instance, e.g., [260], turbulent diffusion and dispersion have also been investigated in terms of the mean-square displacement of a pair of particles relative to their center of mass and relative velocity between pairs.Relative or particle pairs dispersion has been, since the pioneering paper by Richardson [261] and Batchelor [262], extensively studied theoretically, experimentally and numerically, in relation to new stochastic models and DNS calculations [263] or experiments [264,265].Significant reviews on turbulent relative dispersion were those of Sawford [266] on two-particle Lagrangian stochastic models, Salazar and Collins [267] on more recent advances.The reviews are often related to tracer diffusion and HIST, and offer an increased understanding of the Lagrangian properties of turbulent flows.But they as well as other articles not only concern low-inertia particle pairs [268] but also separation of high-inertia particle pairs [269,270] and preferential concentration [271].A new review on turbulent relative dispersion studies would be welcome.(A3) Following a development proposed by Reeks [40], a relation between the particle velocity correlations R P ij (τ) and the fluid velocity correlations R P F,ij (τ) (Relation 2.9, in [20] (p.1900)) is obtained and for long time dispersion, algebraic expressions for the variances of the particle's fluctuating velocity, the particle's integral time scales and the particle's long time dispersion coefficients can be written as:

•
variances of the particle turbulent velocity in the i-th direction: As the previously expressions indicate, the heavy particles' dispersion depends on the inertia parameter St (characterizing the inertia's particle for a given turbulence), the drift velocity parameter γ and the turbulence structure parameter m.
When the heavy particle has a zero drift velocity (absence of external force), the previously analytical expressions take simplified forms: (A10) In the case where the drag law is nonlinear, the above analytical expressions remain valid insofar as the particle Stokesian relaxation time τ s is replaced by the "effective" particle relaxation time: where f is the nonlinear correction factor as proposed by Schiller and Naumann [37] (see Section 2.2.The lateral autocorrelation coefficients (directions perpendicular to the fall direction) can have negative loops, which are more pronounced for larger drift parameters.The lateral autocorrelation coefficients (directions perpendicular to the fall direction) can have negative loops, which are more pronounced for larger drift parameters.

Fluids 2021, 6
are the position and velocity of the particle, respectively, ), t) is the undisturbed fluid velocity at the location of the particle.m f denotes the mass of the fluid displaced by the particle and µ F the dynamic viscosity.The derivative D Dt denotes the time derivative following a fluid element, D time derivative following the moving sphere,

Fluids
less be correlated during a process time EIM T depending on two situations:• the particle (P) leaves the fluid eddy (F) for entering into a new one when the eddy (F) lifetime τ * is elapsed (condition 1),

Figure 1 .
Figure 1.Sketch of the eddy interaction model.Figure 1. Sketch of the eddy interaction model.

Figure 1 .
Figure 1.Sketch of the eddy interaction model.Figure 1. Sketch of the eddy interaction model.
f L , the random fluid velocity is kept constant.If y < ∆t/T f L , a new fluid velocity u f (t) is generated, according to a Gaussian probability-density function.Then a Frenkiel space correlation is used and a Fluids 2021, 6, 145 16 of 48 Deutsch and Simonin[99] proposed the following formulas for the vertical and transverse components, , including CTE and CE but no inertia effect, since T f L is not the fluid Lagrangian time along the discrete particle path without gravity.

Figure 2 .
Figure 2. Construction principle of an inertial particle trajectory in the random walk model (RWM).

Figure 3 .
Figure 3. Matrices A & B used to calculate the turbulent fluid velocity seen by the particle along its path.

Figure 11 .
Figure 11.Normalized variance of particle velocity (g ≠ 0).4.1.4.Conclusion on EIM-DRW Monte-Carlo Methods It is a matter of fact that Lagrangian numerical simulations based on EIM-Monte-Carlo processes cannot predict the dispersion of heavy particles with or without gravity.

Figure 22 .
Figure 22.Particle dispersion coefficient (g = 0).4.3.2.Matrix Method with Gravity (g ≠ 0) Figures 23-25 yield the numerical results obtained with the proposed Lagrangian model.The whole of the numerical results is again compared with the theory developed by Wang and Stock.In Figure 25, the analytical and calculated long time dispersion particle's coefficients ) (∞ i p D

Funding:
This research received no external funding.coefficients have the advantage of integrating both the inertia and the drift velocity of the heavy particle.The fluid Lagrangian integral scales seen by the heavy particle having a drift velocity are then for long time diffusion:

••
long time coefficient of the heavy particle in the i-th direction:D P,11 (∞) = D P,22 (∞) = u 2 particle Lagrangian integral time scale in the i-th direction:T P,11 = T P,22 =
2): f = 1 + 0.15 (Re P ) the heavy particle's diameter, ρ F and µ F are, respectively, the density and the dynamic viscosity of the fluid.In FigureA1are shown the Wang and Stock fluid correlation coefficients for different values of the particle Stokesian relaxation time.These coefficients are all the more so directional dependent as the relaxation time τ s is large.
on the Kolmogorov microscale, • St mE = τ p /T mE based on the Eulerian moving macroscale, • St E = τ p /T E based on the classical Eulerian macroscale, • St L = τ p /T L based on the Lagrangian time scale, • St MF = τ p /T MF based on a time scale of the mean flow T MF = L/U indicating if the particles follow the mean fluid flow.

Table 1 .
Classification of dispersion regimes.

Table 2 .
Particle properties and dimensionless control parameters.