LES and RANS calculations of particle dispersion behind a wall-mounted cubic obstacle

In the present paper, we evaluate the performances of three stochastic models for particle dispersion in the case of a three-dimensional turbulent flow. We consider the flow in a channel with a cubic wall-mounted obstacle, and perform large-eddy simulations (LESs) including passive particles injected behind the obstacle, for cases of low and strong inertial effects. We also perform Reynolds-averaged simulations of the same case, using standard turbulence models, and employ the two discrete stochastic models for particle dispersion implemented in the open-source code OpenFOAM and the continuous Lagrangian stochastic model proposed by Minier et al. (2004). The Lagrangian model is consistent with a Probability Density Function (PDF) model of the exact particle equations, and is based on the modelling of the fluid velocity seen by particles. This approach allows a consistent formulation which eliminates the spurious drifts flawing discrete models and to have the drag force in a closed form. The LES results are used as reference data both for the fluid RANS simulations and particle simulations with dispersion models. The present test case allows to evaluate the performance of dispersion models in highly non-homogeneous flow, and it used in this context for the first time. The continuous stochastic model generally shows a better agreement with the LES than the discrete stochastic models, in particular in the case of particles with higher inertia.


Introduction
Bluff-body flows are of primary importance because present in numerous applications in natural and engineering sciences.In particular, the flow around a surface-mounted blunt obstacle placed in a channel is fundamental to the understanding of the flow in complex two-and three-dimensional geometries, often in relation with the movement of cars, buses, trains or trucks (Simpson, 2001).
A canonical example of such flows is given by a cube with sharp edges mounted over a plate in a channel.This flow has been thoroughly investigated since the pioneering experimental work by Martinuzzi and Tropea (1993), and has become a classical benchmark for separated flows both for direct numerical simulation (DNS) and different models (Rodi, 1997;Iaccarino et al., 2003).
The progress in computational tools together with the practical importance of such flows has motivated several DNS studies of this configuration (Yakhot et al., 1992;Diaz-Daniel et al., 2017) or similar bluff-body flows (Bruno et al., 2014;Cimarelli et al., 2018).However, while DNS remain important to get physical insights, only Reynolds averaged (RANS) models can be routinely used in practice, and largeeddy simulations (LES) have to be further developed as the reference numerical tool for high-Reynolds-number flows.Hence, to develop and assess RANS and LES models remain a crucial research activity for these complex fluid flows.
In many cases of interest for environmental studies and engineering applications, the flow is laden with particles or droplets, as for instance when a car encounters rain or snow.Even considering only the case of small particles heavier than the carrier fluid, these systems are much more complex than their single-phase counterpart and the motivation for the development of reliable models is even stronger because of the complexity of performing relevant experiments (e.g. in cases with snow).In this context, the field of simulating and modelling particleladen flows is yet in its infancy, even in very dilute conditions when the modifications of the turbulence by the solid phase can be neglected and one is interested in the particle dynamics, notably clustering and preferential sampling.https://doi.org/10.1016/j.ijmultiphaseflow.2022.104037Received 22 September 2021; Received in revised form 3 January 2022; Accepted 23 February 2022 Important advances in understanding the basic physics of such flows have been acquired in simple configurations (Balachandar and Eaton, 2010), and combined progress in computational approaches have been made, as recently reviewed in Fox (2012), Elghobashi (2019) and Brandt and Coletti (2022).Nevertheless, as regards to modelling particle-laden complex flows, notably when the geometry is complex and/or recirculation is important, the state-of-the-art is still to use empirical approaches often flawed by errors (Fox, 2014;Minier et al., 2014).The development of accurate engineering models for complex particle-laden flows is therefore one urgent task, to which the present research aims at contributing.
Generally speaking, engineering models for turbulent disperse flows can be developed both in an Eulerian (Simonin, 1996) or Lagrangian framework (Minier and Peirano, 2001;Peirano et al., 2006).In analogy with classical RANS models of single-phase flows, Eulerian models represent the two phases as two fluids.The advantage is the possibility to exploit existing Reynolds stress models for single-phase turbulence.Yet, due to the specific nonlinear character of the particle equations, this approach has to overcome a difficult closure problem, and hence it may be considered relevant only for dense flows (Marchisio and Fox, 2013;Fox, 2014).
The Lagrangian approach appears thus necessary for dilute collisionless conditions, and indeed it has received attention for many years (Stock, 1996).In particular, researchers have tried to develop a rigorous Lagrangian probability density function (PDF) formulation (Minier and Peirano, 2001;Peirano et al., 2006), similarly to the PDF approach to reactive flows (Pope, 2000).This approach consists in modelling the relevant one-point PDF as a diffusion process and solving the corresponding set of stochastic equations, which allows to correct the drawbacks present in previous empirical-based models (Minier et al., 2014).The PDF approach is nowadays well established in reactive flows even in complex geometries (Pope, 2000).It can be considered as the standard approach for dispersion in homogeneous media, like atmospheric flows (Stohl et al., 2005).
Modelling of Particle-laden non-homogeneous flows appears more difficult and therefore it is less developed.Heuristic Lagrangian models have been applied to many free and bounded flows showing yet a modest capability of prediction, notably as to deposition in wallbounded flows (Henry et al., 2012).The PDF approach has shown good performances in several homogeneous and non-homogeneous free-shear flows (Peirano et al., 2006), and most notably when applied to a benchmark bluff-body case (Minier et al., 2004).The Lagrangian PDF approach has been recently extended also to moderately dense homogeneous flows (Innocenti et al., 2019), where it can be a useful complement to Eulerian two-fluid models.
Wall-bounded flows usually ask for a special treatment of the boundary layer and can be considered less developed, even though progresses have been recently achieved.Notably, good results have been obtained in the LES framework (Innocenti et al., 2016), and a first attempt to model moderately dense bounded flows has been made (Innocenti et al., 2021).Concluding this brief description of the state-of-the-art, we can conclude saying that the assessment of present numerical models of particle-laden flows is a crucial ongoing research field (Berrouk et al., 2007;Pozorski and Apte, 2009;Michałek et al., 2013;Alletto and Breuer, 2012;Breuer and Hoppe, 2017;Salehi et al., 2017).Indeed, while the Lagrangian approach is a mature and reliable approach for particle-laden flows in homogeneous and simple non-homogeneous flows, it remains to be analysed and developed to be routinely used in more complex flows.
From a technical point of view, while in principle a fully Lagrangian fluid/particle approach is possible (Minier and Peirano, 2001), in practice a Hybrid Eulerian/Lagrangian approach is always used, as in the present work.This means solving the fluid phase in an Eulerian framework, which may be LES or RANS, and to couple with it a Lagrangian PDF approach for the particles (Peirano et al., 2006).The hybrid approach has many advantages from the computational side and is the only one possible for complex flows.Yet, specific issues arise in this framework, as already pointed out in similar efforts for reactive flows (Muradoglu et al., 2001), which need careful attention and make the benchmark of particle models subtle and difficult.
The main aim of this work is precisely to validate different Lagrangian models on a relevant particle-laden test-case which has not yet been investigated and may represent a standard benchmark for the future improvement of the modelling.
Specifically, we present and discuss LES and RANS results for the cube-mounted benchmark case in the presence of particles.This testcase appears well-suited to represent a challenging yet relevant testcase to evaluate the predictions of different Lagrangian dispersion models.In particular, being the first time this benchmark case is used in presence of particles, LES and RANS are carried out first in the singlephase configuration to be compared with the literature, and then the particle-laden case will be analysed using different fluid and particle models.A remark is in order to clarify the context of the present work.All Lagrangian models for particles are implicitly based on a model for the fluid phase.Furthermore, since particles are coupled to the fluid, fluid statistical moments are used in the particle model in a hybrid approach.Both facts point out to an often overlooked consistency issue, which has been addressed both in single reactive and particle-laden flows (Pope, 1987;Muradoglu et al., 2001;Chibbaro and Minier, 2011;Chibbaro et al., 2014).The consequence is that when studying a particle-laden flow, the model used to tackle the fluid flow should be consistent with that used in the particle model, in the sense indicated by Chibbaro et al. (2014).Otherwise, uncontrolled errors and nonphysical results may arise, most notably in the RANS framework.This also implies that particle models cannot be assessed singularly, independently from the fluid phase models.In particular, it has been shown that the naive idea to couple DNS with standard Lagrangian models, as to assess them directly, leads to huge bias errors (Chibbaro and Minier, 2011).In the present work, we try a first extensive benchmark using several different models, and we hope that these results could be used in the future as a reference to improve the quality of the different models.
The paper is organized as follows: in Section 2, we describe the test case that we considered; in Section 3, we introduce the governing equations of the system and the particle dispersion models; in Section 4, we describe our simulations, separately for the fluid and the particle phases, and we discuss our results; and in Section 5, we summarize our findings and give perspectives on further developments.

Test case
We focus on the case of a cubic obstacle mounted on the lower wall of a channel.This case has been investigated in details by Martinuzzi and Tropea (1993), who employed different experimental techniques to obtain information about the mean fields as well as the second order moments at several positions around the obstacle.The experimental setup employed by Martinuzzi and Tropea (1993) consists of a channel of 390 cm × 60 cm × 5 cm (length × width × height), with a cubic obstacle of side 2.5 cm mounted 52 channel heights far from the inlet.We take the cube side as characteristic length, denoted by   , and the mean velocity at the inlet, denoted by   , as reference velocity.The flow reaches a fully developed state far before the cube.The natural definition of the Reynolds number becomes then: We study the case with  = 50,000, for which high quality experimental data are available online. 1 The same experimental set-up has already been chosen as reference for different benchmark studies (Rodi et al., 1997;Rodi, 1997;Muld et al., 2012).Note that we will present  In our reference system the  axis correspond to the main flow direction, the  axis to the vertical direction and the  axis to the span-wise direction.The domain can be described via three parameters: the length of the channel upstream of the cube   ; the length of the channel downstream the cube   ; the width of the channel   .Due to the fact that the cubic obstacle occupies half of the channel height, the full height is included in the computational domain, while the streamwise and spanwise dimensions of the computational domain may be smaller than the experimental channel to reduce the computational costs.The parameters used here are reported in Table 1 in units of   .As reported in the table, the domain used in present computations is clearly reduced with respect to experiments, and it is in line with previous computations, where it was shown that the choice of such small upstream length was sufficient to get reasonably accurate results for the fluid statistics (Muld et al., 2012).
To give a visual illustration of the most prominent features of the chosen flow configuration, we display vortical structures at an arbitrary time step of our simulation in Fig. 1.It is possible to appreciate the complexity of the turbulent structures in the wake of the obstacles, as well as in the separation regions caused by the sharp edges of the geometry.

Governing equations
In the present work, we consider dilute incompressible gas-particle flows, so that the modulation of the fluid by the particles and particleparticle interactions are neglected.The fluid flow obeys the continuity and Navier-Stokes equations, For the disperse phase, we consider heavy particles (  ≫   ); the equations of motion for such particles can be written as (Gatignol, 1983;Maxey and Riley, 1983): where   = (  (), ) is the fluid velocity at the particle position, i.e. the fluid velocity sampled along the particle trajectory   ().To simplify the analysis and the comparison among different cases, gravity has been neglected.The particle relaxation time is defined as: where the local instantaneous relative velocity is   =   −   and the drag coefficient   is a non-linear function of the particle-based Reynolds number,   =   |  |∕  (Clift et al., 2005).A commonly used empirical relation for the drag coefficient is: (5)

Numerical approach for the fluid phase
The fluid phase has been investigated via LES and RANS.Since the main novelty of the present work resides in the particle statistics, we have simulated the fluid phase with standard models using the standard values for the different model coefficients and parameters.For this reason, we shall not give details here, and we refer to the literature for further details (Pope, 2000).
Specifically, we have performed several LES of the fluid phase (Pope, 2004;Ferziger and Perić, 2002) using the Spallart-Almaras IDDES model (Shur et al., 2008) to estimate the turbulent viscosity .In these simulations, the governing equations and the transport equation for the model variable, denoted by ν, are solved for each cell of the domain and each time step.The length scale for the ν-production term is computed using the cell size far from the wall, the wall distance in the near-wall region, and a blending function in the intermediate regions.
Concerning the RANS, we have used the  − , the realizable  −  and the  −  models.All models are used in their standard implementation in the OpenFOAM package.

Particles tracking
Let us denote variables obtained via a generic numerical simulation with the superscript (⋅)  .The central question in performing particle tracking is how to evaluate   in Eq. (3).
For DNS, whose aim is to solve exactly all the scales of motion so that the uncertainties are due only to numerical errors, we can write for the definition of the fluid seen by particles: Then, to track particles in the simulation amounts to solving the particle equations, together with the fluid ones, using   in Eq. (3); in the DNS spirit, the particle trajectories should be considered as experimental data.
If closure models are employed, they a priori introduce errors related to the unsolved scales.In the following paragraphs, we shall explain the different strategies used in the present work.
When we perform a LES, we have access to the filtered velocity Ũ () = ∫ (, )  () , where the function  indicates the filter.This a coarse-grained version of the full velocity field, yet fluctuating and gathering some detailed information about the large-scale motion, that is for those scales above the filter length.In this case, the first-order approximation is to simply use as for a DNS, and the resulting trajectories provide an approximated picture of the real ones.This approximation is known to lead to important errors (Bianco et al., 2012;Geurts and Kuerten, 2012;Chibbaro et al., 2014), which can be corrected with appropriate stochastic subgrid models, see Innocenti et al. (2016).However, it has been shown that in shear flows the error might be reasonably small with respect to the large-scale statistics (Armenio et al., 1999).Given the complexity of the flow studied in this work, and the lack of previous simulations with particles, in this work we discuss the LES results obtained without introducing a stochastic subgrid model, whose behaviour is left for future studies.
When the fluid phase is solved through a RANS approach, the Reynolds decomposition   =   +  ′  is used and only statistical averages are available.The first and simplest approach to particle tracking is to simply write: This procedure surely allows to close the problem but is a too rough approximation; one therefore needs to add a model for the effects of the turbulent fluctuations.We therefore resort to a statistical approach and formally write: To estimate   (  ()) using only the RANS fields leads to considerable difficulties (Minier and Peirano, 2001).Nevertheless, several methods have been presented to adapt RANS simulations to particle tracking, because of their low computational cost.Generally speaking, this approach has some validity in homogeneous flow (Stock, 1996;Wilson and Sawford, 1996).
The first attempts were based on the assumption that the particle motion can be described via a homogeneous diffusion process, assumption which is justified only if   is Gaussian.In this simple approach, an effective diffusion coefficient for the Eulerian field that represents the particle concentration can be estimated, and   is evaluated as a random variable with Normal distribution  (0, ), and different values at each integration time steps (Dukowicz, 1980).The three components of   are independent and the standard deviation  is linked to the other variables of the turbulent model.An improvement came from the idea of taking directly into account a time-scale when modelling   (Gosman and Ioannides, 1983;Kabanovs et al., 2016).In this approach, the turbulence is still assumed to be isotropic, and the fluctuations are evaluated from the turbulent kinetic energy : The value of   used for any single particle is not selected at each time step, but after a turbulent time interval   .This choice has been inspired by the consideration that the velocity seen by the particle   changes significantly only if the particle moves from a coherent structure to another.Thus, the turbulent time interval is defined as:   = min(  ,   ), where   represents the time needed for the eddy that is carrying the particles to be dissipated; this is estimated as the dissipation length scale   over the magnitude of the randomly generated fluctuation   : and   represents the time needed for the particles to pass through the eddy; this is estimated from a linearized form of Eq. ( 3): In this work, we have used two standard models available in Open-FOAM.More specifically, the two models implemented are • the so-called Stochastic Dispersion model, which we rename Random Dispersion model to avoid confusion in the following, and is partly related to the Gosmann model.In this model, the turbulent time interval,   , is defined as: where   ≈ 0.16 is a model constant.After   elapses, the particle velocity is perturbed by adding a vector of random direction and magnitude generated as a Gaussian distributed random variable of variance: • the so-called Gradient Dispersion model, which adopts the same definition of   to determinate how often the fluctuation is seen by the particle and the same random generated amplitude with variance , but introduce a more refined system to determine the direction of the velocity perturbation: this is chosen parallel to the opposite of the gradient of the turbulent kinetic energy, i.e.: Such models may lead to large errors and nonphysical behaviour in strongly non-homogeneous flows.Namely, they are affected by spurious drifts (Thomson, 1987;Pope, 1987), which may be avoided using rigorous stochastic models (Pope, 1985;Minier et al., 2014).For this reason, we have also implemented and used in our simulations a Lagrangian stochastic model, which is detailed in the next section.

Diffusion stochastic model
The present Lagrangian stochastic model has been developed and extensively validated previously (Minier and Peirano, 2001;Minier et al., 2004;Peirano et al., 2006).We summarize here its main characteristics.
The variables retained in the model are the particle position   , the particle velocity   and the velocity of the fluid seen by the particles   .Therefore, the system will be described by the stochastic process  = (  ,   ,   ).In particular, the model is given directly in terms of the total instantaneous velocity   =   +   , as it has been shown to be a better practice (Minier et al., 2014).On physical analysis ground (Minier et al., 2004), it is a sound choice to model the stochastic process as a generalised Langevin equation, that is to describe   as a diffusion process (Gardiner et al., 1985).The governing equations for the three variables of the stochastic process are as in (Minier et al., 2004): The equations show that the particle position and velocity equations are closed and unaltered with respect to the exact ones in Eq. ( 3).The velocity of the fluid seen by the particles contains a drift term (terms proportional to the differential ) and a diffusion term, which is related to the Wiener process .
For the sake of clarity, we provide below a brief account of the rational behind the model.The starting point is the Langevin equation model used in single-phase PDF modelling (Pope, 2000), since particles reduce to fluid tracers in absence of inertia.The model is hence built so to retrieve the classical simplified Langevin model in the zero-Stokes limit (Pope, 2000;Minier and Peirano, 2001), which assures the model to be free of spurious drift by construction.
When inertia is present, two main effects affect the dispersion of particles: (i) in all cases, the typical integral time scale of the fluid velocity seen by the particles becomes different from the Lagrangian fluid velocity scale.(ii) When the flow is non-homogeneous, particles may drift which induces the Crossing-Trajectory-effect (CTE), that is a significant decorrelation of the fluid velocity fluctuations along the particle path.The first effect is much smaller than the second and is generally neglected in non-homogeneous flows.The crossing-trajectory effect (CTE) is due to the mean shear and has been modelled by changing the timescales in the drift and diffusion terms following Csanady's analysis (Csanady, 1963;Minier and Peirano, 2001).Assuming, for the sake of simplicity, that the mean drift is aligned with the first coordinate axis (the general case is discussed in Minier and Peirano, 2001), the model relations for the timescales are, in the longitudinal direction: and in the transversal directions (axes labelled 2 and 3) In these equations  is the ratio of the Lagrangian and the Eulerian timescales of the fluid,  =   ∕  , and  *  represents the Lagrangian time scale in the absence of mean drifts but accounting for particle inertia.Since in the present work shear effects are dominant, changes induced by the particle inertia are neglected.We therefore assume that  *  =   , and  = 1 for the sake of simplicity.To make the stochastic process consistent, a different kinetic energy has been introduced in the diffusion matrix; with   =   ∕ , , this reads Finally, a mean gradient term is added to the drift model, representing the dispersion of the particles with respect to the fluid trajectory.In this model, there is only one free parameter: the constant  0 .Although this can be considered as a free parameter, useful to calibrate the model (Minier and Pozorski, 1999), we have chosen to fix it to  0 = 2.1, which is a standard value in literature (Pope, 2000).The model Eq. ( 16) constitutes a set of stochastic differential equations (SDEs) in the form of a Langevin Equation (or a diffusion process more rigorously), and therefore they are equivalent to a closed Fokker-Planck equation for the one-point pdf  (,   ,   ,   ) (Gardiner et al., 1985).For this reason, it is possible to refer to it as a PDF model (Pope, 1985;Minier and Peirano, 2001).It is worth noting that the diffusion coefficient is not a constant, but a nonlinear function of space.This means that in non-homogeneous conditions the PDF may be strongly non-Gaussian.(Pope, 2000;Minier and Peirano, 2001;Innocenti et al., 2020) We conclude the description of the dispersion models with few comments about their consistency with fluid models.The discrete models are built on a pure heuristic ground, and strictly speaking cannot be fully consistent with any fluid model.Hence the two standard models implemented in OpenFOAM as well as all other models of the same class of discrete dispersion models, are in principle affected by spurious drifts, that is conservation laws may not be satisfied.The Langevin model is instead fully consistent with the SLM fluid model, which in turn corresponds to the ROTTA Reynolds stress model (Pope, 2000), but numerical studies have shown that the present model can be considered consistent with the standard  −  class models for engineering purposes (Chibbaro and Minier, 2011).

LES computation: numerical setup
First, we show the LES results, which will be used as a reference for the RANS calculations.Concerning the boundary conditions, we use a fit of the experimental data with a second order polynomial function at the position which corresponds to the inlet of the computational domain.Furthermore, we resort to the synthetic turbulent generator ''LEMOS''2 to introduce fluctuations at the inflow, based on empirical parameters lengths and time scales.We have checked that using more complex conditions has no effect on the flow downstream the obstacle.For the inlet condition of the variable in the Spalart-Allmaras model, ν, we have studied different values in the range between ν =  and ν = 10.We have obtained similar results both for the mean values and for the root mean squares, so that we have then always used ν = .The standard wall function implemented in OpenFOAM, which uses the Spalding formula for the law of the wall (Spalding, 1961), limits the turbulent viscosity,   , in the proximity of the wall.Lastly, the boundary condition for the pressure at the inlet has been set to the so called zero-gradient condition, which consists in imposing zero normal gradient at the boundary.At the outlet, a standard outflow condition is imposed for the velocity, which consists of a zero-gradient condition for each velocity component in the absence of backflow, and a fixedvalue condition that set the velocity to 0 in case of backflow (backflow events at the outlet are virtually absent in these simulations).The zerogradient condition is also imposed for turbulent viscosity, and the value of the reference pressure is set to 0. The sides of the computational domain have free-slip conditions whereas no-slip conditions are applied at walls.The pressure field satisfies the zero-gradient condition on the walls.
To show the difficulties inherent to our flow configuration, we display in Fig. 1 a three-dimensional visualization of the vortical structures arising from the walls.Tubes of vorticity of different importance are created and the resulting flow is very highly disordered.It is worth emphasizing that appreciable vorticity can be found also relatively far from the obstacle, which underlines the importance of turbulence fluctuations in this complex bluff-body problem.
Given the complexity of the problem, an extensive grid convergence study has not been possible.Nevertheless, after preliminary analysis, we have focused our attention on two different grids, whose parameters are given in Table 2. To evaluate the mesh quality we employed three quality indicators: non-orthogonality, skewness, and maximum aspect ratio.The non-orthogonality is defined as the angle between the line passing through two adjacent cell centres, and the vector normal to the common face of the two cells.The skewness is defined as the distance between the line passing through two adjacent cell centres and the centre of the common face, normalized with the distance between the cell centres.The cell aspect ratio is defined as minimum between the ratio between the largest and smallest areas of the cell bounding box, and the result of the expression 1∕6 ⋅   ∕ 2∕3  (where   and   denote the surface of the cell bounding box and the cell volume, respectively).Hereafter, we denote the simulation performed on the finer grid as ''LES1'', and that performed on the coarser grid as ''LES2''.In both grids, the distance of the first point from the wall is between  + = 1 and 2, in the region of higher refinement in the proximity of the object, and up to  + = 5 in the region farther from the object where reaching higher accuracy is less relevant ( + is computed based on the mean velocity).Note that the grid LES1 has a higher degree of refinement around the cube, in the region  ∈ (−1, 2.5), than the coarse-grid LES2, while the rest of the domain is discretized similarly.The number of point on the obstacle edges in LES1 and LES2 is 160 and 40, respectively.A boundary-layer mesh is used in both cases to obtain the proper resolution in the wall-normal direction near the surfaces.In the region of turbulent flow around the obstacle and in its wake, the two highest levels of refinement in LES1 result in a value of turbulent viscosity approximately equal to the physical viscosity.In the same regions in LES2, which is designed as LES1 without the third and fourth refinement levels, we measure a mean turbulent viscosity approximately twice as high as the physical viscosity.Note that the turbulent viscosity, denoted by   , is estimated based on the model variable ν, as defined in Spalart and Allmaras (1992).Higher values of   are reached in the region of low refinement farther away from the obstacles, up to approximately ten times the physical viscosity.
A visualization of the results obtained with the two meshes is given in Fig. 2, which displays a snapshot of the  component of the instantaneous velocity around the cube, showing that both grids provide the same qualitative behaviour of the flow.All the statistical results we will show in the following have been obtained via averaging in time for 200 time units, starting from the instant  * = 25.We checked stationarity comparing qualitatively the flow history over 10 time units before and after  * = 25 at selected locations in the obstacle wake.We checked statistical convergence comparing the mean velocity averaged over 100 and 200 time units, which showed lower discrepancies that those between the coarse and fine grids.

LES: comparisons with experiments
Next, we quantitatively analyse the statistical results obtained with the LES approach and the two meshes against the available experimental data.To this end, we consider the mean streamwise velocity component   and the shear Reynolds stress component  ′  ′ , see Fig. 3. Results for the streamwise velocity fluctuations,  ′  ′ , do not add insights and are not shown for the sake of clarity.The locations both on top and behind the cube are chosen as in the literature (Rodi, 1997).The remarks in order are the following: (i) Globally, the results are in reasonable agreement at all locations for both velocity and Reynolds stress.The most important differences are found at 2 behind the obstacle, where the fluctuation peak is underestimated.On the contrary, very good comparison is found on top of the cube both for the average and fluctuating quantities.(ii) Our results are similar to recent LES simulations (Muld et al., 2012) of the same case.(iii) Although the agreement with experiments appear a little better for the results obtained with the fine-grid LES1, notably near the cube, globally the difference between the profiles computed with the two meshes is small.In particular, no significant improvement is found for the fluctuations.Having considered also other meshes (not shown here), the mesh LES2 has been found to be the coarsest for which the results are in acceptable agreement with the experimental data, and convergence is reached.For this reason, the simulations with the particles have been carried out with the mesh LES2.

RANS computations
The boundary conditions can be set consistently for each twoequations model, thanks to the relations among ,  and ; these are: (i) at the inlet, in dimensionless unit:  = 1.5 ⋅ 10 −2 ,  = 2 ⋅ 10 −3 and  = 0.175; (ii) at outlet the zero-gradient condition; (iii) at the lateral patches the slip condition; (iv) the standard wall functions implemented in OpenFOAM are used for solid walls.The wall functions limit the turbulent dissipation, , and the turbulent specific dissipation, , in simulations that use a transport equation for  or , respectively (Menter and Esch, 2001;Popovac and Hanjalic, 2007).We have carried out several simulations for the different RANS models to test grid-convergence, until reaching a satisfying accuracy.The first set of simulations have been performed using a grid roughly corresponding to the mesh used for LES2 without refinement in the cube proximity.A second and third set of simulations, were performed after subsequent refinements, using as initial condition the solution of the previous simulation.The number of cells and quality indicators for the grids employed in the RANS simulations is reported in Table 3.The grid with highest refinement (denoted by RANS1/LES1 in Table 3) is overresolved with respect to the RANS simulations, and all turbulence models gave virtually the same results on RANS2 and RANS1 for the fluid (not shown in the paper).However, using the same resolutions for   Martinuzzi and Tropea (1993).

Table 3
Meshes employed for RANS simulations.The table reports the number of points, cells and four indexes of the quality of the meshes: the percentage of hexahedra among the cell, the maximum non-orthogonality, the maximum skewness, and the maximum aspect ratio of the cells.both sets of simulations avoids additional uncertainty when comparing dispersion models, because the Eulerian grid needs to be used in both cases to sample the Eulerian fields for particle tracking.

N points
The main flow features obtained with LES2 and the RANS are illustrated in Fig. 4. As anticipated, we have used different models, namely standard  −  and  − , and the Realizable  − .
In Fig. 4, a quantitative comparison among the different models is displayed for the average velocity on the centre-plane of the domain.In this figure, the colour map represents the mean velocity in the streamwise direction,  , the black contours denote  = 0, and the red lines are streamlines computed using the mean velocity components on the centre plane.Even though the residual convergence is fulfilled, discrepancy between the LES and the different RANS models are apparent.Overall, the  −  model is able to capture relatively well the length of the recirculation region behind the obstacle, but it fails to describe the separation caused by the obstacle edges.On the other hand, the − model is slightly better around the cube but overpredicts the length of the recirculation region.The Realizable  −  gives the best representation of separation, but also an even longer recirculation region than the standard  −  model.
In Fig. 5, a comparison between mean streamwise velocity and turbulent kinetic energy, , is shown for two vertical profiles extracted behind the obstacles.It is possible to observe that, in the wake region (lower values of ),  is lower for all RANS models than in the LES, which is a result of the longer recirculation bubble.In the region above the wake however, at  ≈ 1.5,  is in better agreement between the LES and the − and Realizable − models, which is a consequence of the larger acceleration above the obstacle in these models due to a more prominent separation.We also observe that the RANS models tend to under-predict  in the wake region and overpredict it above it when

Particle phase
We consider small heavy particles and examine the behaviour of three types of particles.The main dimensionless parameter that characterized infinitesimal particles is the Stokes number, defined here as  =    ∕, where  0 and  0 are the characteristic velocity and length scales of the fluid and   , defined in (4), is the particle relaxation time.The Stokes number is varied changing the particle radius, while keeping constant the density ratio   ∕  = 1000.The Stokes numbers investigated in this study are  = 10 −3 , 10 −1 , and 3.It is worth noting that gravity is not considered in the present study to simplify the problem and focus on turbulence effects.Whether it is actually negligible for a similar test case would depend on the physical dimensions.Similarly, collisions are neglected, assuming that the particle volume fraction is low enough, which is plausible for such a kind of flow.This allows the use of an arbitrary number of particle during the simulation, which makes easier reaching statistical convergence.
We design the particle injection as a theoretical model for cases of complex geometries where dispersion is limited to a certain region of the domain.The injection is continuous, releasing particles from a large number of points spread on a surface behind the cubic obstacle.More specifically, particles are continuously released from 37 2 points equally distributed on the square ( = 5.01, 0.1 <  < 1, 4.5 <  < 5.5) (the cubic obstacle occupies the region (4 <  < 5, 0 <  < 1, 4.5 <  < 5.5)), at a total rate of approximately 10,000 particles per time unit.This means that we handle in total ≃ 10 6 particles for simulations of ≃ 100 time units.
To provide a visual image of the kind of numerical simulations carried out, we present in Fig. 6 an instantaneous visualization of the particle dispersion in the flow as obtained in the LES for two different classes of particles, namely the smallest and the largest ones,  = 10 −3 ,  = 3.The visualizations reveal that large particles follow more closely ballistic paths, and therefore show less important lateral and vertical displacements, although the flow is highly chaotic.Furthermore, small particles display a larger variability in velocity, with very large fluctuations.More importantly, while many small particles explore the region ahead of the injection at the end of the cube, this is true only for a small amount of heavy particles, which anyway remain in the near vicinity of the walls.It turns out that the results of the numerical simulation of particles at  = 10 −1 are practically indistinguishable from those at  = 10 −3 .For this reason, these results obtained for this intermediate Stokes number will not be shown in the following.
As explained previously in Section 3.3, three models have been used to describe the particle phase when the RANS approach is used for the fluid, which have been labelled in the following way: We first discuss the profiles of the average velocity of the particles and of the fluid seen by the particles (relevant for the model accuracy), for two different Stokes numbers.Then, we shall discuss the particle dispersion for each  number.We will compare the average fields computed with the LES against those obtained with the different particle models coupled with two RANS models, namely with the realizable − and the  −  model.Since no experiment is available for the particle dynamics in this configuration, the LES is taken here as the reference result.

Particle velocity
We focus on comparing the particle velocity in the region downstream of the injection.Only a small portion of particle explore the region upstream of the injection, as discussed in the next section, so that in this region a comparison between particle velocity would be less meaningful.
In Fig. 7, we display the mean particle velocity   and the velocity of the fluid seen by particles   at the two locations  = 2,  = 3 behind the cube, for the  = 10 −3 particles.A first comment is that   and   are indistinguishable for this small  number case.Neither the location in the channel nor the RANS model coupled with the particle dynamics have impact on this result, so that for  = 10 −3 , one can safely write   =   in the entire domain, which is physically expected since particles act mainly as tracers in this regime.
Looking at the LES fields, it is possible to appreciate the reattachment at  = 3, while negative velocity is found at  = 2 at low .The profiles appear similar to those computed in the Eulerian LES, Fig. 3, and yet different.In particular, the discrepancy is more important for  > 1.5.The discrepancy between   and  is the result of the non-homogeneous particle distribution due to the injection location.Only particles which are caught by larger coherent structures explore the region above the obstacle wake, thus experience a flow history that does not correspond to all the possible states of the flow which contribute to the Eulerian average  .
Regarding the RANS simulations with dispersion models (denoted hereafter by RANS-PDF) important differences among the various models are nonetheless present.In the cube proximity, the stochastic model (Stoch) (16) gives overall better results when the RANS − is used.In particular,   at higher  and at the first streamwise location examined, closer to the obstacle, is in better agreement with that of the LES, where the other two models underestimate the particle velocity.When the  −  model is used, the OF2 and the Stoch models give similar results and are also in reasonably good agreement with the LES.The  OF1 model, on the other hand, gives the lest accurate predictions in all cases.Comparing the RANS approach to the LES, as expected the  −  model performs better in the separation region  < 1, but worse in the region above the wake, notably for  > 1.5.The discrepancies between the particle velocity estimated using the two different RANS models are more important near the cube, at  = 2.
We next discuss Fig. 8, where the same quantities are shown for the high-St case, i.e.,  = 3.Looking first at the LES, we see that there is an obvious difference not only between   and the Eulerian fluid velocity   (Fig. 3) because of the drift linked to the Lagrangian sampling, but also with the relative velocity   in the case of small inertia, cf.Fig. 7. Important inertial effects are demonstrated by the differences between   and   , which can be traced back to the crossing trajectory effect, typical of such a non-homogeneous flow.It is worth emphasizing that the differences we observe are of the order of 50% at some locations.In particular, the inertia strongly reduces the proportion of particles trapped in the re-circulation region.Indeed it is possible to see that   does not become negative downstream of the cube at  = 2, whereas   does.The profiles become more similar farther from the obstacle where the flow non-homogeneity is less strong.
As regards the RANS-PDF approach, we find the same trend displayed by the results obtained at  = 10 −3 , yet more emphasized at higher Stokes number.The predictions obtained with the OF1 model are largely inaccurate in all cases, although some inertial effects are captured.When the  −  model is used for the fluid phase, the results are satisfactory only with the Stoch model.However, while the qualitative behaviour is correctly captured by this model, the quantitative discrepancies with the LES are found to be important.In particular, the crossing-trajectory-effect (CTE) is exaggerated and the particle velocity is under-estimated.On the contrary, the OF2 model gives reasonable qualitative predictions, notably at  < 1, but largely underestimates the inertial effect displaying   and   profiles that are almost identical, despite the large Stokes number.Furthermore, the same trend is displayed with the  −  model.Unlike the cases with small Stokes number, only the Stoch model is capable to give good predictions while the OF2 model clearly underestimates the particle velocities for  > 1.Moreover, only the Stoch model predicts the important difference between   and   .We can appreciate that, in this case,   is accurately predicted by the Stoch model in the region  > 1.The predicted profile displays a decent agreement also for  < 1 at  = 2, while the re-attachment is not captured at  = 3.The mean particle velocity   is underestimated with respect to the LES results, yet inertial effects due to CTE are overall captured.Fig. 9. Contours of the averaged particles probability distributions on the lateral view of the domain for  = 10 −3 from the LES (black), the stochastic model used in this paper (red) and models OF1 and OF2, (blue and green, respectively).For each case, the two contours enclose 50% and 99% of the particles with age between 1 and 2 dimensionless time units (left column), and 4 and 5 dimensionless time units (right column).For the three RANS-based dispersion models, the mean velocity, turbulent kinetic energy and dissipation are computed using the  −  (top row) and realizable  −  (bottom row) models.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) The comparison between the two different models is less sharp for the higher Stokes number under investigation.As noted before for the smaller Stokes number, the − model appears to give better results in the separation region, with smaller differences than for the low-Stokes particles.Instead, the performance of the Realizable  −  model is largely superior in the central region,  > 1.
Concluding our discussion on the mean particle velocity, the following remarks are in order.
• The OF1 model is not able to describe the particle dispersion neither at small or at large Stokes.This underline the need of considering not only intensity and time scales, but also the flow non-homogeneity in modelling turbulent fluctuations.• The OF2 model includes some effect from the flow nonhomogeneity but does not take into account the CTE in a direct way.As a consequence, the model produces better predictions than the OF1, however it largely underestimates inertial effects because the particle velocity and the fluid velocity seen by the particles are barely different for higher Stokes number.• The Stoch model is the sole to include the CTE, and to fulfil basic consistency rules, notably the model is consistent with scalar diffusion in the limit of tracer dynamics.This is the reason for the superior performance in all cases.• Globally, Realizable  −  + Stoch model is the configuration that performs best.This is probably due to the better agreement of the mean velocity field predicted by the Realizable  − .However, it is interesting to note that the Stoch model is less influenced by the differences between RANS models than OF1 and OF2.

Particle dispersion
In this section we study in more detail particle dispersion, comparing the results obtained with the different models against the LES.We focus first on the small Stokes case, namely  = 10 −3 , for which we have seen that the particles basically behave as tracers.We display in Fig. 9 the particle probability distribution projected on the lateral view of the domain for the LES and the three dispersion models, for particles of age between 1 and 2 dimensionless time units and age between 4 and 5 dimensionless time units, where the particle age is the time spent by the particle in the flow after release.Most of the young particles remain near the cube, and the different models give similar results, as highlighted by the 50% isoline (the line delimiting the region containing 50% of the released particles).More important differences are found when looking at the ''tails'' of the distributions (the regions delimiting 99% of the particles), for which fluctuations have a more significant impact.All dispersion models underestimate the movement of particles upstream of the injection, in particular when the  −  model is used for the fluid (this model performs poorly in the separation regions on the obstacle surfaces).On the other hand, the longer recirculation region predicted by the realizable  −  leads to underestimate the dispersion over the cubic obstacle.We can observe that the OF1 model offers always the worst predictions, in line with what obtained for the velocity statistics.Overall, the best agreement with LES is given by the stochastic approach when coupled to the realizable  −  model, while the stochastic and the OF2 models give similar results when coupled with the  −  model for the fluid model.The results for particles allowed to disperse in the flow for longer times confirm the scenario, even though the stochastic and the OF2 models give closer predictions in all cases investigated here.
In general, the LES shows that particles initially disperse in a rather symmetric way around the obstacle, and that the dispersion becomes more important ahead of the obstacle as time goes on.
The same data are displayed in Fig. 10 in the cross-stream − plane to analyse the dispersion around the sides of the obstacle.For the sake of clarity, we only display particles that spent more time in the flow.Overall, the results confirm what observed above for the dispersion in the - plane.The OF1 model is unable to give physically-sound results, as particles do not disperse in a significant way and remain basically trapped in the cube vicinity.The other two models behave similarly and are in relatively good agreement with the LES data.The CTE is not particularly important in the cross-stream direction, and in particular far from the obstacle, so that both models are able to capture the basic features of the particle dynamics, although dispersion is weakly underestimated by the RANS models, as also seen from the data in the streamwise-wall-normal plane.
Finally, we investigate the case with higher Stokes number,  = 3, that is the case of large-inertia particles.The statistical properties of the dispersion in the - plane are shown in Fig. 11.For inertial particles, the dispersion is initially very small as expected (see data pertaining the younger particles), and therefore data obtained at these early times are not useful to assess the quality of the different models.At later times, Fig. 10.Contours of the averaged particles probability distributions on the frontal view of the domain for  = 10 −3 for LES (black), the stochastic model described in this paper (red), and models OF1 and OF2 (blue and green, respectively).For each case, the inner contours enclose 50% of the particles with age between 4 and 5 dimensionless time units, and the outer contour encloses 99% of the particles.For the three RANS-based dispersion models, the mean velocity, turbulent kinetic energy and dissipation are computed using the  −  (top row) and Realizable  −  (bottom row) models.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 11.Contours of the averaged particles probability distributions on the lateral view of the domain for  = 3 from the LES (black), the stochastic model used in this paper (red) and models OF1 and OF2, (blue and green, respectively).For each case, the two contours enclose 50% and 99% of the particles with age between 1 and 2 dimensionless time units (left column), and 4 and 5 dimensionless time units (right column).For the three RANS-based dispersion models, the mean velocity, turbulent kinetic energy and dissipation are computed using the  −  (top row) and realizable  −  (bottom row) models.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)dispersion is more significant and model comparisons more insightful.Comparing to the case of small-inertia particles, we see that dispersion is less important both upstream and downstream of the obstacle.While the OF1 model still provides non physical results, the OF2 and stochastic diffusion models perform similarly.Overall, the stochastic approach together with the  −  model appears the best, notably looking at the reattachment points before and after the obstacle.Yet, this model seems to overemphasize inertial effects, whereas the OF2 model presents slightly larger fluctuations and therefore dispersion.
The results in the  −  plane confirm these findings, see Fig. 12. Neglecting the OF1 model, the other models provide quite similar results.The predictions are better when using the  −  than the  −  model; the dispersion is underestimated with respect to the LES, suggesting that the reconstruction of the velocity of the fluid seen by the particles is reasonable while the effect of inertia is overestimated.

Computational cost
Four different types of numerical simulations have been carried out in this study: (i) time-dependent simulations of the fluid; (ii) time-dependent simulations of the fluid including particles; (iii) RANS simulations of the fluid; and (iv) particle simulations using dispersion models.The computational cost of these simulations is different because of the different number of time steps or iterations needed to obtain meaningful results and the different cost of each time step or iteration.We report in Table 4 the number of iterations or time steps, as well as the total computational cost for a single simulation.Note that the total computational cost is computed as the wall-clock time required by the simulation, multiplied by the number of cores that are used.Hereafter, for the sake of brevity, we will refer to ''time step'' for all simulations, even though a different set of operation is carried out in each case.
Simulations LES1 and LES2 belong to type (i), representing the time-resolved simulations of the fluid, which have been discussed in Sections 4.1.1 and 4.1.2.The LES with highest resolution is the most expensive simulation by far, with a total computational cost of almost 400,000 CPUh.Note that simulations LES1 and LES2 are run for approximately the same number of dimensionless time units, but a lower time interval is covered in LES1 each time step, explaining the higher number of time steps.The cases denoted by LES2 (P), which are the LES with particles described in Section 4.2, belong to type (ii) and are potentially the most expensive ones, because each time step requires both the solution of the incompressible Navier-Stokes equation and particle tracking.However, the relatively small number of Fig. 12. Contours of the averaged particles probability distributions on the frontal view of the domain for  = 3 for LES (black), the stochastic model described in this paper (red), and models OF1 and OF2 (blue and green, respectively).For each case, the inner contours enclose 50% of the particles with age between 4 and 5 dimensionless time units, and the outer contour encloses 99% of the particles.For the three RANS-based dispersion models, the mean velocity, turbulent kinetic energy and dissipation are computed using the  −  (top row) and realizable  −  (bottom row) models.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 4
Computational cost of the simulations.The table reports the total number of time steps (for time-resolved simulations) or iterations (for averaged simulations) required for convergence, the number of MPI ranks, and the corresponding Total computational cost in CPU hours.Cases LES1, LES2, and RANS1 are simulations of fluid only.LES2 (P) denotes simulations of both fluid and particles.Stoch, OF1, and OF2 denote simulations of particles only, carried out with dispersion models.particles in the domain at each time step, and the fact that both particle collisions and particle-fluid momentum transfer have been neglected, result in a very similar computational cost between LES2 and LES2 (P).Averaged simulations of the fluid only, which belong to type (iii), are the least expensive.In these simulations, each time step corresponds to a iteration of the SIMPLE algorithm and the simulation ends once the residuals of all variables have reached convergence.Particle simulations using dispersion models do not require the solution of the Navier-Stokes equation because only particle tracking is performed, which is a much cheaper operation.However, they need to cover a time interval long enough to reproduce the statistically-stationary condition of the reference simulations, which leads to a relatively large number of time steps.As expected, more complex dispersion models require a higher number of operation for each particle, and are thus more expensive.Note, however, that the continuous stochastic model has been implemented for the first time in OpenFOAM and that the implementation has not been optimized yet.

Conclusions
In this work, we have performed an extensive numerical analysis of a turbulent particle-laden flow in a test-case, which is a standard benchmark for the development of numerical simulations of singlephase fluid mechanics, namely the flow around a wall-mounted cube.To best of our knowledge this test case is examined here for the first time for particle-laden flows, and we have considered cases with both small and large particle inertia.
The purpose of the work is to provide a benchmark test for the development of reliable numerical modelling of turbulent dispersed flows.The test-case is particularly interesting for its engineering relevance and the challenging features to be captured, namely flow separation and stagnation.Indeed, while a serious assessment of the available models has been carried out since long time for the single-phase flows, the same is not true for the particle-laden flows, where simple heuristic models are usually used in applications.
The flow is investigated both with LES and RANS-based models, and a thorough comparison is made.In both cases, a hybrid approach is taken, solving on an Eulerian basis the fluid phase and through Lagrangian tracking the particle phase.
The LES has been performed with standard subgrid models and no specific subgrid model has been added for the particle phase.Different RANS models for the fluid have been considered, and eventually the Realizable  −  and  −  models have been used for the multiphase flows.In the RANS framework, the particle phase has been described with two simple discrete stochastic models already available in the OpenFOAM package, and with a consistent continuous diffusion stochastic model developed previously and implemented in OpenFOAM for this work.It is worth emphasizing here that it is crucial to couple consistent fluid and dispersion models to get accurate predictions, most notably in RANS framework.Therefore, the present class of fluid models appear the best suited for the dispersion models used.To couple them with more sophisticated models like  2 −  might lead to huge errors, even though the fluid model had given better predictions.Moreover, given that this is the first time this test-case is considered for particle-laden flows, no particular calibration of the models have been attempted, in order to emphasize the basic features of the different approaches for engineering-oriented applications.The improvement and the calibration of the models will be the object of future work.
Since neither DNS nor experiments of such a complex flow are currently available for cases including particle dispersion, a well-resolved LES has been carried out and used as reference solution.Of course, errors are to be expected, notably for particle statistics because of the lack of a specific model.However, we have used a fine-grained LES for this test-case such that we can be confident that the reconstruction of the fluid velocity seen by particles may be considered accurate enough at least for the low-order statistics investigated here (Armenio et al., 1999).We therefore consider the present framework as satisfactory with respect to this first benchmark.Nevertheless, significant errors may be encountered looking at the near-to-the-wall region and considering features like particle preferential concentration as obtained in LES, as highlighted in previous studies (Innocenti et al., 2016).To perform a DNS of the present case together with the assessment of a more complex particle LES model is suggested as the subject of future work.
From a physical point of view, particles behave as passive tracers in the small-inertia case and they therefore show the Lagrangian characteristics of the fluid flow.For larger inertia, a significant difference is found between the particle velocity and that of the carrier phase.The particles are still influenced by flow separation and stagnation, but the dynamics is smoother, in particular near to the obstacle.In the far region, particle and fluid velocity profiles become rather similar.The same observations can be made when looking at the geometrical dispersion of particles.
From the point of view of the model performances of the particle simulations based on RANS, the simplest dispersion model (OF1) based on a rough randomization of the turbulence fluctuations seen by the particles is found to give unacceptable results and should be abandoned.Furthermore, it is the most sensitive in respect of the specific RANS model employed for the fluid.The stochastic PDF model (stoch) gives the best agreement with the LES, especially when coupled to the Realizable  −  model for the fluid and in the case of higher particle inertia.The results obtained through the more refined discrete stochastic model (OF2) lie in-between the two other models.This model performs globally worse than the PDF one, yet it is qualitatively comparable to it in many case giving overall reasonable predictions.Since we have not pursued any calibration of the diffusion stochastic model, its formulation free from spurious drifts and taking account for the crossing-trajectory effect appears superior and more promising for future developments.
Overall, the discrepancy between the LES and the RANS/PDF predictions is still important, even in the best case.In fact, we have used the models that are the common choices in engineering applications, and yet it is to be stressed that a significant error on the prediction of the fluid-phase alone is present for such a complex flow.To greatly improve the present predictions, it would be needed to use a more effective RANS model for the fluid.Because of the consistency issues in two-phase flows, that would imply to develop a consistent stochastic model for the particles.While possible, it is not clear whether this kind of complex development is viable for engineering applications.
On a more modest ground, some improvements could be obtained performing a complete calibration of the diffusion stochastic model.The development of effective strategies to perform the calibration will be the objective of future studies.It is important to note that heuristic discrete models cannot be improved by construction, and therefore improvements are only possible by further developing the stochastic model.
Finally, considering the complexity of the problem, and that the flow has not been investigated before, we have simplified the configuration not considering the gravity in the present work.It is known that gravity may have an important crossing-trajectory effect and the assessment of the dispersion models in the general case also constitutes an important future direction.

Fig. 1 .
Fig.1.Visualization of the mid wall-normal plane and turbulent structures identified with the Q criterion(Hunt et al., 1988) at an arbitrary time step of the simulation.The structures are coloured with the average-velocity magnitude, whose range of values is from 0 (dark blue) to 1.8 (dark red).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Grids employed in the simulation, illustrated together with the instantaneous fluid velocity at an arbitrary time step, for (top) LES1 and (bottom) LES2.The grid employed for LES2 is the same as in the RANS simulations.

Fig. 3 .
Fig. 3. (a-b) Vertical profiles of the mean streamwise velocity,  , at different streamwise locations: from left to right,  = 2,  = 3. (c-d) Vertical profiles of the shear Reynolds stress at the same locations.The symbols denote profiles extracted from the dataset created byMartinuzzi and Tropea (1993).

Fig. 4 .
Fig. 4. Streamwise component of the mean velocity for (from top to bottom) LES2,  − ,  − , and Realizable  −  RANS.The black contours denote  = 0, and the red lines are streamlines computed using the mean velocity components on the centre plane.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Particle cloud coloured with the streamwise velocity for (left)  = 10 −3 and (right)  = 3 at an arbitrary time step of the simulations.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 1. OF1: the Gradient Dispersion model available in OpenFOAM.2. OF2: the Random Dispersion model available in OpenFOAM.3. Stoch: the Diffusion Stochastic model we have implemented.

Fig. 7 .
Fig. 7. Mean particle velocity   (thick lines) and the velocity of the fluid seen by particles   (dotted lines) at the two locations  = 2,  = 3 for  = 10 −3 .On the left panels, the particle models are coupled with the  −  model for the fluid.On the right panels, the particle models are coupled to the realizable  −  model.

Fig. 8 .
Fig. 8. Mean particle velocity   (thick lines) and the velocity of the fluid seen by particles   (dotted lines) at the two locations  = 2,  = 3 for  = 3.On the left panels, the particle models are coupled with the  −  model for the fluid.On the right panels, the particle models are coupled the Realizable  −  model.

Table 2
Meshes employed in the LES validation.The table reports the number of points, cells and four indexes of the quality of the meshes: the percentage of hexahedra among the cell, the maximum non-orthogonality, the maximum skewness, and the maximum aspect ratio of the cells.