Direct numerical simulations of sedimenting spherical particles at non-zero Reynolds number

Abstract

We performed direct numerical simulations, using a smoothed profile method to investigate the inertial effects on the static and dynamic properties of a sedimenting suspension over a wide range of volume fractions from 0.01 to 0.4. We found that at Reynolds number Re ≤ 0.5, static and dynamic properties show the typical non-Brownian, Stokes regime characteristics, due to insignificant inertial effects. The microstructure analysis at the high Re revealed that at Re = 1 inertial forces have significant effects and these create a deficiency of particles around a given particle, which is more pronounced in the direction of gravity than in the perpendicular direction. This deficiency decreased the velocity fluctuations and particle diffusion in the vertical direction, whereas both of these properties remain unchanged in the perpendicular direction. Moreover, at Re = 10, strong inertial forces generated a significant deficit of particles in both directions, which decreased velocity fluctuations and particle diffusion in both directions. We also observed that the range of volume fraction affected by inertial forces is increased with the increase of Re. At high volume fraction ϕ ≳ 0.15, intrinsic many-body interactions dominate the phenomena and govern the transport properties thereafter.

---

1. Introduction

In general sedimentation is the settling through fluid of suspended material of various sizes, ranging from large rocks to proteins and peptides. Sediments pose a big problem in industrial drainage, from which they must be continuously removed. Along with the other separation techniques, sedimentation is an essential, cheap and widely used clarification process. Most of the experimental,^1–8 theoretical^9–11 and simulation^12–31 studies on sedimentation are focused in the Stokes regime, where theoretical aspects of the phenomena are well developed, whereas little is known about the behavior of the suspension when fluid inertia is significant. This is because turbulent particle laden flows are difficult to investigate through experiments and simulations, as the former require sophisticated instrumentation and the latter enormous computational resources. Moreover, this apparently simple technique is strongly affected by the host fluid and particle parameters as well as the configuration of the container, which make this phenomena difficult to tackle. Microscopically, it is a typical example of non-equilibrium dynamics, which shows interesting yet perplexing static and dynamic properties.

The first significant theoretical work on sedimentation was performed by Stokes,³² who calculated the settling velocity of a single particle in an incompressible fluid and in the absence of appreciable inertia, known as the Stokes velocity. Later, Richardson and Zaki³³ found the effect of particle volume fraction and fluid back flow on the settling velocity, which make particles undergo fluctuating motions. Researchers characterize this fluctuating motion of particles mainly by calculating velocity fluctuations, their relaxations times and particle diffusion. Since these properties are largely dependent on the particle and fluid parameters, along with the size and shape of the container, the accurate measurement of these properties is important to get an insightful picture of the physical phenomena.

Unlike particles undergoing Brownian motion, sedimenting particles show different behavior in the direction parallel and perpendicular to gravity. Caflisch and Luke³⁴ predicted that the intensity of velocity fluctuations diverges with increasing system size, which is certainly unphysical. Hinch³⁵ solved this paradox by considering a hypothetical blob of fluid of size l, having n_dl^1/3 average number of particles (n_d is the number density of the particles), with statistical fluctuations in particle number. If m is the mass of the particles which balance the buoyancy, then the fluctuations in the weight are . Balancing this weight with the drag force, Hinch came up with the scaling of velocity fluctuations as , where ϕ and a denote the particle volume fraction and radius, respectively.

These predictions stimulated many experimental studies. Batchelor³⁶ investigated the bidispersed suspension, where particles lighter than the fluid rose and those heavier sank. Later, Nicolai performed a series of experiments^1–3 to investigate the effects of particle volume fraction, system size and polydispersity on velocity fluctuations and particle diffusion. These studies quantified the transport properties, but gave little information about the physics. Guazzelli³⁷ investigated the effects of walls and found that these walls act as kinetic traps and reduce the velocity fluctuations of the particle in the near vicinity. Segrè,^4–7 in a series of experiments, also investigated the effects of particle volume fraction ϕ on velocity fluctuations at low volume fraction and found that these velocity fluctuations scale as ϕ^1/3, in contrast to the theoretical prediction of Hinch³⁵ and experimental results of Bernard.³⁸ This contradiction between the theory and experiments was solved by Brenner³⁹ through a theoretical investigation, scaling arguments and numerical simulations. He found that vertical walls play a critical role in this scaling and these reduce the velocity fluctuations. He recalculated the velocity fluctuations of Segrè's work by considering the particles not affected by the side walls and found the ϕ^1/2 scaling, in good agreement with the theoretical prediction.³⁵ The side walls have significant effects in the Stokes regime, whereas Koch et al.⁴⁰ found a weaker logarithmic scaling when particle inertia is significant. This suggests that bounding solid walls which suppress the velocity fluctuations in the Stokes regime, are relatively less pronounced when inertial forces are significant.⁴⁰

The theoretical aspects of sedimentation phenomena in the Stokes regime are well developed, where velocity fluctuations scale as ϕ^1/2 at low volume fraction^21,39,41 and as (L/a)^1/2 with the system size, whereas these properties are not explored extensively at non-zero Reynolds number (Re), where Re is the ratio of inertial forces to viscous forces. For a given particle diameter σ, and fluid density ρ_f and viscosity η, Re is defined as Re = ρ_fV_tσ/η, where V_t is the terminal velocity of a single particle. Koch^40,42–44 explored the physical aspects of sedimentation at non-zero Re and found that at low volume fraction and in the Stokes regime, the flow field produced by a settling sphere has a fore-aft symmetry relative to the horizontal plane through its center, which is broken by the fluid inertia at non-zero Re. This breaking of symmetry generates a wake behind the particle which convects fluid towards its rear and this fluid flux is balanced by the source flow. If a second particle is in the wake of the leading particle then the trailing particle will drift to the leading particle and it also experiences a lift force due to the shear flow in the wake which pushes it outward and cause a particle deficiency in the vicinity of a test particle. If this lift force is small, the two particle will kiss and then tumble into a horizontal orientation and repel each other due to the source flow. This is known as the drafting–kissing–tumbling (DKT) mechanism. Similarly, horizontally oriented particles will repel each other due to the source flow. At high volume fractions, this weak interaction competes with the many-body interactions, which then lead to the isotropic microstructure, similar to that in a hard-sphere distribution. Climent et al.⁴⁵ observed a similar mechanism in their study as well.

Under sedimentation, the microstructure of the suspension is determined by the long-range hydrodynamic interaction (HIs), which are characterized by large time and length scales, features that render the non-equilibrium properties of the colloidal particles difficult to tackle both experimentally and by simulations. Koch⁴⁰ performed lattice Boltzmann simulations to investigate the microstructure changes with the increase in Re and its effects on transport properties at four different volume fractions and system sizes, whereas Climent et al.⁴⁵ used an approximate numerical method, where the particle velocity is calculated by averaging the fluid velocity in a Gaussian envelop to investigate the sedimentation at finite Re. At low Re, where inertia can be ignored, the Navier–Stokes equations (NS) are reduced to a linear set of equations and are easy to tackle computationally, whereas with the introduction of appreciable inertia, the treatment of the nonlinear equations that govern the flow is very difficult. Computational schemes which use unstructured meshes⁴⁶ usually provide good results in the Stokes regime, whereas at non-zero Re, these require very heavy computational resources because of the fast phenomenological changes. Hence, it is believed that the best way to investigate the non-equilibrium properties of a suspension with appreciable inertia is through direct numerical simulations (DNS), with a fixed grid around the particles.

Despite the extensive work on sedimentation as briefly reviewed above, there are still some open problems which need to be investigated, e.g., characterization of correlation length and finite size effects and interplay of thermal and hydrodynamic forces especially at high Peclet number (Pe). Moreover, there is also a lack of conclusive studies on effects of non-zero Re, system and particle size and shape and evolution of the particle microstructure and its effects on transport properties of sedimenting suspension. In our previous studies, we have already investigated the interplay of thermal and hydrodynamic forces^47–49 at Pe ≤ 125, the microstructure, and transport properties of non-Brownian particles,⁴¹ as a function of particle volume fraction and system size and microstructure of colloidal particles.⁵⁰ In this study, we investigated the effects of moderate/non-zero Re on the microstructure, average sedimentation velocity, velocity fluctuations and self-diffusion in a steady state homogeneous suspension of non-Brownian sedimenting particles over a wide range of volume fraction from 0.01 to 0.4. The present study not only quantifies the transport properties of the suspension at non-zero Re, it also gives a better understanding of the phenomenological changes based on a fluid and particle microstructure analysis. We organized this work in such a way that we could test the effect of Re at the low, moderate and high volume fraction regimes.

In this work, we first explained the simulation method and parameters in Section 2 and 3, respectively, and then investigated the DKT mechanism in Section 4. We presented our results in Section 5, where we investigated the hindered settling, fluid behavior around a particle, evolution of particle microstructure, velocity fluctuations, their relaxation times and self-diffusion. We conclude our study in Section 6.

2. Simulation method

In the SP method, the boundary between the colloidal particle and the host fluid is replaced with a continuous interface by assuming a smoothed profile,^51,52 having thickness equal or larger than the grid spacing Δ. This simple modification enabled us to accurately characterize the hydrodynamic interactions, without neglecting many-body interactions. A similar smoothed profile⁵³ was adopted in the previously proposed fluid-particle dynamics method in which particles are modelled by a highly viscous fluid. This large viscosity in turn requires a small time increment. Instead, we treated particles as non-deformable solids such that no additional constraints arose for the time increment. In this section, we briefly explain the salient features of our method.

The colloid's density profile is defined such that it changes smoothly within the interface. Quantities such as the velocity and pressure are defined over the entire computational domain, which includes the colloid as well as the solvent. The motion of the ith colloidal particle is obtained by solving Newton and Euler's equations of motion for a given particle position R_i, translational velocity V_i and rotational velocity Ω_i:

M_iV̇_i = F^H_i + F^c_i + F^ext_i, Ṙ_i = V_i, (1)

I_i·[capital Omega, Greek, dot above]_i = N^H_i + N^ext_i, (2)

where mass and moment of inertia are denoted as M_i and I_i, respectively. The hydrodynamic torque and force exerted by the solvent on the particle are represented as N^H_i and F^H_i, respectively. F^ext_i and N^ext_i are the external force and torque, respectively. Direct inter-particle interactions are denoted by F^c_i and we represent these interactions using the following Weeks–Chandler–Andersen (WCA) type potential,⁵⁴where r_ij = |R_i − R_j|. The parameters σ = 2a and ε = 2.5 denote the length and energy units, respectively.

The fluid motion is obtained by solving the Navier–Stokes equation for a given value of fluid viscosity η and density ρ_f within an Eulerian description as:

(∂_t + v_f·∇)v_f = ρ⁻¹_f∇·(−pI + σ_f) (4)

∇·ν_f = 0 (5)

where σ_f is the stress tensor, ν_f is the fluid velocity and p is the pressure field, under the incompressibility condition (∇·ν_f = 0). The center-of-mass of the suspension is kept constant to avoid an indefinite acceleration of the system. For given values of η and ρ_f, the remaining units of mass, time, pressure and energy are respectively defined as ρ_fΔ³, ρ_fΔ²/η, η²/ρ_fΔ² and η²Δ/ρ_f. The basic idea of the SP method is to solve the modified NS equation over the entire domain, by treating the colloids as fluid particles, where particles are represented using a smooth phase field 0 ≤ φ (x, t) ≤ 1, which removes the troublesome boundary conditions at the particle surface. Here φ = 0 stands for the fluid, 0 < φ < 1 describes the interface and φ = 1 the particle domain. The rigidity of the particles can be maintained by introducing a body force φf_p (ref. 51,55) in the NS equation, which is computed assuming momentum conservation between the fluid and the particles.

(∂_t + v·∇)v = ρ⁻¹_f∇·(−pI + σ_f) + φf_p, (6)

v = (1 − φ)v_f + φv_p, (7)

The detailed formulas, algorithm and applicability of the SP method with numerous test cases can be found in previous publications.^55–57 Recently, SP method has also been extended to self-propelled swimmers,⁵⁸ non-spherical rigid bodies⁵⁹ and compressible fluids.⁶⁰

A rigorous study⁶¹ by Luo et al. revealed that compared to a high resolution DNS⁶² based on high-order spectral/hp element discretization on hybrid grids, the SP method is computationally less expensive and advantageous for simulating moving particles because it avoids the complex discretizations around the particles. This study simulated many test cases and found that the SP method is accurately resolving the lubrication force and far and near field flows. In addition, the computational cost for this method depends upon the number of grid points (system size), however it is insensitive to the number of particles. This enabled us to study dense colloidal dispersion, which is difficult to do with other methods.

3. Simulation parameters

In order to study the inertial effects, we varied Re from 0.05 to 10 by increasing the gravity. In principal, we can go beyond Re = 10, but we have to reduce the time step to a much smaller value and need to refine the mesh as well, which requires enormous computer resources and is beyond the scope of this study. A cubic periodic box of dimension L/a = 32 is used, where the particle radius a is kept constant at 4Δ and the interface thickness is set to unity in the unit of grid spacing Δ in all simulations. These particles are non-Brownian, having Pe ≈ ∞, in accordance with our previous study⁴⁷ which shows that particles having Pe ≥ 29 show complete non-Brownian characteristics. Gravity was introduced in the −z direction. We set the particle to fluid density ratio to 5 and fluid density equal to one. Direct interparticle interactions are introduced by a Weeks–Chandler–Andersen (WCA) type potential with powers of 36 : 18, as defined in eqn (3). Simulations are run for 1000–1500t_t; sufficient time to yield statistically meaningful data for this analysis, where t_t denotes the terminal time (t_t = a/V_t).

In order to characterize the accuracy of our method, we calculated the drag coefficient of an isolated sphere and compared it with an empirical relation,⁶³ as shown in Fig. 1. The drag coefficient C_D which expresses the resistance of fluid on a moving object is calculated from the empirical expression.

Fig. 1 (Color online) Drag coefficient of an isolated sphere as a function of particle Reynolds number.

In order to simulate this test case, we used the same scheme as adopted by Koch et al.⁴⁰ Since our input parameter in the simulations is gravity (g), we first fix the particle Re and calculated the respective g. We ran simulations of an isolated particle for a long time until the steady state terminal velocity is achieved. From this terminal velocity, we calculated the Re and then C_D from eqn (9). Thereafter, we compared the C_D from the simulations with the empirical relation in Fig. 1. It shows a good agreement with the empirical relation up to Re = 14, whereas a deviation is evident at Re > 14. In principal, we can enhance the accuracy at high Re as well, but this require enormous computer resources. Using the simulation parameters for a single particle, we simulated multiple spheres in a periodic cubic cell.

4. Drafting–kissing–tumbling mechanism

Before discussing our main results, we investigated the DKT mechanism for particle pair interactions by simulating two interacting particles. We place particles one above the other, making sure that the trailing particle is in the wake of the leading particle. We calculated the distance between the particles in all three directions Δr_β, where β ∈ x, y, z and plotted it as a function of the terminal time for three different Re, as shown in Fig. 2. The vertical distance between the particles decreases initially, showing the drafting of the trailing particle, as it experiences a low pressure due to the wake of the leading particle. After drafting, the two particles interact with each other, which is the kissing phase, which can be seen through the plateau in all three figures at a vertical distance equal to σ. After the kissing phase, the lift force or the source flow sweep the particle to a horizontal orientation, which is the tumbling phase. The horizontal distance between the particles increases in the tumbling phase, whereas the vertical distance falls to zero. We have divided these three figures in three sections to highlight the DKT mechanism, where “D” shows the drafting, “K” shows the kissing and “T” shows the tumbling phase. We can also see that the time required to undergo a DKT type interaction decreases tremendously as Re is increased, which shows the fast phenomenological changes at high Re. In order to capture such a fast change, one has to use very sophisticated instruments in experiments and a very small time step and a refined mesh in simulations. This is why most of the experimental and simulation studies on sedimentation are in the Stokes regime.

Fig. 2 (Color online) Distance between two particles in all three directions Δr_β, where β ∈ x, y, z, as a function of terminal time for three different Re. Each figure is divided into three sections, where “D” shows the drafting, “K” shows the kissing and “T” shows the tumbling phase of the DKT mechanism.

5. Results and discussion

5.1 Hindered settling

A particle settling in the vicinity of other particles experiences a hindered settling due to the drag force induced by the fluid back flow and the particle-particle interactions. Consequently, it reduces the average sedimentation velocity of the suspension with respect to the terminal velocity of an isolated particle. The effect of fluid back flow and particle particle interaction is very small at low volume fractions, whereas both of these effects increase with increasing volume fraction. We computed the average sedimentation velocity V_sed of the suspension and compared it with the theoretical,^{33,36,64–67} experimental¹ and simulation results,^21,45 as shown in Fig. 3 and 4, for different Re, where V_sed = 〈V_iz〉. Fig. 3(a) shows that the average sedimentation velocity normalized by the terminal velocity, decreases with the increase of volume fraction due to this hindered settling. The decay of V_sed as a function of volume fraction at Re = 0.05 and 0.1 is same both qualitatively and quantitatively. It also shows a good quantitative agreement with other studies in the Stokes regime, depicting weak inertial forces effects, as shown in Fig. 3(b). In contrast, at Re ≥ 0.5, a decrease in V_sed is evident with the increase of Re due to enhanced inertial effects, which are more pronounced at low volume fractions. This decrease in V_sed is attributed to two factors: (1) fast DKT mechanism at high Re, where a particle spends less time in the wake of another particle as compared to low Re, which decreases the average sedimentation velocity of the particles; (2) particles are cross stream oriented for most of the time and feel the full effects of fluid back flow at high Re, which also decreases the average sedimentation velocity. At high volume fractions, the inertial effects are suppressed by the increasing number of particles and the average sedimentation velocity at high volume fraction is independent of Re, as shown in Fig. 3(a). We compared our results for Re ≤ 0.1 with other studies in Fig. 3(b). Batchelor³⁶ was the first who found the effect of the particle volume fraction on the average sedimentation velocity by assuming a uniform distribution in the separation of pairs of spheres and ignoring the fluid back flow effects, which is valid only for low volume fraction, as evident in Fig. 3(b). This hindered settling in the Stokes regime is well described by the Richardson and Zaki³³(RZ) power law, V_sed/V_t = (1 − ϕ)ⁿ, where the power law exponent n varies from 4.7 to 6.55 depending on studies.^{1,21,22,45,65,66} In our study the best suited exponent is 5.3, well within these two extremes for Re = 0.05 and 0.1, as shown in Fig. 3(b). A small deviation from the theoretical prediction for the Stokes regime is evident at Re = 0.5 for low volume fractions, showing the onset of inertial effects.

Fig. 3 (Color online) Average sedimentation velocity V_sed of particles normalized by the terminal velocity V_t of an isolated sphere, as a function of volume fraction for different Re. Theoretical predictions^{33,36,40,64,65} for Stokes regime are represented with lines, whereas points show the experimental¹ and simulation data.^21,45 (a) shows the comparison of V_sed of the present study for all Re, whereas (b) compares the results for Re ≤ 0.1 (Stokes regime) with other studies to highlight the Stokes regime behavior. Legends of our data in (a) are also valid for (b).

Fig. 4 (Color online) Average sedimentation velocity V_sed of particles normalized by the terminal velocity V_t of an isolated sphere, as a function of volume fraction for Re = 0.5, 1 and 10. Simulation results are fitted with the suggested hindered function of eqn (13) and are also compared with the other theoretical predictions^40,65 and simulation results.⁴⁵ Legends are valid for all three figures.

Fig. 4(a)–(c) show that the theoretical predictions for the Stokes regime are no longer applicable when inertial effects are relevant, rather a modified expression⁶⁵ of the form.

V_sed/V_t = k(1 − ϕ)ⁿ, (10)

has been suggested, where the terminal velocity at a given Re can be calculated as, V_t = Reη/ρ_fσ. The value of the prefactor k and the exponent n varies with the particle Re and Koch et al.,⁴⁰ found that the k for Re = 1 and 10 is 0.92 ± 0.03 and 0.86 ± 0.04, respectively, and suggested a quadratic polynomial expression for the value of n of the form,

n = 4.23 − 0.0526Re + 0.00111Re². (11)

Similary, Di Felice⁶⁵ found a value of k ≈ 0.81 at intermediate Re in their experimental study and suggested an expression for n of the form.

The value of k is found by fitting the data and it ranges from 0.8 to 0.95 in the literature.^{40,65,68–70} Di Felice also suggested that the value of the exponent n at low ϕ and Re > 1 is 1.5 times the corresponding value from eqn (12), caused by the DKT type interactions.

Both Koch and Di Felice provided ample information about the mechanism which caused the initial fast decay of the hindered settling function at low volume fraction but the mathematical expression they suggested fails to provide the necessary information at ϕ ≈ 0. In order to provide this missing information, we suggested a modified hindered settling expression of the form.

V_sed/V_t = k(1 − ϕ)ⁿ + (1 − k)exp(−ϕ/0.008), (13)

where values of k and n can be found as,and

This simple fitting function gives the complete information for a given value of Re and volume fraction. It reduces to RZ law at Re ≤ 0.1, whereas at 0.1 < Re ≤ 10 second term gives the information for the fast decay of the expression, as shown in Fig. 4.

Moreover, we also found that the rapid decay of average sedimentation velocity at low ϕ and Re = 10 is attributed to the strong influence of the anisotropic microstructure⁴⁰ which will be explained later, having pronounced effects at high Re. At low ϕ and high Re, a particle is deficit of neighboring particles, due to the DKT mechanism.⁴⁰ This mechanism makes the spheres well separated, having a cross-stream orientation. These well separated particles feel the full effect of the fluid back flow, which in turn causes a rapid decay in the sedimentation velocity. At low Re, a particle spends a significant time in the wake of the leading particle, which causes the increase in the average sedimentation velocity, whereas at Re = 10, Fig. 2 suggests that a particle spends very little time in the wake and most of its time being spent in the cross stream orientation, feeling the full back flow effect, which causes a rapid decrease in the average sedimentation velocity at low volume fraction. As the concentration of particles increases, this mechanism competes with the many body interactions and as the latter overpower the phenomena, the rate of decrease slows down.

5.2 Fluid behavior

Before discussing the particle properties, we would like to look at the fluid behavior around a particle at different Re. At low Re in the Stokes regime, it is well established that the fluid has a fore-aft symmetry around a particle, whereas at high Re, fluid inertia breaks this symmetric behavior.

We calculate the fluid velocity as a function of position and time denoted as ν_f (x, t). In order to look at the fluid around the particle in cylindrical coordinates, we shifted the origin at the particle position as, u_f (r_cyl, z) = 〈ν_f (r)〉, and averaged over time, number of particles and angular coordinate θ. Here, r = x − R_i, is the position vector of the fluid from the particle in which x and R_i show the positions of fluid and particle, respectively. If r^x, r^y and r^z are the components of r in the Cartesian coordinates, then , θ = arc tan(r^y/r^x) and z = r^z. In order to characterize the asymmetric behavior of fluid before (B) and after (A) the particle, we plotted |Δu^AB_f(r_cyl, z)| = |u^A_f(r_cyl, z)−u^B_f(r_cyl, z)|/V_sed, in Fig. 5. The fluid velocities before and after the particle are defined in term of their components as, and , respectively, for z ≥ 0. Fig. 5(a)–(d) show this difference at ϕ = 0.025, 0.04, 0.07 and 0.1, respectively, for Re = 0.1, 1, 10. By comparing these plots from left to right for different Re at a particular ϕ, we found that Re = 0.1 shows a symmetric behavior, whereas with increasing Re the difference between the velocities after and before the particle increases, showing the effects of the fluid inertia. If we compare these figures from top to bottom for Re = 1 and 10 for different ϕ, it is evident that the inertial effects are suppressed with the increase of particle volume fraction, caused by the presence of more number of particles. Fig. 5 clarifies the asymmetric behavior of the fluid before and after the particle and to get a more insightful picture of the fluid behavior, we also plotted the difference of velocities at Re = 1 and 10 with respect to the corresponding velocity at Re = 0.1 as, , in Fig. 6 for different ϕ, where γ ∈ 1, 10. Comparing Fig. 6(a) and (b) from left to right reveals that this difference decreases with the increase of particle volume fraction, as the increased packing of particles suppresses the effects of fluid inertia. Similarly, comparing figures from top to bottom reveal that this difference increases with the increase of Re. Furthermore, this difference is stronger before than after the particle at Re = 10. If we look at the rear end of the particle at Re = 10 in Fig. 6(b) this difference is almost negligible, whereas at the radial distance σ to 2σ the difference in the velocities is stronger because of the strong inertia of the fluid.

Fig. 5 (Color online) The absolute value of the difference in the fluid velocity before (B) and after (A) the particle |Δu^AB_f(r_cyl,z)| = |u^A_f(r_cyl,z) − u^B_f(r_cyl,z)|/V_sed, in cylindrical coordinates. Figures “(a)”, “(b)”, “(c)” and “(d)” show this difference at ϕ = 0.025, 0.04, 0.07 and 0.1, respectively, for different Re.

Fig. 6 (Color online) The absolute value of the difference in the fluid velocity at Re = 1 and 10 with the Re = 0.1 as, , for different ϕ, where γ ∈ 1, 10.

The fluid behavior around a particle has a great importance in the determination of the particle configuration and their transport properties. As it is shown that at high Re, normalized velocities are stronger, hence this strong behavior will push away the nearby particles. The effects of this asymmetric behavior of the fluid before and after the particle on the particle configuration is discussed in the next section.

5.3 Microstructure

The microstructure of the particles strongly affects the transport properties of the suspension. We explored the evolution of the microstructure with Re and volume fraction and its effects on transport properties. In order to have an idea of the average accumulation of the particle pairs, we calculated the pair distribution function (PDF) as.in which N is the total number of particles, r_ij = R_i − R_j, 〈…〉 denotes an ensemble average, and the summation is taken over all particle pairs. We plotted the spherically averaged PDF, known as the radial distribution function (RDF), for different volume fractions at different Re, as shown in Fig. 7. For ϕ = 0.01 and Re ≤ 0.1, the RDF shows a distribution typical of hard spheres, whereas at Re ≥ 0.5, the RDF shows the deficiency of neighbouring particles and this deficiency increases with increasing Re, in good agreement with Koch.⁴⁰ As ϕ is increased, the deficit is decreased and at ϕ = 0.03 and Re ≤ 1, particles already show the typical hard sphere distribution, whereas at Re = 10, the deficit of particle pairs is still apparent, as shown in Fig. 7(b). With a further increase in volume fraction, particles start showing the RDF of a hard sphere distribution, as shown in Fig. 7(c) and (d). The increase in particle inertia at Re ≥ 1 not only induced the significant deficit of the particle pairs but it also shifted the peak and this shift decreases with the increase of volume fraction. The RDF provides useful information of the particle deficiency and the shift in the peak with increasing Re, but it does not give any information about the particular direction where this deficiency exists. In order to investigate this preference, we compute the PDF in cylindrical coordinates g(r_cyl,z) averaged over the angular coordinate θ, for three different Re, as shown in Fig. 8, where lines are drawn at σ, , and 2σ to visualize the long range ordering. Fig. 8 shows the non-Brownian distribution, where particles show the preference to orient themselves in a horizontal position with respect to a test particle and this preference decreases with increasing volume fraction. At high volume fraction, the microstructure is isotropic. At Re = 1, PDF in cylindrical coordinates shows the deficit of particle pairs and this deficit is more pronounced in the vertical direction, as shown in Fig. 8(b).

Fig. 7 (Color online) Radial distribution function as a function of radial distance for different Re at four different volume fractions.

Fig. 8 (Color online) Pair distribution function of particles in cylindrical coordinates (g(r_cyl, z)) for different ϕ at three different Reynolds number.

At Re = 10, PDF in cylindrical coordinates shows the deficit of particle pairs and this deficiency is stronger than the one at Re = 1, in both directions. Compared to Re = 1, PDF in cylindrical coordinates at Re = 10 shows that the preference of particles in the horizontal direction is stronger and a shift in the peak is also evident. At high volume fraction ϕ ≳ 0.15, an isotropic microstructure is visible. The comparative analysis of these maps revealed that with the increase of Re, the DKT mechanism grows stronger and its effects is felt over a larger volume fraction range. It is also evident from this microstructure analysis that the deficit of particles is stronger in the vertical direction than the horizontal direction at all Re. This is because, a particle in the wake of a leading particle experiences a DKT type interaction, which leads the particles to a horizontal orientation. Moreover, a horizontally orientated particle faces the source flow, which keeps the particle in the horizontal orientation as well. Hence, the preferred position of particle pairs is horizontal, which make the deficit stronger in the vertical direction, until many-body interactions suppress this DKT type interaction and make the microstructure isotropic. The effects of this microstructure on the dynamic properties will be discussed in the next section.

The change in the microstructure of the particles at high Re is in accordance with the fluid behavior at high Re. The large difference of the fluid velocity at high Re with respect to Re = 0.1, pushed the particle away from a given particle and created a deficiency near the given particle and DKT type interactions orient the particle pairs horizontally. Since this difference in fluid velocity is increased with the increase of Re, which pushed the particle further away from a given particle and created the deficiency of the particle pairs in the near vicinity. The range of radial distance over which this deficiency is evident increased with the increasing Re due to the strong fluid velocity. Moreover, the range of volume fraction which showed this difference in fluid velocity is also increased with the increase of Re, which in turn increase the range of volume fraction which shows the deficiency of nearby particles. If we look closely at ϕ = 0.02 in Fig. 8, small high density regions are evident above and below the particle for all Re, whereas these high density regions are not visible at ϕ = 0.05. These regions are caused by the finite size effects. If we compare the microstructure of this study for the Stokes regime of system size L/a = 32 with our previous study⁴¹ of system size L/a = 64, we can not spot these high density regions at ϕ = 0.02 for L/a = 64. This shows that our simulations are affected by the finite size effects at low volume fractions. These finite size effects are suppressed with increasing volume fractions and Re due to presence of large number of particles and the formation of wakes and eddies, respectively. Similarly, other properties like velocity fluctuations, their relaxation times and self diffusion coefficients are also affected by the finite system size at low volume fractions. More details on these finite size effects can be seen in our previous publications.^41,47

5.4 Velocity fluctuations

The particle velocity fluctuations induced by HIs are quantified both in the direction parallel and perpendicular to gravity as, and , respectively. These velocity fluctuations are plotted in Fig. 9, as a function of particle volume fraction for different Re. For Re ≤ 0.5, velocity fluctuations in both directions show the Stokes regime characteristics: they scale as ϕ^1/2 at low volume fraction regime (ϕ ≲ 0.04), remain unchanged at moderate volume fractions (0.04 < ϕ ≲ 0.12) and decrease sharply at high volume fraction regimes, in good agreement with other studies.^1,5,21,45 At Re = 1, vertical velocity fluctuations are decreased, whereas horizontal velocity fluctuations still show the Stokes regime behavior. The behavior of velocity fluctuations in both directions is in accordance with the microstructure changes, where a deficit of particles is observed in the vertical direction, whereas there was no significant deficiency of particles in the horizontal direction.

Fig. 9 (Color online) Scaling of the hydrodynamic velocity fluctuations normalized by the terminal velocity with respect to volume fraction for different Re. Figures “(a)” and “(b)” show the velocity fluctuations in z and x directions, respectively, for Re = 0.1, which shows the typical Stokes regime characteristics. These results are also compared with the previous studies. Figures “(c)” and “(d)” show the velocity fluctuations in z and x directions, respectively, for different Re. Lines of slope 1/2 are drawn as guide to eye.

At Re = 10, with the enhanced inertial force effects, both vertical and horizontal velocity fluctuations are decreased, due to a pronounced DKT mechanism, which caused a significant deficit of the particles in both directions. Apart from the decrease in the velocity fluctuations, the range of volume fraction which show the decrease in velocity fluctuations also increased with increasing Re. This increase in the range is caused by the deficit of the neighbouring particles in the vicinity of a test particle, which increases with increasing Re, before many body interactions suppress the effects of inertial forces and the microstructure is similar to the Stokes regime microstructure.

The vertical and horizontal relaxation times decrease with the increase of volume fraction for all Re. These relaxation times are calculated as, , where β ∈ x, z and C^β(t) is the velocity fluctuations autocorrelation function. The integration is performed from zero to the time at which these correlations start to oscillate around zero. The horizontal relaxation times at low volume fraction increase with the increase of Re. At high volume fractions and Re ≤ 1 horizontal relaxation times show the same quantitative behavior, whereas for Re = 10, these are increased. A jump at ϕ ≈ 0.3 is attributed to the ordering of the particles, as shown in Fig. 10(a). In contrast, the vertical relaxation time decreases with the increase of Re, at low volume fraction. At high volume fractions and Re ≤ 1 vertical relaxation times also show the same quantitative behavior, as shown in Fig. 10(b). The increase in the relaxation times at high Re and volume fraction remains an open question.

Fig. 10 (Color online) Decay of relaxation times normalized by the terminal time, as a function of volume fraction for different Re. Figures “(a)” and “(b)” show the x and z directions relaxation times, respectively.

5.5 Self-diffusion

In the absence of an external force, particles exhibit an isotropic self-diffusion, whereas an external force such as gravity breaks this isotropic nature. We have seen the effects of the increase of inertial forces on the velocity fluctuations and their relaxation times. Both velocity fluctuations and their relaxation times dictate the diffusion of the particle (D^β ≈ (ΔV^β)²τ^β). In order to quantify this dispersive motion, we calculated the self-diffusion from the slope of mean square displacement in the x and z directions, as defined in eqn (17) and (18), respectively.

The saturated value of these diffusion coefficients normalized by the product of the terminal velocity and particle radius are plotted in Fig. 11, as a function of volume fraction, for different Re. It shows that at Re ≤ 0.5 the vertical diffusion coefficient shows the typical Stokes regime ϕ^1/2 scaling at low volume fraction regime, decreasing at moderate volume fraction regime and then falls sharply at high volume fraction. The horizontal diffusion coefficient also increases as ϕ^1/2 at low volume fractions, remains unchanged at moderate volume fraction and decreases sharply at high volume fractions, in good agreement with other studies.^1,21,22

Fig. 11 (Color online) A comparison of the long time steady-state self-diffusion coefficients D^β(β ∈ x, z), normalized by the product of particle radius a and terminal velocity V_t for different Re. Figures “(a)” and “(b)” show the horizontal and vertical directions, respectively.

The vertical diffusion coefficient decreases with the increase of Re at Re ≥ 1, whereas at high volume fraction this behavior becomes independent of Re, due to the dominance of many body interactions. Similar to the vertical velocity fluctuations, the range of volume fraction which show the decrease in the diffusion coefficients is increased with the increase of Re, which is attributed to change in microstructure and fluid behavior. Similar to horizontal velocity fluctuations, the diffusion coefficient in the horizontal direction remains unchanged at Re = 1, due to small change in microstructure in horizontal direction, whereas at Re = 10 a quantitatively large decrease in diffusion coefficient is evident up to ϕ ≲ 0.12, caused by the pronounced changes in the microstructure and fluid behavior.

6. Concluding remarks

We performed simulations of sedimentation of spherical particles at non-zero Reynolds number for a volume fraction ranging from 0.01 to 0.4. Inertial forces are increased with increasing Re, giving rise to a drafting–kissing–tumbling (DKT) mechanism.⁴⁰ With increasing Re, this mechanism creates a deficit of neighbouring particles which brings significant phenomenological changes and affects the transport properties of the suspension.

We observed that at Re = 0.1 the fluid shows a symmetric behavior before and after the particle, whereas at Re = 1 and 10 an asymmetric behavior is observed. We also found that the difference in the fluid velocities normalized by their respective average sedimentation velocities at Re = 1 and 10 with respect to Re = 0.1 increased with the increase of Re. Moreover the range of volume fraction which show this difference also increased with the increase of Re, before the large number of particles suppress the inertial effects. The microstructure analysis of the particles revealed that at Re ≤ 0.5, the microstructure shows the typical Stokes regime, non-Brownian particles characteristics, where at low volume fraction particles prefer to orient themselves in the horizontal direction and this preference decreases with increasing volume fraction and at high volume fraction ϕ ≳ 0.15, the microstructure is isotropic due to dominance of many-body interactions. With the increase of inertial forces and the enhanced effects of the DKT mechanism, the microstructure shows the deficit of neighbouring particles around a test particle and this deficiency is stronger in the vertical direction than in the horizontal direction at Re = 1. At Re = 10, with the pronounced effects of inertial forces, both directions are affected by this deficiency and we also observed a shift in the radial distribution function's first peak. At high volume fractions, where many-body interactions dominate, microstructure resembles the typical Stokes regime's microstructure. Apart from the deficiency of the particles, the range of volume fraction affected by this deficiency is also increased with increasing Re. These changes in the microstructure of the particles are in accordance with the change in the fluid velocities observed at high Re.

The fluid and microstructural changes at high Re, strongly affected the transport properties of the suspension and at Re = 1, vertical velocity fluctuations and diffusion are decreased. In contrast, horizontal velocity fluctuations and diffusion remain unchanged at Re = 1, due to less pronounced effects of inertial forces in the horizontal direction. Furthermore, at Re = 10, due to strong effects of inertial forces, which caused the significant deficit of neighbouring particles in both directions, which in turn decreased both the vertical and horizontal velocity fluctuations and diffusion. Moreover, the vertical relaxation times decreased with the increase of Re at low volume fraction, whereas horizontal relaxation time increased with the increase of Re. Our analysis also revealed that at low volume fraction, where particles can interact with fewer neighbouring particles, microstructure governs the transport properties of the particles. With increasing Re, the rise of the DKT mechanism creates the deficiency of the neighbouring particles, which means particles still interact with the fewer other particles even at higher volume fraction range, causing the microstructure to control the transport properties of the suspension to the higher volume fraction range. At high volume fractions ϕ ≳ 0.15, intrinsic many-body interactions suppress the DKT mechanism and take over control of the transport properties. The present study has been performed using spherical particles, research is under way using non-spherical rigid bodies.

Acknowledgements

This work was supported by KAKENHI 23244087 and the Global COE Program, “International Center for Integrated Research and Advanced Education in Materials Science” from the Ministry of Education, Culture, Sports, Science and Technology of Japan. The author would like to extend his gratitude to Dr Takahiro Murashima and Dr Hideyuki Mizuno for useful discussions and suggestions.

References

1. H. Nicolai, B. Herzhaft, E. J. Hinch, L. Oger and E. Guazzelli, Phys. Fluids, 1995, 7, 12–23.
2. H. Nicolai and E. Guazzelli, Phys. Fluids, 1995, 7, 3–5.
3. H. Nicolai, Y. Peysson and E. Guazzelli, Phys. Fluids, 1996, 8, 855–862.
4. P. N. Segrè, O. P. Behrend and P. N. Pusey, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 5070–5083.
5. P. N. Segrè, E. Herbolzheimer and P. M. Chaikin, Phys. Rev. Lett., 1997, 79, 2574–2577.
6. P. Segrè, F. Liu, P. Umbanhowar and D. A. Weitz, Nature, 2001, 409, 594–597.
7. P. N. Segrè and J. McClymer, J. Phys.: Condens. Matter, 2004, 16, S4219–S4230.
8. S.-Y. Tee, P. J. Mucha, L. Cipelletti, S. Manley, M. P. Brenner, P. N. Segrè and D. A. Weitz, Phys. Rev. Lett., 2002, 89, 3–6.
9. D. L. Koch and E. Shaqfeh, J. Fluid Mech., 1991, 224, 275–303.
10. A. Levine, S. Ramaswamy, E. Frey and R. Bruinsma, Phys. Rev. Lett., 1998, 81, 5944–5947.
11. D. L. Koch, Phys. Fluids, 1994, 6, 2894–2900.
12. A. J. C. Ladd, Phys. Fluids A, 1993, 5, 299–310.
13. A. J. C. Ladd, Phys. Rev. Lett., 1996, 76, 1392–1395.
14. A. J. C. Ladd and R. Verberg, J. Stat. Phys., 2001, 104, 1191–1251.
15. A. J. C. Ladd, Phys. Rev. Lett., 2002, 88, 048301.
16. A. J. C. Ladd, Phys. Fluids, 1997, 9, 491–499.
17. R. Verberg and A. J. C. Ladd, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 65, 016701.
18. N.-Q. Nguyen and A. J. C. Ladd, J. Fluid Mech., 2005, 525, 73–104.
19. J. T. Padding and A. Louis, Phys. Rev. Lett., 2004, 93, 220601.
20. J. T. Padding and A. Louis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 031402.
21. J. T. Padding and A. A. Louis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 011402.
22. F. R. Cunha, G. C. Abade, A. J. Sousa and E. J. Hinch, J. Fluids Eng., 2002, 124, 957.
23. S.-Y. Tee, P. J. Mucha, M. P. Brenner and D. A. Weitz, Phys. Fluids, 2012, 113304, 113304.
24. M. Miguel and R. Pastor-Satorras, Europhys. Lett., 2001, 54, 45–50.
25. E. Kuusela, K. A. I. Höfler and S. Schwarzer, J. Eng. Math., 2001, 41, 221–235.
26. E. Kuusela and T. Ala-Nissila, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 061505.
27. E. Kuusela, J. M. Lahtinen and T. Ala-Nissila, Phys. Rev. Lett., 2003, 90, 094502.
28. E. Kuusela, J. M. Lahtinen and T. Ala-Nissila, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 066310.
29. E. Kuusela, J. M. Lahtinen and T. Ala-Nissila, Europhys. Lett., 2004, 65, 13–19.
30. P. Mucha and M. Brenner, Phys. Fluids, 2003, 15, 1305–1313.
31. P. J. Mucha, S.-Y. Tee, D. A. Weitz, B. I. Shraiman and M. P. Brenner, J. Fluid Mech., 2004, 501, 71–104.
32. G. G. Stokes, Mathematical and physical papers, Johnson Reprint Corporation, New York, 1996.
33. R. F. Richardson and W. N. Zaki, Trans. Inst. Chem. Eng., 1954, 32, 35–53.
34. R. E. Caflisch and J. H. C. Luke, Phys. Fluids, 1985, 28, 759.
35. E. J. Hinch, Disorder and Mixing, Kluwer Academic publisher, Netherland, 1988.
36. G. K. Batchelor, J. Fluid Mech., 1972, 52, 245.
37. Y. Peysson and E. Guazzelli, Phys. Fluids, 1998, 10, 44.
38. G. Bernard-Michel, A. Monavon, D. Lhuillier, D. Abdo and H. Simon, Phys. Fluids, 2002, 14, 2339.
39. M. P. Brenner, Phys. Fluids, 1999, 11, 754.
40. X. Yin and D. L. Koch, Phys. Fluids, 2007, 19, 093302.
41. A. Hamid, J. J. Molina and R. Yamamoto, Soft Matter, 2013, 9, 10056–10068.
42. D. L. Koch, Phys. Fluids A, 1993, 5, 1141.
43. X. Yin and D. L. Koch, Phys. Fluids, 2008, 20, 043305.
44. M. Shin, D. L. Koch and G. Subramanian, Phys. Fluids, 2009, 21, 123304.
45. E. Climent and M. Maxey, Int. J. Multiphase Flow, 2003, 29, 579–601.
46. A. A. Johnson and T. E. Tezduyar, Comput. Meth. Appl. Mech. Eng., 1995, 134, 351–373.
47. A. Hamid and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 022310.
48. A. Hamid and R. Yamamoto, J. Phys. Soc. Jpn., 2013, 82, 024004.
49. A. Hamid and R. Yamamoto, AIP Conf. Proc., 2013, 1518, 444–447.
50. A. Hamid, J. J. Molina and R. Yamamoto, Mol. Simul., 2014 DOI:10.1080/08927022.2014.929124.
51. Y. Nakayama and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 036707.
52. Y. Nakayama, K. Kim and R. Yamamoto, Adv. Powder Technol., 2010, 21, 206–211.
53. H. Tanaka and T. Araki, Phys. Rev. Lett., 2000, 85, 1338–1341.
54. A. Ahmed and R. J. Sadus, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 061101.
55. Y. Nakayama, K. Kim and R. Yamamoto, Eur. Phys. J. E: Soft Matter Biol. Phys., 2008, 368, 361–368.
56. R. Yamamoto, Y. Nakayama and K. Kim, Comput. Phys. Commun., 2005, 169, 301–304.
57. R. Yamamoto, Y. Nakayama and K. Kim, Int. J. Mod. Phys. C, 2009, 20, 1457–1465.
58. J. J. Molina, Y. Nakayama and R. Yamamoto, Soft Matter, 2013, 9, 4923–4936.
59. J. J. Molina and R. Yamamoto, J. Chem. Phys., 2013, 139, 234105.
60. R. Tatsumi and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 066704.
61. X. Luo, M. R. Maxey and G. E. Karniadakis, J. Comput. Phys., 2009, 228, 1750–1769.
62. R. M. Kirby, T. C. Warburton, I. Lomtev and G. E. Karniadakis, Appl. Numer. Math., 2000, 33, 393–405.
63. R. Clift, J. R. Grace and M. E. Weber, Bubbles, Drops and Particles, Academic Press, London, 1978.
64. H. Hayakawa and K. Ichiki, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, R3815–R3818.
65. R. Di Felice, Int. J. Multiphase Flow, 1999, 25, 559–574.
66. J. Garside and M. Al-Dibouni, Ind. Eng. Chem. Process Des. Dev., 1977, 16, 206–2014.
67. J. F. Brady and L. J. Durlofsky, Phys. Fluids, 1988, 31, 717.
68. R. Di Felice and E. Parodi, AIChE J., 1996, 42, 927–931.
69. R. Di Felice and R. Kehlenbeck, Chem. Eng. Technol., 2000, 23, 1123–1126.
70. Y. Chong, D. Ratkowsky and N. Epstein, Powder Technol., 1979, 23, 55–66.

---