W.P. Breugem
Please Note
54 records found
1
We studied the hydrodynamic stability of a viscous liquid jet enclosed by a much less viscous fluid in a narrow vertical tube. In literature, this flow pattern is also known as perfect core-annular flow. The main objective is to unravel the competition between capillary and shear-driven instability mechanisms acting on the flow. The temporal stability of the flow was tested under laminar conditions for a small axisymmetric sinusoidal perturbation of the interface. To this purpose, numerical simulations were conducted using a finite-volume two-phase flow solver combined with a geometric Volume-of-Fluid method to capture the interface between the immiscible fluids. The simulation results are interpreted using linear stability theory for thin liquid jets in free space. The main conclusion is that perfect core-annular flow is hydrodynamically unstable, either through a capillary or a shear-driven instability. The competition between the two instability mechanisms is characterized by the Weber number based on the annular layer thickness, Wea. For Wea≪1, the flow is prone to a capillary instability, while for Wea≫1, the liquid jet may undergo atomization. Evidence is also found for a reduced growth rate of capillary instabilities in the presence of strong shear at high Wea.
Newtonian die-swell phenomenon revisited
Theory and simulations
Characteristics of slurry transport regimes
Insights from experiments and interface-resolved Direct Numerical Simulations
A common way to transport solids in large quantities is by using a carrier fluid to transport the solids as a concentrated solid/liquid mixture or slurry through a pipeline. Typical examples are found in dredging, mining and drilling applications. Dependent on the slurry properties and flow conditions, horizontal slurry pipe flow is either in the fixed-bed, sliding-bed or fully-suspended regime. In terms of non-dimensional numbers, the flow is fully characterized by the bulk liquid Reynolds number (Re), the Galileo number (Ga, a measure for the tendency of particles to settle under gravity), the solid bulk concentration (ϕb), the particle/fluid density ratio (ρp/ρf), the particle/pipe diameter ratio (Dp/Dpipe), and parameters related to direct particle interactions such as the Coulomb coefficient of sliding friction (μc). To further our fundamental understanding of the flow dynamics, we performed experiments and interface-resolved Direct Numerical Simulations (DNS) of slurry flow in a horizontal pipe. The experiments were performed in a transparent flow loop with Dpipe=4 cm. We measured the pressure drop along the pipeline, the spatial solid concentration distribution in the cross-flow plane through Electrical Resistance Tomography (ERT), and used a high-speed camera for flow visualization. The slurry consisted of polystyrene beads in water with Dp=2mm, ρp/ρf=1.02, Ga between 40–45 and ϕb between 0.26–0.33. The different flow regimes were studied by varying the flow rate, with Re varying from 3272 till 13830. The simulations were performed for the same flow parameters as in the experiments. Taking the experimental uncertainty into account, the results from the DNS and the experiments are in reasonably good agreement. The results for the pressure drop agree also fairly well with popular empirical models from literature. In addition, we performed a parametric DNS study in which we solely varied Re and Ga. In all flow regimes, a secondary flow of Prandtl's second kind is present, ascribed to the presence of internal flow corners and a ridge of densely packed particles at the pipe bottom during transition towards the fully-suspended regime. In the bulk of the turbulent flow above the bed, secondary flow transport of streamwise momentum dominates over turbulent diffusion in regions where the secondary flow is strong and vice versa where it is weak. The transition between flow regimes appears to be governed by the competition between the net gravity force on the particles and shear-induced particle migration from particle–particle interactions. This competition can be expressed by the Shields number, θ. For θ≲0.75, gravity is dominant and the flow is in the fixed-bed regime. For θ≳0.75, shear-induced migration becomes progressively more important for increasing θ. Low-concentration zones flanking the sliding bed start to form at the top corners of the bed, and gradually expand downwards along the pipe wall till the pipe bottom is reached. For θ≳1.5, shear-induced migration is responsible for lifting the particle bed away from the wall, associated with the onset of the suspended regime. For θ≫1, gravity is of minor importance and the mean flow eventually reaches axi-symmetry with a high-concentration particle core at the pipe center and negligible secondary flow.
Interparticle friction in sheared dense suspensions
Comparison of the viscous and frictional rheology descriptions
In the literature, two different frameworks exist for describing the rheology of solid/liquid suspensions: (1) the “viscous” framework in terms of the relative suspension viscosity, ηr, as a function of the reduced solid volume fraction, f=fm, with fm the maximum flowable packing fraction, and (2) the “frictional” framework in terms of a macroscopic friction coefficient, μ, as a function of the viscous number, Iv, defined as the ratio of the viscous shear to the wall-normal particle stress. Our goal is to compare the two different frameworks, focusing on the effect of friction between particles. We have conducted a particle-resolved direct numerical simulation study of a dense non-Brownian suspension of neutrally buoyant spheres in slow plane Couette flow. We varied the bulk solid volume fraction from fb ¼ 0:1 to 0.6 and considered three different Coulomb friction coefficients: μc ¼ 0, 0.2, and 0.39. We find that ηr scales well with f=fm, with fm obtained from fitting the Maron–Pierce correlation. We also find that μ scales well with Iv. Furthermore, we find a monotonic relation between f=fm and Iv, which depends only weakly on μc. Since ηr ¼ μ=Iv, we thus find that the two frameworks are largely equivalent and that both account implicitly for Coulomb friction. However, we find that the normal particle stress differences, N1 and N2, when normalized with the total shear stress and plotted against either f=fm or Iv, remain explicitly dependent on μc in a manner that is not yet fully understood.
Results from particle-resolved Direct numerical simulations are presented for dense suspensions of frictional non-colloidal spheres in viscous pressure-driven channel flow. The bulk solid volume fraction varies between ϕb=0.2 and 0.6, and the Coulomb friction coefficient is either μc=0 or 0.5. The main objectives are to unravel the influence of (1) ϕb and μc on the flow development time and of (2) heterogeneous shear on the steady-state suspension rheology. Starting from an initially homogeneous distribution, the particles show shear-induced migration toward the core until equilibrium is reached. The flow development time decays exponentially with increasing ϕb/ΦR, where ΦR is a friction-dependent reference bulk concentration beyond which particle contacts cause a rapid increase in the particle stress. The steady-state rheology is studied by means of the ‘viscous’ and ‘frictional’ rheology frameworks. Excluding the central core and wall regions, the data for the local relative suspension viscosity collapse onto a single curve as function of the normalized local concentration ϕ¯/ϕm, where ϕm is the friction-dependent maximum flowable packing fraction. The frictional rheology shows ‘subyielding’ at low viscous number Iv in the core region, where the macroscopic friction coefficient μ drops below the minimal value found for homogeneous shear flows. A modified frictional rheology model is presented that captures subyielding. Finally, a model is presented for ϕ¯/ϕmp as function of Iv, where ϕmp is a modified maximum flowable packing fraction. It captures both ‘overcompaction’ in the core beyond ϕm at high ϕb and maximum core concentrations below ϕm at lower ϕb.
Controlling the breakup of spiralling jets
Results from experiments, nonlinear simulations and linear stability analysis
From nearly homogeneous to core-peaking suspensions
Insight in suspension pipe flows using MRI and DNS
Magnetic resonance imaging (MRI) experiments have been performed in conjunction with direct numerical simulations (DNS) to study neutrally buoyant particle-laden pipe flows. The flows are characterized by the suspension liquid Reynolds number (Res), based on the bulk liquid velocity and suspension viscosity obtained from Eilers' correlation, the bulk solid volume fraction (φb), and the particle-to-pipe diameter ratio (d/D). Six different cases have been studied, each with a unique combination of Res and φ, while d/D is kept constant at 0.058. The selected cases ensure that the comparison is performed across different flow regimes, each exhibiting characteristic behavior. In general, an excellent agreement is found between experiment and simulation for the average liquid velocity and solid volume fraction profiles. Root-mean-square errors as low as 1.7% and 5.3% are found for the velocity and volume fraction profiles, respectively. This study presents accurate and quantitative velocity and volume fraction profiles of semidilute up to dense suspension flows using both experimental and numerical methods. Three different flow regimes are identified, based on the experimental and numerical solid volume fraction profiles. These profiles explain observations in the drag change. For low bulk solid volume fractions a drag increase (with respect to an equal Res single-phase case) is observed. For moderate volume fraction distributions the drag is found to decrease, due to particle accumulation at the pipe center. For high volume fractions the drag is found to decrease further. For solid volume fractions of 0.4 a drag reduction higher than 25% is found. This drag reduction is linked to the strong viscosity gradient in the radial direction, where the relatively low viscosity near the pipe wall acts as a lubrication layer between the pipe wall and the dense core.
Particle-resolved Direct Numerical Simulations have been performed on the gravitational settling of mono-disperse solid spheres in a viscous fluid and triply periodic domain. In a comprehensive study, the bulk solid volume concentration was varied from ϕ=0.5 to 30%. To study the effect of inertia, three different Galileo numbers were considered in the inertial regime, Ga=144, 178 and 210, for which a single settling sphere exhibits distinctly different wake and path characteristics. The particle/fluid mass density ratio was fixed at 1.5. We find that for ϕ=2−30% the suspension microstructure and dynamics depend predominantly on the bulk concentration. In qualitative agreement with previous studies in literature, three different sedimentation regimes can be distinguished: (1) the dilute concentration regime for ϕ≲2% with preferential settling of particles in vertical trains, (2) the moderate concentration regime for 2%≲ϕ≲10% with preferential settling of particles in horizontal pairs with an interparticle distance of ∼ 1.5 particle diameters, and (3) the dense concentration regime for ϕ≳10% with a nearly random (“hard-sphere”) distribution of the particles in space. The clustering of particles is dictated by, respectively, trapping of particles in the wake of other particles, a drafting–kissing–tumbling (DKT) instability by which two vertically aligned particles quickly reorient themselves into a horizontally aligned particle pair, and short-range multiparticle interactions through viscous lubrication and to a lesser extent collisions between particles. In all cases, hindered settling at a reduced speed is observed as compared to a single settling sphere. The well-known Richardson–Zaki relation for the mean sedimentation velocity appears valid only for the dense concentration regime. We provide ample evidence that in the dense regime the characteristic velocity and time scales of particle motion are proportional to gDp and Dp/g, respectively, with g the gravitational acceleration and Dp the particle diameter. We also observe an ω−3 scaling of the particle velocity spectra for ωDp/g≳0.4 and we propose a model to explain this scaling behavior, based on the inertial response of the particles to small-scale flow perturbations. Kinematic waves, i.e., vertically propagating plane waves in the local concentration field, are observed in all cases, though unrelated particle motions are responsible for significant loss of the spatio-temporal coherence of the waves. The wave speed was determined from repeated space–time autocorrelations of the local concentration field and appears in reasonable agreement with Kynch sedimentation theory using the Richardson–Zaki relation. The passage of kinematic waves causes perturbations in the particle velocity at a frequency that matches well with peak frequencies in the particle velocity spectra for concentrations up to ϕ≈10%. The time-lagged cross-correlation of the vertical and horizontal particle velocity suggests that kinematic waves may trigger DKT instabilities, while conversely DKT instabilities may be responsible for the onset of kinematic waves. Finally, we suggest that obstruction and perturbation of the particle wake by neighboring particles could offer an explanation for the small influence of the Galileo number on the suspension behavior for ϕ=2−30%.
When two spherical particles submerged in a viscous fluid are subjected to an oscillatory flow, they align themselves perpendicular to the direction of the flow leaving a small gap between them. The formation of this compact structure is attributed to a nonzero residual flow known as steady streaming. We have performed direct numerical simulations of a fully resolved, oscillating flow in which the pair of particles is modeled using an immersed boundary method. Our simulations show that the particles oscillate both parallel and perpendicular to the oscillating flow in elongated figure-8 trajectories. In absence of bottom friction, the mean gap between the particles depends only on the normalized Stokes boundary layer thickness δ∗, and on the normalized, streamwise excursion length of the particles relative to the fluid Ar∗ (equivalent to the Keulegan-Carpenter number). For Ar∗≲1, viscous effects dominate and the mean particle separation only depends on δ∗. For larger Ar∗ values, advection becomes important and the gap widens. Overall, the normalized mean gap between the particles scales as L∗≈3.0δ∗1.5+0.03Ar∗3, which also agrees well with previous experimental results. The two regimes are also observed in the magnitude of the oscillations of the gap perpendicular to the flow, which increases in the viscous regime and decreases in the advective regime. When bottom friction is considered, particle rotation increases and the gap widens. Our results stress the importance of simulating the particle motion with all its degrees of freedom to accurately model the system and reproduce experimental results. The insights of the particle pairs provide an important step towards understanding denser and more complex systems.
Path instabilities of a sphere rising or falling in a quiescent Newtonian fluid have been studied experimentally. The rich palette of possible instabilities is dependent upon two dimensionless quantities, namely the Galileo number (Ga) and the particle/fluid mass density ratio (ρ¯). In recent literature, several (Ga,ρ¯) regime maps have been proposed to characterize path instabilities, based on both numerical and experimental studies, with substantial disagreements among them. The present study attempts to shed light on path instabilities for which previous studies disagree. A detailed experimental investigation has been conducted for 219 different combinations of Ga and ρ¯, grouped around four values of ρ¯ (∼ 0.87, 1.12, 3.19 and 3.9) and Ga in the range of ∼ 100 to 700. Our results agree well with literature for the low Ga range in which a particle takes a steady vertical or steady oblique path and for which all previous studies agree with each other. For the higher and more controversial Ga range, we discuss consensus and disagreements with previous studies. Some regimes, which were only recently observed in numerical simulations, have been observed experimentally for the first time. Also, intriguing bi-stable regimes (i.e., coexistence of two stable asymptotic states) have been observed. For all four investigated density ratios, an update of the regime map is proposed. Finally, for both the rising and falling spheres, the drag coefficient as function of terminal settling Reynolds number has been determined, which for the investigated density ratios does not differ significantly from that of flow past a fixed sphere.
In this work, a sensor is developed to obtain the hydraulic state gradient, based on Archimedes’ principle. A cylindrical heavy object is submerged in the fluidised bed and lowered gradually while its weight is measured at various heights using a sensitive force measuring device.
Based on accurate fluidisation experiments with calcite grains, the voidage is determined and a straightforward empirical model is developed to estimate the particle size as a function of superficial fluid velocity, kinematic viscosity, suspension density, voidage and particle density. The surface area and specific space velocity can be estimated accordingly, which represent key performance indicators regarding the hydraulic state of the fluidised bed reactor. The prediction error for voidage is 5 ± 2 % and for particle size 9 ± 4 %.
The newly developed soft sensor is a more time-effective method for obtaining the hydraulic state in full-scale liquid-solid fluidised bed reactors. ...
In this work, a sensor is developed to obtain the hydraulic state gradient, based on Archimedes’ principle. A cylindrical heavy object is submerged in the fluidised bed and lowered gradually while its weight is measured at various heights using a sensitive force measuring device.
Based on accurate fluidisation experiments with calcite grains, the voidage is determined and a straightforward empirical model is developed to estimate the particle size as a function of superficial fluid velocity, kinematic viscosity, suspension density, voidage and particle density. The surface area and specific space velocity can be estimated accordingly, which represent key performance indicators regarding the hydraulic state of the fluidised bed reactor. The prediction error for voidage is 5 ± 2 % and for particle size 9 ± 4 %.
The newly developed soft sensor is a more time-effective method for obtaining the hydraulic state in full-scale liquid-solid fluidised bed reactors.
Using a slender-jet approach, we numerically investigate the control of jet breakup using temperature modulation at the nozzle with a specified frequency and amplitude. Our results show that temperature modulation does lead to instability through capillary and Marangoni stresses, providing control of the droplet formation in terms of intact length and resultant drop size distribution, which is otherwise irregular due to inevitable presence of background noise. For understanding the mechanisms underpinning the breakup of a thermally modulated jet in the presence of noise, it is useful to decompose the surface tension forces into a contribution from curvature-gradient forces and a contribution from surface tension-gradient forces, associated with axial variations in the jet curvature and the temperature-dependent surface tension coefficient, respectively. We show that in the limit of slow axial heat diffusion and slow cooling to the ambient, as considered here, the breakup of a thermally modulated jet is governed by the ratio of the surface tension-gradient force to the imposed random perturbation force at nozzle exit. This so-called “thermal modulation strength number” depends on the amplitude and frequency of the thermal modulation, the sensitivity of the surface tension coefficient to variations in temperature, the Weber number, and the strength of the Gaussian white noise added to the nozzle exit velocity. We show that the thermal modulation strength number governs the shift in breakup characteristics from forward to rear pinchoff for increasing modulation strength as well as the nature of the instability. When thermal modulation is weak, the surface tension-gradient forces act only as a trigger, and curvature-gradient forces soon take over and grow exponentially downstream from the jet due to inertio-capillary growth. When thermal modulation is strong, the surface tension-gradient forces not only act as a trigger, but remain significant until breakup. The thermal modulation strength number is thus useful to the design of thermal modulation in practical applications as a possible alternative to often-used mechanical excitation mechanisms to control jet breakup.
The effect of drag reducing riblets on the flow structure was examined experimentally for a turbulent boundary layer at Reθ = 9890 and riblet spacing s+ = 13.4. Trapezoidal riblets were used, which were attached to the water tunnel wall as a coating. Force measurements were performed to quantify the amount of drag reduction. Then, the mechanism underlying this reduction was investigated by stereo-PIV measurements in the cross-stream plane. To determine the effect of the drag reducing riblets, the results were compared with the smooth flat plate. Time-averaged turbulent statistics such as turbulent kinetic energy and Reynolds shear stress were found to be lower over the riblets compared to the flat surface. Two-point correlations of the fluctuating velocity components were calculated to elucidate the average flow structure size and strength, where riblets significantly suppressed the turbulent structures. Quadrant analysis of the Reynolds shear stress was performed to assess the change in ejection and sweep events and the results were found to be in correspondence with previous works.
This work aims to analyse and explain the different causes of this spread. To this end, terminal settling experiments were conducted in a quiescent fluid with particles varying in density, size, and shape. For the settling experiments, opaque and transparent spherical polydisperse and monodisperse glass beads were selected. In this study, we also examined drinking-water-related particles, like calcite pellets and crushed calcite seeding material grains, which are both applied in drinking water softening. Polydisperse calcite pellets were sieved and separated to acquire more uniformly dispersed samples. In addition, a wide variety of grains with different densities, sizes, and shapes were investigated for their terminal settling velocity and behaviour. The derived drag coefficient was compared with well-known models such as the one of Brown and Lawler (2003).
A sensitivity analysis showed that the spread is caused, to a lesser extent, by variations in fluid properties, measurement errors, and wall effects. Natural variations in specific particle density, path trajectory instabilities, and distinctive multi-particle settling behaviour caused a slightly larger degree of the spread. In contrast, a greater spread is caused by variations in particle size, shape, and orientation.
In terms of robust process designs and adequate process optimisation for fluidisation and sedimentation of natural granules, it is therefore crucial to take into consideration the influence of the natural variations in the settling velocity when using predictive models of round spheres. ...
This work aims to analyse and explain the different causes of this spread. To this end, terminal settling experiments were conducted in a quiescent fluid with particles varying in density, size, and shape. For the settling experiments, opaque and transparent spherical polydisperse and monodisperse glass beads were selected. In this study, we also examined drinking-water-related particles, like calcite pellets and crushed calcite seeding material grains, which are both applied in drinking water softening. Polydisperse calcite pellets were sieved and separated to acquire more uniformly dispersed samples. In addition, a wide variety of grains with different densities, sizes, and shapes were investigated for their terminal settling velocity and behaviour. The derived drag coefficient was compared with well-known models such as the one of Brown and Lawler (2003).
A sensitivity analysis showed that the spread is caused, to a lesser extent, by variations in fluid properties, measurement errors, and wall effects. Natural variations in specific particle density, path trajectory instabilities, and distinctive multi-particle settling behaviour caused a slightly larger degree of the spread. In contrast, a greater spread is caused by variations in particle size, shape, and orientation.
In terms of robust process designs and adequate process optimisation for fluidisation and sedimentation of natural granules, it is therefore crucial to take into consideration the influence of the natural variations in the settling velocity when using predictive models of round spheres.
We use direct numerical simulations (DNS) to investigate the turbulent modulation due to the presence of bubbles in vertical channels flowing downward. The Reynolds number for single-phase flow based on half channel height h* and friction velocity is Reτ= 180. A density and viscosity ratio of ρd*/ρc*=0.01 and μd*/μc*=0.018 is chosen for two void fractions of ϵ=1.2% and ϵ=2.4%. For each void fraction three different bubble sizes are simulated: D/h=0.2130, 0.2684 and 0.3382, where D denotes the diameter of the bubbles. Numerical simulations are based on multiple markers Coupled Level-Set/Volume-of-Fluid (CLSVOF) method. To improve the efficiency of this method, a fast pressure-correction method is used in order to enable the simulation to exploit a constant coefficient Poisson equation which can be solved with FFT-based technique. Extensive verification and validation were performed and perfect accuracy and agreement are obtained. In all the simulations performed in this work, the new Poisson solver showed a minimum speedup of 22 times. Accumulation of bubbles in the core region of the channel for all cases is observed, which forms a bubble-free layer in the near-wall region. The presence of bubbles resulted in considerable modification in the mean velocity profile compared to single-phase flow. Another common observation is that all the components of velocity fluctuations in the near-wall region decrease with increasing void fraction and decreasing wall layer thickness. The opposite happens in the core region, where the presence of bubbles favours turbulence. With respect to the bubble size, the wall-normal and spanwise velocity fluctuations decrease in the near-wall region for smaller bubbles, however, the streamwise velocity fluctuations remained almost unaffected. The investigation of turbulent kinetic budgets shows that, unlike single-phase flow, the dissipation terms rises to large values in the core region of the channel. This behaviour is referred to the presence of bubbles and hence enhancement of turbulent kinetic energy in the core region.
The purpose of this work is to experimentally establish the combined influence on the flow and thermal resistance of an exhaust pipe wall formed by a porous, compliant layer with overlying discrete roughness elements exposed to the pulsating exhaust gas flow of a combustion engine. Through measuring the streamwise pressure drop over and radial temperature differences in different pipe samples for a range of flow states with different Reynolds numbers and non-dimensional pulsation frequencies, the effects were discerned. The configurations of the sample walls covered a range of mesh pitches, compliant-layer densities, and compliant-layer compression ratios. The (non-sinusoidally) pulsating exhaust gas flow spanned the following range: Reb (= ubD/νb) = 1⋅ 104 - 3⋅ 104, Tb = 500 - 800 ∘C, ω+(= ωνb/uτ2) = 0.003 - 0.040. The friction factors were found to be effectively constant with Reynolds number and non-dimensional pulsation frequency while the variation with insulation density/compression was not significant. Additionally, for both mesh pitches, the measured friction factors were in line with those reported in literature for similar geometries with steady flow and solid walls. Together this indicates that neither compliance nor the pulsations in the exhaust gas flow significantly affect the friction for this configuration. Comparison of the samples based on the derived thermal resistance showed a similar influence of the fluid-wall interface as for the friction. Additionally a distinct influence of compression, independent of the insulation density, was observed that increases with increasing temperature. It was concluded that the increased resistance was due to additional radiation resistance because of fibre reorientation due to compression.
Macroscopic properties of sedimenting suspensions have been studied extensively and can be characterized using the Galileo number (Ga), solid-to-fluid density ratio (πp) and mean solid volume concentration (ϕ¯). However, the particle–particle and particle–fluid interactions that dictate these macroscopic trends have been challenging to study. We examine the effect of concentration on the structure and dynamics of sedimenting suspensions by performing direct numerical simulation based on an Immersed Boundary Method of monodisperse sedimenting suspensions of spherical particles at fixed Ga= 144 , πp= 1.5 , and concentrations ranging from ϕ¯ = 0.5 to ϕ¯ = 30 %. The corresponding particle terminal Reynolds number for a single settling particle is ReT= 186. Our simulations reproduce the macroscopic trends observed in experiments and are in good agreement with semi-empirical correlations in literature. From our studies, we observe, first, a change in trend in the mean settling velocities, the dispersive time scales and the structural arrangement of particles in the sedimenting suspension at different concentrations, indicating a gradual transition from a dilute regime (ϕ¯ ≲ 2 %) to a dense regime (ϕ¯ ≳ 10 %). Second, we observe the vertical propagation of kinematic waves as fluctuations in the local horizontally-averaged concentration of the sedimenting suspension in the dense regime.
During the last decade, many approaches for resolved-particle simulation (RPS) have been developed for numerical studies of finite-size particle-laden turbulent flows. In this paper, three RPS approaches are compared for a particle-laden decaying turbulence case. These methods are, the Volume-of-Fluid Lagrangian method, based on the viscosity penalty method (VoF-Lag); a direct forcing Immersed Boundary Method, based on a regularized delta function approach for the fluid/solid coupling (IBM); and the Bounce Back scheme developed for Lattice Boltzmann method (LBM-BB). The physics and the numerical performances of the methods are analyzed. Modulation of turbulence is observed for all the methods, with a faster decay of turbulent kinetic energy compared to the single-phase case. Lagrangian particle statistics, such as the velocity probability density function and the velocity autocorrelation function, show minor differences among the three methods. However, major differences between the codes are observed in the evolution of the particle kinetic energy. These differences are related to the treatment of the initial condition when the particles are inserted in an initially single-phase turbulence. The averaged particle/fluid slip velocity is also analyzed, showing similar behavior as compared to the results referred in the literature. The computational performances of the different methods differ significantly. The VoF-Lag method appears to be computationally most expensive. Indeed, this method is not adapted to turbulent cases. The IBM and LBM-BB implementations show very good scaling.