Influence of Concentration on Sedimentation of a Dense Suspension in a Viscous Fluid

Macroscopic properties of sedimenting suspensions have been studied extensively and can be characterized using the Galileo number (Ga), solid-to-fluid density ratio (πp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi _p$$\end{document}) and mean solid volume concentration (ϕ¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{\phi }}$$\end{document}). 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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Ga=144$$\end{document}, πp=1.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi _p=1.5$$\end{document}, and concentrations ranging from ϕ¯=0.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{\phi }}=0.5$$\end{document} to ϕ¯=30%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{\phi }}=30\%$$\end{document}. The corresponding particle terminal Reynolds number for a single settling particle is ReT=186\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re_T = 186$$\end{document}. 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%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{\phi }} \lesssim 2\%$$\end{document}) to a dense regime (ϕ¯≳10%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{\phi }} \gtrsim 10\%$$\end{document}). 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.


Introduction
Sedimentation is the collective settling of particles suspended in a viscous fluid under the influence of gravity. It plays a major role in a multitude of industrial processes such as slurry transport, water treatment plants, land reclamation projects and fluidized beds. It 1 3 also plays a prominent role in the dynamics of environmental processes such as volcanic eruptions, sediment transport in rivers and rain.
Researchers have employed a multitude of techniques to study the dynamics of sedimenting suspensions, an overview is provided by Guazzelli and Morris (2011) and Davis and Acrivos (1985). Sedimentation of monodisperse suspensions of non-colloidal particles can be fully characterized using the Galileo number Ga = √ ( p − 1)gD 3 p ∕ 2 f , solid-to-fluid density ratio p = p ∕ f and mean solid volume concentration ̄ , where is the gravitational acceleration, p and f are the density of the solid and fluid, respectively, f is the dynamic viscosity of the fluid and D p is the diameter of the particle. For concentrations ̄> O(1%) , Davis and Birdsell (1988) describe that the settling of particles is hindered due to the increased vertical hydrostatic pressure gradient over the suspension, a net upward motion of the fluid and the nearby presence of other neighboring particles in their vicinity. As a result the average settling velocity of dense suspensions decreases for increasing concentration. In literature many semi-empirical laws have been proposed for the average settling velocity of a suspension ( V S ) as function of the terminal settling velocity of a single particle ( V T ) and bulk concentration described by Garside and Al-Dibouni (1977). A popular correlation is the one proposed by Richardson and Zaki (1954a, b) which predicts the average settling velocity of dense suspensions with good accuracy. Based on their experiments they demonstrated that the average settling velocity of a suspension exhibits a power-law relation in the bulk void fraction: where the exponent n is a function of the terminal particle Reynolds number Re T = V T D p ∕ f . Richardson and Zaki (1954a) determined n and V T by fitting Eq. 1 to their experimental results for various concentrations in the dense regime (typically ̄> 10% ). However, in later studies it was found that Eq. 1 underestimates the settling velocity for ̄< 10% . The value of V T obtained from fitting Eq. 1 to the dense regime appears to underestimate the actual value of V T for a single settling particle. Yin and Koch (2007) therefore proposed to multiply Eq. 1 with a certain correction factor of typically 0.86 − 0.92 when the actual value of V T is used. Furthermore, Garside and Al-Dibouni (1977) proposed an empirical and implicit logistic curve equation for the (relative) settling velocity that agrees well with experimental data over the entire range of ̄ , including the regime of ̄< 10%.
While Richardson and Zaki consistently found that the settling velocity of dense suspensions is always smaller than the settling velocity of a single sphere, it was recently found that very dilute suspensions may show enhanced settling velocities for certain Ga and p . This has been connected to stable vertical columnar structures within which particles settle at enhanced velocities. Heitkam et al. (2013) observed in their experiments at very low Reynolds number an increase of the sedimentation velocity of particles confined in circular capillaries. Experimental observation of columnar structures in dilute particle suspensions in homogenous turbulence was reported by Nishino and Matsuhita (2004). These columnar structures were also observed recently in a Direct Numerical Simulation (DNS) study by Uhlmann and Doychev (2014). They found that at Ga = 178 and p = 1.5 a dilute suspension with ̄= 0.5% settles faster than a single particle at the same Ga and p values. A single settling particle experiences instabilities in its wake as a function of Ga and p . As a result the path the particle travels may not be rectilinear and aligned with the direction of gravity. Path instabilities experienced by a single settling particle as a function of Ga and p have recently been studied in detail by Uhlmann and Dušek (2014) and Jenny et al. (2004). At Ga ≈ 178 and p ≈ 1.5 a particle has a (1) 1 3 relatively long planar wake and settles at a slightly oblique trajectory. Uhlmann and Doychev suggested that the lateral motion and the elongated wakes of particles settling in close proximity to one another may lead to one particle being captured in the wake of another particle, which subsequently may result in the formation of columnar structures in which particles settle at enhanced velocities. The exact mechanism of this phenomenon is not fully understood yet. An enhanced settling velocity was also found later in experiments by Huisman et al. (2016), who investigated concentrations up to ̄= 0.1% at Ga varying between 110 and 310 and p = 2.5 . However, in their experiments vertical clustering was observed even at Ga and p values where in case of a single particle no path instabilities were expected.
At dense concentrations, suspensions exhibit instabilities in the form of vertically propagating, horizontally oriented kinematic waves and shocks. Kynch theory provides an expression for the velocities of the kinematic shocks and waves as a function of concentration and sedimentation volume flux (Kynch, 1952). The development and propagation of these waves have been researched extensively in the context of fluidized beds by Kytömaa andBrennnen (1991), El-Kaissy andHomsy (1976), Batchelor (1988) and Jackson (1963). However, it has been challenging to relate the development of these instabilities to individual particle-particle and particle-fluid interactions.
The development of non-intrusive experimental techniques mentioned in Shih et al. (1986) and Williams et al. (1990), without the need for optical access has been able to provide insight into the dynamics of sedimenting suspensions. In addition, with the availability of computational resources and efficient numerical methods, it has been possible to provide a detailed account of the behavior of individual particles in these suspensions. In-depth descriptions of computational methods used to study particle laden flows are provided by Prosperetti and Tryggvason (2009) and more recently by Maxey (2017).
In this work, we aim to provide a description of gravity-driven monodisperse sedimentation of dense suspensions in a viscous fluid. The main questions we would like to address are: (1) how are the macroscopic properties of a dense sedimenting suspension related to particle-particle and particle-fluid interactions? and (2) how is this influenced by the concentration? We have performed DNS of gravity-driven sedimenting suspensions of non-colloidal spherical particles in a triply periodic computational domain, with the solid volume concentration varying from ̄= 0.5% to ̄ = 30% at a fixed Ga = 144 and p = 1.5 . The choice of Ga and p was motivated by a numerical study of Uhlmann and Dušek (2014) of a single settling particle in which they considered the same Ga and p values. At a comparable Ga = 121 , in the DNS performed by Uhlmann and Doychev (2014) no significant particle clustering was reported. This was confirmed by the DNS performed by Fornari et al. (2016) at Ga = 144 and concentrations ̄= 0.5% and ̄= 1% . While the previous works of Uhlmann and Doychev and Fornari et al. focused on the dynamics of very dilute suspensions ̄< O(1%) , we focus our analysis on the influence of concentration on the dynamics of sedimenting suspensions with an emphasis on dense suspensions for which ̄> O(1%).
We use an interface-resolved DNS based on an Immersed Boundary Method for the fluid/solid coupling described by Breugem (2012). In addition, a soft-sphere collision model described by Costa et al. (2015) is used for frictional particle collisions. We first describe the computational method and provide validation of the method. Next, we study the particle-particle and particle-fluid interactions by investigating the mean structural configurations and the average flow field around a particle as a function of concentration. Lastly, we focus our attention on investigating macroscopic trends in the average settling velocity, dispersion of particles within the suspension and development of kinematic waves.

Governing Equations
Fully resolved DNS is carried out in a triply periodic rectangular domain filled with a viscous fluid in which immersed non-colloidal spherical particles are allowed to settle under gravity. The two phases in the simulation (fluid and particulate) are treated independently and coupled through a no-slip boundary condition enforced on the surface of the particle. The solution to the fluid phase is computed on a fixed Eulerian mesh and the moving surface of the particle is represented using a Lagrangian mesh that translates with the particle. For the fluid/solid coupling the Immersed Boundary Method of Breugem (2012) is used. In this method the no-slip/ no-penetration condition at the particle/fluid interface is approximately enforced by (1) interpolating the provisional velocity in the fractional-step integration scheme from the Eulerian to the Lagrangian mesh, (2) calculating the force ( f IBM ) required to correct the fluid velocity to the local particle velocity, and (3) spreading this force from the Lagrangian to the Eulerian mesh. This approach is a blend of the regularized -function approach proposed by Peskin (1972) and the direct-forcing approach of Fadlun et al. (2000) initially described by Uhlmann (2005) and later improved to obtain second-order accuracy by Breugem (2012).
The fluid phase is governed by the incompressible Navier-Stokes equations and the par- where is the velocity, p is the modified pressure (total pressure -p h ). Here, ∇p h is the contribution to the hydrostatic pressure gradient from the submerged weight of the particles. For a homogenous suspension with concentration ̄ , ∇p h = ( p − 1)̄̃ and the gravitational acceleration g is non-dimensionalized as ̃ = ∕| | . The penultimate term in Eq. 3 becomes singular for p = 1 , but note that this corresponds to the case in which particles are not settling at all. The translational and angular velocities of monodisperse spherical particles are determined from the Newton-Euler equations, given by: where p and p are the translational and rotational velocity of the particle, is the stress tensor for a Newtonian fluid with I the unit tensor, is the position vector with respect to the particle centroid, is the outward normal directed from the surface ( V ) of the particle and c and c are the force and torque acting on the particle from inter-particle collisions. We consider non-colloidal particles in our simulations and hence the inter-particular interactions exclude electrostatic and Van der Waals forces. Brownian motion is neglected as well.

Collision Model and Lubrication Correction
A soft-sphere collision model described by Costa et al. (2015) is used to model frictional particle collisions. The collision model simulates a spring-damper interaction by allowing partial overlap between colliding entities. The collision force consists of a normal and tangential component. The normal repulsive component is represented by a spring-dashpot model, where k n and n represent the stiffness and damping coefficient, respectively, and n and u n are the overlap distance and the relative velocity between the particles in the normal direction respectively. The tangential force component is modeled with a spring-dashpot model in the stick regime and Coulomb's friction model in the slip regime, where k t , t and c are the stiffness, damping coefficient in the tangential direction and the coefficient of sliding friction respectively, and t and u t are the overlap distance and the relative velocity between the particles in the tangential direction. k n,t and n,t are determined from the reduced mass of the particles, the dry coefficient of restitution and a preset collision time described by Costa et al. (2015). In our simulations e n,dry = 0.97 , e t,dry = 0.10 and c = 0.15 based on experimental data for oblique glass particle-wall collisions in an aqueous glycerine solution discussed by Costa et al. (2015) and Joseph and Hunt (2004).
Lubrication effects are automatically accounted for in our DNS, although underresolved at inter-particle distances smaller than a grid cell. When prior to a collision the distance between the particles is lower than a threshold of the order of the mesh size and the particles are not colliding yet, a lubrication force correction is added to the rhs of Eq. 4. The (dimensional) lubrication correction is given by, where = 2 n ∕D p and x is the threshold gap between two particles, and is the Stokes amplification factor given by Jeffrey (1982). The collision model has been validated against several benchmark experiments and the results show a good quantitative agreement, see Costa et al. (2015) for details.

Numerical Method
The Navier-Stokes equations are solved by a fractional step approach and a three-step Runge-Kutta scheme is used for integration in time. The spatial discretization uses secondorder central finite differences on a uniform, staggered and isotropic grid. The Eulerian 1 3 mesh employs a Cartesian coordinate system where the y-axis is aligned with the vertical direction. Gravity is the only external force acting on the system and it is directed vertically downwards in the negative y-direction. The domain size was chosen to minimize the effect of periodic boundary conditions. The time taken for fluctuations in particle velocities to decorrelate was determined by computing the auto-correlation of particle velocities in each component direction, shown in Fig. 9. A vertical decorrelation distance was calculated as a product of the vertical decorrelation time scale and the mean settling velocity of the suspension. The evolution of the settling velocity of suspensions at different concentrations is shown in Fig. 1. We have checked a posteriori that the domain size in each direction was several times larger than the decorrelation distance, see fourth and last column in Table 1. The domain size for all the simulations (except ̄= 0.5% ) was fixed at 25D p × 100D p × 25D p . The domain size for the case ̄= 0.5% was chosen to be 37.5D p × 200D p × 37.5D p in order to account for a longer decorrelation distance and to increase the number of particles for statistical convergence. A grid resolution of D p ∕ x = 16 was chosen for our simulations and validation for it is provided in the next section. An illustration of the computational domain with particles is shown in Fig. 2. The  Table 1 Physical parameters in the DNS: global solid volume fraction ̄ , size of the computational domain L (in the coordinate direction x ), the number of particles N p , T obs is the duration of the observation interval over which statistics were calculated and the product V s L is the decorrelation distance. The Galileo number Ga = 144 and density ratio p = 1.5 in all simulations domain size, mean solid volume concentration, number of particles and observation time for calculating the statistics is provided in Table 1. At the start of the simulation, particles in the domain are initialized at random locations within the domain with zero velocity and allowed to settle under gravity. The fluid is initialized over the entire domain with zero velocity as well. After an initial transient of about t = 50 √ D p ∕g the particles settle at the mean settling velocity, after which statistics are collected. At every time step the hydrostatic pressure gradient in Eq. 3 is implicitly imposed from the requirement that the overall bulk flow (particle and fluid) has to be zero. This mimics the presence of a bottom wall in a batch sedimentation process.

Validation
The present numerical code is validated against experimental results by Ten Cate et al.
for settling of a sphere in a viscous fluid at Re T = 11.6 . The experiment was performed in a container of dimensions 100 × 100 × 160 mm 3 . The fluid used in the experiment is silicon oil with a density of 962 kg∕m 3 and dynamic viscosity 0.113 kg/ms and the solid used is a nylon sphere with density 1120 kg∕m 3 and diameter of 15 mm, which corresponds to Ga = 19.85 or Re T = 11.6 and p = 1.16 . The numerical simulation was set up to be similar to the experiment. A no-slip, no-penetration boundary condition is imposed on all 6 walls of the container. In Fig. 3, the solid line represents the computed result from the DNS and the blue dots represent the experimental data. The numerical results are found to be in good agreement with the experimental data. In addition, we performed a DNS of a single settling particle at Ga = 144 and p = 1.5 . A moving frame of reference with inflow/outflow conditions was used in the vertical direction and periodic boundary conditions were imposed in the horizontal directions. The  (2014) for a similar case in a cylindrical domain with a diameter of 5.34D p . Our terminal settling velocity was close to the one of Uhlmann and Dušek with an error of 0.45% . The particle was found to settle in a rectilinear fashion as expected. The vertical velocity field relative to the particle is shown in Fig. 6a. The terminal Reynolds number, Re T = V T D p ∕ f , is 186. This is also very close to the expected value of 184.5 obtained from Abraham's correlation for the drag coefficient (Abraham, 1970), indicating a negligible effect of the lateral domain size. The DNS has also been validated against several different flows previously by Breugem (2012) and Costa et al. (2015).

Snapshots of Instantaneous Particle Distribution and Velocity
From our simulations, we observe a distinct change in the structure and dynamics of sedimenting suspensions at different concentrations. We support our observations with a number of statistical correlations that demonstrate this change. Instantaneous snapshots of the computational domain at three different solid volume concentrations are shown in Fig. 4. We observe different structural arrangements of particles as the solid volume concentration is increased. At ̄= 2% we observe a tendency for vertical aggregation of particles. This can be observed in the trains of particles that settle significantly faster compared to the average settling velocity of the suspension, indicated by the red colored particles at the flanks of the domain. On moving to the denser solid volume concentrations of ̄= 10% and ̄= 20% the particles exhibit a seemingly random distribution, however it is hard to discern any structural trends from visual inspection alone. At ̄= 10% and ̄= 20% we observe a trace amount of particles traveling upwards (i.e., V∕V s < 0 , V > 0 and V s < 0 ). The distribution of colors (velocities), indicate the presence of fairly large scale structures in the sedimenting suspensions at ̄= 10% and ̄= 20%.

Particle-Conditioned Average Particle Distribution
The local average distribution of particles in the vicinity of each particle is studied by computing the particle-conditioned average. A solid phase-indicator function is defined over the entire domain. This function is defined to have a value of 1 within solid particles and 0 elsewhere. By averaging the solid phase-indicator function in the vicinity of each particle, the local particle-conditioned average is obtained.
The particle-conditioned average in a vertical plane passing through the center of the particle for different solid volume concentrations is shown in Fig. 5. The yellow colors in the plot indicate regions of higher than average concentration and the blue regions indicate regions of lower than average concentration. In Fig. 5, we observe that the regions away from the center of the particle show a concentration equal to the mean solid volume concentration. For ̄= 0.5% we notice that there is a tendency for vertical aggregation of particles indicated by the cone-shaped profiles in the vertical direction. We also observe a significantly lower than average concentration of particles in the regions adjacent to the vertical columns. In the case of ̄= 10% , we notice that there is a higher concentration of particles adjacent to the reference particle in the horizontal direction, but this anisotropy in particle arrangement is restricted to distances less than 2D p . This suggests Fig. 4 Instantaneous snapshot of a thin slab of the computational domain at solid volume concentrations a ̄= 2% , b ̄= 10% and c ̄= 20% . The slab thickness is 6.25 particle diameters. Particles are colored discreetly by their vertical velocity scaled with the mean sedimenting velocity ( V s < 0 ) of the suspension. The vertical trains of particles in the instantaneous snapshot for the case ̄= 2% are outlined using rectangles in a. Positive values of V∕V s indicate that particles are settling along the direction of gravity, i.e vertically downwards the increased probability of a particle to settle adjacent to a neighboring particle in the horizontal direction.
In the case of ̄= 30% , we notice first, spherical contours around the particle at radii of 1.16D p and 2.09D p and second, a slightly higher concentration in the horizontal direction. However, this effect is much weaker as compared to ̄= 10% . The concentric circles near the reference particle can be explained by the kinematic constraint that particles cannot overlap with the reference particle; this effect fades away within a few particle diameters. This results in a largely uniform hard sphere like distribution. A similar arrangement was also observed in the dense regime of sedimenting spherical particles at Re T ≪ 1 , shown using a radial distribution function by Guazzelli and Morris (2011), and using a pair probability distribution function by Yin and Koch (2007) at ̄= 20% for Re T = 10.

Average Flow Field Around a Particle
The average vertical flow field around a particle relative to the mean settling velocity of the suspensions is determined by computing the average vertical fluid flow field in two mutually perpendicular vertical planes centered around each particle. The mean settling velocity of the suspension is subtracted from the flow field to obtain the velocities relative to the particle's frame of reference. The mean is computed over all particles in the suspension and at 5 different instants in time over the course of the simulation. The averaged flow fields for the concentrations ̄= 2% and ̄= 30% as shown in Fig. 6b and c respectively. For comparison, the instantaneous flow field around a single particle is shown in Fig. 6a using a simulation with inflow/outflow conditions described in Sect. 2.4. Comparing Fig. 6b and c, we observe that the extent of the wake at ̄= 30% is much smaller than at ̄= 2% , while the latter is similar to the wake structure of a single particle. The presence of other neighboring particles in close proximity to the reference particle disrupts the formation of elongated wakes at ̄= 30% . This limits the influence of particle-fluid interactions and suggests a dominance of particle-particle interactions by lubrication and collisions, while the opposite holds for ̄= 2% as the wake in this case is similar to the case of a single settling particle.

Particle Velocity Statistics
From literature, it is known that the settling velocity of a suspension decreases with the increasing solid volume concentration, with the exception of very dilute suspensions at specific Ga and p such as mentioned in the introduction. For dense suspensions the settling velocity can be described by Eq. 1 suggested by Richardson and Zaki (1954a). Our simulations reproduce a similar trend for the solid volume concentrations simulated. We observe that Eq. 1 is valid for suspensions at solid volume concentrations ̄> O(10%) , while Richardson and Zaki underpredicts settling velocity at lower concentrations; the Fig. 6 a Relative vertical flow field around a settling particle for an isolated single particle (the red contour marks the location where the vertical fluid velocity V f = 0 , which indicates the extent of the recirculation region), (b) averaged flow field relative to a particle for a suspension with ̄= 2% , (c) Idem for ̄= 30% . Colors and contour lines in each figure span from 0 to the terminal settling velocity V T for the single particle and the mean settling velocity V s for concentrations ̄= 2% and ̄= 30% Fig. 7 Double-logarithmic plot of settling velocities of suspensions at different ̄ . The red dot indicates the terminal settling velocity V T of a single settling particle and the blue dot indicates the terminal settling velocity V T of a single settling particle computed from fitting the power relation proposed by Richardson and Zaki (1954a) and given in Eq. 1 to the settling velocities for ̄≥ 10% deviation increases for lower concentration in agreement with previous studies mentioned in the introduction section. In our DNS suspensions at lower solid volume concentrations show a deviation from this trend and settle faster than that predicted by the power-law relation. This can be seen in the double logarithmic plot of the settling velocity as a function of concentration, shown in Fig. 7.
The red dot in the figure is the terminal settling velocity V T = 0.91 √ gD p of a single particle settling under gravity. Richardson and Zaki (1954a) determined V T,ext of a particle settling under gravity by extrapolating Eq. 1 to a suspension at infinite dilution i.e. ̄= 0% . We compute V T,ext = 0.76 √ gD p (indicated in blue in Fig. 7) in our DNS using the same approach and the corresponding Re T = 154.3 . Note that V T,ext computed using this approach is different from the real V T of a single settling particle (indicated by the red dot in Fig. 7). From fitting Eq. 1 to the results for ̄≥ 10% we find that n = 3.0 . Richardson and Zaki (1954a) performed their experiments in 2 different pipes and found a clear dependency of the exponent n on the ratio of the particle to tube diameter. For Re T = 154.3 their proposed correlation for n varies between n = 2.69 for D pipe ∕D p = ∞ and n = 3.12 for D pipe ∕D p = 25 . In our simulations the ratio of particle diameter to the lateral extent of the domain is equal to 25. Note, however, that we make use of periodic boundary conditions, so we expect that the value from our DNS is in between these limits. This is indeed the case as our measured value for n = 3 is in between the 2 limits mentioned before : 2.69 < n < 3.12.

Particle Dispersion Statistics
The dispersion of the particles in the suspension is measured from the mean square displacement of the particle as a function of time. The expression to compute the mean square displacement is given by: where, Y is the displacement, is time interval over which the displacement is measured and t is the simulation time over which the mean is computed. The term −V s is a correction for the mean vertical displacement over a time interval . The overline represents an average over t, the time over which statistics were computed and over the displacements of all the particles in the vertical direction. The expression for the mean square displacement provided in Eq. 9 is analogous for the lateral directions, but without −V s term. Einstein (1956) predicted ballistic and diffusive transport at short and long times, respectively, for Brownian motion of small particles. Similarly, for sedimenting suspensions the particles are in the ballistic regime for short times, where the mean square displacement scales quadratically with time, while for long times they are in diffusive regime, where the mean square displacement scales linearly with time. In turbulent flows dispersion of a passive scalar also displays ballistic and diffusive regime as shown by Taylor (1921) similar to Einstein's theory of Brownian motion. The time scale = L at which the two regimes (indicated by the linear and quadratic fits) intersect is the integral time scale. This is indicated by the red line in Fig. 8 for the case of ̄= 2% . The integral time scale marks the transition from the ballistic to the dispersive regime. The integral time scale is a measure of the typical time over which the particle velocity decorrelates with itself or, alternatively, is a typical diffusive time scale. The integral time scale associated with each solid volume concentration is computed and expressed as a function of ̄ in Fig. 9. It can be observed from this plot that the rate of diffusion increases with increasing concentration up to ̄∼ 6% and 10% in the horizontal and vertical direction, respectively, while it remains nearly constant for higher concentrations. The dispersive time scale for the vertical direction is larger than for the lateral directions. Because of symmetry, the dispersive time scales for the x and z directions should be identical. This can indeed be observed, except for some discrepancies in the smallest concentrations related to some lack of statistical convergence.

Kinematic Waves
Kynch theory of sedimentation makes use of two assumptions (Kynch, 1952). First, the concentration of a sedimenting suspension is assumed to be uniform in the lateral directions. This enables the distribution of the local concentration of the sedimenting suspension to be expressed as a function of vertical position and time (y, t) . Second, the settling velocity of the suspension is assumed to be quasi-steady only dependent on the local concentration. Using these assumptions, mass conservation requires that:  To study this, we calculate the plane-averaged solid volume concentration as function of the vertical height and time. This is shown in Fig. 11 for the case ̄= 20%.
From this space/time plot, we observe local fluctuations in the vertical concentration profile around the global solid volume concentration (indicated by the blue and yellow bands) that show a wave-like pattern that propagates downwards, though at a lower velocity than the average settling velocity.
In order to measure the velocity of the wave-like structures in the sedimenting suspension, we performed an autocorrelation of the plane averaged solid volume concentration. The autocorrelation is given by: where � (y, t) = (y, t) −̄ is the fluctuation in the plane-averaged solid volume concentration. The overline represents an average over y and t. The correlations over large displacements in space and time are expected to correspond to kinematic waves. By means of a linear fit through the correlation peaks over large displacements in y and t, the velocity of the correlated structure (i.e. kinematic wave) can be determined. The velocity of the kinematic waves V KW corresponding to the cases ̄= 10% to ̄= 30% are computed in a similar manner. From Kynch sedimentation theory, the velocity of the kinematic waves V KW at a particular concentration can be determined by computing the slope of the sedimentation flux curve in Fig. 10, given by: We compare the kinematic wave velocities at different concentrations from the DNS and theory, shown in Fig. 12. The measured results from the DNS follow a similar trend as predicted by the theory. Deviations between the DNS data and Eq. 13 may be related to dispersive effects not accounted for in Eq. 10, as discussed by Jackson (1963).
(12) ( y, t) = � (y, t) � (y + y, t + t) Comparison of kinematic wave velocity of a sedimenting suspension as a function of concentration between theory (black dashed line) with the measured sedimentation flux from DNS (red dots)

Conclusions and Discussion
Particle-particle and particle-fluid interactions influence the macroscopic properties of a sedimenting suspension and are a strong function of the global particle volume concentration. From our results we observe for increasing concentration a gradual transition from a dilute ( ̄≲ 2% ) to a dense regime ( ̄≳ 10% ). From the instantaneous snapshots of the suspensions at different concentrations and the particle-conditioned concentration, we notice that structural arrangement of particles is different in each regime. From the conditionally averaged particle concentrations, we observe with increasing ̄ a gradual transition from a preference for vertical particle aggregation in the dilute regime to a uniform distribution in the dense regime with a slightly higher preference of a particle to settle adjacent to another particle. A particle settling in the wake of another particle experiences less drag and hence tends to draft towards that particle. The trailing particle can either spend an extended duration of time drafting towards the leading particle or can come into contact (kissing) with the leading particle and tumble into a more stable horizontal configuration. The latter is known as the drafting, kissing and tumbling (DKT) mechanism. The former is expected to play a significant role in the dilute regime where the wake of a settling particle is elongated and relatively undisturbed by the presence of other particles as compared to the dense regime. This leads to vertical aggregation of particles which we observed at the most dilute concentrations considered, i.e. ̄= 0.5% and ̄= 2% . At ̄> 2% , the slight increase in the horizontal concentration close to the reference particle observed in the particle-conditioned concentration can be explained by the DKT mechanism. The negative velocities of a few particles observed in the dense regime could be due to the presence of particle-rich and particle-poor regions that set up local convection of fluid causing a few particles in the latter region to travel upwards. Lubrication and collisions likely play an important role in the dynamics of the dense regime, though a more quantitative analysis of their influence is required.
In tandem with this change, we observe first, a deviation of the mean settling velocity from the power-law relation suggested by Richardson and Zaki (1954a) in the dilute regime and second, the dispersive time scales were found to decrease rapidly up to ̄∼ 10% , while they remained relatively unaltered in the dense regime. This highlights the change in dynamics from the dilute to the dense regime. Yin and Koch (2007) suggested that the power-law scaling of the settling velocity is associated with a hard sphere like (random) particle distribution. In our simulation we study the structural arrangement of particles in the suspension by means of the particle-conditioned concentrations and we can confirm that the power-law scaling is indeed associated with such a distribution as observed in the dense regime. The statistically homogeneous distribution of particles could be a possible reason for the good agreement of our results with the power-law relation provided by Richardson and Zaki, while the vertical arrangement of particles in the dilute regime promotes a settling velocity higher than predicted by the power-law relation for a homogenous suspension.
We observed the presence of kinematic waves in the dense regime and we computed their velocities by means of an autocorrelation of the fluctuations in local vertical concentration of the sedimenting suspension. The results are found to be in good agreement with that predicted by Kynch theory. The origin of this agreement can be attributed to two factors: an agreement with the assumptions of the theory and a trigger mechanism to initiate the propagation of kinematic waves. First, in the dense regime, particles exhibit a hard sphere distribution which is largely uniform and any local anisotropy in the particle concentration is limited to distances less than 2D p . Second, the settling velocity of the sedimenting suspensions are dependent on the local-concentration as indicated by the good agreement with the power-law relation presented by Richardson and Zaki, Eq. 1. We speculate that the slight increase in local horizontal concentration of particles could be the trigger that initiates the propagation of kinematic waves. Though the agreement with Kynch theory is good, the theory does not account for dispersion of particles within the suspension. It would therefore be interesting to extend this investigation to a more detailed analysis such as proposed by Jackson (1963).