Study of Polydisperse Particulate Systems with a ‘Direct-Forcing/Fictitious Domain’ Method

Abstract

Natural sediments responsible for the morphodynamic of the estuaries and coast are of different sizes and densities. Some are cohesive and some are non-cohesive. The transport in suspension and their sedimentation of such a polydisperse suspension are different than the ones for a monodisperse suspension. A fully resolved model based on the Direct-Forcing/Fictitious Domain method (DF/FD) was developed and applied to simulate settling of monodisperse particles in a water column. The behaviour of the suspension corresponds qualitatively to experimental results and average settling velocities follow a Richardson-Zaki type law. Then the model is applied to the sedimentation of suspension composed of particles of three diameters. The segregation of the bed is obtained naturally. The excess pore pressure is drawn and compared with the theory.

1. Introduction

Estuaries are places where silting-up events usually occur. Different processes are the root of particle sedimentation on the bed, which contributes to the filling process of estuaries. Flocculation processes (aggregation and break up), erosion, deposition and compaction are major phenomena that must be considered in order to predict the long term behaviour of estuarine sedimentary dynamic. So, in sedimentary transport, as well as other problems with particulate flows, the complex dynamics is not well described [1].

Most of the study had been realised at macro-scales, due to the huge number of particles that constitute these particulate systems, where the fluid and the particles can be considered as a mixture. So it is of single-phase models at large scale without [2] or with sedimentation consolidation model [3]. They can simulate the behaviour of large-scale systems such as estuaries or ports. At smaller scale, Euler–Euler two-phase flow models have been used with success [4, 5, 6, 7, 8, 9, 10]. The predictive ability and accuracy of theses codes depend essentially on the quality of the closure equations that model phenomena at smallest scales. Processes such as drag or lift forces acting on a suspension are not well established, and there is clearly a need to better understand smallest-scale phenomena in order to better formulate these constitutive relations.

Direct numerical simulation models for particulate flows have received a great attention for 20 years. In these methods, the fluid phase is governed by the Navier–Stokes equations, and the rigid inclusions are governed by Newton’s laws. The flow field around each particle is resolved, and hydrodynamic interactions are results of simulations. Such a model was proposed by Glowinski et al. [11], Peskin [12], Yu and Shao [13], Wachs [14].

Natural sediments responsible for the morphodynamic of the estuaries and coast are of different sizes and densities. Some are cohesive and some are non-cohesive. The transport in suspension and their sedimentation of such a polydisperse suspension are different than the ones for a monodisperse suspension [15]. Thus, the complex behaviour of such particulate systems is not well understood and not really explored by numerical investigations.

We propose here to examine the sedimentation of a polydisperse suspension by the use of a two-dimensional fully resolved model that we developed and validated with monodisperse suspensions in water column. Section 2 briefly describes the mathematical and numerical background of the model. Section 3 presents a validation of the model. Simulation results on the sedimentation of monodisperse and polydisperse particle systems are discussed in Section 4. Finally, Section 5 presents the conclusions and perspectives of the work.

2. Mathematical and numerical background

In the direct numerical simulation models for particulate flows, the fluid phase is governed by the Navier–Stokes equations, and the rigid inclusions are governed by Newton’s laws. The flow field around each particle is resolved, and hydrodynamics interactions are results of simulation. The model presented here is based on the ‘Direct-Forcing/Fictitious Domain’ method of Yu and Shao [13] and the Direct Numerical code of Guillou and Makhloufi [16]. Herein we briefly recall the basic equations with special closure laws. More details can be found in Verjus [17] and Verjus et al. [18].

2.1 Fictitious domain formulation

Let consider a Newtonian fluid and one rigid cylindrical particle. Let P(t) represent the domain of the Particule and Ω the whole domain including the fluid and the particle. Suppose that the particle density, volume and moment of inertia, translational velocity and angular velocity are ρ_s, V_s, J_s, U and ω, respectively. The fluid density and viscosity are ρ_f and μ. Let us introduce scale parameters: L_c for length, U_c for velocity, L_c/U_c for time, ρ_fU_c² for the pressure and ρ_fU_c²/L_c for the pseudo body force. Then, the dimensionless fictitious domain equations can be written as Yu and Shao [13]:

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}.\nabla \mathbf{u}=\frac{{\nabla }^2\mathbf{u}}{Re}-\nabla p+\mathbf{\lambda }\text{in}\mathrm{\Omega },$$E1

$$\mathbf{u}=\mathbf{U}+\mathbf{\omega }\times \mathbf{r}\text{in}P\left(t\right),$$E2

$$\nabla \mathbf{.}\mathbf{u}=0\text{in}\Omega ,$$E3

$$\left({\rho }_r-1\right)V_d\left(\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}-\frac{1}{\mathit{Fr}}\frac{\mathbf{g}}{g}\right)=-\underset{P}{\int }\mathbf{\lambda }\mathrm{d}x,$$E4

$$\left({\rho }_r-1\right)\frac{\mathrm{d}\left({\mathbf{J}}_{\mathbf{d}}\mathbf{.}\mathbf{\omega }\right)}{\mathrm{d}t}=-\underset{P}{\int }\mathbf{r}\times \mathbf{\lambda }\mathrm{d}x,$$E5

where u represents the fluid velocity, p the fluid pressure, λ the pseudo body force, U and ω are respectively the velocity and angular velocity of the particle, r is the position vector with respect to the particle mass centre, Re is the Reynolds number defined by ρ_fU_cL_c/μ, Fr is the Froude number defined by U_c²/gL_c, ρ_r is the particle-fluid density ratio defined by ρ_s/ρ_f, J_d is the dimensionless moment of inertia defined by J_s/ρ_sL_c⁴ and V_d the dimensionless volume defined by V_s/L_c².

2.2 Numerical method

A fractional-step method is used in order to decouple the problem described by the Eqs. (1)–(5). For a time instant, three major steps are written: one for the prediction of the fluid velocity (and pressure) (Eqs. (6) and (7)); one for the calculation of the particle’s motion and the interaction with the fluid (Eqs. (7)–(9)); and the last one to correct the fluid’s velocity (Eq. (10)). The present model is an extension of the SUDRES code [16] for the DNS of the fluid flow in two dimensions. A Finite Difference Method on a staggered grid and a projection technique are used to solve the fluid problem. Finally, the time marching-algorithm used is as follows:

• Step 1: Calculation of the fluid velocity and pressure with a projection method

$$\frac{{\mathbf{u}}^{\ast }-{\mathbf{u}}^n}{\mathrm{\Delta }t}=-{\nabla }^{\ast }p-\frac{1}{2}\left(3{\left(\mathbf{u}.\nabla \mathbf{u}\right)}^n-{\left(\mathbf{u}.\nabla \mathbf{u}\right)}^{\left(n-1\right)}\right)+\frac{1}{2}\left\{\frac{1}{Re}{\nabla }^2{\mathbf{u}}^n+\frac{1}{Re}{\nabla }^2{\mathbf{u}}^{\ast }\right\}+{\mathbf{\lambda }}^n$$E6

$$\nabla \cdot {\mathbf{u}}^{\ast }=0$$E7

• Step 2: Calculation of the particle subproblem for U^n + 1, ω^n + 1, u^n + 1 and the body force λ^n + 1, with:

◦ Step 2.1: Calculation of U^n + 1, ω^n + 1 depending on u^* and λⁿ

$${\rho }_r{V}_d\frac{{\mathbf{U}}^{n+1}}{\mathrm{\Delta }t}=\left({\rho }_r-1\right)V_d\left(\frac{{\mathbf{U}}^n}{\mathrm{\Delta }t}-\frac{1}{\mathit{Fr}}\frac{\mathbf{g}}{g}\right)+\underset{P}{\int }\left(\frac{{\mathbf{u}}^{\ast }}{\mathrm{\Delta }t}-{\mathbf{\lambda }}^n\right)\mathit{dx}$$E8

$${\rho }_r\frac{{\mathbf{J}}_{\mathbf{d}}\mathbf{.}{\mathbf{\omega }}^{n+1}}{\mathrm{\Delta }t}=\left({\rho }_r-1\right)\frac{{\mathbf{J}}_{\mathbf{d}}\mathbf{.}{\mathbf{\omega }}^n}{\mathrm{\Delta }t}+\underset{P}{\int }\mathbf{r}\times \left(\frac{{\mathbf{u}}^{\ast }}{\mathrm{\Delta }t}-{\mathbf{\lambda }}^n\right)\mathit{dx}$$E9

◦ Step 2.2: Calculation of λ^n + 1 depending on U^n + 1, ω^n + 1, u^* and λⁿ

$${\mathbf{\lambda }}^{n+1}=\frac{{\mathbf{U}}^{n+1}+{\mathbf{\omega }}^{n+1}\times \mathbf{r}-{\mathbf{u}}^{\ast }}{\mathrm{\Delta }t}+{\mathbf{\lambda }}^n$$E10

• Step 3: Correction of the fluid velocity u^n + 1 regarding λ^n + 1 and λⁿ

$${\mathbf{u}}^{n+1}={\mathbf{u}}^{\ast }+\mathrm{\Delta }t\left({\mathbf{\lambda }}^{n+1}-{\mathbf{\lambda }}^n\right)$$E11

The discretization of the particle (Lagrangian mesh associated to the particle) is still an open question. In the Distributed Lagrangian Multiplier Method, different techniques are used to discretize the particle [19], but the Collocation Element Method is often used. In the present study, a structured meshing strategy is used to discretize the particle (Figure 1). So for the solid problem, the particle is dicretized with the Collocation Point Method. The arrangement of the Lagrangian points is presented in Yu and Shao [13]. Special attention is dedicated to the space between particle nodes: it is greater than the fluid nodes space as suggested by Glowinski et al. [11]. Bilinear interpolations are used to interpolate quantities defined on the fluid’s mesh (Eulerian) on the particle’s mesh (Lagrangian) and vice versa. A trapezoidal rule is used to perform integrals. The collision strategy employed in this code is based on a normal repulsive force acting when two particles are too close each other [11].

Figure 1.

Meshes of the fluid and solid domains.

3. Validation

The numerical model is validated by three cases: the simulation of the flow around a fixed particle, the sedimentation of one particle and the sedimentation of two particles.

3.1 Flow around a fixed particle

The first test case is the fixed cylinder (circular particle) in Poiseuille flow. A cylinder is placed at the centre of a channel of width L = 4D, where D is the diameter of the cylinder. The length of the channel is L = 25D. Several meshes were tested for h = D/16, D/32, D/64. The time step varies from 0.0005 to 0.005 depending on the Reynolds number and the CFL condition. The particular Reynolds number varies from 0.5 to 100. The stream lines and the vorticity contours for different Re_p are presented in Figure 2. For Re_p = 0.5, the flow is symmetric. For higher Re_p, the inertia becomes important until the apparition of eddies, and at last the periodic detachment of it. Table 1 presents the drag coefficients for different Re_p and meshes. The values are very close to the ones found by Yu and Shao [13]. The results are better for finer meshes. So the flow dynamic around the cylinder and the interaction with the cylinder are well reproduced.

Figure 2.

Flow past a circular cylinder (vorticity on the left and stream line on the right) in a channel (width/diameter = 4) for Re_p = (0.5, 20, 40, 100).

Rep	0.5	1	5	10	20	0	100
Present h = D/16	93.00	46.68	10.28	6.05	3.96	2.84	2.16
Present h = D/32	90.34	46.68	10.01	5.89	3.85	2.75	2.08
Present h = D/64	89.11	44.73	9.88	5.83	3.80	2.72	2.06
Yu and Shao [13]	87.95	44.16	9.95	5.75	3.76	2.68

Table 1.

Drag coefficients for several Reynolds numbers and for several mesh sizes (h = D/16, D/32, D/64).

3.2 Sedimentation of one particle

The second test case concerns the sedimentation of circular particle in a fluid at rest. The gravity is along the channel axis. The domain and numerical parameters are identical to the ones used for the first case. The Froude and Reynolds numbers are based on an estimation of the settling velocity U_c given by Happel and Brenner [20] at very low Reynolds number (Eq. (12)). It depends on the cylinder diameter, the viscosity, the gravity constant and the density of the fluid and particles and a constant K given by (Eq. (13)) with L* = L/D. Simulations for Re_p = 0.1, L* = 4 and ρ_p/ρ_f = 1.1 were performed and presented in Figure 3. The particle accelerates first and then reaches a constant velocity very close to the theoretical value of U_c (Re_p = 0.1) given by Eq. (12).

Figure 3.

Sedimentation of a circular particle (Re_p = 0.1): (a) Isocontour of vorticity at t = 0.25; (b) time development of the settling velocity in a vertical channel.

$$U_C=\frac{D^2}{16\mathit{K\mu }}\left({\rho }_p-{\rho }_f\right)g$$E12

$$K=\frac{1}{ln\left(L^{\ast }\right)-0.9157+1.7244{\left(L^{\ast }\right)}^{-2}-1.7302{\left(L^{\ast }\right)}^{-4}+2.4056{\left(L^{\ast }\right)}^{-6}-4.5913{\left(L^{\ast }\right)}^{-8}}$$E13

3.3 Sedimentation of two particles

This test concerns the sedimentation of two particles in a close domain. The simulations for Re_p = 1, L = 8D, and ρ_p/ρ_f = 1.1 were performed and presented on Figure 4. The cylinders are located at 2D from the axis and their centres are distant of 2D on the z-axis. The domain is of size 8D*160D covers with a uniform mesh of size h = D/16.

Figure 4.

Sedimentation of two particles: (a) trajectories of the particles; (b) equilibrium positions of the particles; (c) time evolution of the dimensionless rotation velocity.

This configuration corresponds to the experiment of Jayaweera and Mason [21]. As described by Feng et al. [22], the trailing cylinder accelerates in the wake of the leading one and then turns around this one until the centre of the cylinders is horizontally aligned. Afterwards, the two particles move apart on the horizontal direction. This phenomenon is found in the simulation. As Feng et al. [22] showed the trailing particle oscillates around the channel axis, whereas the leading one oscillates with smaller amplitude around an axis close to wall. Figure 4c shows that the rotation of the particles is accounted for. The rotational velocity is quite small but not nil. More details about the chaotic sedimentation of two particles at low Reynolds number were presented in Verjus et al. [17].

4. Numerical study

4.1 Sedimentation of a mono-disperse suspension

In this section, the sedimentation of a great number of particles is considered. The settling velocity depends on the number of particles. The more there are, smaller is the velocity. In such a case, the terminal velocity formulation is close to the one proposed by Richardson-Zaki [23]: $U_c=U_{\mathit{ct}}{\left(1-\phi \right)}^5$, where U_ct is the terminal setting velocity of one particle and is the solid volume fraction.

The problem concerns the sedimentation of N_t particles in a close domain. The particles are the same as the previous test (ρ_p/ρ_f = 1.1, Re_p = 0.1). The width is as L = 22D, and the height on which is the initial homogeneous suspension is H₀ = 88D (for the first simulations, then it will be 176D). The mesh size is h = D/16. The time evolution of suspension is presented in Figure 5. After some instants the initial structure is destabilised, three domains are visible: the first one conserves the initial concentration, the second one which is close to the bottom has a solid volume fraction close to the saturation one (about 0.7), and then in the third part intermediate solid volume fraction is present (clearly visible here).

Figure 5.

Snapshot of the sedimentation of an initially homogeneous suspension of 785 particles (L = 22D, H0 = 88D, the initial solid volume fraction is ϕ = 0.32).

A comparison with macroscopic parameter such as the settling velocity of Richardson-Zaki can be made. But a control surface (volume in 3D) on which meaning operation will be made must be defined. So the width of this surface is the length L (width of the domain) and the height is H. Each control volume contains N particles. Spatial meaning and temporal operators can be defined as (T is the period on which the control volume has a solid volume fraction close to the initial value):

$$\tilde{U}\left(t\right)=\frac{1}{N}\sum _1^N{U}_i\left(t\right)\overline{U}\left(t\right)=\frac{1}{T}{\int }_0^T\tilde{U}\left(t\right)\mathrm{d}t$$E14

In Figure 6 (left), the clear water interface is drawn for several solid volume fractions. At the beginning of the motion, the interface settles with a constant velocity and then reaches a regime for which the evolution slows down. This is hindered settling regime. After a while, there is no motion. It corresponds to the gel point. The hindering regime depends on the number of particles, so its time duration is longer for the higher initial solid volume fraction. The right picture presents the evolution of the dimensionless settling velocity versus the solid volume fraction (0.026, 0.058, 0.103, 0.136, 0.233 and 0.32 which correspond to 136, 306, 544, 850, 1204 and 1666 particles) for H0 = 176D. The control volume is placed at the centre of the channel. Qualitatively, the settling velocities resulting from the simulations follow the curve provided by the Richardson-Zaki formula with an exponent of 5. But the best fit of the numerical results is obtained for an exponent of 7.2 ± 0.6. That difference can be due to the 2D-simulation, whereas the Richardson-Zaki formula is made for spherical particles.

Figure 6.

Left: Time evolution of the clear water-suspension interface (settling curve) for three initial solid volume fractions 0.32, 0.162, 0.058 (L = 22D). Right: evolution of the ratio of the calculated settling velocity (W_s) by the one of one particle (W_t) versus the solid volume fraction. The dash line corresponds to the settling velocity calculated with the Richardson-Zaki formula and an exponent of 5.

4.2 Sedimentation of a poly-disperse suspension

In this section, we consider the sedimentation of a poly-disperse suspension. The simulations are two dimensional, and particles are represented by discs with three different diameters: d₁ = d₅₀ 2/3, d₂ = d₅₀ and d₃ = d₅₀ 4/3. The mean diameter is noted d₅₀. All the particles have the same ratio of density ρ_p/ρ_f = 1.5. The computational domain is closed and of sizes 28 d₅₀ × 84d₅₀ respectively in the horizontal and vertical directions. Initially, the fluid and the particles are at rest with a homogeneous distribution as in Figure 7. The particulate Reynolds number Re_p = W_s.d₅₀/ν is equal to 0.1. Re_p is based on the settling velocity obtained with the Richardson-Zaki [23] formula. The initial solid volume fraction is 0.107 which corresponds to 300 particles. The time step is fixed at 0.001, and the fluid mesh is 673 × 2017 nodes which corresponds to a ratio of h = d₅₀/16, where h is the fluid mesh size.

Figure 7.

Snapshots of the sedimentation of tri-disperse particles in a closed box (grey particles for diameter d₁, black ones for diameter d₂ and blue ones for diameter d₃). In this simulation the parameters are: Re_p = 0.1; width/d₅₀ = 28 and the initial volume fraction of 0.107. The contour presents the net pressure (unity: 10⁻² Pa).

Figure 7 presents the motions of particles and the fluid flow for various instants. At the beginning of simulation, the suspension is homogeneous in the domain. After few seconds, most of smallest particles are expulsed upwards, whereas the two other kinds of particles fall on the bottom. However, the biggest particles fall faster than medium ones and then reach the bottom first. In the suspension, two different areas appear: one which is essentially composed of smallest particles in the upper part and a second one constituted with a non-homogeneous mixture of small and medium particles. The bed is constituted in majority of big particles. At the end of simulation, particles are fully segregated, and the bed is composed of three layers filled mainly with respectively one class of particle.

Figure 8 presents the vertical evolution of the net pressure at different instants. This pressure corresponds to the total pressure less the hydrostatic pressure in that case. At the beginning of the simulation, the pressure gradient is constant. Two areas appear at instant t = 25.905 s. Surprisingly, these two areas seem to be linear, and we note that the highest-pressure gradient is in the region near the lower wall. At the end of the simulation, the pressure gradient is constant. In mixture theory (e.g. [24]), the net pressure gradient is calculated with:

Figure 8.

Net pressure versus height at different instants. Right figure is the simulation with the parameters presented in the article. Left figure, the simulation case differs by the number particles by species (114 particles of diameter d₁, 112 particles of diameter d₂ and 74 particles of diameter d₃) which leads to d₅₀ = 1.33.

$$-\frac{\mathit{dp}}{\mathit{dz}}\approx \left({\rho }_p-{\rho }_f\right)\mathrm{\Phi }g$$E15

Where Φ is the solid volume fraction. The pressure gradient obtained with the code for the first time instant is −0.56 Pa/cm and that calculated with (Eq. (15)) is −0.52 Pa/cm. At the end of simulation, when the particles are approximately fixed, the pressure gradient tends to −3.21 Pa/cm. With Eq. (15) and the maximal volume fraction for discs (0.7), the pressure gradient is equal to −3.43 Pa/cm. Results of simulations, with other initial volume fraction confirm this good agreement and seem independent of the polydispersivity. The small difference between Eq. (15) and the results obtained with the code is due to the fact that the volume fraction we calculated depends on the total volume we choose to calculate it. For example, at the beginning of the simulation, the volume fraction is calculated with the total volume of the domain. If only the volume where the particles are located (which give Φ = 0.112) is considered, the accordance with Eq. (15) is better (pressure gradient is equal to −0.55 Pa/cm). Simulations with more particles would reduce the error we made on the calculation of the volume fraction. Moreover, at the end of simulation we have taken the maximal volume fraction. In our cases, the total blockage never appends. Particles in close contact never touch, due to the collision strategy. The flow in the pore remains possible (Cf. Figure 9).

Figure 9.

Particles arrangement in the bed at the end of the simulation. Velocity vectors and contours are plotted (unity: cm/s).

5. Conclusions

A fully resolved model based on the Direct-Forcing/Fictitious Domain method (DF/FD) was developed to study the sedimentation of particles. It was validated on simple cases and applied to the sedimentation of mono-disperse particles. The macroscopic behaviour of the solution corresponds to experimental results and average settling velocities follow a Richardson-Zaki type law. The application of the model to a poly-disperse suspension leads to the apparition of a natural segregation of the particles on the bed. As the pressure is computed by the model, the pressure between the grains can be extract. So the excess pore pressure along the vertical is extract. It appears that at the end of the simulation, but not of the rearrangement of the bed, the pressure profile approaches the law proposed by Guazzelli and Morris [24].

The results presented here are limited to low Reynolds numbers, so we will explore other range of Reynolds numbers. The simulation of sedimentation of particles with complex shapes is possible, but a supplementary work is necessary concerning the collision strategy.

Acknowledgments

We would like to thank the Syndicat Mixte du Cotentin for financing the computational resources.