Sedimentation of finite-size spheres in quiescent and turbulent environments
UUnder consideration for publication in J. Fluid Mech. Sedimentation of finite-size spheres inquiescent and turbulent environments
Walter Fornari † , Francesco Picano and Luca Brandt Linn´e Flow Centre and Swedish e-Science Research Centre (SeRC), KTH Mechanics,SE-10044 Stockholm, Sweden Department of Industrial Engineering, University of Padova, Via Venezia 1, 35131 Padua,Italy(Received ?; revised ?; accepted ?. - To be entered by editorial office)
Sedimentation of a dispersed solid phase is widely encountered in applications and envi-ronmental flows, yet little is known about the behavior of finite-size particles in homo-geneous isotropic turbulence. To fill this gap, we perform Direct Numerical Simulationsof sedimentation in quiescent and turbulent environments using an Immersed BoundaryMethod to account for the dispersed rigid spherical particles. The solid volume fractionsconsidered are φ = 0 . − ρ p /ρ f = 1 .
02. The par-ticle radius is chosen to be approximately 6 Komlogorov lengthscales. The results showthat the mean settling velocity is lower in an already turbulent flow than in a quiescentfluid. The reduction with respect to a single particle in quiescent fluid is about 12% and14% for the two volume fractions investigated. The probability density function of theparticle velocity is almost Gaussian in a turbulent flow, whereas it displays large positivetails in quiescent fluid. These tails are associated to the intermittent fast sedimentationof particle pairs in drafting-kissing-tumbling motions. The particle lateral dispersion ishigher in a turbulent flow, whereas the vertical one is, surprisingly, of comparable mag-nitude as a consequence of the highly intermittent behavior observed in the quiescentfluid. Using the concept of mean relative velocity we estimate the mean drag coefficientfrom empirical formulas and show that non stationary effects, related to vortex shedding,explain the increased reduction in mean settling velocity in a turbulent environment.
Key words:
1. Introduction
The gravity-driven motion of solid particles in a viscous fluid is a relevant process in awide number of scientific and engineering applications (Guazzelli & Morris 2012). Amongthese we recall fluvial geomorphology, chemical engineering systems, as well as pollutanttransport in underground water and settling of micro-organisms such as plankton.The general problem of sedimentation is very complex due to the high number of factorsfrom which it depends. Sedimentation involves large numbers of particles settling in dif-ferent environments. The fluid in which the particles are suspended may be quiescent orturbulent. Particles may differ in size, shape, density and stiffness. The range of spatialand temporal scales involved is wide and the global properties of these suspensions canbe substantially different from one case to another. Because of these complexities, our † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] N ov W. Fornari, F. Picano and L. Brandt general understanding of the problem is still incomplete.1.1.
Settling in a quiescent fluid
One of the earliest investigations on the subject at hand is Stokes’ analysis of the sedi-mentation of a single rigid sphere through an unbounded quiescent viscous fluid at zeroReynolds number. This led to the well-known formula that links the settling velocityto the sphere radius, the solid to fluid density ratio and the viscosity of the fluid thatbears his name. Later, the problem was studied both theoretically and experimentally.Hasimoto (1959) obtained expressions for the drag force exerted by the fluid on threedifferent cubic arrays of rigid spheres. These relate the drag force only to the solid volumefraction, but were derived under the assumption of very dilute suspensions and Stokesflow. The formulae have later been revisited by Sangani & Acrivos (1982). A differentapproach was instead pursued by Batchelor (1972) who found a relation between themean settling velocity and the solid volume fraction by using conditional probabilityarguments. When the Reynolds number of the settling particles ( Re t ) becomes finite,the assumption of Stokes flow is less acceptable (especially for Re t > φ of about 25%) and for low Reynolds numbers. Subsequentinvestigations improved the formula so that it could also be applied in the intermediateReynolds numbers regime (Garside & Al-Dibouni 1977; Di Felice 1999).Efficient algorithms and sufficient computational power have become available only rel-atively recently and since then many different numerical methods have been used toimprove our understanding of the problem (Prosperetti 2015). Among others we recallthe dynamical simulations performed by Ladd (1993), the finite-elements simulationsof Johnson & Tezduyar (1996), the force-coupling method simulations by Climent &Maxey (2003), the Lattice-Boltzmann simulations of Yin & Koch (2007), the Oseenletsimulations by Pignatel et al. (2011), and the Immersed Boundary simulations of Kempe& Fr¨ohlich (2012) and Uhlmann & Doychev (2014). Thanks to the most recent tech-niques it has become feasible to gain more insight on the interactions among the differentphases and the resulting microstructure of the sedimenting suspension (Yin & Koch 2007;Uhlmann & Doychev 2014). Uhlmann & Doychev (2014), most recently, simulated thesettling of dilute suspensions with particle Reynolds numbers in the range 140 −
260 andstudied the effect of the Archimedes number (namely the ratio between gravitational andviscous forces) on the microscopic and macroscopic properties of the suspension. Theseauthors observe an increase of the settling velocity at higher Archimedes, owing to parti-cle clustering in a regime where the flow undergoes a steady bifurcation to an asymmetricwake. Settling in stratified environments has also been investigated experimentally, i.e.by Bush et al. (2003), and numerically, i.e. by Doostmohammadi & Ardekani (2015). edimentation of finite-size spheres
Sedimentation in an already turbulent flow
The investigations previously reported consider the settling of particles in quiescent oruniform flows. There are many situations though, where the ambient fluid is in factnonuniform or turbulent. As in the previous case, the first approach to this problem wasanalytical. In the late 40’s and 50’s Tchen (1947) and later Corrsin & Lumley (1956)proposed an equation for the motion of a small rigid sphere settling in a nonuniformflow. In the derivation, they assumed the particle Reynolds number to be very low sothat the viscous Stokes drag for a sphere could be applied. The added mass and theaugmented viscous drag due to a Basset history term were also included. Maxey & Riley(1983) corrected these equations including also the Faxen forces due to the unsteadyStokes flow.In a turbulent flow many different spatial and temporal scales are active. Therefore thebehaviour and motion of one single particle does not depend only on its dimensions andcharacteristic response time, but also on the ratios among these and the characteristicturbulent length and time scales. The turbulent quantities usually considered are the Kol-mogorov length and time scales which are related to the smallest eddies. Alternatively,the integral lengthscale and the eddy turnover time can also be used. It is clear that aparticle smaller than the Kolmogorov legthscale will behave differently than a particle ofsize comparable to the energetic flow structures. A sufficiently large particle with a char-acteristic time scale larger than the timescale of the velocity fluctuations will definitelybe affected by and affect the turbulence. A smaller particle with a shorter relaxation timewill more closely follow the turbulent fluctuations. When particle suspensions are con-sidered, the situation becomes even more complicated. If the particles are solid, smallerthan the Kolmogorov lengthscale and dilute, the turbulent flow field is unaltered (i.e.,one-way coupling). Interestingly, the turbulent dynamics is instead altered by microbub-bles. The presence of these microbubbles leads to relevant drag reduction in boundarylayers and shears flows (e.g. Taylor-Couette flow)(Sugiyama et al. et al. et al. et al. et al. et al.
W. Fornari, F. Picano and L. Brandt in order to confirm these results and to study the turbulence modulation due to thepresence of particles (Hwang & Eaton 2006).The results on the mean settling velocities of particles of the order or larger than theKolmogorov scale are not conclusive. Good et al. (2014) studied particles smaller thanthe Kolmogorov scale and with density ratio O (1000), whereas Variano (experiments;private communication) and Byron (2015) studied finite-size particles at density ratioscomparable to ours. Good et al. (2014) found that the mean settling velocity is reducedonly when nonlinear drag corrections are considered in a one-way coupling approachwhen particles have a long relaxation time (a linear drag force would always lead to asettling velocity enhancement). For finite-size almost neutrally-buoyant particles, Vari-ano and Byron (2015) observe instead that the mean settling velocity is smaller than ina quiescent fluid. In relative terms, the settling velocity decreases more and more as theratio between the turbulence fluctuations and the terminal velocity of a single particlein a quiescent fluid increases. It is generally believed that the reduction of settling speedis due to the non-linear relation between the particle drag and the Reynolds number.Nonetheless, unsteady and history effects may also play a key role (Olivieri et al. et al. et al. (1995)tried to motivate these findings in terms of the relative motion between the fluid and theparticles. When the period of the fluid velocity fluctuations is smaller than the particleresponse time, a significant relative motion is generated between the two phases. Dueto the drag non-linearity, appreciable upward forces can be produced on the particlesthereby reducing the mean settling velocity.Unsteady effects may become important when considering suspensions with moderateparticle-fluid density ratios, as suggested by Mordant & Pinton (2000) and Sobral,Oliveira & Cunha (2007). The former studied experimentally the motion of a solid spheresettling in a quiescent fluid and explain the transitory oscillations of the settling velocityfound at Re ≈ O (100) by the presence of a transient vortex shedding in the particlewake. The latter, instead, analyzed an equation similar to that proposed by Maxey &Riley (1983), and suggested that unsteady hydrodynamic drags might become importantwhen the density ratio approaches unity.1.3. Fully resolved simulations
As already mentioned, most of the numerical studies of settling in turbulent flows used ei-ther one or two-way coupling algorithms. In order to properly understand the microscop-ical phenomena at play, it would be ideal to use fully resolved simulations. An algorithmoften used to accomplish this is the Immersed Boundary Method with direct forcing forthe simulation of particulate flows originally developed by Uhlmann (2005). The codewas later used to study the clustering of non-colloidal particles settling in a quiescentenvironment (Uhlmann & Doychev 2014). With a similar method Lucci et al. (2010)studied the modulation of isotropic turbulence by particles of Taylor length-scale size.Recently, Homann et al. (2013) used an Immersed Boundary Fourier-spectral methodto study finite-size effects on the dynamics of single particles in turbulent flows. Theseauthors found that the drag force on a particle suspended in a turbulent flow increases asa function of the turbulent intensity and the particle Reynolds number. We recently useda similar algorithm to examine turbulent channel flows of particle suspensions (Picano,Breugem & Brandt 2015).The aim of the present study is to simulate the sedimentation of a suspension of parti-cles larger than the Kolmogorov lengthscale in homogeneous isotropic turbulence with a edimentation of finite-size spheres ρ p /ρ f = 1 .
02) and investigate particle and fluid velocity statistics,non-linear and unsteady contributions to the overall drag and turbulence modulation.The suspensions considered in this work are dilute ( φ = 0 . − pdf of the settlingvelocities, we observe that the pdf is well approximated by a Gaussian function centeredaround the mean in the turbulent cases. In the laminar case instead, the pdf showsa smaller variance and a larger skewness, indicating that it is more probable to findparticles settling more rapidly than the mean value rather than more slowly. Theseevents are associated to particle-particle interactions, in particular to the drifting-kissing-tumbling motion of particle pairs. We also calculate mean relative velocity fields andnotice that vortex shedding occurs around each particle in a turbulent environment.Using the concept of mean relative velocity we calculate a local Reynolds number andthe mean drag coefficient from empirical formulas to quantify the importance of unsteadyand history effects on the overall drag, thereby explaining the reduction in mean settlingvelocity. In fact, these terms become important only in a turbulent environment.
2. Methodology
Numerical Algorithm
Different methods have been proposed in the last years to perform Direct Numerical Sim-ulations of multiphase flows. The Lagrangian-Eulerian algorithms are believed to be themost appropriate for solid-fluid suspensions (Ladd & Verberg 2001; Zhang & Prosperetti2010; Lucci et al. ∇· u f = 0 (2.1) ∂ u f ∂t + u f · ∇ u f = − ρ f ∇ p + ν ∇ u f + f (2.2)where u f , ρ f and ν = µ/ρ f are the fluid velocity, density and kinematic viscosity respec-tively ( µ is the dynamic viscosity), while p and f are the pressure and the force field usedto mantain turbulence and model the presence of particles. The particles centroid linearand angular velocities, u p and ω p are instead governed by the Newton-Euler Lagrangian W. Fornari, F. Picano and L. Brandt equations, ρ p V p d u p dt = ρ f (cid:73) ∂ V p τ · n dS + ( ρ p − ρ f ) V p g (2.3) I p d ω p dt = ρ f (cid:73) ∂ V p r × τ · n dS (2.4)where V p = 4 π a / I p = 2 ρ p V p a / a the particle radius; g is the gravitational acceleration; τ = − p I + 2 µ E is the fluidstress, with I the identity matrix and E = (cid:16) ∇ u f + ∇ u Tf (cid:17) / r is the distance vector from the center of the sphere while n is the unit vector normal tothe particle surface ∂ V p . Dirichlet boundary conditions for the fluid phase are enforcedon the particle surfaces as u f | ∂ V p = u p + ω p × r .In the numerical code the coupling between the solid and fluid phases is obtained byusing an Immersed Boundary Method. The boundary condition at the moving particlesurface (i.e. u f | ∂ V p = u p + ω p × r ) is modeled by adding a force field on the right-handside of the Navier-Stokes equations. The problem of re-meshing is therefore avoided andthe fluid phase is evolved in the whole computational domain using a second order finitedifference scheme on a staggered mesh. The time integration is performed by a third orderRunge-Kutta scheme combined with a pressure-correction method at each sub-step. Thesame integration scheme is also used for the Lagrangian evolution of eqs. (2.3) and (2.4).The forces exchanged by the fluid and the particles are imposed on N L Lagrangianpoints uniformly distributed on the particle surface. The force F l acting on the l − th Lagrangian point is related to the Eulerian force field f by the expression f ( x ) = (cid:80) N L l =1 F l δ d ( x − X l )∆ V l . In the latter ∆ V l represents the volume of the cell containing the l − th Lagrangian point while δ d is the Dirac delta. This force field is obtained throughan iterative algorithm that mantains second order global accuracy in space. Using thisIBM force field eqs. (2.3) and (2.4) are rearranged as follows to maintain accuracy, ρ p V p d u p dt = − ρ f N l (cid:88) l =1 F l ∆ V l + ρ f ddt (cid:90) V p u f dV + ( ρ p − ρ f ) V p g (2.5) I p d ω p dt = − ρ f N l (cid:88) l =1 r l × F l ∆ V l + ρ f ddt (cid:90) V p r × u f dV (2.6)where the second terms on the right-hand sides are corrections to account for the inertiaof the fictitious fluid contained within the particle volume. In eqs. (2.5),(2.6) r l is simplythe distance from the center of a particle. Particle-particle interactions are also consid-ered. When the gap distance between two particles is smaller than twice the mesh size,lubrication models based on Brenner’s asymptotic solution (Brenner 1961) are used tocorrectly reproduce the interaction between the particles. A soft-collision model is usedto account for collisions among particles with an almost elastic rebound (the restitutioncoefficient is 0 . et al. (2013), Lashgari et al. (2014) and Picano et al. (2015).2.2. Parameter setting
Sedimentation of dilute particle suspensions is considered in an unbounded computationaldomain with periodic boundary conditions in the x , y and z directions. Gravity is chosento act in the positive z direction. A zero mass flux is imposed in the simulations. A edimentation of finite-size spheres a . Non-Brownian rigid spheres are considered with solid to fluid density ratio ρ p /ρ f = 1 .
02. Hence, we consider particles slightly heavier than the suspending fluid. Twodifferent solid volume fractions φ = 0 .
5% and 1% are considered. In addition to the solidto fluid density ratio ρ p /ρ f and the solid volume fraction φ , it is necessary to introduceanother nondimensional number. This is the Archimedes number (or alternatively theGalileo number Ga = √ Ar ), Ar = (cid:16) ρ p ρ f − (cid:17) g (2 a ) ν (2.7)a nondimensional number that quantifies the importance of the gravitational forces actingon the particle with respect to viscous forces. In the present case the Archimedes number Ar = 21000. Using the particle terminal velocity v t we define the Reynolds number Re t = 2 av t /ν . This can be related by empirical relations to the drag coefficient of anisolated sphere when varying the Archimedes number, Ar . Different versions of theseempirical relations giving the drag coefficient as a function of Re t and Ar have beenproposed. As Yin & Koch (2007) we will use the following relations, C D = (cid:40) Re t (cid:104) . Re (0 . − . log Re t ) t (cid:105) , . < Re t (cid:54) Re t (cid:2) . Re . t (cid:3) , < Re t <
260 (2.8)since C D = 4 Ar/ (3 Re t ) (Yin & Koch 2007) we finally write, Ar = (cid:40) Re t (cid:104) . Re (0 . − . log Re t ) t (cid:105) , . < Re t (cid:54) Re t (cid:2) . Re . t (cid:3) , < Re t <
260 (2.9)The Reynolds number calculated from eq. (2.9) is approximately 188 for Ar = 21000.In order to generate and sustain an isotropic and homogeneous turbulent flow field, arandom forcing is applied to the first shell of wave vectors. We consider a δ -correlatedin time forcing of fixed amplitude ˆ f (Vincent & Meneguzzi 1991; Zhan et al. Re λ = λu (cid:48) /ν , where u (cid:48) is the fluctuating velocity and λ = (cid:112) νu (cid:48) /(cid:15) is thetransverse Taylor length scale. This is about 90 in our simulations. The ratio betweenthe grid spacing and the Kolmogorov lengthscale η = ( ν /(cid:15) ) / (where (cid:15) is the energydissipation) is approximately 1 . η . Finally, theratio among the expected mean settling velocity and the turbulent velocity fluctuationsis v t /u (cid:48) = 3 .
4. The parameters of the turbulent flow field are summarized in table 1. Forthe definition of these parameters, the reader is referred to Pope (2000).2.3.
Validation
To check the validity of our approach we performed simulations of a single sphere settlingin a cubic lattice in boxes of different sizes. Since this is equivalent to changing the solidvolume fraction, we compared our results to the analytical formula derived by Hasimoto(1959) and Sangani & Acrivos (1982), | V t | = | v t || V s | = 1 − . φ / (2.10)where | V s | = 29 (cid:16) ρ p ρ f − (cid:17) ga ν (2.11) W. Fornari, F. Picano and L. Brandt z V t V t = 1−1.7601*(2a/L)*( π /6) (1/3) Computed (Re t = 1)Yin & Koch, Re t = 1 Figure 1.
Terminal velocity of a periodic array of spheres. Present results Re t = 1 denoted byred circles are compared with the ones obtained by Yin and Koch by black squares at the same Re t and with the analytical solution for Re t = 0 of eq. (2.10). η/ (2 a ) u (cid:48) k λ/ (2 a ) Re λ (cid:15) T e Re L Table 1.
Turbulent flow parameters in particle units, where k is the turbulent kinetic energy, λ is the Taylor microscale, T e = k/(cid:15) is the eddy turnover time and Re L is the Reynolds numberbased on the integral lengthscale L = k / /(cid:15) . is the Stokes settling velocity. In terms of the size of the computational domain, l z ,eq. (2.10) can also be rewritten as V t = 1 − . a/l z )( π / / . In figure 1 we showthe results obtained with our code for Re t = 1 together with the data by Yin & Koch(2007) and the analytical solution. Although the analytical solution was derived with theassumption of vanishing Reynolds numbers, we find good agreement among the variousresults.The actual problem arises when considering particles settling at relatively high Reynoldsnumbers. If the computational box is not sufficiently long in the gravity direction, aparticle would fall inside its own wake (due to periodic boundary conditions), therebyaccelerating unrealistically. Various simulations of a single particle falling in boxes of dif-ferent size were preliminarily carried out, in particular 48 a × a × a , 32 a × a × a and 32 a × a × a . The first two boxes turn out to be unsuitable for our purposes.We find a terminal Reynolds number Re t = 200 in the longest domain considered, whichcorresponds to a difference of about 6% with respect to the value obtained from theempirical relations (2.9). As reference velocity we use the value obtained from simula-tions performed in the largest box at a solid volume fraction two orders of magnitudesmaller than the cases under investigation, φ = 5 · − (as in Uhlmann & Doychev2014), corresponding to a terminal velocity such that Re t = 195, 4% larger than thevalue from the empirical relations (2.9). Further increasing the length in the z directionwould make the simulations prohibitive. Note also that simulations in a turbulent envi-ronment turn out to be less demanding as turbulence disrupts and decorrelates the flowstructures induced by the particles. The final choice was therefore a computational boxof size 32 a × a × a with 256 × × φ = 0 . φ = 1%. In all cases, the particles are initially distributed randomlyin the computational volume with zero initial velocity and rotation. edimentation of finite-size spheres Figure 2.
Instantaneous snapshot of the velocity component perpendicular to gravity on dif-ferent orthogonal planes, together with the corresponding particle position for φ = 0 . A snapshot of the suspension flow for φ = 0 .
5% is shown in figure 2. The instantaneousvelocity component perpendicular to gravity is shown on different orthogonal planes.The simulations were run on a Cray XE6 system at the PDC Center for High Perfor-mance Computing at the KTH, Royal Institue of Technology. The fluid phase is evolvedfor approximately 6 eddy turnover times before adding the solid phase. The simulationsfor each solid volume fraction were performed for both quiescent fluid and turbulent flowcases in order to compare the results. Statistics are collected after an initial transientphase of about 4 eddy turnover times for the turbulent case and 15 relaxation times( T p = 2 ρ p a / (9 ρ f ν )) for the quiescent case. Defining as reference time the time it takesfor an isolated particle to fall over a distance equal to its diameter, 2 a/v t , the initialtransient corresponds to approximately 170 units. Statistics are collected over a timeinterval of 500 and 300 in units of 2 a/v t for the quiescent and turbulent cases, respec-tively. Differences between the statistics presented here and those computed from halfthe samples is below 1% for the first and second moments.
3. Results
We investigate and compare the behavior of a suspension of buoyant particles in aquiescent and turbulent environment. The behaviour of a suspension in a turbulent flowdepends on both the particles and turbulence characteristic time- and length-scales: ho-mogeneous and isotropic turbulence is defined by the dissipative, Taylor and integralscales, whereas the particles are characterized by their settling velocity v t and by theirStokesian relaxation time T p = 2 ρ p a / (9 ρ f ν ) = 11 . a/v t ) through-out the paper). A comparison between characteristic time scales is given by the Stokesnumber, i.e. the ratio between the particle relaxation time and a typical flow time scale St f = T p /T f . In the present cases, the Stokes number based on the dissipative scales(time and velocity) is St η = T p /T K = 8 . W. Fornari, F. Picano and L. Brandt p, τ /v t pd f quies. φ =0.5%quies. φ =1.0%turb. φ =0.5%turb. φ =1.0% 0 1 210 −2 −4 a) −1 −0.5 0 0.5 102468 V p,n /v t pd f quies. φ =0.5%quies. φ =1.0%turb. φ =0.5%turb. φ =1.0% −1 0 110 −2 −4 b) Figure 3.
Probability density functions of the particles velocities along the directions a) parallel V p,τ and b) perpendicular V p,n to gravity for φ = 0 .
5% and 1%. The velocities are normalizedby v t . Quiescent fluid cases denoted by blue curve with squares for φ = 0 .
5% and green curvewith triangles for φ = 1%, turbulent flow data by red curve with circles for φ = 0 .
5% and blackcurve with downward-pointing triangles for φ = 1%. In the insets we show the pdf s in lin-logscale. number St T e = T p /T e = 0 .
24. This value of St T e reflects the fact that the particles areabout 20 times smaller than the integral length scale L . The strong coupling betweenparticle dynamics and turbulent flow field occurs at scales of the order of the Taylorscale for the present cases. Indeed, the Taylor Stokes number is St λ = T p /T λ = 2 . T λ = λ/u (cid:48) . It should be noted that the Taylor length is slightly larger than theparticle size, 3 . a , while particles fall around 3.4 times faster than typical fluid velocityfluctuations, v t /u (cid:48) = 3 .
4. Hence particles are strongly influenced by the fluid fluctuationsoccurring at scales of the order of λ .3.1. Particle statistics
We start by comparing the single-point flow and particle velocity statistics for the twocases studied, i.e. quiescent and turbulent flow. The results are collected when a statisti-cally steady state is reached. Due to the axial symmetry with respect to the direction ofgravity, we consider only two velocity components for both phases, the components par-allel and perpendicular to gravity, V α,τ and V α,n respectively, where α = f, p indicatesthe solid and fluid phases.In figure 3 we report the probability density function of the particle velocities for bothvolume fraction investigated here, φ = 0 .
5% and 1%; the moments extracted from thesedistribution are summarized in table 2. The data in 3a) show the pdf for the componentof the velocity aligned with gravity V p,τ normalized by the settling velocity for φ → v t . This is extracted from the simulation of a very dilutesuspension discussed in section 2.3.In the quiescent cases, the mean settling velocity slightly reduces when increasing thevolume fraction φ , in agreement with the findings of Richardson & Zaki (1954) andDi Felice (1999) among others. The sedimentation velocity decreases to 0.96 at φ = 0 . φ = 1%. Di Felice (1999) investigated experimentally the settling velocity ofdilute suspensions of spheres ( φ = 0 . .
01 to 1000). Following the empiricalfit proposed in Di Felice (1999), we obtain (cid:104) V f,τ (cid:105) (cid:39) .
98 at φ = 0 . edimentation of finite-size spheres Re t < φ = 0 .
5% is instead close to theestimated value from eq. (2.9). This is in agreement with what reported in Uhlmann &Doychev (2014). These authors found that the particle mean settling velocity of a dilutesuspension increases above the reference value only when the Archimedes number Ar islarger than approximately 24000 and clustering occurs. In our case Ar (cid:39) (cid:104) V p,τ (cid:105) is 12% at φ = 0 .
5% and 14% at φ = 1%, see table 2.This result unequivocally shows that the turbulence reduces the settling velocity of asuspension of finite-size buoyant particles, in agreement with the experimental findingsby Variano and co-workers (Variano, private communication) and Byron (2015). We alsonote that the reduction of the settling velocity with the volume fraction is less importantfor the turbulent cases.Looking more carefully at the pdf s we note that fluctuations are, as expected, largerin a turbulent environment. In addition, the vertical particle velocity fluctuations arethe largest component in a quiescent fluid, whereas in a turbulent flow the fluctuationsare largest in the horizontal directions, as summarized in table 2. In the quiescent casethe rms of the tangential velocity σ V p,τ is 0.15 and 0.17 for φ = 0 .
5% and φ = 1%respectively, while it increases up to 0 .
23 in the corresponding turbulent cases. Thedifference in the width of the pdf is particularly large in the directions normal to gravitywhere the rms of the variance σ V p,n is 0.3 for both turbulent cases, while it is 5 and4 times smaller for the quiescent flows at φ = 0 .
5% and φ = 1%. We believe that theinteractions among the particle wakes, mainly occurring in the settling direction, promotethe higher vertical velocity fluctuations found in the quiescent cases. The shape of the pdf is essentially Gaussian for the turbulent cases, showing vanishing skewness and normalflatness. Interestingly an intermittent and skewed behaviour is exhibited in the quiescentcases. The flatness K is around 9 for both components at φ = 0 .
5% and slightly reducesto 5 . φ = 1%. The settling velocity of the quiescent cases is skewed towards intensefluctuations in the direction of the gravity. The skewness S is higher for the more dilutecase, being 1 .
26 at φ = 0 .
5% and 0 . φ = 1%.We interpret the intermittent behavior suggested by values of K > pdf s shown in figure 3a)are indeed associated to a specific behavior: as particles fall, they tend to be acceleratedby the wakes of other particles, before showing drafting-kissing-tumbling behavior (Fortes et al. (cid:104) V p,τ (cid:105) . In thequiescent case the fluid is still, the wakes are the only perturbation present in the fieldand are long and intense so their effect can be felt far away from the reference particle.The more dilute the suspension the more intermittent the particle velocities are. On thecontrary, when the flow is turbulent the wakes are disrupted quickly and therefore fewerparticles feel the presence of a wake. The velocity disturbances are now mainly due toturbulent eddies of different size that interact with the particles to increase the varianceof the velocity homogeneously along all directions leading to the almost perfect normaldistributions shown above, with variances similar to those of the turbulent fluctuations.The features of the particle wakes will be further discussed in this manuscript to supportthe present explanation.2 W. Fornari, F. Picano and L. Brandt
Quiescent φ = 0 .
5% Turbulent φ = 0 .
5% Quiescent φ = 1% Turbulent φ = 1% (cid:104) V p,τ (cid:105) +0 .
96 +0 .
88 +0 .
93 +0 . σ V p,τ +0 .
15 +0 .
23 +0 .
17 +0 . S V p,τ +1 .
26 +0 .
01 +0 .
70 +0 . K V p,τ +9 .
65 +2 .
92 +6 .
01 +3 . (cid:104) V p,n (cid:105) +0 . · − − . · − − . · − − . · − σ V p,n +0 .
06 +0 .
31 +0 .
08 +0 . S V p,n − . · − +0 . − . · − +0 . K V p,n +8 .
95 +2 .
78 +5 .
55 +2 . Table 2.
First four central moments of the probability density functions of V p,τ and V p,n normalized by v t . S V p,τ ( S V p,n ) and K V p,τ ( K V p,n ) are respectively the skewness and the flatnessof the pdf s. We also note that the sedimenting speed in the quiescent fluid is determined by twoopposite contributions. The excluded volume effect that contributes to a reduction of themean settling velocity with respect to an isolated sphere and the pairwise interactions(the drafting-kissing-tumbling), increasing the mean velocity of the two particles involvedin the encounter. To try to identify the importance of the drafting-kissing-tumbling effect,we fit the left part (where no intermittent behaviour is found) of the pdf pertaining thequiescent case at φ = 0 .
5% with a gaussian function. The mean of V p,τ is reduced toabout 0 .
93 (value due only to the hindrance effect) instead of 0.96 in the full simulation;thereby the increment in mean settling velocity due to drafting-kissing-tumbling can beestimated to be of about 3%.The probability density function of the particle angular velocities are also different inquiescent and turbulent flows. These are shown in figures 5a) and b) for rotations aboutan axis parallel to gravity, | ω p,τ | = | ω p,z | , and orthogonal to it, (cid:113) ω x + ω y . In the settlingdirection the peak of the pdf is always at | ω p,z | = 0. As for the translational velocities,the pdf s are broader in the turbulent cases. Due to the interaction with turbulent eddies,particles tend also to spin faster around axis perpendicular to gravity. From figure 5b)we see that the modal value slightly increases in the quiescent cases when increasingthe volume fraction. In the turbulent cases, the modal value is more than 3 times thevalues of the quiescent cases and the variance is also increased. Unlike the quiescentcases, the curves almost perfectly overlap for the two different φ , meaning that turbulentfluctuations dominate the particle dynamics. Turbulence hinders particle hydrodynamicinteractions.Figure 6 shows the temporal correlations of the particle velocity fluctuations, R v τ v τ (∆ t ) = (cid:104) V (cid:48) p,τ ( p, t ) V (cid:48) p,τ ( p, t + ∆ t ) (cid:105) σ V p,τ (3.1) R v n v n (∆ t ) = (cid:104) V (cid:48) p,n ( p, t ) V (cid:48) p,n ( p, t + ∆ t ) (cid:105) σ V p,n (3.2)for the turbulent and quiescent cases at φ = 0 .
5% and φ = 1%. Focusing on the dataat the lower volume fraction, we observe that the particle settling velocity decorrelatesmuch faster in the turbulent environment, within ∆ t ∼
50, while it takes around oneorder of magnitude longer in a quiescent fluid. Falling particles may encounter intensevortical structures that change their settling velocity. The turbulence strongly alters edimentation of finite-size spheres a) t=1789b) t=1794c) t=1799d) t=1801 Figure 4.
Drafting-kissing-tumbling behavior among two spherical particles in the quiescentcase with φ = 0 . W. Fornari, F. Picano and L. Brandt ω p,z | pd f quies. φ =0.5%quies. φ =1.0%turb. φ =0.5%turb. φ =1.0% a) ω p,x2 + ω p,y2 ) pd f quies. φ =0.5%quies. φ =1.0%turb. φ =0.5%turb. φ =1.0% b) Figure 5.
Probability density functions of a) | ω p,z | and b) (cid:112) ω p,x + ω p,y for φ = 0 .
5% and1%. The angular velocities are normalized by v t / (2 a ), the settling velocity of a single particlein a quiescent environment and its diameter. Quiescent fluid cases denoted by blue curve withsquares for φ = 0 .
5% and green curve with triangles for φ = 1%, turbulent flow data by redcurve with circles for φ = 0 .
5% and black curve with downward-pointing triangles for φ = 1%. ∆ t v t /(2a) R v n v n , R v τ v τ R v n v n R v τ v τ a) ∆ t v t /(2a) R v n v n , R v τ v τ R v n v n R v τ v τ b) ∆ t v t /(2a) R v n v n , R v τ v τ R v n v n R v τ v τ c) ∆ t v t /(2a) R v n v n , R v τ v τ R v n v n R v τ v τ d) Figure 6.
Time correlations of the particle velocity fluctuations. The blue curve with circlesrepresents the correlation of V (cid:48) p,τ ( R v τ v τ ), while the red curve with diamonds is used for thecorrelation of V (cid:48) p,n ( R v n v n ). a) Quiescent case with φ = 0 . φ = 0 . φ = 1 .
0% and d) Turbulent case with φ = 1 . the fluid velocity field seen by the particles, which in the quiescent environment is onlyconstituted by coherent long particle wakes. This results in a faster decorrelation of thevelocity fluctuations along the settling direction in the turbulent environment. Moreover, R v n v n crosses the null value earlier than for the settling velocity component. This resultis not surprising since the particle wakes develop only in the settling direction. edimentation of finite-size spheres R v n v n of the turbulent case oscillates around zerobefore vanishing at longer times. We attribute this to the presence of the large-scaleturbulent eddies. As a settling particle encounters sufficiently strong and large eddies,its trajectory is swept on planes normal to gravity in an oscillatory way. To provide anapproximate estimate of this effect, we consider as a first approximation the turbulentflow seen by the particles as frozen since the particles fall at a higher velocity thanthe turbulent fluctuations (3.4 times). Since the strongest eddies are of the order ofthe transversal integral scale we can presume that these structures are responsible forthis behavior. In particular, the transversal integral scale L T = L / (cid:39) u (cid:48) rms =0 .
3, so we expect a typical period of t = L T /u (cid:48) rms (cid:39)
26 which is of the order of theoscillations found for both turbulent cases, i.e. t ≈
20. Note that a similar behavior hasbeen observed by Wang & Maxey (1993) for sufficiently small and heavy particles, termedthe preferential sweeping phenomenon.The same process can be interpreted in terms of crossing trajectories and continuityeffects as described by Csanady (1963). An inertial particle falling in a turbulent envi-ronment changes continuously its fluid-particle neighborhood. It will fall out from theeddy where it was at an earlier instant and will therefore rapidly decorrelate from theflow. In order to accommodate the back-flow necessary to satisfy continuity, the normalcorrelations must then contain negative loops (as those seen in figures 6b and d). Follow-ing Csanady (1963) we define the period of oscillation of the fluctuations as the ratio ofthe typical eddy diameter in the direction of gravity (i.e. the longitudinal integral scale L ), and the particle terminal velocity v t obtaining t = L /v t (cid:39)
16. This value is similarto the period of oscillations in the correlations in figure 6.As shown in figure 6c), the quiescent environment presents a peculiar behavior of thesettling velocity correlation at φ = 1%. In particular we observe oscillations around zeroof long period, T = 160 − φ = 0 .
5% and 1 . (cid:104) V p,τ (cid:105) t , is subtracted from theinstantaneous displacement in the settling direction ∆ z ( t ) to highlight the fluctuationswith respect to the mean motion. For all cases we found initially a quadratic scaling intime ( (cid:104) ∆ x i (cid:105) ∼ t ) typical of correlated motions while the linear diffusive behavior takesover at longer times ( (cid:104) ∆ x i (cid:105) ∼ D e t with D e the diffusion coefficient).The turbulent cases show a similar behavior for both volume fractions. The crossovertime when the initial quadratic scaling is lost and the linear one takes place is about t (cid:39)
10 and t (cid:39)
50 for the normal and tangential component, respectively. This differenceis consistent with the correlation timescales previously discussed. The dispersion rates aresimilar in all directions in a turbulent environment. The quiescent cases present differentfeatures. First of all, dispersion is much more effective in the settling direction than inthe normal one. The dispersion rate is smaller than the turbulent cases in the horizontaldirections, while, surprisingly, the mean square displacement in the settling direction issimilar to that of the turbulent cases being even higher at φ = 0 . t (cid:39) W. Fornari, F. Picano and L. Brandt −1 −4 −2 t v t /(2a) < ∆ x i > turb. < ∆ x >turb. <( ∆ z−
Mean square particle displacement in the directions parallel and perpendicular togravity for both quiescent and turbulent cases for a) φ = 0 .
5% and b) φ = 1%. The mean particledisplacement (cid:104) V p,τ (cid:105) t in the settling direction is subtracted from the instantaneous displacementwhen computing the statistics. In the inset, we show (cid:104) ∆ x i (cid:105) /t as a function of time as well asthe diffusion coefficients D e for the turbulent cases. In the insets of figure 7a) and b) we show (cid:104) ∆ x i (cid:105) /t as a function of time. In all casesexcept the quiescent one at φ = 0 . D e = (cid:104) ∆ x i (cid:105) / (2 t ). For the turbulent case with φ = 0 . D e = 0 .
52 and 0 .
28 in the directions parallel and perpendicular to gravitywhile for φ = 1% we obtain D e = 0 .
57 and 0 .
32. In the quiescent cases the diffusioncoefficients in the horizontal directions are approximately 0 .
03, whereas the coefficientin the gravity direction at φ = 1% is about 0 .
40. Csanady (1963) proposed a theoreticalestimate of the diffusion coefficients for pointwise particles. Using these estimates, weobtain approximately D e = 1 . . . edimentation of finite-size spheres Quiescent φ = 0 .
5% Turbulent φ = 0 .
5% Quiescent φ = 1% Turbulent φ = 1 σ V f,τ +0 .
18 +0 .
28 +0 .
25 +0 . σ V f,n +0 .
04 +0 .
27 +0 .
06 +0 . Table 3.
Fluctuation rms of the fluid velocities parallel σ V f,τ and perpendicular σ V f,n togravity. The turbulent fluid velocity undisturbed rms is ∼ . Fluid statistics
Table 3 reports the fluctuation intensities of the fluid velocities for all cases considered.These are calculated by excluding the volume occupied by the spheres at each time stepand averaging over the number of samples associated with the fluid phase volume. Asexpected, the fluid velocity fluctuations are smaller in the quiescent cases than in theturbulent regime. In the quiescent case the rms of the velocity fluctuations is about 50%larger in the settling direction than in the normal direction because of the long rangedisturbance induced by the particles wakes. The increase of the volume fraction enhancesthe fluctuations in both directions. Fluctuations are always larger in the turbulent case,with the most significant differences compared to the quiescent cases in the normal direc-tion, where the presence of the buoyant solid phase brakes the isotropy of the turbulentvelocity fluctuations.Hence, the solid phase clearly affects the turbulent flow field. Although the presentstudy focuses on the settling dynamics, it is interesting to briefly discuss how turbulence ismodulated. Modulation of isotropic turbulence by neutrally buoyant particles is examinedin Lucci et al. (2010); however the results change due to buoyancy as investigated here.Typical turbulent quantities are reported in table 4 where they are compared with theunladen case at φ = 0. The energy dissipation (cid:15) increases with φ becoming almost doubleat φ = 1%. This behavior is expected since the buoyant particles inject energy in thesystem that is transformed into kinetic energy of the fluid phase that has to be dissipated.The higher energy flux, i.e. dissipation, is reflected in a reduction of the Kolmogorovlength η . The particles reduce the velocity fluctuations, decreasing the turbulent kineticenergy level. The combined effect on k and (cid:15) result in a decrease of the Taylor microscale λ and of Re λ ; likewise the integral length L and Re L also decrease. The reductionof the large and small turbulence scales is associated to the additional energy injectionfrom the settling particles. Energy injection occurs at the size of the particles, which isbelow the unperturbed integral scale L explaining the lowering of the effective integral L and of Taylor λ length-scales. This additional energy is transferred to the bulk flowin the particle wake. Associated to this energy input there is a new mechanism fordissipation that is the interaction of the flow with the no-slip surface of the particles.The mean energy dissipation field in the particle reference frame for the turbulent casewith φ = 0 .
5% is therefore shown in figure 8. After a statistically steady state is reached,the norm of the symmetric part of the velocity gradient tensor E ij and the dissipation (cid:15) = 2 νE ij E ij are calculated at each time step on a cubic mesh centered around eachparticle; the dissipation is calculated on the grid points outside the particle volume. Thedata presented have been averaged over all particles and time to get the mean dissipationfield displayed in the figure. The maximum (cid:104) ε (cid:105) is found around the particle surface withmaximum values in the front; the mean dissipation drops down to the values found inthe rest of the domain on the particle rear. The overall energy dissipation is thereforemade up of two parts: the first associated to the dissipative eddies far from the particlesurfaces and the second associated to the mean and fluctuating flow field near the particle8 W. Fornari, F. Picano and L. Brandt φ η/ (2 a ) k λ/ (2 a ) Re λ (cid:15) T e Re L Table 4.
Turbulent flow parameters in particle units for φ = 0, φ = 0 .
5% and φ = 1%. z/(2a) y / ( ) −2 −1 0 1 2−2−1012 0.050.10.150.20.250.30.35 Figure 8.
Contours of the mean fluid dissipation, (cid:104) ε (cid:105) ( v t / (2 a )), averaged in the particle frameof reference. surface. To conclude, the settling strongly alters the typical turbulence features via ananisotropic energy injection and dissipation, thus breaking the isotropy of the unladenturbulent flow. The energy is injected by the fluctuations in the particle wake whereasstronger energy dissipation occurs in the front of each particle. As a consequence, thefluid velocity fluctuations change in the directions parallel and perpendicular to gravityas shown in table 3. 3.3. Relative velocity
An important quantity to understand and model the settling dynamics is the particle tofluid relative motion. Although it is still unclear how to properly calculate the slip velocitybetween the two phases, we consider spherical shells around each particle, centered onthe particles centroids, inspired by the works of Bellani & Variano (2012) and Cisse et al. (2013). We calculate the mean difference between the particle and fluid velocities in eachshell as (cid:104) U rel (cid:105) x ,t,NP = (cid:42) u p − (cid:90) Ω(∆) u f d V (cid:43) t,NP where Ω(∆) is the volume of a shell of inner radius ∆. A parametric study on the slipvelocity is performed by changing the inner radii of these spherical shells from ∆ = 0 . . δ constant and equal to 0.063 in units of 2 a . In figure 9 we report the component of U rel parallel to gravity as a function of the shells inner radii ∆ / (2 a ). As the shell inner radiiincrease, | U rel,τ | tends exponentially to an asymptotic value which corresponds to themean particle velocity in the same direction, (cid:104) V p,τ (cid:105) . This is expected since the correlationbetween the fluid and particle velocity goes to zero at large distances. The quiescent casesstill show a 0 .
5% difference between | U rel,τ | and (cid:104) V p,τ (cid:105) at ∆ / (2 a ) = 5. This difference isagain attributed to the long coherent wakes of the particles. edimentation of finite-size spheres ∆ /(2a) < U r e l , τ > / v t Figure 9.
The component of U rel aligned with gravity as a function of the shell inner radii∆ / (2 a ) for the four cases studied. The dashed lines represent instead the mean particles velocitiesin the direction of gravity, (cid:104) V p,τ (cid:105) .∆ / (2 a ) (cid:104) U rel,τ (cid:105) σ U rel,τ S U rel,τ K U rel,τ .
75 +0 .
81 0 .
06 +8 .
040 123 . .
00 +0 .
88 0 .
09 +5 .
655 68 . .
25 +0 .
92 0 .
10 +4 .
547 47 . .
50 +0 .
93 0 .
11 +3 .
912 37 . .
75 +0 .
94 0 .
12 +3 .
482 31 . .
50 +0 .
95 0 .
13 +2 .
919 24 . .
00 +0 .
95 0 .
14 +2 .
18 16 . .
00 +0 .
96 0 .
14 +2 .
22 17 . Table 5.
Moments of the pdf ( U rel,τ ) for φ = 0 .
5% in the quiescent case. The thickness of theshell used to compute the slip velocity is δ/ (2 a ) = 0 . / (2 a ) (cid:104) U rel,τ (cid:105) σ U rel,τ S U rel,τ K U rel,τ .
75 +0 .
76 0 .
07 +0 .
072 3 . .
00 +0 .
83 0 .
09 +0 .
010 4 . .
25 +0 .
85 0 .
10 +0 .
087 4 . .
50 +0 .
86 0 .
11 +0 .
128 4 . .
75 +0 .
87 0 .
11 +0 .
117 4 . .
50 +0 .
87 0 .
13 +0 .
054 3 . .
00 +0 .
88 0 .
16 +0 .
046 3 . .
00 +0 .
88 0 .
18 +0 .
047 3 . Table 6.
Moments of the pdf ( U rel,τ ) for φ = 0 .
5% in a turbulent environment. The thicknessof the shell used to compute the slip velocity is δ/ (2 a ) = 0 . The probability density function of these relative velocities, U rel,τ , and their first fourcentral moments are computed and reported in tables 5 and 6 as a function of ∆ / (2 a )for the quiescent fluid and turbulent cases at φ = 0 . W. Fornari, F. Picano and L. Brandt rel, τ /v t pd f quies. φ =0.5%quies. φ =1.0%turb. φ =0.5%turb. φ =1.0%
50 100 150 200 250 300 0246Re p a) rel, τ /v t pd f quies. φ =0.5%quies. φ =1.0%turb. φ =0.5%turb. φ =1.0%
50 100 150 200 250 300 012345Re p b) Figure 10.
Comparison of the probability density functions of tangential relative velocity andcorresponding particle Reynolds number Re p . Colors and symbols are the same as previousfigures. In Panel a) the cases for ∆ / (2 a ) = 1 . / (2 a ) = 2 . velocities. The pdf s pertaining the four cases considered, calculated in spherical shellswith an inner radius of 1 . . Re p based on U rel,τ , is alsodisplayed in each figure. In the former case, ∆ / (2 a ) = 1 .
5, the shell radius is of theorder of the Taylor scale to highlight the particle dynamics, while the relative velocity isapproaching the asymptotic values for the larger shell. The pdf s of the relative velocityappear narrower than those of the particle absolute velocity, indicating that the particlestend to be transported by the large-scale motions, filtering the smallest scales.The distributions pertaining the simulations in a turbulent environment are nearlyGaussian with modal values well below one. The quiescent cases show skewed distribu-tions with long tails at high velocities, as observed for the particle velocities in figure 3a).The particles settle on average with a velocity close to that of a single particle, with occa-sional events of higher velocity due to the drafting-kissing-tumbling dynamics. The lowerthe volume fraction the more intermittent is the dynamics.Knowing the tangential fluid velocity (cid:104) V f,τ (cid:105) r ( p, t ), averaged in each shell and at eachtime step, and the corresponding tangential particle velocities V p,τ ( p, t ), it is then possibleto find their joint probability distribution P ( V f,τ , V p,τ ) (for sake of simplicity we write (cid:104) V f,τ (cid:105) r as V f,τ ). These are evaluated in shells at ∆ / (2 a ) = 1 .
75 for each case studied andreported in figure 11. In each case, the integral of P ( V f,τ , V p,τ ), (cid:104) V p,τ | V f,τ (cid:105) = (cid:90) ∞−∞ P ( V f,τ , V p,τ ) V f,τ dV f,τ , (3.3)is also reported (continuous lines). This represents the most probable particle velocity V p,τ given a certain fluid velocity V f,τ , or, equivalently, the most probable fluid velocitysurrounding a particle settling with velocity V p,τ . In the turbulent cases these integralsare well approximated by straight lines (displayed with dashed lines in the figure) V p,τ = C V f,τ + C . (3.4)In both cases C is approximately 1 while C is about 0 .
86 for φ = 0 .
5% and 0 .
84 for φ = 1%. These values are in agreement with the values found for the average relativevelocities of shells at ∆ / (2 a ) = 1 .
75. In a quiescent flow, conversely, the integral ineq. (3.3) gives a curved line and no best-fit is therefore reported. In these cases, thejoint probability distribution is broader, particularly in the region of higher particle edimentation of finite-size spheres V f, τ /v t V p , τ / v t −0.5 0 0.50.20.40.60.811.21.41.6 a) V f, τ /v t V p , τ / v t −0.5 0 0.500.511.5 b) V f, τ /v t V p , τ / v t −0.5 0 0.50.20.40.60.811.21.41.6 c) V f, τ /v t V p , τ / v t −0.5 0 0.500.511.5 d) Figure 11.
Joint probability distributions of V f,τ and V p,τ evaluated in spherical shells locatedat 1 .
75 particle diameters from each particle. We report in panels a) and c) the quiescent casesfor φ = 0 .
5% and 1%, while the respective turbulent cases are reported in panels b) and d). Thecontinuous lines represent the integrals of the joint probability distributions. In the turbulentcases the dashed lines represent the best-fit of these integrals. velocities, V p,τ . This is again due to the intense particle interactions and the drafting-kissing-tumbling behaviour described in figure 4, which confirms the high flatness of theprobability density functions of the relative particle velocities.Further insight can be obtained by plotting isocontours of the average particle relativevelocities and their fluctuations, in both quiescent and turbulent flows. To this end, wefollow the approach by Garcia-Villalba et al. (2012). We place an uniform and structuredrectangular mesh around each particle, with origin at the particle center. By means oftrilinear interpolations we find the fluid and relative velocities on this local mesh andaverage over time and the number of particles to obtain the mean relative velocity fieldand its fluctuations.The vertical velocity of the single sphere settling in a quiescent fluid is reported infigure 12. The relative normal and tangential velocities and their fluctuations are insteaddisplayed in figure 13 for suspensions with φ = 0 .
5% in both quiescent and turbulentenvironment. In the quiescent fluid simulations, a long wake forms behind the represen-tative particle and, as seen from the single point particle velocity correlations, it takes along time for this velocity fluctuations to decorrelate. In the turbulent case instead, thewakes are disrupted by the background fluctuations.Interesting observations can be drawn from the relative velocity fluctuation fields.Comparing figures 13c) and 13d) we note that intense vortex shedding occurs around theparticles in the turbulent case, with important fluctuations of U rel,τ . From figures 13e)2 W. Fornari, F. Picano and L. Brandt z/(2a) y / ( ) −6 −4 −2 0 2−2−1012 00.20.40.60.811.2 Figure 12.
Contour plot of the velocity component in the direction of gravity for the singlesphere settling in a quiescent fluid. and 13f) we also see that the relative velocity fluctuations are drastically increased in thehorizontal directions in a turbulent environment. Noteworthy, vortex shedding occurs atparticle Reynolds numbers below the critical value above which this is usually observed(Bouchet et al.
Drag analysis
As in Maxey & Riley (1983), we write the balance of the forces acting on a single spheresettling through a turbulent flow. The equation of motion for a spherical particle reads43 πa ρ p d V p dt = 43 πa ( ρ p − ρ f ) g + (cid:73) ∂ V p τ · n dS (3.5)where the integral is over the surface of the sphere ∂ V p , n is the outward normal and τ = − p I + 2 µ E is the fluid stress. As commonly done in aerodynamics, we replace thelast term of equation (3.5) with a term depending on the relative velocity U rel and adrag coefficient C D πa ρ p d V p dt = 43 πa ( ρ p − ρ f ) g − ρ f πa | U rel | U rel C D (3.6)with πa the reference area. Generally, the drag coefficient C D is a function of a Reynoldsnumber and a Strouhal number which accounts for unsteady effects. In the present case,we consider it to be a function of the Reynolds number based on the relative velocity Re p = 2 a | U rel | /ν (in a turbulent field it is proper to define Re p in terms of the relativevelocity between particle and fluid) and of the Strouhal number defined as St = d V p dt (2 a ) | U rel | . (3.7)The drag on the particle depends on both nonlinear and unsteady effects (such as theBasset history force and the added mass contribution) through these two non-dimensionalnumbers.Commonly, the unsteady contribution is neglected and C D is assumed to depend only edimentation of finite-size spheres z/(2a) y / ( ) −6 −4 −2 0 2−2−1012 00.20.40.60.81 a) z/(2a) y / ( ) −6 −4 −2 0 2−2−1012 00.20.40.60.81 b) z/(2a) y / ( ) −6 −4 −2 0 2−2−1012 00.050.10.150.20.25 c) z/(2a) y / ( ) −6 −4 −2 0 2−2−1012 00.050.10.150.20.25 d) z/(2a) y / ( ) −6 −4 −2 0 2−2−1012 0.020.040.060.08 e) z/(2a) y / ( ) −6 −4 −2 0 2−2−1012 00.10.20.30.4 f) Figure 13.
Fields of (cid:104) U rel,τ (cid:105) , U (cid:48) rel,τ , U (cid:48) rel,n for the quiescent (left column) and turbulent cases(right column). on the Reynolds number. Since we want to investigate both nonlinear and unsteady ef-fects, we decide to express the drag coefficient as C D ( Re p , St ) = C D ( Re p ) [1 + ψ ( Re p , St )],yielding43 πa ρ p d V p dt = 43 πa ( ρ p − ρ f ) g − ρ f πa | U rel | U rel C D ( Re p ) [1 + ψ ( Re p , St )] , (3.8)where ψ = 0 ∀ Re p if St = 0 (steady motion). We can therefore identify a quasi-steadyterm and the extra term which accounts for unsteady phenomena.By ensemble averaging equation (3.8) over time and the number of particles, andassuming the settling process to be at statistically steady state, we can find the mostimportant contributions to the overall drag. The steady-state average equation reads0 = 43 πa ( ρ p − ρ f ) g − ρ f πa (cid:104)| U rel | U rel C D ( Re p ) [1 + ψ ( Re p , St )] (cid:105) . (3.9)Denoting the entire time dependent term simply as Ψ( t ) and rearranging, we obtain the4 W. Fornari, F. Picano and L. Brandt following balance43 πa (cid:18) ρ p ρ f − (cid:19) g = 12 πa (cid:104)| U rel | U rel C D ( Re p ) (cid:105) + Ψ( t ) . (3.10)The term on the left-hand side is known whereas the time-dependent term Ψ( t ) is ofdifficult evaluation. The nonlinear steady term can be calculated using the approach de-scribed in section 3.2. At each time step we calculate the relative velocity in the sphericalshells surrounding each particle. From these we compute the particle Reynolds number Re p = 2 a | U rel | /ν and using equation (2.8) (where we replace the terminal Reynoldsnumber with the new particle Reynolds number) we obtain the drag coefficient. The firstterm on the right-hand side can therefore be evaluated by averaging over the numberof particles and time steps considered. Finally we divide everything by the buoyancyacceleration to estimate the relative importance of each term100% = S ( Re p ) + Ψ( t ) (3.11)where S ( Re p ) represents the nonlinear steady term while the unsteady term has beendenoted again as Ψ( t ) for simplicity.This approach is applied to the results of the simulations of a single sphere and to thequiescent and turbulent cases with φ = 0 .
5% and 1%. The inner radius of the samplingshells is chosen to be 5 particle radii and the results obtained are reported in table 7.The single sphere simulation provides an estimate of the error of the method. Since ourterminal Reynolds number is smaller than the critical Reynolds number above whichunsteady effects become important and our velocity field is indeed steady, the term Ψ( t )should be negligible. The critical Reynolds number Re c has been found to be approxi-mately 274 by Bouchet et al. (2006) while we recall that our terminal Reynolds numberfor the isolated particle is 200. The nonlinear steady term provides in this case, however,an overestimated value of the drag, with a percentage error of about +3% with respectto the buoyancy term. The possible causes of this error are the long wake and the factthat equation (2.8) is empirical. The data from the single-particle simulations are usedto correct the results from the other runs, i.e. the data are normalized with the totaldrag obtained in this case. In the quiescent case at φ = 0 . t ) adds up to approximately 10% of the total at φ = 0 .
5% and to about12% of the total at φ = 1%.Note that one can write the steady drag as mean and fluctuating component S ( Re p ) = (cid:104)| U rel | U rel C D ( Re p ) (cid:105) = (cid:104) U rel (cid:105) C D ( (cid:104) Re p (cid:105) ) + S (cid:48) ( (cid:104) Re p (cid:105) ) . The fluctuations S (cid:48) ( (cid:104) Re p (cid:105) ) would be responsible of the reduction of the settling velocityif this were to be attributed to nonlinear drag effects only (see also Wang & Maxey1993). We verified that for our results, the total and mean component differ by about2%, (cid:104)| U rel | U rel C D ( Re p ) (cid:105) ≈ (cid:104) U rel (cid:105) C D ( (cid:104) Re p (cid:105) ). This leads us to the conclusion thatthe main contribution to the overall drag is due to the steady term but the reductionof the mean settling velocity in a turbulent environment is almost entirely due to thevarious unsteady effects. These can be related to unsteady vortex shedding, see figure 13,as in the experiments of a single sphere by Mordant & Pinton (2000) These observationsare also in agreements with the results in Homann et al. (2013). These authors observethat the enhancement of the drag of a sphere towed in a turbulentt environment can beexplained by the modification of the mean velocity and pressure profile, and thus of theboundary layer around the sphere, by the turbulent fluctuations. edimentation of finite-size spheres case S ( Re p ) Ψ( t )Quiescent φ = 0 .
5% 99 .
5% 0 . φ = 0 .
5% 90 .
4% 9 . φ = 1 .
0% 94 .
1% 5 . φ = 1 .
0% 87 .
8% 12 . Table 7.
Percentage contribution of the non-linear and unsteady terms for the quiescent andturbulent case with φ = 0 .
4. Final remarks
We report numerical simulations of a suspension of rigid spherical slightly-heavy par-ticles in a quiescent and turbulent environment using a direct-forcing immersed bound-ary method to capture the fluid-structure interactions. Two dilute volume fractions, φ = 0 .
5% and 1%, are investigated in quiescent fluid and homogeneous isotropic turbu-lence at Re λ = 90. The particle diameter is of the order of the Taylor length scale andabout 12 times the dissipative Kolmogorov scale. The ratio between the sedimentationvelocity and the turbulent fluctuations is about 3.4, so that the strongest fluid-particleinteractions occur at approximately the Taylor scale.The choice of the parameters is inspired by the reduction in sedimentation velocityobserved experimentally in a turbulent flow by Byron (2015) and in the group of Prof.Variano at UC Berkeley. In the experiment, the isotropic homogeneous turbulence isgenerated in a tank of dimensions of several integral lengthscales by means of two facingrandomly-actuated synthetic jet arrays (driven stochastically). The Taylor microscaleReynolds number of the experiment is Re λ = 260. Particle Image Velocimetry usingRefractive-Index-Matched Hydrogel particles is used to measure the fluid velocity andthe linear and angular velocities of finite-size particles of diameter of about 1 . ρ p /ρ f = 1 .
02, 1 .
006 and 1 . v t and the turbulence fluctuating velocity u (cid:48) is about1, higher than in our simulations where this ratio is 3 .
3. Byron (2015) observes reductionsof the slip velocity between 40% and 60 % when varying the shape and density of theparticles. As suggested by Byron (2015), the larger reduction in settling velocity observedin the experiments is most likely explained by the larger turbulence intensity.The new findings reported here can be summarized as follows: i) the reduction of set-tling velocity in a quiescent flow due to the hindering effect is reduced at finite inertiaby pair-interactions, e.g. drafting-kissing-tumbling. ii) Owing to these particle-particleinteractions, sedimentation in quiescent environment presents therefore significant inter-mittency. iii) The particle settling velocity is further reduced in a turbulent environmentdue to unsteady drag effects. iv) Vortex shedding and wake disruption is served also insubcritical conditions in an already turbulent flow.In a quiescent environment, the mean settling velocity slightly decreases from thereference value pertaining few isolated particle when the volume fraction φ = 0 .
5% and φ = 1%. This limited reduction of the settling velocity with the volume fraction is inagreement with previous experimental findings in inertia-less and inertial flows. TheArchimedes number of our particle is 21000, in the steady vertical regime before theoccurrence of a first bifurcation to an asymmetric wake. In this regime, Uhlmann &Doychev (2014) observe no significant particle clustering, which is is confirmed by thepresent data.The skewness and flatness of the particle velocity reveal large positive values in a6 W. Fornari, F. Picano and L. Brandt quiescent fluid, and accordingly the velocity probability distribution functions displayevident positive tails. This indicates a highly intermittent behavior. In particular, it ismost likely to see particles sedimenting at velocity significantly higher than the mean:this is caused by the close interactions between particle pairs (more seldom triplets).Particles approaching each other draft-kiss-tumble while falling faster than the average.In a turbulent flow, the mean sedimentation velocity further reduces, to 0.88 and 0.86at φ = 0 .
5% and φ = 1%. The variance of both the linear and angular velocity increases ina turbulent environment and the single-particle time correlations decay faster due to theturbulence mixing. The velocity probability distribution function are almost symmetricand tend towards a Gaussian of corresponding variance. The particle lateral dispersion is,as expected, higher in a turbulent flow, whereas the vertical one is, surprisingly, of com-parable magnitude in the two regimes; this can be explained by the highly intermittentbehavior observed in the quiescent fluid.We compute the averaged relative velocity in the particle reference frame and thefluctuations around the mean. We show that the wake behind each particle is on averagesignificantly reduced in the turbulent flow and large-amplitude unsteady motions appearon the side of the sphere in the regions of minimum pressure where vortex shedding istypically observed. The effect of a turbulent flow on the damping of the wake behinda rigid sphere has been discussed for example by Bagchi & Balachandar (2003), whilethe case of a spherical bubble has been investigated by Merle et al. (2005). Using theslip velocity between the particle and the fluid surrounding it, we estimate the nonlineardrag on each particle from empirical formulas and quantify the relevance of non-stationaryeffects on the particle sedimentation. We show that these become relevant in the turbulentregime, amount to about 10-12% of the total drag, and are responsible for the reductionof settling velocity with the respect to the quiescent flow. This can be compared withthe simulations in Good et al. (2014) who attribute the reduction of the sedimentationvelocity of small (2 a < η ) heavy ( ρ p /ρ f ≈ Flowingmatter . REFERENCESAliseda, A, Cartellier, Alain, Hainaux, F & Lasheras, Juan C
Journal of Fluid Mechanics , 77–105.
Bagchi, Prosenjit & Balachandar, S
Physics of Fluids (1994-present) (11), 3496–3513. Balachandar, S & Eaton, John K
Annual Reviewof Fluid Mechanics , 111–133. edimentation of finite-size spheres Batchelor, GK
Journal of fluid mechanics (02), 245–268. Bec, J´er´emie, Homann, Holger & Ray, Samriddhi Sankar
Physical Review Letters (18),184501.
Bellani, Gabriele & Variano, Evan A
New Journal of Physics (12), 125009. Bergougnoux, Laurence, Bouchet, Gilles, Lopez, Diego & Guazzelli, Elisabeth
Physics of Fluids (1994-present) (9), 093302. Bosse, Thorsten, Kleiser, Leonhard & Meiburg, Eckart
Physics of Fluids(1994-present) (2), 027102. Bouchet, G, Mebarek, M & Duˇsek, J
European Journal of Mechanics-B/Fluids (3), 321–336. Brenner, Howard
Chemical Engineering Science (3), 242–251. Breugem, Wim-Paul
Journal of Computational Physics (13),4469–4498.
Bush, John WM, Thurber, BA & Blanchette, F
Journal of Fluid Mechanics , 29–54.
Byron, Margaret L arXiv preprint arXiv:1506.00478 . Ceccio, Steven L
Annual Review of Fluid Mechanics , 183–203. Cisse, Mamadou, Homann, Holger & Bec, J´er´emie
Journal of Fluid Mechanics , R1.
Climent, E & Maxey, MR
International journal of multiphase flow (4), 579–601. Corrsin, Se & Lumley, J
Applied Scientific Research (2), 114–116. Csanady, GT
Journal of theAtmospheric Sciences (3), 201–208. Di Felice, R
International journal of multiphase flow (4), 559–574. Doostmohammadi, A & Ardekani, AM
Physics of Fluids (1994-present) (2), 023302. Elgobashi, S
Appl.Sci. Res (3-4), 301–314. Fortes, Antonio F, Joseph, Daniel D & Lundgren, Thomas S
Journal of Fluid Mechanics , 467–483.
Garcia-Villalba, Manuel, Kidanemariam, Aman G & Uhlmann, Markus
International Journal of Multiphase Flow , 54–74. Garside, John & Al-Dibouni, Maan R
Industrial & engineering chemistry process designand development (2), 206–214. Good, GH, Ireland, PJ, Bewley, GP, Bodenschatz, E, Collins, LR & Warhaft, Z
Journal of Fluid Mechanics , R3.
Guazzelli, Elisabeth & Morris, Jeffrey F
A physical introduction to suspensiondynamics , , vol. 45. Cambridge University Press.
Gustavsson, K, Vajedi, S & Mehlig, B
Physical Review Letters (21), 214501.
Hasimoto, H W. Fornari, F. Picano and L. Brandt application to viscous flow past a cubic array of spheres.
Journal of Fluid Mechanics (02),317–328. Homann, Holger, Bec, J´er´emie & Grauer, Rainer
Journal of FluidMechanics , 155–179.
Hwang, Wontae & Eaton, John K st ∼
50) particles.
Journal of Fluid Mechanics , 361–393.
Johnson, Andrew A & Tezduyar, Tayfun E
Computer Methods in Applied Mechanics and Engineering (3),351–373.
Kempe, Tobias & Fr¨ohlich, Jochen
Journal of Computational Physics (9), 3663–3684.
Ladd, AJC & Verberg, R
Journal of Statistical Physics (5-6), 1191–1251.
Ladd, Anthony JC
Physics of Fluids A:Fluid Dynamics (1989-1993) (2), 299–310. Lambert, Ruth A, Picano, Francesco, Breugem, Wim-Paul & Brandt, Luca
Journal of Fluid Me-chanics , 528–557.
Lashgari, Iman, Picano, Francesco, Breugem, Wim-Paul & Brandt, Luca
Physical Review Letters (25), 254502.
Lucci, Francesco, Ferrante, Antonino & Elghobashi, Said
Journal of Fluid Mechanics , 5–55.
Maxey, Martin R & Riley, James J
Physics of Fluids (1958-1988) (4), 883–889. Merle, Axel, Legendre, Dominique & Magnaudet, Jacques
Journal of Fluid Mechanics ,53–62.
Mordant, N & Pinton, J-F
The EuropeanPhysical Journal B-Condensed Matter and Complex Systems (2), 343–352. Olivieri, Stefano, Picano, Francesco, Sardina, Gaetano, Iudicone, Daniele &Brandt, Luca
Physics of Fluids (1994-present) (4), 041704. Picano, Francesco, Breugem, Wim-Paul & Brandt, Luca
Journal of Fluid Mechanics , 463–487.
Pignatel, Florent, Nicolas, Maxime & Guazzelli, Elisabeth
Journal of Fluid Mechanics , 34–51.
Pope, Stephen B
Turbulent flows . Cambridge university press.
Prosperetti, Andrea
Journal of Fluid Mechanics , 1–4.
Richardson, JF & Zaki, WN
Chemical Engineering Science (2), 65–73. Sangani, AS & Acrivos, A
International journal of Multiphase flow (3), 193–206. Siewert, C, Kunnen, RPJ & Schr¨oder, W
Journal of Fluid Mechanics , 686–701.
Sobral, YD, Oliveira, TF & Cunha, FR
Powder Technology (2), 129–141.
Stout, JE, Arya, SP & Genikhovich, EL
Journal of the atmospheric sciences (22), 3836–3848. Sugiyama, Kazuyasu, Calzavarini, Enrico & Lohse, Detlef
Journal of Fluid Mechanics ,21–41.
Tchen, Chan-Mou edimentation of finite-size spheres Toschi, Federico & Bodenschatz, Eberhard
Annual Review of Fluid Mechanics , 375–404. Tunstall, EB & Houghton, G
Chemical Engineering Science (9), 1067–1081. Uhlmann, Markus
Journal of Computational Physics (2), 448–476.
Uhlmann, Markus & Doychev, Todor
Journal of Fluid Mechanics , 310–348.
Vincent, A & Meneguzzi, M
Journal of Fluid Mechanics , 1–20.
Wang, Lian-Ping & Maxey, Martin R
Journal of Fluid Mechanics ,27–68.
Yin, Xiaolong & Koch, Donald L
Physics of Fluids (1994-present) (9), 093302. Zhan, Caijuan, Sardina, Gaetano, Lushi, Enkeleida & Brandt, Luca
Journal of Fluid Mechanics , 22–36.
Zhang, Quan & Prosperetti, Andrea
Physics of Fluids (1994-present)22