Stochastic Large Eddy Simulation of an Axisymmetrical Confined-Bluff-Body Particle-Laden Turbulent Flow

Flow turbulence modulation by dispersed solid particles in a bluff body was studied using two-way-coupled stochastic large eddy simulat ion. Point-force scheme was used to model the inertial particle back effects on the fluid motion. The fluid velocity field seen by inertial part icles was stochastically constructed based on the filtered flow field obtained from well resolved large eddy simulations. For that purpose a Langevin-type diffusion process was used with the necessary modifications to account for particle inert ia, cross-trajectory effects and the two-way coupling. The numerical results regarding mean and turbulence statistics for both phases show a very good agreement with the experimental findings for both light and heavy mass loadings (22% and 110% respectively). This numerical investigation demonstrates also the ability of the stochastic-LES-part icle approach to predict turbulence modification by inert ial particles.


Introduction
In turbulent dispersed flows, particles are driven by the carrier phase and as a consequence this phase experiences reaction forces from these particles. The contaminating character of particles locally mod ifies the turbulence, since inside the particle there is no flow field and no eddies. This kind of local modificat ion is expected to be negligib le if the particle d iameter is much smaller than the Kolmogorov scale [1]. Ho wever, when higher particle loadings are considered (volume fraction > 10 -6 ), particles significantly affect the turbulent flow field and the two-way coupling between the phases has to be taken into account.
The presence of particles in a turbulent flo w can modulate the turbulence in several ways. Particles can distort streamlines or modify velocity gradients leading to a change in turbulence generation.Large particles can generate wakes yielding a turbulence enhancement while s mall part icles cause turbulence damping by the drag forces on them since kinetic energy is converted into heat. This has been the general trend shown by the experimental data [2][3][4][5] based on which many demarcation criteria fo r attenuation or enhancement of gas turbulence in the presence of part icles were proposed [6][7][8][9].
With the steady increase in computing power, there have been numerous efforts to nu merically quantify turbulence modulation by inertial particles. Ho wever, h ighly resolving the flow around thousands to millions of particles to get an accurate particle/turbulence interaction has been prohibited by the number of grid points required. Eaton and Segura [10] showed that one million grid points is needed to mesh a local spherical grid of diameter equal to 25d p to keep errors, in comparison with high-resolution experiments, under 1%. Thus, physical models have been developed and "plugged" to well-resolved numerical simu lations to render prediction of turbulence modulation tractable. A co mmonly applied model is the point-force model that was first used by Squires and Eaton [11] and subsequently improved bySundaram and Collins [12] and Lo mholtet al. [13]. The point-force scheme was originally an adaptation of the Particle-source-in-cell (PSIC) [14]. As stated by Eaton and Segura [10] , the point-force model, when added into a Direct Nu merical Simu lation (DNS), is believed to capture the effects of particles on the energy-containing eddies while being incapable of capturing the ext ra v iscous dissipation associated with the small-scale-turbulence/particle interaction. It is then expected that this fact will be exacerbated when the point-force model is used in the context of Large Eddy Simu lation (LES) where the small-scale or the sub-grid-scale motion is discarded by the spatial filtering operation inherent in LES calculat ions. Segura et al. [15] perfo rmed a well-resolved LES of a channel flow using the point-force coupling model. The single-phase flow and particle motions for light loading cases (mass loading ≈ 20%) were accurately predicted compared to the experimental results of Paris and Eaton [16] where as the turbulence attenuation was grossly under-predi cted. There are some reasons behind the failure of the point force scheme to capture the high levels of local dissipation around the particles. One reason is that most of the Lagrangian models relate the particle drag to the undisturbed flu id velocity wh ich is not availab le in a two-way coupled LES. The undisturbed flu id velocity field is defined as the velocity field that would exist if the particle was not present. Another reason is the under-prediction of particle drag by Lagrangian tracking models because the models do not account for the effects of subgrid scale turbulence that is missing in the LES field. Thus, the representation of the ext ra dissipation that occurs at the sub-grid scales may yield amelioration of the point-force LES predict ions for Lagrangian two -phase flows.
In recent reviews on two-way coupled turbulence simu lat ions of gas-particle flows, Eaton [17] and Balachandar and Eaton [18] acknowledged the shortcomings of the point-force model in accurately capturing particle-turbulence interaction. They suggested the development of LES with a sub-grid scale model that is locally modified by particles since high-resolution numerical simulat ions of particle-laden flows are still beyond hope. This suggestion is in line with the outcomes of Elghobashi [19] work who pointed out that in case of a significant two-way coupling, the sub-grid scale turbulence model might need modificat ion specifically if the particulate phase couples strongly with the small scales turbulence (matching between the particle and the sub-grid scale turbulence time scales). Thus, the two-way turbulence modeling should be carried in a rigorous way that account for the particle presence. However, it is well known that there is no precise 'geo metrical' knowledge on the structure of turbulence in the presence of a suspension of discrete particles and this makes it extremely difficu lt to modify the one-way turbulence modeling as well as the theory of Kolmogorov on wh ich many turbulence closures for two-phase turbulent flows are based. In the case of two-way coupling, it is assumed that the presence of solid particles changes the transfer rates of energy and energy dissipation while the nature and structure of turbulence remain the same [20]. This imp lies that one-way turbulence models can be used for two-phase turbulent flows after taking into account the fluid-particle coupling through the void fraction. However, Boiv in et al. [21] indicated that the presence of particles in isotropic turbulence yield a non-uniform distortion of the energy spectrum. Squires and Eaton [11] demonstrated that the non-uniform distortion is dependent upon the particle relaxation time. In a different work by Squires and Eaton [22], the non-uniform d istortion of the turbulence energy spectrum by particles was attributed to preferential concentration. Also, Elghobashi and Truesdell [23] found that the coupling between particles and flu id yielded an increase in small-scale energy. The relative increase in the energy of the high-wave-number co mponents of the velocity field resulted in a larger turbulence dissipation rate. They also stated that the effect of gravity caused an anisotropic modulation of the turbulence and an enhancement of turbulence energy levels in the direction aligned with grav ity. A ll these findings point towards the fact that the nature and the structure of the energy transfer mechanis ms of turbulence are mod ified by the presence of particles. In this case, modeling of the energy dissipation has to be altered in a manner that takes into account such a modification in turbulence structure. However, the emp irical input required for both Reynolds-averaged and LES approaches at the present time makes it difficult to obtain a rigorous two-phase turbulence modeling as suggested by Eaton [17].
Due to this severe limitation, LES modeling of two-waycoupled particle-laden turbulent flow using the point-force scheme should ideally include the sub-grid turbulence motion as seen by the inertial part icles. The stochastic approach based on Langevin equation emerges as a promising candidate in this regard. Previous numerical investigations by Berrouk et al. [24,25] and Berrouk and Laurence [26] in the context of one-way coupled gas-particle turbulent flows showed its promising potential in accounting for the effects of sub-grid turbulence motion on inert ial particle dispersion and deposition. In the present work, a Langevin-type diffusion process is used in conjunction with Lagrange LES as a tentative model to compute mean, turbulence statistics and turbulence modulation as predicted by the experimental wo rk of Borée et al. [27] on gas-particle turbulent flow in a confined bluff body. The experiment of Borée et al. [27] presents a challenging validation case for the stochastic Lag rangian particle-LES model since it contains mu ltip le co mplex flow features such as strong recirculating zones created by dump geometry and mult iple stagnation points beside the high Reynolds number. Th is configuration is typical of an industrial applicat ion where the objective is to control the mixing of a fuel (pulverized coal) with the air. Minier et al. [41] showed the potential of RANS/PDF approach through the simulat ion of the experiment of Borée et al. [27]. Riber et al. [28] simu lated the experiment of Borée et al. [27] for the light mass loading case (M=22%) using both hybrid and Lagrange LES. They investigated the effects of the numerical schemes, gas and particle boundary conditions, the grid refinement and the sub-grid scale models on the particle-LES accuracy and they confirmed the potential of LES approaches for two-phase turbulent flows and their relative insensitivity to the details of the nu merical solver. Moreover, they recommended accounting for the effects of subgrid flu id turbulence on the particle dynamics which is particularly crucial for the heavy mass loading case.

Governing Equati ons for Flui d Fl ow
The equations for the carrier flu id flow are solved using the LES methodology. The LES model is based on spatial filtering of the velocity field : plays the role of a low-pass filter that eliminates scales smaller than the filter width ∆. This results in smoothing the signal u which is usually random and unpredictable. The velocity can be decomposed into resolved and sub-grid scale components: Inserting this definit ion into the Nav ier-Stokes equations, the filtered continuity and mo mentum equations are obtained: The resulting filtered Navier-Stokes equations contain an unknown sub-grid scale (SGS) stress term that must be defined using a closure approximation: In this wo rk, the subgrid scale (SGS) stress tensor is modelled using the algebraic eddy viscosity model proposed by Smagorinsky [29]: the resolved rate-of-strain tensor. The value of the constant s C is evaluated based on the dynamic procedure [30,31].
The filter width, ∆ is taken equal to The two-way coupling effects are incorporated in the form of a reaction force exerted by the particles on the fluid. This reaction force is equal and opposite to the sum of fluid forces, only drag force in the present simulation, acting on the particles. Therefore the source term i r S in the mo mentum equation is given by: particles, respectively. p u is the velocity of a g iven particle p and s u is the fluid velocity seen by the part icle p along its trajectory. It is computed based on the resolved velocity using a stochastic diffusion process as we shall detail in the next section.

Governing Equati on for Particle Motion
The particle motion is calculated in the Lagrangian frame of reference. This approach requires interpolation of the surrounding carrier flu id velocity ū i onto the particle position. This was done using a tri-linear interpolation scheme. The particle mot ion in the Lagrangian frame of reference is described by the follo wing set of equations [32] (8) Here p x is the particle position, g is the gravity force per unit of mass, The effect of particle-part icle collision is neglected for these simu lations because it is deemed that the mo mentum exchange due to p-p collisions is much smaller than the mo mentu m exchange due to the forces exerted by the gas since the volume loading ration is still below values usually characterizing contact-dominated flows. We shall discuss this point later in mo re details.
The unknown in the system of Eqs. (8) is the fluid velocity s u seen by inertial part icles along their trajectories as they move through the turbulent field. The Eulerian velocity field computed using Equations of section 2.1 and interpolated to the particle positions, contains only part of the in formation of the velocity field that inert ial particle should see. This informat ion is linked to the filtered velocity of the large scales. Any information about the sub-grid scale mot ion is lost because of the filtering operation. Th is information on the turbulence at the subgrid scale level is crucial for the transport of inertial particles with response times smaller than the smallest LES-resolved time scales. The Langevin model can be used to reconstruct the Lagrangian instantaneous fluid velocity increment as seen by inertial part icles based on the LES filtered velocity. The model should account for the inertial character of the particles and the presence of a body force. The general form of the Langevin model chosen for the velocity of the flu id seen by particles is: and delta-correlated in the time do main [33].
The theoretical and nu merical formu lations of the Langevin model have been extensively discussed in the framework of particle-laden RA NS [34,35]. The use of the Langevin model was extended by Berrouk et al. [24] for the modeling of time increment of the flu id velocity seen by inertial part icles in LES framework. A detailed derivation of the different terms of the Lagrangian model is provided by Berrouk et al. [24,26] for the case of one-way coupled LES. Hereafter, we shall discuss only the evaluation of Langevin equation's terms in the case of inhomogeneous and anisotropic turbulence and how to account for particle inert ia, cross-trajectory effects (CT), and the two-way coupling.
For t wo-way coupled situations, the Langevin equation reads:  (11) St is Stokes nu mber based on the Eulerian SGS time scale and β is the ratio between the Lagrangian and the Eulerian time scales. It is assumed that β keeps the same value across the different scales of turbulence: It was shown that its value is Reynolds number dependent [37] and varies considerably in the literature. For this study, it is expected that its influence on the model predictions is very small, since small universal scales are modelled. In this study, β is chosen to be 0.356 [36]). Eq.
(11) was developed for particles interacting with homogeneous and isotropic turbulence. Its use in the present context to account for inertia effect on Lagrangian subgrid time scale is more appropriate compared to its use to include inertia effect on Lagrangian t ime scale in the framewo rk of RANS/SM where the construction of a wide spectrum of anisotropic turbulence fluctuations is sought through the stochastic modeling. For LES, the Lagrangian time scale for the sub-grid fluctuations SGS L T , is co mputed using the sub-grid kinetic energy SGS k and its dissipation rate SGS ε .
The last two quantities are evaluated according to the SGS model used to take into account sub-grid effects on the large scales. It reads following Hein z [38]: In the case of the Smagorinsky model [29], if equilibriu m is assumed at the cut-off, the production balances the dissipation. Thus, the SGS kinetic energy and its dissipation rate can be evaluated as: [39].
The directional dependence of the flu id Lagrangian SGS time scales SGS L T , is neglected since sub-grid scales are assumed to be ho mogeneous and isotropic. To account for the cross trajectory effect due to the presence of a body force, the Lagrangian time scale is exp ressed in the case of inert ial particle as function of the instantaneous relative velocity between the fluid and the inertial part icles [20]. If we assume direction (1) is the one aligned with the direction of the mean relative velocity and (2) and (3) are the transversal ones, we can use Csanady formu las [40] to compute the different anisotropic time scales: is the mean slip velocity between fluid and inertial particles. k is the resolved turbulent kinetic energy. In fact, Csanady formulas [40] also take into account the continuity effect. The continuity effect postulates that the inertial particle d ispersion in a direction perpendicular to the mean drift is twice as faster as inertial particle dispersion in a direction parallel to the mean drift.
The diffusion coefficient d C is evaluated according to the following formu lation [20]: Where SGS k is the modified SGS kinetic energy: Here i u' is the fluid fluctuating velocity and The mean statistical behavior of the cloud of particles that is present in one cell at time t should be taken into account. This results in an additional term in the equation of the drift [20]. In the fluid case, the drift term entering the stochastic differential equation is chosen to be a sum of a filtered term and a subgrid fluctuation term that is characterized by a time scale i SGS T , . Th is form is retained for thetwo-phase flow case with a modification of both altered and fluctuating terms to account for the inertia and cross trajectory effects.
In case of two-way coupling, it is assumed that solid particles impact the turbulent kinetic energy transfer rate and dissipation with a little effect on the nature and structure of turbulence. Based on that assumption, Min ier et al. [20,41] proposed a model fo r the seen fluid velocity in case of two-way coupling in RANS framework. They added a source term to the mean part of the drift term to account for the two-way coupling situations. This idea is retained for the LES calcu lations and depicted by adding the extra acceleration term in Eq. (10).

Numerical Simulation Overview
A flow solver fro m the R&D section of Electricite de France named Code_Saturne (www.code-saturne.org) was used to simulate Borée et al. [27] work. The discretization in Code_Saturne is based on the collocated finite-vo lu me approach. It allows solving Navier-Stokes and scalar equations on hybrid and non-conform unstructured grids. Velocity and pressure coupling is ensured by a prediction/correction method with a SIMPLEC algorith m. The collocated discretization requires a Rhie and Chow interpolation in the correction step to avoid oscillatory solutions. A second order centered scheme (in space and time) is used. The flow solver has been extensively tested for LES of single-phase flows [42,43]. The stochastic differential equations (SDE) system that comprises Eqs. (8) and (10) is integrated using an appropriate weak second-order integration scheme [35]  . When these scales go to zero, a hierarchy of stochastic differential systems is obtained. The Langevin-type model used in this study degenerates to a stochastic model for turbulent diffusion when approaches 0, that is, the inertial particles behave like flu id particles. The weak second-order integration scheme consists of a prediction step and a correction step. The prediction step is a weak first-order integration scheme (Euler scheme). By freezing the coefficients on the integration intervals and by resorting to Ito's calculus, it can be shown that the SDE system (Eqs. 8 and 10) has an analytical solution [26,35]. As test case for the developed model, an axisymmetrical confined-bluff-body particle-laden turbulent flow is simu lated. It consists of the experimental work of Bo rée et al. [27] perfo rmed using the flow loop Hercule at EDF R&D. In this experimental work, a lo w Reynolds number inner jet (Re ≈ 4500) laden with polydispersed glass particles interact with an outer jet characterized by a higher Reynolds number (Re ≈ 40000) creating therefore a zone of strong recircu lation (Fig. (1)). The outlet velocity of the inner jet was chosen low enough in order to obtain two stagnation points on the axis of the single-phase flow. The stagnation point configuration obtained is therefore interesting as inertia p roperties of the particles and flu id/particle coupling in the inner jet is expected to play a dominant role. Fu ll description of the test facility and of the working parameters is given by Borée et al. [27].
The geometry with mesh is shown in Figs. (2) and (3) where an overview of the grid resolutions in longitudinal and front cutting planes is presented respectively. The combustion chamber is 1.5m long with a 0.1m injection section placed in front. The dimensions of the bluff body are summarized in Tab. (1). Unstructured mesh consisting of 1.930,000 hexahedral cells is used for better resolution distribution as shown in Figs. (2) and (3). The non-conforming embedded refinement allows to concentrate more gridpoints in regions characterized by steep gradients and small energy-containing eddies like the near wall region and recirculation zones (0<z<0.25m and 0<R<0.075). In this region the y + is kept under 2 (less than 1 for 0<z<0.15m and 0<R<0.035). The rat ios ( + / + and + / + ) are kept under 5 to avoid oscillat ions in the gradient computations. The velocity profile at the inlet (z=-0.1 m) is prescribed by imposing the experimental mean and turbulence fluctuations measured at z=0 m. Non-slip conditions are imposed at the walls wh ile the outlet is purely convective.
A polydispersed glass particles are considered inBorée work [27] with material density ( = 2470 kg/m 3 ) and diameter that covers a wide range of size classes from to . Figure (4) shows the particle d istribution in mass and in number. Two mass loadings are considered in this experiment and they are both high enough to give rise to a t wo-way coupling between the carrier and the dispersed phases. The experimental findings show that for the highest mass loading (M=110%), the turbulence is modulated by the particles to an extent that it suppressed the two stagnation points that were observed for the smaller mass loading case (M =22%) and part icle-free case. The particle mean and fluctuation profiles are imposed at the inner pipe in let (z=-0.1 m, R= 0 m) and correspond to the ones measured experimentally at z=0 m. The time step used to solve the continuous phase is ∆t f = 10 -4 while the one used for the particulate phase is one order of magnitude smaller: ∆t p = 10 -5 . Injection conditions concerning the position and initial velocity of fluid and particle were documented in the experimental work [27]. For both phases, the instantaneous velocity vector was constructed by superimposing to the mean value a fluctuating value multip lied by a random Gaussian number.
The initial particle position is chosen from a random Gaussian distribution in a cross section of 100mm in diameter (inject ion tube diameter) at the centre of the chamber. The flu id velocity seen is set equal to the particle velocity.

Results and Discussion
In this section, the numerical results of the stochastic LES are co mpared to the experimental results of Borée et al. [27]. Experimental measurements are obtained with a Phase Doppler Anemo meter and experimental data for the continuous phase of the particle-free case (M=0%), continuous and dispersed phases for the particle-laden cases (M=22% and M=110%) are available for rad ial profiles of different statistical quantities at diffe rent axial distances downstreamof the injection. These quantities include the mean axial and radial velocit ies for both fluid and particle phases. Axial profiles along the axis of symmetry for these quantities are also provided.
Before detailing these results, the well-resolvedness of the present LES, justification of the use of the stochastic model,and neglect of the particle collisions are discussed. To ensure the large eddy simu lation well-resolvedness, the averaged amount of the turbulent kinetic energy discarded by the grid resolution and the LES modeling shouldnot exceed 20% of the total turbulent kinetic energy according to Celik et al. [44] and Pope [45].
A-posteriori estimation of the filtered out kinetic energy based on Eqs. (14) demonstrates that the present LES is adequate according to the LES index of quality [44,45]. Indeed, the ratio of the SGS kinetic energy to the total kinetic energy (KER) is g lobally well under 20% almost everywhere in the domain and consequently the present LES is adequate. This fact is depicted by Figs. (5a) and (6) for the streamwise and radial d irections respectively. Although the present LES is well resolved, there might be turbulent eddies discarded by the applied exp licit filter and which can affect the turbulent transport of all or part of the particle classes injected. This can be checked by comparing the particle response time to the SGS turbulence time scale. As a general rule, solid particles do sense turbulence with time scales smaller than their response times. This can be   Table (2) shows the response time of all the particle classes injected in the present simulat ions. These results assume Stokes flo w regime around particles with particle Reynolds number in the vicin ity of 1. In fact, nu merical Results (not included herein) show that the particle-and time-averaged particle Reynolds number varies between 0.2 and 10 depending on the particle class and its position in the flow. The last equation in the system of equations (8) that computes the particle response time as function of the particleReynolds number shows that a Part icle Reynolds number in the vicinity of 10 reduces the particle response time by 40%.Thus, particle response times presented in Tab.
(2) are only the upper limits and can be, depending on the particles' positionsin the flow field, much smaller and thus particles associated with them can see more resolved or discarded turbulence.  The time scale of the turbulence discarded by the LES filtering operation or the SGS turbulence t ime scale SGS T , is estimated using Eqs. (11)(12)(13) and depicted by Figs.(5b) and (7) for the streamwise and radial directions respectively.  Fig.  (4). The transport of these classes of particles should be affected by the discarded subgrid scale SGS turbulence. Hence the effect of latter should be modelled which justifies the use of the stochastic model described above.
Particle-part icle collision is usually taken into account for cases where the particle volu me fraction p α exceeds 10 -3 [18]. Present numerical results on the maximu m part icle volume fract ion (Fig. (10a)) show that the latter has a time-averaged value of around for M =22% mass loading and for the M=110% mass loading. Figure (10b) shows that these maximu m part icle vo lu me fractions are found for a limited nu mber of cells equalling 0.05% of the total cell number for M=22% mass loading and 0.5% for the M=110% mass loading. Thus, an algorithm handling particle collisions was deemed not necessary. The decision of not including a time-expensive p-p collision algorithm in the present simu lation is also justified by the fact that the time between successive inter-particle collisions, , is larger than the largest response time of all part icle classes present in the chamber (see Tab. (2)), whereby the fluid dynamic transport of the particles is the dominant effect.  The overall results regarding both carrier and part icle phases obtained using the stochastic LES model and the point-force scheme to account for the t wo-way coupling are in good agreement with the measurements of Bo rée et al. [27]. In the next sections, numerical results for continuous and dispersed phases are co mpared to experimental findings and the effects of the mass loading on both mean and turbulence statistics of the continuous phase are discussed.

Part icle-Free Case: M=0%
In the single-phase configuration, the experimental observations of Borée at al. [27] allow to schemat ically decompose the flow into three longitudinal regions as depicted by Fig. (1). Region A ranges fro m the jet nozzle to the first stagnation point S1. In this reg ion, the jet is rapid ly stopped by the recircu lating flow. The central jet is surrounded by a recirculation upward flow which feeds both the initial entrain ment in the jet and the annular shear layer developing at the edge of the bluff body. The positive pressure gradient due to the section increase should play an important role in the particle dispersion. Reg ion B ranges fro m S1 to S2 and represents the recirculation region. An intense upward flow on the axis is detected. Profiles measured further downstream of the second stagnation point S2 show that region C downstream S2 corresponds to the development of a wake flow. In two-phase configuration, this description remains valid fo r lo w mass-loading only. It was shown that no mean stagnation points on the flow axis is detected for the M=110% case as a result of a strong two-way coupling in the dense central jet.  (11a) and (11b) d isplay the streamwise profiles of flu id mean and RMS streamwise velocities for the particle-free case (M=0%). The numerical p redictions are in a very good agreement with the measurements in particular for the mean streamwise velocity (Fig. (11a)): very good predictions of the jet penetration, the position of the t wo stagnation points, and the width of the recirculation zone. The level of the turbulent intensity (Fig. (11b)) is also well predicted although a small discrepancy between the numerical and experimental pred ictions is observed in the zone z < 0.1m.  Radial profiles for fluid mean streamwise and rad ial velocities are depicted by Figs 12 and 13 respectively at six stations along the axis. Very good agreement with the experiment is obtained for the mean streamwise velocity while some departure is noticed for the radial velocity for the zone z≈0.2m and R<0.05m The radial profile for the fluid turbulent kinetic energy is shown in Fig. (14). The agreement of the nu merical results with measurements is also good despite some discrepancies near the zone z=0.16m and R<0.05 Overall, the numerical set-up (grid, boundary conditions, numerical discretization, and turbulence modeling) adopted for this b luff-body configuration gives close predictions of the experimental findings. The accuracy of the fluid solver is good enough to test the dispersed phase with reasonable confidence.

Part icle-Laden Case: M =22%
This section discusses the results for the low particle mass loading: M=22%. Figures (15a) and (15b) display the axial profiles of fluid mean and RMS streamwise velocities respectively for both particle-free (M=0%) and particle-laden (M=22%) cases. Both numerical results and experimental observations show that the impact of the dispersed phase on the gas phase at this low mass loading ratio is very negligib le outside the circu lation zone wh ile limited inside it : the central jet penetrates slightly further in the chamber, also slightly modifying the location of the recircu lation zone (Fig. (15a)). The slight discrepancy between the physical and the numerical results is noticed for the position and the value of the pick of the fluid RMS streamwise velocity for both single-phase and two-phase predictions (Fig. (15b)). This can be attributed to inadequate resolution in a region (z < 0.25m) characterized by high shear and strong recirculat ion. Radial pro files for the flu id mean streamwise velocity is depicted by Fig. (16) at six stations along the axis. Both experimental and numerical findings predict very minor modification of the flu id mean streamwise velocity by solid particles. Very good matching between the physical and numerical predictions is noticed.
Same remarks apply also for the radial pro file fo r the fluid mean radial velocity. Indeed, Fig. (17) shows that slight modulation of the fluid mean radial velocity by solid particles is predicted by both experimental and nu merical findings. Also, good matching between the experimental and numerical results is noticed except in a large part of the recircu lation zone (0.1m<z<0.24m & r<0.05m).
The radial profile for the fluid turbulent kinetic energy for single-and two-phase calculations is shown in Fig. (18) atsix stations along the axis. Both experimental and numerical results show that the modification of the gas turbulence by particles is minor except for the region 0.04m<R<0.08m downstream z=0.16m. Differences between the numerical results and the experiment are noticed in the recirculation zone but the overall agreement remains reasonable  The high mass loading case (M=110%) is more numerically challenging since the two-way coupling between the two-phases is expected to be strong. Figure (19) displays the axial profiles of mean streamwise velocity and turbulent kinetic energy for both single-phase calculations (M=0%) and two-phase calculations (M=110%). Nu merical calculations are directly co mpared with measurements. It is clear that the presence of the dispersed phase at a high mass loading (M =110%) yielded a significant modification of the jet penetration and a total suppression of the recircu lationzone as it is equally predicted by both experimental and nu merical findings. The discrepancy between the physicaland the numerical experiments results is observed for the extent of the jet penetration whereas very good matching was obtained regarding the min imu m value of the mean axial velocity as depicted by Fig. (19a).For the flu id RMS streamwise velocity, Fig. (19b) shows that differencesbetween the numerical and the experimental results are more significant in the zone (0.8m < z < 0.16m) which is the same zone for wh ich differences was also noticed between experimental and results regarding the mean axial velocity. In this narro w zone the turbulence intensity isunderestimated whereas the turbulence enhancement outside this zone was very well captured. Regarding the radial profile of turbulent kinetic energy that is depicted by Fig. (20), the numerical results capture well the damping of turbulence predicted by the experiment however, there is a noticeable d iscrepancy near the center of the chamber at z=160 mm. Nu merical results (not presented here) show that the highest part icle vo lu me fraction was reco rded in the reg ion (0.8m < z < 0.16m) withparticle concentrations approaching 10 -3 . Thus, collision between particles in this narrow zone might affect theparticle distribution in this region and hence the inter-phase coupling. This shortcoming can be compounded by the fact that the point-force scheme performs usually poorly for very high particle volu me fractions. This can exp lain the differences between the experimental and nu merical predictions that were noticed in this narrow zone.
Radial profiles for mean streamwise and mean rad ial velocities are depicted by Figs. (21) and (22) respectively. Both experimental and numerical findings predict minor modification by particles. Very good matching between the two types of results is noticed.

ParticlePhase
Overall, Nu merical predict ions of the particle-phase mean and RMS velocities in both axial and radial d irections are in good agreement with the experimental findings. As expected, the use of the stochastic LES model yields imp roved particulate results compared to results obtained using LES flow field (The latter are not shown here). Nu merical results for the lo w mass-loading case (M=22%) show better agreement with measurements compared to the high mass loading case (M=110%). Th is is expected since the point-force scheme is known for its decreased performance for particle void fraction close to 10 −3 .
It is worth to mention that in each control volume for the flu id phase, the different velocit ies of all the particles inside are mass-averaged following Fig. 4 (23) displays the streamwise profiles of mean, RMS axial and RMS rad ial particle velocities along the axis (R=0) for the low mass-loading case (M=22%). The results of the LES using the stochastic model are directly co mpared with measurements. Very good agreement between physical and numerical results is obtained except for the recirculation zone where minor d iscrepancies are noticed. In this criticalzone, particles decelerate, accu mulate near the stagnation points, and stop before changing direction to escape from this zone by its sides. Thus, this zone is difficu lt to predict accurately for the particulate phase which exp lain the differences between the numerical results and the experiment.
Radial profiles of the mean streamwise and mean radial particle velocities are depicted by Figs. (24) and (25) respectively at seven stations. Very good matching between the experimental and numerical results is noticed.    Radial profiles of the RMS streamwise and RMS rad ial particle velocities are shown in Figs. (26) and (27) respectively at seven stations.Reasonable agreement between the numerical andexperimental and numerical results is noticed in particular for the RMS streamwise particle velocity. For the RMS radial particle velocity, mo re noticeable discrepancies are observed for the region around z=0.24m that corresponds to the second stagnation point.

Mass loading: M=110%
Figure (28) displays the streamwise profiles of mean, RMS axial and RMS rad ial particle velocities along the axis (R=0) for the lo w mass-loading case (M=110%). The results of the LES using the stochastic model are directly co mpared with measurements. Very good agreement between physical and numerical results is obtained except for the RMS radial particle velocity downstream of z=0.15m    Radial profiles of particle RMS radial velocity for particle-laden (M=110%) configuration. Circle: Experiment; solid line: Numerical simulation. (a) z=0.003m; (b) z=0.08m; (c) z=0.16m; (d) z=0.20m; (e) z=0.24m; (f) z=0.32m ; (g) z=0.40m Radial profiles of the mean streamwise and mean radial particle velocities are depicted by Figs. (29) and (30) respectively at seven stations. Good matching between the experimental and numerical results is noticed.
Radial profiles of the RMS streamwise and RMS rad ial particle velocities are shown in Figs. (31) and (32) respectively at seven stations. Reasonable agreement between the numerical and experimental results is noticed in particular for the RMS streamwise particle velocity.
For the RM S radial part icle velocity, mo re noticeable discrepancies are observed for the region around z=0.16m.

Conclusions
This work has presented the stochastic large eddy simu lation approach to gas flow modificat ion by solid particles. The bluff-body case studied experimentally by Borée et al. [27] was chosen as a validation case because it contains mult iple co mplex flo w features which are typical of combustion chambers.
The point-force scheme was used to model the coupling between phases while Langevin-type diffusion process was deployed to account for the SGS turbulence discarded by filtering in LES calculations. The experimental findings show that for the high mass-loading case (Μ = 110%), the turbulence is modulated by the particles to an extent that it suppressed the two stagnation points that were observed for the smaller mass loading case (Μ = 22%) and part icle-free case.
Furthermore, the existence of a recirculation zone where particles interact with negative axial flu id velocities constitutes a much more challenging test case for simu lations compared to cases where the fluid and the particle mean velocities are of the same sign. This test case was an opportunity to assess the handling of two-way coupling situations in Code-Saturne and to see to what extent the stochastic model described earlier lends itself to LES of dispersed non-equilibriu m turbulent shear flows in a two-way coupling situation.
Nu merical results on mean streamwisevelocities, mean radial velocities, and turbulent kinetic energy in radial and axial directions show very good agreement with the experimental findings. Some d iscrepancies were noticed, in particular for the high mass loading case (Μ=110%). This is due to the limitation of the point-fo rce scheme which should be replaced by a more accurate approach to the inter-phase coupling modeling in future works.