Open access peer-reviewed chapter

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

Written By

Romuald Verjus and Sylvain S. Guillou

Submitted: 12 February 2022 Reviewed: 24 March 2022 Published: 27 April 2022

DOI: 10.5772/intechopen.104654

From the Edited Volume

Modeling of Sediment Transport

Edited by Davide Pasquali

Chapter metrics overview

163 Chapter Downloads

View Full Metrics


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.


  • polydisperse particulate flows
  • sedimentation
  • direct numerical simulation
  • modelling

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, Vs, Js, U and ω, respectively. The fluid density and viscosity are ρf and μ. Let us introduce scale parameters: Lc for length, Uc for velocity, Lc/Uc for time, ρfUc2 for the pressure and ρfUc2/Lc for the pseudo body force. Then, the dimensionless fictitious domain equations can be written as Yu and Shao [13]:


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 ρfUcLc, Fr is the Froude number defined by Uc2/gLc, ρr is the particle-fluid density ratio defined by ρsf, Jd is the dimensionless moment of inertia defined by JssLc4 and Vd the dimensionless volume defined by Vs/Lc2.

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


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

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


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


• Step 3: Correction of the fluid velocity un + 1 regarding λn + 1 and λn


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 Rep are presented in Figure 2. For Rep = 0.5, the flow is symmetric. For higher Rep, 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 Rep 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 Rep = (0.5, 20, 40, 100).

Present h = D/1693.0046.6810.286.053.962.842.16
Present h = D/3290.3446.6810.015.893.852.752.08
Present h = D/6489.1144.739.885.833.802.722.06
Yu and Shao [13]87.9544.169.955.753.762.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 Uc 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 Rep = 0.1, L* = 4 and ρpf = 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 Uc (Rep = 0.1) given by Eq. (12).

Figure 3.

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


3.3 Sedimentation of two particles

This test concerns the sedimentation of two particles in a close domain. The simulations for Rep = 1, L = 8D, and ρpf = 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]: Uc=Uct1φ5, where Uct is the terminal setting velocity of one particle and is the solid volume fraction.

The problem concerns the sedimentation of Nt particles in a close domain. The particles are the same as the previous test (ρpf = 1.1, Rep = 0.1). The width is as L = 22D, and the height on which is the initial homogeneous suspension is H0 = 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, H= 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):


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 H= 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 (Ws) by the one of one particle (Wt) 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: d1 = d50 2/3, d2 = d50 and d3 = d50 4/3. The mean diameter is noted d50. All the particles have the same ratio of density ρpf = 1.5. The computational domain is closed and of sizes 28 d50 × 84d50 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 Rep = Ws.d50 is equal to 0.1. Rep 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 = d50/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 d1, black ones for diameter d2 and blue ones for diameter d3). In this simulation the parameters are: Rep = 0.1; width/d50 = 28 and the initial volume fraction of 0.107. The contour presents the net pressure (unity: 10−2 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 d1, 112 particles of diameter d2 and 74 particles of diameter d3) which leads to d50 = 1.33.


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.



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


  1. 1. Guillou S, Thiébot J, Chauchat J, Verjus R, Besq A, Nguyen DH, Pouv KS. The filling dynamics of an estuary: From the process to the modelling.“Sediment Transport in Aquatic Environments”.Ed. A. Manning. London, UK: IntechOpen; 2011. Chap. 6, p. 125-146. ISBN 978-953-307-586-0.
  2. 2. Rtimi R, Sottolichio A, Tassi P. Hydrodynamics of a hyper-tidal estuary influenced by the world’s second largest tidal power station (Rance estuary, France). Estuarine, Coastal and Shelf Science. 2021;250:107143. DOI: 10.1016/j.ecss.2020.107143
  3. 3. Thiébot J, Guillou S, Brun-Cottan JC. An optimisation method for determining permeability and effective stress relationships of consolidating cohesive sediment deposits. Continental Shelf Research. 2011;31(10) Supplement 1:S117-S123. DOI: 10.1016/j.csr.2010.12.004
  4. 4. Hsu T, Jenkins JT, Liu PLF. On two-phase sediment transport: Dilute flow. Journal of Geophysical Research. 2003;108(C3):3057. DOI: 10.1029/2001JC001276
  5. 5. Amoudry L, Hsu TJ, Liu PLF. Two-phase model for sand transport in sheet flow regime. Journal of Geophysical Research. 2008;113:C03011. DOI: 10.1029/2007JC004179
  6. 6. Nguyen KD, Guillou S, Chauchat J, Barbry N. A two-phase numerical model for suspended-sediment transport in estuaries. Advances in Water Resources. 2009;32:1187-1196. DOI: 10.1016/j.advwatres.2009.04.001
  7. 7. Chauchat J, Guillou S. On turbulence closures for two-phase sediment-laden flow models. Journal of Geophysical Research. 2008;113:C11017. DOI: 10.1029/2007JC004708
  8. 8. Chauchat J, Médale M. A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow. Computer Methods in Applied Mechanics and Engineering. 2010;199:439-449. DOI: 10.1016/j.cma.2009.07.007
  9. 9. Nguyen DH, Levy F, Pham Van Bang D, Guillou S, Nguyen KD, Chauchat J. Simulation of dredged sediment releases into homogeneous water using a two-phase model. Advances in Water Resources. 2012;48:102-112. DOI: 10.1016/j.advwatres.2012.03.009
  10. 10. Khaled F, Guillou SS, Méar Y, Ferhat H. Impact of blockage ratio on the transport of sediments in the presence of a hydrokinetic turbine: Numerical modelling of the interaction sediments-turbine. International Journal of Sediment Research. 2021;36:696-710. DOI: 10.1016/j.ijsrc.2021.02.003
  11. 11. Glowinski R, Pan T-W, Hesla TI, Joseph DD. A distributed lagrange multiplier/fictitious domain for particulate flows. The International Journal of Multiphase Flow. 1999;25:755-794. DOI: 10.1016/S0301-9322(98)00048-2
  12. 12. Peskin CS. The immersed boundary method. Acta Numerica. 2002;11:479-517. DOI: 10.1017/S0962492902000077
  13. 13. Yu Z, Shao X. A direct-forcing fictitious domain method for particulate flows. Journal of Computational Physics. 2007;227:292-314. DOI: 10.1016/
  14. 14. Wachs A. A DEM-DLM/FD method for direct numerical simulation of particulate flows: Sedimentation of polygonal isometric particles in a Newtonian fluid with collisions. Computers & Fluids. 2010;38(8):1608-1628. DOI: 10.1016/j.compfluid.2009.01.005
  15. 15. Davis RH. In: Tory EM, editor. Velocities of Sedimenting Particles in Suspension, in: Sedimentation of Small Particles in Viscous Fluid, Advances in Fluid Mechanics. Southampton: Computational Mechanics Publications; 1996. pp. 161-198
  16. 16. Guillou S, Makhloufi R. Effect of a shear-thickening rheological behaviour on the friction coefficient in a plane channel flow: A study by direct numerical simulation. Journal for Non-Newtonian Fluid Mechanics. 2007;144:73-86. DOI: 10.1016/J.JNNFM.2007.03.008
  17. 17. Verjus R, Guillou SS, Ezersky A, Angilella JR. Chaotic sedimentation of particle pairs in a vertical channel at low Reynolds number: Multiple states and routes to chaos. Physics of Fluids. 2016;28:123303. DOI: 10.1063/1.4968559
  18. 18. Verjus R. Contribution à la modélisation des processus de sédimentation-consolidation dans les ports et les estuaires: Etude bi-échelle. [thesis]. University of Nancy; 2015
  19. 19. Yu Z, Phan-Thien N, Fan Y, Tanner RI. Viscoelastic mobility problem of a system of particles. Journal for Non-Newtonian Fluid Mechanics. 2002;104:87-124. DOI: 10.1016/S0377-0257(02)00014-9
  20. 20. Happel J, Brenner H. Low Reynolds Number Hydrodynamics, 552P. The Hague: Martinus Nijhof; 1973
  21. 21. Jayaweera KOLF, Mason BJ. The behaviour of freely falling cylinders and cones in a viscous fluid. Journal of Fluid Mechanics. 1965;22:709-720. DOI: 10.1017/S002211206500109X
  22. 22. Feng J, Hu HH, Joseph DD. Direct simulation of initial-value problems for the motion of solid bodies in a Newtonian fluid. 1. Sedimentation. Journal of Fluid Mechanics. 1994;261:95-134. DOI: 10.1017/S0022112094000285
  23. 23. Richardson JF, Zaki WN. Sedimentation and fluidization: Part I. Transactions of the Institution of Chemical Engineers. 1954;32:35-53. DOI: 10.1016/S0263-8762(97)80006-8
  24. 24. Guazzelli E, Morris JF. A Physical Introduction to Suspension Dynamics. Cambridge: Cambridge University Press; 2012

Written By

Romuald Verjus and Sylvain S. Guillou

Submitted: 12 February 2022 Reviewed: 24 March 2022 Published: 27 April 2022