Dedico aquest treball, el més gran que he fet mai, als meus pares, Ana i Juan Ignacio, i a la Vera, el meu amor
I hereby declare that except where specific reference is made to the work of others, the contents of this dissertation are original and have not been submitted in whole or in part for consideration for any other degree or qualification in this, or any other university. This dissertation is my own work and contains nothing which is the outcome of work done in collaboration with others, except as specified in the text and Acknowledgements.
Writing this doctoral dissertation often felt like writing a book, perhaps a long one. A somewhat lonely activity, a fight to keep a focus for a long time chasing a quality standard that keeps getting away. And as so often writers explain, this is best achieved in solitude: it is a lonely job.
On the other hand, one is not ready to write it the first day. I arrived at CIMNE knowing very little about programming, or physics, or about how to do research. I must thank Eugenio Oñate for, in spite of this, giving me the opportunity to enter the world of numerical methods and science in general. Eugenio has taken care of my career and has cultivated my confidence in a very generous way over these years. I am especially grateful for his welcoming me a second time after I came back from my failed entrepreneurial adventures.
A year into my doctorate, Riccardo Rossi kindly took me in as his candidate as codirector. He has been the one to listen to my worries and patiently help me every time I got stuck.
An advantage of working in CIMNE is that there is no shortage of smart researchers, always ready to discuss your problems. And I certainly took advantage of this during the course of this work: First of all, thank you Jordi Cotela. I owe most of my work on the extension of the finite element code to you (and also my understanding of subscale stabilization). You provided me with code, lessons and fruitful discussions. I am sorry I could not finish the work on turbulence we were collaborating on in timeI must also thank my other labmates: Julio Marti, Pavel Ryzhakov and Ricardo Reyes for their many discussions and advice.
Even before I began my doctorate, I had entered the DEMPack team. I owe a lot to my colleagues there too: Miguel Ángel Celigueta, with whom I have collaborated in many developments and held numerous stimulating discussions on programming and on the discrete element method; to Salvador Latorre, with whom I have shared many jobs; and Ferran Arrufat, who helped me with crucial developments at the end of this work. I am indebted to these people, as well as to Merce López and Joaquin Irazábal for covering my back so many times when I was working on this manuscript. It has been a pleasure sharing all those good times during the lunch breaks with these friends.
I must thank several other researchers who have offered their help. I do not want to forget any of them: Roberto Flores, a friend and a such a talented scientist, who helped me sort out my ideas in several occasions; Enrique Ortega, who gave me valuable ideas on how to use polynomial smoothers; Ramón Codina and Joan Baiges, who helped me understand concepts on subscales stabilization and other questions; Prashanth Nadukandi, who helped me with matrix manipulations and with whom I enjoyed discussing possible research topics; and Pooyan Dadvant and Carlos Roig, who gave me programming solutions I would have taken an eternity to find.
I am convinced that the best research is done in collaboration. I must thank the very fruitful (and pleasant) collaborations I have had with Alex Ferrer, who played a crucial role in the optimization process described in the third chapter and Ignasi Pouplana, with whom we worked on particle impact drilling.
On occasions, I encountered difficulties when reading crucial sections of the literature. Thankfully, most times the authors were kind to answer my questions. I want to thank Benoit Pouliot, Ksenia Guseva and Andrew Bragg for their patient and elaborate answers.
I spent two extremely productive months at UP Berkeley thanks to Tarek Zohdi, who kindly hosted me two times and collaborated with me in writing a paper. I thank him and Eugenio for making that exchange possible.
My research would have not been possible without financial support from Generalitat de Catalunya, under Doctorat Industrial program 2013 (DI 024). I also thank the support received from the project PRECISE  Numerical methods for PREdicting the behaviour of CIvil StructurEs under water natural hazards (MINECO  BIA201783805R  01/01/2018 – 31/12/2020)
Finally, I want to thank my family for their extreme patience and constant love and support.
In this work we study the numerical simulation of particleladen fluids, with an emphasis on Newtonian fluids and spherical, rigid particles.
Our general strategy consists in using the discrete element method (DEM) to model the particles and the finite element method (FEM) to discretize the continuous phase, such that the fluid is not resolved around the particles, but rather averaged over them. The effect of the particles on the fluid is taken into account by averaging (filtering) their individual volumes and particlefluid interaction forces.
In the first part of the work we study the Maxey–Riley equation of motion for an isolated particle in a nonuniform flow; the equation used to calculate the trajectory of the DEM particles. In particular, we perform a detailed theoretical study of its range of applicability, reviewing the initial effects of breaking its fundamental hypotheses, such as small Reynolds number, sphericity of the particle, isolation etc. The output of this study is a set of tables containing orderofmagnitude inequalities to assess the validity of the method in practice.
The second part of the work deals with the numerical discretization of the MRE and, in particular, the study of different techniques for the treatment of the historydependent term, which is difficult to calculate efficiently. We provide improvements on an existing method, proposed by van Hinsberg et al. (2011), and demonstrate its accuracy and efficiency in a sequence of tests of increasing complexity.
In the final part of the work we give three application examples representative of different regimes that may be encountered in the industry, demonstrating the versatility of our numerical tool. For that, we describe necessary generalizations to the MRE to cover problems outside its range of applicability. Furthermore, we give a detailed account of the stabilized FEM algorithm used to discretize the fluid phase and compare several derivative recovery tools necessary to calculate some of the interphase coupling terms. Finally, we generalize the algorithm to include the backwardcoupling effects according to the theory of multicomponent continua, allowing the code to deal with arbitrarily dense flow regimes.
Roman Symbols
 particle's radius
 particle crosssectional area
 body (continuum mechanics)
 body belonging to component (continuum mechanics)
 drag coefficient in Shah's (2007) model
 local sound speed
 Stokes–Einstein–Smoluchowski diffusivity
 threedimensional affine Euclidean space
 friction factor (boundary layer)
 added mass force
 characteristic scale of
 total body work
 sum of contact forces
 steady drag force
 characteristic scale of
 Boussinesq–Basset history force
 characteristic scale of
 hydrodynamic force
 Momentum exchange between component and the rest of components
 Herron lift force
 Rubinow and Keller lift force
 characteristic value of
 Saffman lift force
 Froude number
 total surface work
 undisturbedflow force
 characteristic scale of
 static submerged weight
 gravity acceleration
 space of functions with square integrable derivatives in
 space of functions in vanishing on
 space of vectors in fulfilling the Dirichlet boundary conditions
 spheric particle's moment of inertia
 consistency index
 Boltzmann's constant
 Boussinesq–Basset kernel function
 Knudsen number
 characteristic length scale of the smallest flow structures
 space of square integrable functions in
 characteristic filter length scale
 total momentum
 displaced fluid's mass
 particle mass
 equivalent particle mass
 behaviour index
 total moment of momentum
 number of space dimensions
 number density
 total number of mesh nodes
 number of time steps between quadrature steps
 number of space dimensions
 neighbourhood function
 power set of set
 fluid pressure
 Péclet number associated to the settling velocity
 heat flux
 space of pressure scalars
 discrete counterpart of
 particle position
 field of the real numbers
 strictly positive real numbers
 internal energy source/sink
 Reynolds number
 rotational Reynolds number
 particle Reynolds number
 Reynolds number for powerlaw fluids
 particle Reynolds number for powerlaw fluids (Metzner and Reed, 1955)
 Shearrelated particle Reynolds number
 Cartesian product involving copies of
 strain rate tensor
 settling coefficient
 Brinkman over Oseen length scales squared
 Schmidt number
 sphere centred at zero with radius equal to
 Stokes number
 Stokes number for bubble rebound
 collisional Stokes number
 rotational Stokes number
 scaledependent Stokes number
 surface traction on
 characteristic time scale of the smallest flow structures
 mass balancerelated stabilization parameter
 momentum balancerelated stabilization parameter
 matrix of stabilization parameters
 total torque due to body actions
 sum of contact torques
 discrete time variable for the particles
 characteristic filter time scale
 hydrodynamic torque
 total torque due to surface actions
 fluid velocity
 characteristic velocity scale of the smallest flow structures
 internal energy
 total energy
 particle velocity
 space of velocity vectors
 discrete counterpart of
 particle volume
 slip velocity
 power spent by the body forces
 power spent by the surface forces
 space of solutions
 space of subscales
 space of subscales that vanish at the boundary
 discrete counterpart of
 dimensionless distance for powerlaw fluids
Greek Symbols
 particle mass fraction
 Basset slip coefficient
 configuration (continuum mechanics)
 motion of a body (continuum mechanics)
 motion of a body of component (continuum mechanics)
 configuration at time (continuum mechanics)
 configuration of component at time (continuum mechanics)
 continuous phase time step
 disperse phase time step
 Levi–Civita symbol
 amplitude of surface irregularities
 radial distribution function
 local strain rate
 subset of where Dirichlet boundary conditions apply
 subset of where Neumann boundary conditions apply
 mean free path of the fluid molecules
 fluid's dynamic viscosity
 yield viscosity
 fluid's kinematic viscosity
 molecular collision frequency
 continuum phase domain
 particle's angular velocity
 characteristic value of
 characteristic scale of the vorticity in pureshear flow
 moment exchange of component with the rest of components
 relative density of the disperse component
 fluid's density
 Cauchy's stress tensor
 tangential accommodation coefficient
 deviator stress tensor
 sonic time scale
 molecular diffusion time scale
 viscous diffusion time scale
 rotational particle relaxation time
 particle response time
 shear stress at the wall
 dimensionless frequency normalized by
 absolute temperature
 inverse of the particle mobility
 activation coefficient for term
Other symbols
 filtering operator over the component
 filtering operator per unit component volume over the component
 material derivative of the fluid continuum
 Laplacian operator
 boundary of set
 inner product
 integral of the product of two functions
 component variable over volume fraction
Acronyms / Abbreviations
ALE  arbitrary LagrangianEulerian method
ALU  arithmetic logic unit
ASGS  algebraic subgrid scales
BEM  boundary element method
CFD  computational fluid dynamics
CFL  Courant–Friedrichs–Lewy number
COR  coefficient of normal restitution
DEM  discrete element method
DNS  direct numerical simulation
DOF  degree of freedom
EE  Euler–Euler (method)
EL  Euler–Lagrange (method)
FDE  fractional differential equations
FEM  finite element method
FFC  recovery method of Pouliot et al. (2012)
FFT  fast Fourier transform
FLOP  floating point operations
FPU  floating point unit
FSI  fluidstructure interaction
FVM  finite volume method
LBM  lattice Boltzmann method
LES  large eddy simulation
MAE  method of approximation by exponentials
MFP  mean free path
MPM  material point method
MRE  Maxey–Riley equation
MSD  mean square displacement
OSS  orthogonal subscales
PFEM  particle finite element method
PIC  particleincell
PID  particle impact drilling
PPC  particles per cell
PPR  recovery method of Zhang and Naga (2005)
PSM  pseudospectral method
RDF  radial distribution function
RVE  representative elemental volume
SC  standard case
SM  streaming multiprocessors
VMS  Variational multiscale (method)
How is a moving particle affected by its surrounding fluid? This question has intrigued scientists since the Scientific Revolution and earlier, going back to Galileo [325], Leonardo Da Vinci [144], and even Aristotle [293]; perhaps because it refers to a simple system, ideal for thought experiments. The problem attracted attention for pragmatical reasons early on as well, as it is central to the study of ballistics (the science of predicting the trajectory of projectiles), which began to be systematized as a discipline by Tartaglia in the sixteenth century. Despite his ignorance of Newtonian dynamics, Tartaglia was able to provide useful range tables and instructions, although his understanding of the problem was still very rudimentary. For instance, he considered the trajectory of projectiles to be formed by the succession straightcircumference archstraight, influenced by earlier tradition [335]. About a century later, Galileo provided the parabolic solution for a projectile in void; noting, in 1638, that he expected his solution to be valid only for slowmoving projectiles, while faster moving trajectories should make the beginning of the parabola less tilted and curved than its end. Thus, for a very long time, the intuitive notion that a resistance existed proved too difficult to translate into scientific knowledge.
In 1671 Newton wrote the following text in passing, as an analogy to describe the motion of his light corpuscles, when explaining the breakdown of a white ray into its different colors [259]^{1}:
…I remembred that I had often seen a Tennisball struck with an oblique Racket describe such a curve line. for a circular as well as a progressive motion being communicated to it by that stroak, its parts on that side where the motions conspire must presse & beat the contiguous air more violently then on the other, & there excite a reluctancy & reaction of the air proportionally greater.
The quotation reflects Newton's very early insight into the fundamental mechanisms of fluidparticle interaction, made possible only thanks to the idea of force, the notion of continued action (as opposed to discrete) and its necessity to have the ball deviate from a straight line (law of inertia), and the actionreaction postulate; all of which were recent discoveries at the time. Moreover, it exemplifies the reductionist approach, the process of continued zooming in that would ultimately give the scientific method its modern power and that particularly resonates with the common theme of numerical methods: 1) break down into simpler entities, 2) model the interactions, 3) aggregate the result of all these interactions.
By running experiments on pendulums with spherical bobs, Newton himself was able to establish that the air resistance was approximately proportional to the square of the velocity in his Principia [260], a result still standing today (for fastmoving particles, in the socalled Newtondrag regime; see e.g., [218]). Together with his laws of motion, this finding was to set the science of ballistics on a much more solid ground^{2}.
Still, little progress could be made from the theoretical point of view for many years, due to limited understanding of fluid dynamics. Famously, in 1752 D'Alembert proved the paradoxical result that inviscid flows, for which Euler had correctly given their equations of motion, yielded zero resistance on submerged objects. Almost a century later, SaintVenant (who had actually derived the Navier–Stokes equations before Stokes) took the first steps toward its resolution by suggesting that it was necessary to take the fluid viscosity into account, as it would surely introduce some tangential resistance [7].
But it was not until the work of Stokes [326] that the first (correct) quantitative solution to the problem could be derived from first principles, for the case of very small Reynolds number^{3}. Stokes calculated the drag force on a small sphere moving at steady speed in an infinite Newtonian fluid using a simplified version of his own discovery, the Navier–Stokes equations. The force turned out to depend linearly on the velocity in this regime, in contrast to the theory of Newton, whose pendulums operated at higher Reynolds numbers. Insightfully, Stokes observed that
Since in the case of minute globules falling with their terminal velocity the part of the resistance depending upon the square of the velocity, as determined by the common theory, is quite insignificant compared with the part which depends on the internal friction of the air, it follows that were the pressure equal in all directions in air in the state of motion, the quantity of water which would remain suspended in the state of cloud would be enormously diminished. The pendulum thus, in addition to its other uses, affords us some interesting information relating to the department of meteorology.
Where the reference to a pendulum is due to Stokes being in fact primarily interested in solving that problem, from which he derived his famous drag relation as the limiting case of an infinite oscillation period. As it turns out, the problem of cloud formation is still a subject of current studies, most often based on his equation [70,150].
Indeed, Stokes' result became (and remains) very important to fundamental as well as applied science, having lead to at least three Nobel prizes, according to Dusenbery [111]. For instance, Einstein used it to calculate the Brownian diffusion rate (see [290]), a result that provided the first strong evidence of the existence of atoms (i.e., strictly speaking, water molecules) as well as a practical means to calculate their mass!
Most analytical successes in studying this problem have been achieved in the realm of the tiny to the very small; in a world where viscosity counteracts all external forces and where inertia is, to all effects, negligible (or at most a small correction). Neglecting inertia renders the NavierStokes equations linear, which removes the tremendous difficulties associated with nonlinear effects, including the analytical intractability of problems and the onset of turbulent chaos. For large particles, when the Reynolds number grows much past one, we have not been able to go too far past the empirical solution of Newton.
But even under low Reynolds number conditions, the analytical difficulties do appear and it took many years until a complete picture was available. Successive advances where made over several decades by many prominent scientists, such as Boussinesq (1885), who derived the unsteady terms that arise when the particle decelerates in otherwise the same conditions as Stokes; Oseen (1910, 1913), who correctly found the first effects of inertia on the steady drag (although not entirely satisfactory, see 2.2.1); Faxén, who accounted for the nonuniformity of the background flow and Saffman (1965), who calculated the lift force in the presence of a linear shear flow (see the historical account by [248]). The equation of motion that we will study in the following two chapters, the Maxey–Riley equation, that generalizes the works of Boussinesq and Faxén was not derived until the year 1983.
With this brief historical account we hope to motivate our interest in studying the motion of a single particle in a fluid. One of the lessons learned during the course of these years of research is to what extent this is a fundamental theoretical piece towards understanding particleladen flows. Although initially unplanned, such endeavour has ended up taking a large part of our work (2,3).
(^{1}) Specifically, Newton is addressing the origins of the lateral force induced by simultaneous rotation and translation of a spherical particle, causing it to drift, curving its trajectory sideways. This phenomenon would later be rediscovered and empirically studied by Robins in the eighteenth century (although the phenomenon is commonly known as the 'Magnus effect', after Magnus, who studied it much later and who, in fact, cites Robins in his work [229]) [219,325].
(^{2}) Although Newton was unable to improve upon Galileo's parabolic solution by introducing the quadratic air resistance into the equations. The first accurate ballistics tables where produced around mideighteenth century by Robin, who obtained them empirically, and by Euler who used Robin's data to fix the required empirical constants and a numerical method to integrate the nonlinear equation describing the trajectories based on Newton's squaredvelocity drag law [325].
(^{3}) This still did not resolve d'Alembert's paradox, since the fact that the NavierStokes equations coincide with Euler's equations in the limit of very large Reynolds number (large velocities over viscosity ratio), which yielded zero drag, seemed inconsistent with the known empirical observation that the drag did not actually tend to zero with viscosity at high velocities. A resolution of this problem did not come until 1904 with the work Prandtl, see [6], who introduced the concept of the boundary layer that provided a sound transition between both conflicting models.
With the growing computational power available and the development of multiphase solvers, one may be tempted to drop the analytical route altogether and simply calculate the forces on individual particles by solving the Navier–Stokes equations around them. A particleladen flow system is, after all, a set of solid bodies submerged in a fluid, and could in principle be simulated by using their surface to impose the corresponding boundary conditions on the fluid as in any fluidstructure interaction problems.
In fact, researchers are increasingly following this route, producing fullyresolved simulations to study systems of submerged particles (the number of particles in each study is indicated between parentheses): In Johnson and Tezduyar [186] the settling of a sphere in a container was studied with a bodyfitted mesh (); in [273] the fluidization of a bed of spheres was studied numerically and experimentally (); in [138] the latticeBoltzmann method was used to simulate spheres setling in a stationary fluid (); while in [349] an embeddedbody approach was used to study particle clustering in turbulent flows (). These methods are particularly helpful in basic research, as they can be used to calibrate empirical models [384], to study qualitative aspects of the flow, including turbulent effects [344] and to study complex systems including nonNewontian flow models and particles interaction in exquisite detail [107,250].
But the numbers in the previous paragraph are still clearly way too small for the vast majority of particleladen flows of interest: sand grains or pebbles in a river bed, soil slurries, avalanches and pyroclastic flows; sedimentation tanks and vessels used to clarify water or separate small particles by size; pneumatic conveying systems carrying sugar, flour, coal seeds, nuts and conglomerate pellets and in general solid matter in gas or liquid ducts; fluidized beds, used in the chemical industries to enhance chemical reactions and homogenize and dry particulates; bubbles in nuclear reactors, oil, rock and gas bubbles in oil wells; microscopic biological systems like blood cells and platelets in small vesicles or suspended particles in the bronchioles, such as contaminants or vaporised sprays; plankton and organic matter in the ocean [231]; drops in a forming cloud, snow crystals and solid contaminant, like rubber particles, in the atmosphere; and even rocky particles in lowdensity atmospheres, relevant to the study of planet formation [358].
Thus, one soon realizes that this approach remains hopeless as a tool to study most systems as a whole; and it is expected to remain so for a long time, unless some unexpected, revolutionary discovery leads to a jump of several orders of magnitude in computational power.
On the other hand, note that if we are able to give up the resolution of the flow around the particles and instead modelling the hydrodynamic interactions using the coarse scale description of the flow we are able to immediately generate a radical cutdown in the computational cost.
A quick calculation illustrates this point very convincingly: First, take a suspension of particles that we will consider to be of uniform size for simplicity. Let us assume that discretization size required to run a fully resolved simulation (where all the scales present in the fluid motion are considered) to be one tenth of the diameter of the particles. Note that it is unreasonable to think that a substantially coarser size should be sufficient to capture the dynamics of the flow around the particles in sufficient detail.
On the other hand, consider the same problem solved by an alternative, coarsegrained approach, where the hydrodynamic interactions are modelled as a distributed effect on the fluid, and where only the details of the flow much larger than the particles (say, ten times as big) are numerically resolved. Note that with such approach it would hardly make sense to consider finer resolutions at all because, at scales comparable to the particle size, the real flow would be locally distorted by the presence of the individual particles that are being averaged over. But without the information about the individual distortions the solution would necessarily be very poorly modelled at such fine scales. In other words, it would just not pay off to consider that level of detail.
Comparing both situations, one has that the number of computational cells required scales as in the first case and in the second case. That is a factor 1x10^{6} in number of cells between the two approaches, which amounts to a significantly greater number in terms of computational resources^{1}.
Note that such strategy corresponds to a mixedscale method, where one phase (the particles) is described to a finer level of detail (individual particles) than the other (the fluid), which is described to a coarser level. This work is concerned with precisely this type of mixed methods.
Clearly, one can continue the process of coarsening of the description level further by homogenizing the particles phase as well, giving rise to a multicomponent continuum description of the particulate system. The method most consistent with this philosophy is the socalled Euler–Euler (EE) method [180], also called the EulerianEulerian or twofluid model; as opposed to the Euler–Lagrange (EL) method, which refers to the mixed method discussed above. In the EE method the motion of the particles is determined by a set of conservation differential equations, analogous to those of the fluid phase. Thus, both the particles and the fluid phases have their associated systems of conservation equations. The equations that correspond to each phase are coupled to each other through the averaged volume fractions (volumetric proportions), whose unknown values vary pointwise, as well as by momentum exchange terms (see Appendix H).
The equations involved can be derived following a purely formal procedure that relates the different averaged (largescale) variables of interest. Several averaging methods exist, including time, volume and ensemble averages, but they all lead to similar equations [109]. As a consequence, and in a completely analogous situation to that encountered in turbulence modelling, the resulting equations include a number of terms that depend on finerscale details, and that must be closed (expressed as a function of the averaged variables) to obtain a wellposed problem. This closure is normally based on the introduction of new physical models that most often include unknown parameters to be calibrated from experiments [182]. These methods have been applied to the description of particleladen flows for many years; especially in the chemical [332,139,192] and the nuclear [196,158] industries.
In fact, EEtype methods remain the only realistic approach in many industrial settings. The reason is found again in the large numbers involved, which penalize the explicit description of each and every particle in the domain. This is a limitation common to all particlebased methods, including the discrete element method (DEM), since most industrial problems involve a much larger number of particles than what is computationally feasible.
The latter point is illustrated in Fig. 1, where a bunch of industrial applications are identified in a domain sizenumber of particles graph. As a reference, simulations involving several million particles (green curve) are representative of what can typically be achieved in a personal computer and a few billion (red curve) particles that corresponding to what can ordinarily be achieved in powerful computer clusters, using massive parallelization techniques. The violet 2x10^{12}curve is representative of the numbers of particles included in today's largest cosmological simulations [281] and can be interpreted as an extremely optimistic upper bound. By mentally locating one's preferred application in the plot of Fig. 1 one quickly arrives at an intuitive grasp of the limitations of particlebased methods as a numerical tool to simulate real systems with realistic numbers of particles.
Figure 1: Number of particles associated to cubic domains for a fixed proportion of volume occupied by the particles (1x10^{1}) as a function of the equivalent side length of the domain. The curves corresponding to different numbers of particles are included, as well as the upper limit sizes for clay, silt and sand, as representative granular materials. The scattered points correspond to different industrial application examples: coal, salt and nuts pneumatic conveying (volume per linear meter for typical pipe diameters); fluidized bed experiments of Geldart [143], showing the different Geldart categories; and catalytic cracking (FCC) for the production of gasoline. 
Nonetheless, it must be mentioned at this point that there exists the possibility of using a coarsegrained description, where a smallerthanreal number of (larger) particles is used instead of the onetoone description [125,304,271] in analogy to what is common practice in molecular dynamics [295,176]. And, while the level of maturity of this approach is still low, with several problems still unresolved [200,258], the technology seems to be yielding promising results with spectacular computational gains [225].
Furthermore, even if we ignore the possibility of coarsegraining, hybrid methods and, in particular, CFDDEM methods retain an interest due to their many advantages over exclusively continuumbased methods:
So, although the number of realworld applications where hybrid methods can be applied is in practice limited, it is important to recognize their value as validation tools that operate at an intermediate scale between fully resolved models and continuumonly models. Therefore, in addition to their value as direct simulation tools for the industry, they are useful to feed the empirical data needed to close EEtype equations and to understand the fullsize objects of study through the analysis of simplified, smallsize systems [132].
(^{1}) In computational fluid dynamics, the computational cost is mostly determined by the resolution of the system of equations resulting from the spatial discretization, whose size scales as the number of unknowns. The computational cost of the resolution of the system can at the very best be linear for the latestgeneration iterative solvers [146] (grow proportionally to the number of unknowns), although such scaling in practice is difficult to achieve. Furthermore, the smallest time step also needs to be modified, typically decreased proportionally to the element size, in order to preserve the accuracy increase desired; see the discussion around the Courant–Friedrichs–Lewy condition in Section 4.7.7
This dissertation has been completed with financial support from the Doctorat Industrial grant from the Generalitat de Catalunya, that supports research work performed as part of the work of a company in collaboration with a research center and with a focus on applied research. In our case, the research center role was played by the International Center for Numerical Methods in Engineering (CIMNE) and the role of the company was played by Computational and Information Technologies S.A. (CITECHSA).
The main objective of CITECHSA was to start the development of a computational fluid dynamics (CFD) tool able to simulate a wide range of particleladen flows. Our realization at the time of beginning the work that a new DEM application tool was at that point at an early stage of development by the DEM group of Kratos (see next section) motivated our choice of the DEM for the simulation of the particles' phase. This choice also created a strong dependence between the DEM developments and those concerning the coupling, that explained the author's enrolment into the DEM group in CIMNE at around that time.
As a consequence, DEMrelated tasks have occupied a substantial part of the research time. However most of the resulting algorithmic developments do not go beyond a state of the art, which is sufficiently documented. Furthermore, the application examples and analyses of interest are too detached from the main thread of this work to include them. Thus only a very brief account of some DEMrelated elements are described in Appendix A. Nonetheless the reader should keep in mind the DEM background of the author, and the (perhaps) clear emphasis given to the particles' side of the research.
The objectives of the research reported in this document are the following:
Our research activities have consisted in a combination of bibliographical investigation, numerical analysis (designing and running numerical experiments using existing or new code, mainly on a personal computer and sometimes on a small cluster) and programming.
The bulk of the programming work has been carried out within the framework Kratos Multiphysics [88], or Kratos for short. Kratos is based on C++
, with a Pythonbased external layer. It is organized according to a marked objectoriented philosophy that, at its coarsest level, can be modelled as a central core (Kratos Core) connected to a list of independent applications (FluidDynamicsApplication, StructuralApplication, DEMApplication, etc.).
The core contains the definition of the fundamental abstract classes common to a large part of the implementations. It provides a common language and a set generic protocols and algorithmic tools like linear solvers, search algorithms or input/output utilities.
On the other hand, the developments related to specific problems, such as solving the Navier–Stokes equations or performing structural analyses, are coded in the corresponding specialized applications (in this case perhaps FluidDynamicsApplication and StructuralApplication). Moreover, an application can relate any subset of other, existing applications to solve a more complex problem. In order to fulfil our particular goals, we have created SwimmingDEMApplication, that couples FluidDynamicApplication [82,311]^{1} to DEMApplication [64,177].
FluidDynamicsApplication is an Euleriandescription, finite element methodbased CFD simulator and it is in charge of solving the fluid equations of motion. DEMApplication is a generalpurpose discrete element method simulation suite that is responsible of tracking the particles, computing contacts between them and evolving their motion in time.
SwimmingDEMApplication contains all the interphase coupling tools, the hydrodynamic calculation algorithms, the modified fluid finite elements necessary to solve the modified backwardcoupled flow equations (see Chapter 4) and several utilities that were developed where no similar tool was found in the core, but that have not yet reached a stage of maturity to make them suitable for their implementation in the Kratos core (for instance, the derivative recovery tools, see Section 4.4).
The whole of Kratos, including all the applications mentioned above and many others, is freely available for download^{2}. FluidDynamicApplication has an interface usable from the Pre/Postprocessor GiD [244], that, by default, is included in Kratos. Moreover, DEMApplication powers any of the four free packages grouped under the name DEMPack [65], also based on the GiD interface and also included in Kratos. Among these, FDEMPack gives access to most of the capabilities of SwimmingDEMApplication, developed during the course of the research reported in this work.
(^{1}) Our application has been designed so as to facilitate changing the fluidsolving application. Incidentally, it has already been coupled to PFEMApplication, a Lagrangian description fluid solver based on the particle finite element method (PFEM) [268,66].
(^{2}) The code can be retrieved from https://github.com/KratosMultiphysics/Kratos (retrieved on June 15, 2018).
The core of this work is contained in Chapters 2 to 4. Of these, Chapter 2 is the most theoretical in character. It is devoted to the analysis of the Maxey–Riley equation, as a fundamental model for the description of the motion of individual particles in a fluid, which at this point is assumed to be described by a known field. Section 2.2 systematically explores the range of validity of the model with respect to several criteria, first in terms of the nondimensional values that appear in the equation itself and later in terms of additional variables involving a selection of simplifications introduced a priori in the development of the theory. In Section 2.3 we apply a scaling analysis to the different terms in the model, providing estimates of their relative magnitude that may be applied in practice to simplify the basic model by neglecting the less important terms. Section 2.4 contains a summary of the main results of the chapter, including Tables 2, 4 and 5, that list the most important numerical estimates.
Chapter 3 shifts from the analytic point of view of Chapter 2 to the study of different numerical techniques to solve the Maxey–Riley equation, providing a detailed description of the algorithms involved. We emphasize the treatment of the history force term, whose numerical solution remains challenging, and that is often ignored for this reason. In Section 3.2 we provide a state of the art concerning the numerical treatment of this term, comparing the accuracies of several methods of quadrature. Section 3.3 contains the description of our method of choice for the quadrature, proposed by Hinsberg et al. [353], and of the extensions and improvements we have developed based on it. Section 3.4 contains the description of the overall algorithm for the time integration of the equation of motion including the history term. In Section 3.5 we make a stop to connect the general theory of fractional calculus to the Maxey–Riley equation, which had only been very superficially sketched before. This section appears at this point to take advantage of the terminology and concepts introduced in the previous developments. Section 3.6 is devoted to a systematic analysis of the accuracy and efficiency of the numerical method of solution, applying it to a sequence of benchmarks of increasing complexity. We close the chapter with a summary of the most important results and developments.
Chapter 4 moves even further toward practice, being eminently applied in nature. We consider a series of applications representative to different families of industrial problems, describing several developments involved in their solution as we go along. Section 4.2 discusses how the equation of motion considered in the previous chapters can be modified using empirical relations to extend its range of applicability beyond the limits studied in Chapter 2, for this will be necessary in the industrial applications that are considered later on. In Section 4.3 we introduce the fluid phase as a problem to be solved for the first time, for which we make use of a wellestablished stabilized finite element formulation. We describe it in sufficient detail to clarify the terminology and set the stage for the generalizations introduced in Section 4.3. In Section 4.4 we discuss the problem of derivative recovery, as a step necessary to obtain accurate estimates of the derivatives of the fluid field, once a solution is produced by the fluid solver. We give a brief stateoftheart and compare several alternatives compatible to the finite element method. The description of the particlesfluid coupling is described very briefly in Section 4.5 where we basically bring together the different algorithmic parts involved. We quickly turn to the first application example, which is the subject of Section 4.6. It consists on the numerical simulation of the phenomenon of air bubble trapping in Tshaped pipe junctions, which was only recently first studied by Vigolo et al. [356]. This test exemplifies the use of the oneway coupled strategy, with no interparticle interactions, representative of lowReynolds number internal flows with lowdensity suspensions of small particles. We continue on to the study of the next test application in Section 4.7, where we study particle impact drilling, a technology used in the oil and gas industries to produce high penetration rate drilling systems. Here we again make use of a oneway coupled strategy for the fluidparticles coupling, although this time we consider the interparticle contact as well. This is by far the most detailed of the examples provided and corresponds to consultancy work made during the Doctorat Industrial. Before moving on to the final example, the theory related to the backwardcoupled flow must be introduced first, based on the theory of multicomponent flows (its rudiments can be found in Appendix H) and the general finite element formulation introduced in Section 4.3, which we specialize for this problem. Following this theoretical interlude, we present the last application example, a fluidized bed of GeldardD particles in Section 4.9. We close with a summary of the chapter in Section 4.10.
We end the work with Chapter 5, where we give a schematic summary of our developments. In Section 5.2 we discuss some of our current and future work. There is indeed a considerable number of research lines that have been opened with this work, but that we have not followed as far as we would have liked due to the breadth of our scope.
We provide a series of appendices that carry substantial content. Our intent in doing so has been more to lighten the main text than to expand it with additional content. Perhaps the idea of modularity, associated to the programming principle, has contaminated our writing style, but we hope it helps comprehension. For instance Appendix A contains most of the little text devoted to the discrete element method while Appendix H is devoted to the description of the continuum theory related to the backwardfluid flow. The other appendices include formulas and numerical data of interest mainly to the interested reader willing to program the associated algorithms.
Finally, the Nomenclature section provides a quick guide for the symbols used throughout. Only the symbols that are exclusively used in a single context, close to their definitions, have been left out of the list, to avoid overpopulating it.
The Maxey–Riley Equation (MRE) [238] describes the motion of a small, rigid, spherical particle immersed in an incompressible, Newtonian fluid. It is an expression of Newton's Second Law, relating the acceleration of the particle, taken as a pointmass, to the sum of the external forces (usually just its weight) plus the total force applied on it by its containing fluid. The smallness of the particle justifies the two fundamental assumptions made in its derivation: 1) that the relative flow can be described by the Stokes equations near the particle; and 2) that the flow is well represented by its secondorder Taylor expansion about the particle's center (near the particle). The equation is thus applicable to the study of disperse particulate flows, for systems where the low volume fractions justify the consideration of each particle as an isolated body. Examples of such systems can be found in the study of warmcloud rain initiation [122,314,363], turbulent dispersion of microscopic organisms [275,365], the dynamics of suspended sediments [324,155] and also in more fundamental areas, such as the study of turbulent dispersion [324,27,53].
The hypotheses involved in the derivation of the MRE are very often too stringent for the applications of interest. For such cases, there exist a huge variety of models applicable outside its strict range of applicability. These models, which usually involve a high degree of empiricism, are nonetheless very often constructed on the basis of the MRE (e.g.they are asymptotically convergent to it), showing similar mathematical properties and physical significance. Examples of such extensions apply to higher particle Reynolds numbers, nonspherical shape or finitesized particles; see Section 4.2. The analysis and numerical treatment of the MRE is therefore of major importance from both the applied and the fundamental perspectives. This chapter is devoted to the study of several aspects related to the MRE, more specifically to its range of applicability and scaling properties. Our goal is to provide a solid ground on which to base engineering decisions concerning its use.
Let us recall the exact formulation of the MRE. Let define the background fluid velocity field, where is some spatial open domain contained in the fluid and is the ambient spatial dimension (either two or three). Let be the trajectory of a particle through this fluid, and its velocity. The MRE states:

(2.1) 
where is the mass of the particle, is the mass of the displaced fluid volume, is the particle radius, is the density of the fluid, is the dynamic viscosity of the fluid ( is the kinematic viscosity) and the acceleration due to gravity. The notation refers to the material derivative of the fluid. Eq.~2.1, together with and the initial conditions and form an initial value problem that must be solved to obtain the trajectory of the particle.
The terms on the right hand side of Eq.~2.1 have distinct physical interpretations and can be identified as: the force applied to the volume displaced by the particle in the undisturbed flow (), the added mass or virtual mass term (), the Stokes drag term (), the Boussinesq–Basset history term () and the term due to the weight of the particle minus its (hydrostatic) buoyancy ().
The terms bearing second order derivatives are the second order corrections, known as Faxén corrections, that take into account the finite size of the particle with respect to the secondorder spatial variations of the flow field around it, see Section 2.2.2.
The exact formulation of the MRE as was originally given by Maxey [238] differs slightly from the version presented here, which adopts the Auton et al. [14] form of the added mass force. This formulation contains the material acceleration of the fluid continuum (i.e.the derivative following the fluid), evaluated at the current position of the particle. The original formulation presented instead the total derivative (i.e.the rate of change of the fluid velocity when following the particle). The former form, derived under the assumption of inviscid flow in [14] is however valid to the same order as the other terms derived by Maxey and Riley within their assumptions, while remaining accurate at higher Reynolds numbers [361].
Another difference with respect to the original formulation involves the Boussinesq–Basset term, on which the time derivative appears outside the integral sign in Eq.~2.1, as opposed to inside, as in the original Maxey and Riley paper. As pointed out in [89] this is a more general formulation, which allows for a nonzero initial relative velocity (and is also equivalent to the generalization later given by Maxey [240], see [222]; see also Appendix C).
Let us rewrite the MRE in dimensionless form. This will serve to introduce some fundamental dimensionless numbers. For the sake of notational simplicity, the same symbols have been overloaded to designate the dependent and independent variables, although here they represent nondimensional quantities.

(2.2) 
where

(2.3) 
where and , and are the dimensional scalars by which (and ), and have been normalized. These quantities are conventionally defined such that is a characteristic length scale of the flow and is defined such that is the characteristic magnitude of the gradient of the unperturbed velocity near the particle location, while is defined such that . By the term unperturbedwe refer to the flow resulting from subtracting the (Stokes) flow produced by the presence of the particle under the same farfield boundary conditions from the actual flow. The application of Buckingham's Pi Theorem to this equation reveals that the set of dimensionless parameters that describe it is in fact minimal. By defining , we obtain the relation

(2.4) 
The quantity is commonly referred to as the particle's response time [216]. It is equal to the time it takes for a static particle to accelerate to (where ) parts of the surrounding fluid velocity, under the action of a constant, uniform flow (and neglecting ). The Stokes number can thus be interpreted as the ratio of the particle's characteristic response time to the fluid's characteristic time, and it measures the dynamical importance of the particle's inertia. It is useful to additionally define

(2.5) 
where is the characteristic magnitude of the relative (or slip) velocity . The quantity is known as the particle's Reynolds number and it characterizes the importance of inertial effects versus viscous effects in the flow produced by the relative motion between the particle and the background flow.
Eq.~2.1 is derived from the full Navier–Stokes problem as an asymptotic relation, valid as [238]

(2.6) 
In practice, these asymptotic relations are interpreted as requiring the lefthand sides to be much smaller than one ^{1}, that is

Note that the LHS of Eq.~2.7 can be understood as a sort of downscaled version of the ambient to the size of the particle, since the velocity variation seen by the particle scales as . This condition, along with Eq.~2.8 can be used to simplify the Navier–Stokes equations to the Stokes equation in the derivation of the MRE by resorting to scaling analysis of the full Navier–Stokes equations [238]. The approximation affects all the forces except , which is treated independently and only relies on Eq.~2.9 (and a sufficient degree of smoothness of the flow).
The condition Eq.~2.9 is necessary to ensure that the disturbance flow caused by the presence of the sphere is well approximated at the surface by its second order Taylor expansion around its center. This approximation remains good as long as the mean flow around the particle is approximately quadratic, and the error is expected to drop very fast with decreasing (as its third power). Note that, in practice, an accurate representation of the fluid velocity field will require a sufficient resolution of whichever numerical method is employed in its computation. The calculation of the spatial derivatives requires special attention. This issue will be addressed in Chapter 4.
Our main goal in this section is to gain a more complete understanding of the limits of applicability of the MRE than that directly offered by Eqs.~2.72.9. The original motivation for such an endeavour was to provide a bounded parametric space for the subsequent discussions about the relative importance of the different terms in the MRE, so as to make the analysis simpler and more precise. But we have come to recognize that such analysis should also be useful in itself, as a reference for the engineer or investigator concerned about the suitability of the equation for their particular problems. We are thus interested in answering questions such as
The survey we present is inevitably incomplete, though this kind of analysis is best realized through an iterative, cumulative process to which we hope to contribute a part. There have certainly been some notable efforts in this direction, the most important of which is the systematic analysis that E. Loth has undertaken over the last fifteen years [216,219,217,218,221]. However, his interest was mainly focussed outside the strict range of applicability of the MRE, exploring the ways in which it should be extended rather than its limits of applicability. Other relevant works are the best practice guideline [334] and the book [85].
Let us start by deriving some bounds to the dimensionless parameters in Eq.~2.3. First, the relative density can be extremely large for gas flows, so that in practice it can be considered as any positive number. For liquid flows, however, it becomes bounded of order 1x10^{1} for any commonly encountered materials (for a suspension of osmium particles in gasoline one has , most likely an upper value for ordinary applications). Note that for atmospheric air instead of water this value would be about a thousand times greater.
At the other extreme, the relative density of air bubbles in water can be practically taken to be zero with a minimal change in the value of (typically, small bubbles and low Reynolds numbers lead to spherical bubbles; see [47] for an indepth discussion on the circumstances under which bubbles can be modelled as rigid and spherical).
Moving on to the Stokes number, we find that it can also be bounded above by relating it to the particle Reynolds number as

(2.10) 
where we have interpreted Eq.~2.9 as and used the fact that for very large the relative velocities will be of the order or larger than the absolute velocities of the smallest scales of the flow. We will come back to this issue in Section 2.3. Using Eq.~2.10 and the above argument for liquids allows us to bound by the particle Reynolds number and in many cases (e.g.water suspended mineral particles) by an additional order of magnitude, so that if , then , speaking in orderofmagnitude terms. In atmospheric air, a similar argument would lead to bound below two to three orders of magnitude above . Using the more conservative would lead to lower upper bounds. We summarize the situation in Table 1, that includes some estimates for the upper bounds corresponding to different representative material combinations.
Continuous phase  Dispersed phase  
air  water  sand  copper  
air  1  10  50  100 
water  1  1  1  5 
Apart from Eqs.~2.72.9 and Eq.~2.10, there are additional assumptions involved in the derivation of Eq.~2.1 that we want to visit. To begin with, the problem is posed for a single sphere in an infinite domain. That is, the presence of nearby particles and of solid boundaries is not contemplated, though these elements are most often present in the applications. Mere qualitative restrictions are of no practical use to the engineer, other than as a reminder of them being a source of error. Thus, there is a strong need to provide predictive quantitative measures too. Unfortunately, to the best of our knowledge, these are still open issues, and not much more than a few rough rules of thumb can be found in the literature.
It is also of interest to assess the range of validity of the different terms in the MRE, one at a time. The reason is that each term describes a distinct physical effect, with its own characteristic response to variations in the dimensionless parameters that characterize the flow. Indeed, in a given situation it might be of outermost importance to correctly calculate the steady drag force and not, say, the added mass force. Take for example the study of the deposition rate of microscopic particles in a container. Due to their tiny inertia the particles reach their terminal velocity very quickly. After the short transient the added mass force cancels, having a negligible impact on the much larger deposition time. In such cases, it is common practice to simply neglect the added mass force altogether, making to all effects the range of applicability of its particular formulation irrelevant. This sort of argument has allowed to justify the application of the MRE to a remarkably wide range of situations, that has been stretched even further by tweaking specific terms. A paradigmatic example of this practice is the use of empirical drag force laws combined with the usual formulation for the other terms (and usually neglecting the history force). We will revisit this type of approaches in Chapter 4, where we will need to move outside the range of applicability of the MRE.
It is important to realize that the use of this kind of empirical extensions almost invariably relies on the validity of the same additive structure of the different forces present in the MRE. Encouragingly, there has been an important body of research supporting such additive division both from a theoretical [223,14] and an empirical [222] points of view in a variety of situations and well beyond the range of applicability of the MRE.
All these questions are the subject of the following sections. We will look at the different directions in which the hypotheses involved in the derivation of the MRE can fail, in a systematic way. Since we are interested in the first effects, we will assume that it is adequate to look at each term independently. Furthermore, the additive nature of the different forces suggests that each one can be treated independently. We will therefore do so when possible, making an effort to set specific numerical bounds to the applicability ranges.
(^{1}) The imprecision in this requirement must be supplemented by experience or by some conventional rule adapted to the particular fields of application. A common criterion is to interpret as meaning at least two orders of magnitude smaller than one or [197].
The well known Stokes solution of the steady, low Reynolds number flow past a sphere (see Eq.~2.31) is obtained by completely neglecting fluid inertia. By applying the noslip boundary conditions on the particle surface and the farfield velocity conditions at infinity and expanding the stream function in powers of , the hydrodynamic force on the particle can be calculated, to leading order, yielding the drag force of the MRE without the Faxén terms (this is the problem solved by Stokes in 1851). However, if one wishes to increase the order of this approximation following the same procedure, one soon realizes that there is no way to fulfil the far field boundary conditions in this case, since the higher order contributions to the perturbation caused by the particle do not vanish at infinity. This phenomenon is known as Whitehead's paradox [248].
Its resolution was possible due to ideas from Oseen, who noted that the assumption by which inertial terms are disregarded (under Eq.~2.8) is only valid near the particle, where viscous effects dominate. But, far from the particle (in particular at a distance such that , see [284]), the approximation breaks down. This implies that there is an inherent inconsistency in requiring that the Stokes solution be valid in an infinite domain, and that it is necessary to consider inertia far from the particle in order to calculate the higher order corrections to the drag force. The final word on the issue was nonetheless given by Proudman and Pearson [284], who, using the technique of matched asymptotic expansions, rederived Oseen's inertial, firstorder correction to the steady drag force, corrected a flaw in Oseen's original reasoning and added an additional term to the expansion:

(2.11) 
This way, the contribution from the outer region becomes a correction of the steady drag force, whose leading term, given by , could be used to provide the order of magnitude of this correction.
In order to illustrate how it could be used in practice, let us assume that we wish to establish the upper limit to the particle Reynolds number to ensure that the correction is smaller than 5 as a convention to fix the range of validity of the MRE. Then the leading term in Eq.~2.11 predicts this to happen for , and by this point the error in this expression is only , taking as a reference the empirical drag coefficient by Clift et al. (1978) (which has a root mean square error of with respect to the empirical data gathered by Brown and Lawler [50] covering the range ). This means that its range of validity is large enough to estimate the firsteffects of inertia for error tolerances lower than 5%. In Fig. 2 all these different approximations to the drag force are compared.Figure 2: Drag force predicted by different approximations, as measured by the drag coefficient . 
The above expression of the force is only valid for longterm steady motions. Lovalenti and Brady [223] studied the motion of a small sphere immerse in an arbitrary space and time varying flow, although considering the particle small enough, so that the undisturbed flow can be considered uniform over the diameter of the particle. The problem is thus the generalization of the MRE (excluding Faxén terms, as these authors assume flow uniformity; see Section 2.2.1) for small, though nonzero particle Reynolds number. Their analysis is based on a different approach compared to the asymptotic expansion matching used by Proudman and Pearson [284]). In fact, they find an expression for a uniformly valid flow over the whole domain to first order in . The hydrodynamic force on the particle is then computed to first order in by approximating the resulting integrals of the stress tensor. The result is an expression for the steady and unsteady forces at finite Reynolds number. The inertial effects include the steady Oseen correction to the drag, as well as corrections to the unsteady forces like the Basset Force and the added mass force. These corrections confirm the observations from fullyresolved numerical simulations [243] that higher rate of convergence of the particle velocity to terminal conditions (i.e.it varies but is for smoothly accelerating motions where the slip velocity tends to increase) than that predicted by the MRE equation, which is , as predicted by the Boussinesq–Basset kernel. According to Lovalenti and Brady [223]: This fact may explain why experimentalists have measured a steady terminal velocity when the length of their apparatus would not have permitted this if the Basset force was correct..
This change in the decay regime of unsteady forces only applies at long times. For shorttimescale motions the unsteady force from the MRE is accurate. The boundary between both regimes is marked by the time it takes for vorticity to diffuse away from the particle up to the socalled Oseen distance, where inertial and viscous contributions become of the same order. At that point, convection transport takes over as the dominant physical mechanism for the transport of vorticity, explaining the change in regime. From the application point of view it is important to realize whether the flow is dominated by high frequency flows, so that the motion is mainly explained by the MRE unsteady forces, or else the characteristic frequencies of the flow are low and thus the (substantially more complicated) expressions provided by Lovalenty and Brady apply. The key is to compare those characteristic frequencies of the particle's motion with the quantity , the viscous diffusion time, which corresponds to the order of magnitude of the time it takes vorticity to be diffused away to the Oseen length. In other words, the MRE equation applies whenever [224]

(2.12) 
After a little algebra, this condition can be rewritten as a function of the usual nondimensional parameters as

(2.13) 
As the value of increases, the condition above may cease to be fulfilled and the form of the history term must be corrected using, for example, the formulation given by Lovalenti and Brady [223]; see also Lovalenti and Brady [224]. In turbulent dispersion, the importance of this correction will not be critical for most flows such that , except perhaps for solid particles in gas. The same conclusion was reached by [211] based on the empirical formulation for the history force of Mei and Adrian [243].
The effect of vorticity
As stated, the derivations in Lovalenti and Brady [223] assume a uniform flow at the scale of the particle. If one allows for some nonuniformity to be present in the flow, additional effects appear. In particular, this allows for a break in the symmetry of the flow on the plane orthogonal to the direction of , which gives rise to lateral forces on the particle that are known as lift forces. Famously, Saffman [299] (with the corrections in [300]) gave the correct expression, to first order in , for the lift force that an isolated, nonrotating particle experiments under constant shear (linear nonuniformity of the flow).

(2.14) 
where is the slip velocity. This result is valid under the restrictions , and ; where the shear Reynolds number is defined by

(2.15) 
that is, the particle Reynolds number multiplied by a dimensionless shear rate gradient.
Unfortunately, the typical values for in turbulence are often smaller than those of [364], invalidating the application of the formulation. However, in this case the expression is in fact overestimating the lift force [364] and so, without loss of generality, we may apply the formulation as an upper bound for the remaining discussion.
But the fluid vorticity can arise due to the solid rotation of the fluid as well, also giving rise to a lift. Its low analytical formula was provided by Herron et al. [164] and reads

(2.16) 
which has a coefficient about 2x10^{1} larger compared to . This expression is subject to the restriction . Note that and are incompatible and are only strictly valid for their corresponding ideal fluid motions. In fact, Candelier and Angilella [56] have proved analytically that for a particle settling in a solidbody rotation fluid, the lift force can even take the opposite sign as that indicated by Eq.~2.14, which is valid only for pure shear, steady flows. In this case, furthermore, the relative motion is not stationary, due to the migration of the particle in the radial directions.
Indeed, as with rectilinear motion, the relative flow unsteadiness also generates historydependent contributions. Those have been studied in [60,55] for pure shear (generalizing Saffman's result) and solidbody rotation motions respectively. Nonetheless, the order of magnitude of the Saffmann lift force, compared to the steady drag force gives an idea of the magnitude of this correction. An analogue estimate can be obtained if fluid motion is closer to that assumed in Eq.~2.16.

(2.17) 
We will thus consider the lift force to be vary small within the range of applicability of the MRE. Nonetheless, one must keep in mind that the correction is orthogonal to the drag force and thus its effect can introduce important systematic changes in the particles' trajectories, especially in flows with a predominant direction. Since the largest values of the shear rate are usually close to the boundaries of the domain, this is very often the case anyway, for instance in internal flows. When considering the flow of particles along ducts, it is therefore important to consider the possibility that the inertial lift force might play a role, except for extremely small particles. In any case, we leave the study of the effect of boundaries for future work, see Section 2.2.6, and thus we will not elaborate this point any further here.
The effect of particle rotation
The linearity of Stokes flow means that, for a spherical particle, rotation and translation are decoupled. Indeed, the combined motion is the result of adding the effects of one and the other, and it is clear that, by symmetry, the relative rotation of the particle alone cannot produce a force in any particular direction. It is however necessary to examine how fast this rotation can be until the first inertial effects appear, invalidating the MRE.
The firsteffects of the force modification due to nonzero particle angular velocity were calculated by Rubinow and Keller [294], to zerothorder accuracy in . This is enough for our arguments, since the MRE assumes this effect to be very small. This formulation predicts that due to the rotation of the particle, there arises a lift force given by

(2.18) 
where is the angular velocity of the particle. In other words: the Reynolds number defined by the slip velocity produced by rotation should be small enough. The lift force that arises due to the rotation of the particle relative to the fluid is generically called Magnus effect [219].
The condition of having a negligible contribution to the force could then be expressed as a function of the particle Reynolds number associated to the rotational motion, given by

(2.19) 
where . Note that , and so this is assumed to be very small within the theory of the MRE.
Maxey [238] argues that in turbulent flows the particle will tend to acquire a rotational velocity of the same order as the local shear rate. Taking into account that the firstorder rotationinduced effect of inertia is a lift force with a magnitude of order [294], the effect results of order ; more precisely, we have

(2.20) 
where give the estimation of local shear rate; and being the characteristic velocity and distance of the flow, which is certainly a small quantity by Eq.~2.7.
There are circumstances in which the particle angular velocity might be substantially higher. For instance, a collision against a wall or another particle might transform a good part of the translational energy into rotational energies, leading to high angular velocities, especially for the smaller particles ^{1}. The fact that this force is perpendicular to the direction of translation makes this effect even more important. Nevertheless, collisions are not expected to be too frequent in this regime (see Section 2.2.4), substantially alleviating this effect.
In this chapter we ignore the rotational degrees of freedom of the particle, except for the present discussion, with the argument that these effects are, to first order of approximation, decoupled. Nonetheless, it is worth to consider momentarily what the rotational dynamics of a particle within the range of applicability of the Maxey–Riley regime looks like. Feuillebois and Lasek [128] derived the expression of the instationary motion of a small, rigid sphere spinning in a viscous, Newtonian fluid, providing historydependent terms that decay much faster than the translational analogues. In the steadystate limit, the equation reads

(2.21) 
This equation allows us to calculate the rotational relaxation time as

(2.22) 
from which a rotational Stokes number can be defined (assuming the steadystate forces are dominant)

(2.23) 
which shows that this Stokes number is of the same order as the translational Stokes number, under the assumption that the fluid time scales are similar. Note that for large Stokes numbers the assumption above that the particles rotation will be of the order of the fluid vorticity most of the time (in fact equal to the local angular velocity or one half of the vorticity in the limit of zero inertia) can be violated, as the particles may not have time to accommodate to the fluid vorticity. In such cases the analysis becomes more involved and it is likely that finite Reynolds number effects must be taken into account for the angular dynamics. Nonetheless, for this orderofmagnitude analysis the present formulation will suffice, especially since the rate at which the steadystate formulation becomes inaccurate is quite slow, and it still yields reasonable results at well past unity, see [219].
(^{1}) Indeed, a fixed percentage of energy conversion from translation to rotation leads to higher angular velocities for smaller particles. Specifically, the rotational kinetic energy is and the translational kinetic energy is . Thus, after a collision we have , for light and/or small particles. Furthermore, this gives ; and the ratio can become quite large.
As the particle radius grows with respect to the characteristic length scale of the unperturbed flow the Faxén corrections become increasingly important. Eventually, even these secondorder corrections become insufficient to accurately characterize the surrounding flow and the MRE breaks down as an appropriate model. The rate at which this happens can be analysed by taking one more term in the series expansion that leads to the appearance of the Faxén corrections [222] (see also Eq. 8.175 in [91]). It is useful to consider the expression of the drag force as a surface integral to measure the error. The Faxén corrected drag force becomes, after taking one more term in the same expansion considered in [222]

(2.24) 
Where the is the th successive composition of the vector Laplacian.
Now, for low , the fluid is expected to present only longwave length variations across the particle diameter. Furthermore, since our interest lies on the integrated value of such variations, it appears reasonable to assume that the characteristic length scale of the new term is also the particle radius. This allows us to compare its order of magnitude to that of the usual Faxén correction. Let us therefore assume that the effect of the new term is of the same order as that produced by quadratic variation of across the particle diameter. Then, an adequate orderofmagnitude estimate of its value is given by (see, e.g. [245] for a discussion on the estimation of order of magnitude scales)

(2.25) 
So that the next order correction will be significantly smaller than the classical Faxén correction up to the point when is not small any more (in the limit of quasisteady Stokes flow it holds exactly). This means that in practice the small requirement will guarantee that the Faxén corrections will make a good enough job and thus the restriction Eq.~2.9, that is, by itself does not appear to be restrictive, but rather indirectly through the violation of . A similar argument can be made with the other forces bearing their corresponding Faxén correction. A significant study supporting this conclusion in the context of isotropic turbulence was provided by Homann and Bec [171], who tested the performance of the Faxén terms for a neutrally buoyant particle in direct numerical simulation (DNS) of isotropic turbulence. Their conclusion was that the first effects of the finite size of the particles were well captured by these terms up to (where is the Kolmogorov microscale). From then on, inertial (finite ) effects quickly kick in.
No particle is perfectly spherical in reality and thus it is necessary to acknowledge this fact in assessing the reliability of the MRE for predicting the trajectories of real particles. Departure from the ideal spherical shape greatly complicates rigorous analysis and, most often, precludes it completely. The situation is somewhat similar to the consideration of nonsphericity in granular flows using DEM. It is known that its effects are important although only very rough approximations to them can be captured when using only spheres. And while the possibility of using nonspherical shapes also exists, it is greatly limited in practice due to increased numerical costs and difficulties in achieving a good characterization of realistic shapes.
Nevertheless, there has been substantial theoretical progress in characterizing the hydrodynamic interaction of particles of nonspherical shape at low ; see [219] and references therein, [114,59,58]. Gavze [142] derived a generalized version of Eq.~2.1 for a body of arbitrary shape written in terms of its viscous resistance tensors [48]. Within this model translation and rotation become coupled.
Of particular relevance to our arguments is the study of Zhang and Stone [389], who derived the first order asymptotic correction to the unsteady equation of motion of a sphere in quiescent fluid in an amplitude parameter, , of the deviation from an ideal spherical shape using the reciprocal theorem (see, e.g. [152]). Their equations particularize the work in [142] to the limit of weak nonsphericity, in which the coupling between rotation and translation comes in as a higherorder correction and can thus be neglected to firstorder approximation.
Specifically, let the particle's surface be defined in polar coordinates around its center as the set of points with coordinates , and such that

(2.26) 
where is the radius of the equivalentvolume sphere.
Their equation for the hydrodynamic force on the translating particle reads:

(2.27) 
where

(2.28) 
and the integral is over the unit sphere , after applying a uniform scaling by its radius over the whole space. Note that represent the exterior unit normals to the unit sphere. The surface of the particle is by construction (although this is not emphasized in [389]) assumed to be starshaped with respect to its center of mass^{1}. We can then use these results as a conservative estimate of the first perturbations due to nonsphericity to the MRE.
Eq.~2.27 allows us to estimate an upper bound to the magnitude of the corrections to the MRE that are needed to take into account generic (weak) nonsphericity. Let us consider the induced matrix norm on matrices denoted by , based on the usual Euclidean norm for vectors, denoted with the same symbol. Then, for example, the corrected drag force fulfils

(2.29) 
Similar estimates can be worked out for the other hydrodynamic effects, yielding similar relative orders in their corrections. This bound is related, although a bit less sharp than the one provided by the application of a theorem by Hill and Power [166], which establishes that the magnitude of the steady drag (it only refers to the steady drag force) on the particle lies between that obtained on an inscribed and a circumscribed spherical particles. That is, the center of the circumscribed and inscribed spheres do not have to be the center of gravity. In this case therefore, for steadystate motions, the relative error made by taking the Stokes expression with the average diameter is smaller than

(2.30) 
In the limit of small corrections to the spherical shape, however, all inscribed and circumscribed spheres tend to collapse to the ones around the center of mass and the result then coincides, for the drag force, with the one provided by the asymptotic technique by Zhang and Stone [389].
(^{1}) In geometry, a body is said to be starshaped if there exists a point within it such that, for any other point in the body, the straight segment that joins both points is itself fully contained in the body. This appears to be a reasonably weak assumption for approximately spherical particles. In particular, it is much less restrictive than convexity. See [156] for an analytic correction to the added mass force for a (nonstar shaped) rough sphere.
The MRE is based on the assumption that the particle is isolated, far away from any neighbouring particles or any other boundary. But, for many applications, such assumption can become too restrictive. Even in the absence of interactions such as electrical potentials or van der Waals forces between the particles, there are still two types of interparticle influences:
Starting with the longrange hydrodynamic interactions, it is clear that they depend on the configuration of neighbours and their velocities, requiring a statistical approach for a general analysis. Under the lowReynolds number hypothesis, such influences are linear on the particles' separations (and velocities, including angular velocities), meaning that their spatial decay is slow. Every particle in a suspension produces a disturbance flow around it which, under the hypothesis of a low particle Reynolds number, is accurately represented by the creeping flow equations solution, at least up to a distance comparable to the Oseen distance () [248] from the particle's center, i.e.

(2.31) 
where the velocity components are Cartesian (with the first component being aligned with the relative flow) but have been parametrized in spherical coordinates, with the radial coordinate (normalized by the sphere's radius) and the azimuthal angle. Note that, by symmetry, the field must be independent of the polar angle, as it is. Eq.~2.31 confirms that the decay of the disturbance intensity is dominated by a linear behaviour, which means that its influence can still be appreciable at significant distances as shown in Fig. 3. In Fig. 4 the average distance to the nearest th neighbour, as a function of , is shown for reference; for in a random array of spheres. The calculation of this distance is done according to the formulation in [340], that takes into account the space taken up by the finitesize particles and is robust for all values of up to close packing.
Figure 4: Average distance (normalized by the radius of the sphere) to th closest neighbour in a random array of monodisperse spheres. 
A description of the flow around a given particle, under the influence of its neighbours, would require the solution of the Stokes equations around all of them simultaneously. And since the solution depends on the velocities of the particles, which are in turn instantaneously determined by the fluid velocity around them, all interactions are completely coupled, making the problem fully multibody in general. In the dilute limit that we are interested in (we are studying the firsteffects of the presence of neighbours), it is, however, possible to give the interactions a pairwise treatment. This means that they can be calculated as the sum of the interaction between the individual pairs of neighbours with only a small error [23]. Still, the main difficulty lies in describing a representative configuration of relative positions and velocities for a representative ensemble of particles, valid for the general case.
Switching to the shortrange interactions, these include two types of forces: the shortestrange hydrodynamic interactions, termed lubrication effects, and the direct contact interactions, that appear once the fluid has been completely squeezed out of the gap between pairs of particles coming into contact. Here too, a statistical analysis seems unavoidable. Fortunately, these effects are more adaptable to the general DEM approach, as these are by definition pairwise, shortrange forces that can be naturally included in a contact model. Nonetheless, in some situations it may still be interesting to neglect them completely as:
We wish to characterize the conditions for which these effects start to become important, causing the MRE to start loosing accuracy in modelling particleladen flows. Such characterization has been attempted before [116,118] in the context of general particleladen flows. Specifically, Elghobashi [116] establishes limits on the global solid fraction (the proportion of volume occupied by the particles in the whole domain; see Section 4.8 and Appendix H) for which threeway interactions become important is

(2.32) 
where is the local solid fraction. This corresponds to an average nearestneighbour separation of about 1x10^{1} particle radii. This is an oftcited limit [232,159,206]
We will attempt to enrich the general picture here, surveying the most relevant results and the different aspects to take into account with the ultimate goal of producing a useful guide for practical applications. Once again, we focus on turbulent flows to provide the orderofmagnitude estimates.
Before proceeding, we must make one further remark concerning the scope of the discussion. In some particleladen systems, the particles may interact forming aggregates, coalescing or breaking. We do not consider these cases here, but rather the case in which the total number of particles and their shape is preserved throughout the simulation.
Hindered settling
In a random arrangement of particles, the added effect of all particles around any one of them can be very substantial, due to the longrange character of the interactions that we discussed above. In fact, the naive summation of the moduli of the disturbances would not work at all for our orderofmagnitude purposes, because of the slow decay of the interactions with distance. Indeed, take the sum of the norms of all the pairwise interactions between the particles inside a ball centred at the target particle and the particle itself, which bounds the contribution of the particles considered above. By considering a sequence of balls of growing radii, one obtains a sequence of sums that, due to the linear decay of the forces with distance, does not lead to an absolutely convergent series, because the number of terms grows with the third power of the radius in a statistically homogeneous distribution. Thus, such straightforward strategy fails to produce any bound at all.
This difficulty was overcome by Batchelor [23], who was able to rigorously calculate the (ensemble) average, to firstorder in the solid fraction, of the effect of a stable, random uniformly distributed suspension of identical rigid spheres on the settling velocity of any given particle. This author used an ingenious method in which known averaged quantities were used to cancel out the slow varying contributions exactly, reducing the remaining terms to rapidly varying quantities that could be calculated, to first order in the concentration, based on the interaction of the two closest particles only. The result of his analysis can be written as

(2.33) 
where is the terminal velocity of the target particle when isolated, the actual settling velocity of the particle in the suspension, is the solid volume fraction and is the settling coefficient, which under the particular assumptions of the theorem is ; and . This result is based on the further assumption that the probability distribution of the distance between pairs of particles is uniform. A similar expression has been derived for polydisperse suspensions [24,93], yielding smaller values for . For instance, for both small spheres surrounded by a suspension of larger spheres or the opposite situation (See [100]). Therefore, one may use Eq.~2.33 as a conservative order of magnitude of the first effect of the presence of neighbouring particles in sedimentation of statistically homogeneous suspensions. Note that in this case ( is negative) the presence of other particles tends to slow down the target particle as it sediments. This is true specifically because of the assumption of having a uniform suspension that can be treated as unbounded, since nonuniformity often leads to a reduced resistance and thus an increased settling velocity with respect to an isolated particle, as will be shown below.
Drag modification
Another relevant theoretical result can be taken as a reference to study the effects of the presence of a random array of particles around a target particle for the calculation of the drag force, which is enough to assess firstorder effects due to the surrounding suspension. The theory is based in the asymptotic analysis of the hydrodynamic forces in the low Reynolds number and solid fraction limit, which is relevant to the present discussion. It was put forward by Kaneda [189] and experimentally verified up to solid fractions of around 4x10^{2} for small Reynolds numbers based on Brinkman's screening length, , as the simulations for higher values of were found to be prohibitively expensive, in [167].
The derivations by Kaneda were performed for simplicity under the hypothesis , although the resulting expression is shown to have the correct asymptotic behaviour at the opposite limits and . In this work it is argued that the first effects of inertia (see Section 2.2.1) and those of the nonvanishing solid fraction cannot be treated independently as both are inextricably linked for all but essentially zero Reynolds numbers. The expression for the modified drag is

(2.34) 
with

(2.35) 
where ^{1} is a nondimensional constant such that (when ) and where is the usual Stokes drag; and

(2.36) 
Note that

(2.37) 
Which, as stated, recover the first order correction due to Oseen for the drag force at (and thus ), see Eq.~2.11, and also the classic expression by Brinkman [49], at .
Figure 5 shows the magnitude of the correction coefficient for a wide range of solid volume fractions, for different values of . Note that the correction to the drag is significant, even for very small Reynolds numbers if the solid volume fraction is greater than 1x10^{3}, in accordance with the rule of thumb of Eq.~2.32, and it still appreciable past the 1x10^{4} mark.Although Kaneda's result applies to arrays of fixed spheres under a uniform background flow, it clearly reflects the inextricable relationship between the lowReynolds number and lowconcentration hypotheses through the dimensionless number . This number can be understood as the quotient between two length scales:

(2.38) 
where is the Oseen length (distance from the particle at which inertia becomes as important as viscosity) and (distance from the particle at which the Stokesian disturbance caused by it switches from a linear to a cubicorder decay due to the presence of neighbours). The physical interpretation of the two regimes in Eq.~2.36 goes as follows: when is very small (), inertial effects enter much before the Stokesian disturbance has been screened by other particles and thus their influence is dictated by finite effects. On the contrary, if , the influence of distant particles is screened before it can reach the target particle, and only the effect of Stokesian flow disturbances is felt.
Finally, note that the concentration dependence of the drag force is proportional to the square root of the solid fraction, to first order, under the assumption of having a small Reynolds number based on the system size. This contrasts with the result by Batchelor discussed in Section 2.2.4, under which the resistance depends linearly on the solid fraction under the same hypothesis. The reason is the hypothesis of a fixed bed used by Kaneda, which constraints the possibilities that the particles have to accommodate to the resistance, causing it to be more important. The hypothesis of a free bed adopted by Batchelor is more relevant to the behaviour of a cloud of particles suspended in a fluid [92]. Nonetheless, the description of the two regimes studied by Kaneda does remind us that the assumption of having a small Reynolds number used by Batchelor is actually tied to the condition of having a large enough solid fraction to make sure that , a condition that was not explicitly stated by this author.
Recently, Pignatel et al. [278] have experimentally studied the transition between regimes of a falling (finite) cloud of particles. In this case, the Oseen length scale is compared to the clouds initial diameter. Furthermore, a Reynolds number associated with the cloud's length scale is also considered. Several regimes are identified, but what is most interesting to the present discussion is their identification of the Stokes regime (Batchelor's hypothesis) for , . We will come back to considering sedimenting clouds in the proceeding paragraphs, which address inhomogeneous suspensions.
Clustering
In reality, the assumption of homogeneity in the spatial distribution of a suspension is rarely fulfilled. Instead, it is often the case that the particles approach each other forming more or less defined clumps, a phenomenon known as clustering. This is observed very prominently in turbulence, where inertial particles tend to concentrate in filamentlike structures causing the solid fraction to fluctuate strongly within the domain [127,305,350,262].
A number of theories have been proposed to explain these inhomogeneities, and it seems that, rather than a single mechanism, there are several [309,44,157], and their relative importance may depend on several factors, including the particle Reynolds number and, especially, the Stokes number. The most wellknown of these mechanisms consists in the progressive expulsion of heavy particles from highvorticity regions into strainingflow regions [112] as a result of their inertia. This is known as the centrifuge mechanism and it has been well documented both from simulations [239] and experiments [305]. In other words, the inhomogeneities arise because the particles spend more time in highstrain regions than in high vorticity regions, so they tend to concentrate there. Such oversampling of regions of the flow with specific properties is generally known as preferential concentration [323,145]^{2} in the literature.
While the preferential concentration effect is the most predominant at small Stokes numbers [45], other mechanisms become important as grows. For instance, the sweepstick mechanism [149] consists in the tendency of particles to move away from regions of high accelerations and stick to the lowacceleration zones. The sweepstick mechanism becomes important at around and higher (inertial range).
Another effect that becomes important at finite Stokes numbers is the socalled sling effect [122] which develops in intense turbulence. This effect is caused by rare events where very large gradients develop, producing jets of particles that are ejected from the trajectories defined by the streamlines of the vortices from which they come, just as a sling does with a stone. The ejected particles can enter regions where the local suspended particles have significantly different velocities, rendering the particles velocity field (the hypothetical field formed by the trajectories of an infinitude of particles taking up all space) multivalued, a phenomenon known as caustics [367]. These caustics are responsible for an important increase of the collision rate of particles, and is currently believed this plays a crucial role in explaining the fast rain initiation times in turbulent warm clouds [122,123].
Other mechanisms for the appearance of inhomogeneities are turbophoresis [257] and clustering related to boundary layers [234], although we will stick to homogeneous, isotropic turbulence in our discussion for the sake of simplicity. An important and perhaps surprising fact about the clustering effects is that they are known to continue below the smallest scales of the flow, the Kolmogorov microscales in turbulence. Indeed, inertial particles tend to cluster forming multifractal structures [26] at subKolmogorov scales, in a process that is still not fully understood. Furthermore, this phenomenon is not exclusive of turbulent flows, and can even be produced by random, uncorrelated fields that contain a smallest lengthscale, as shown by Bec [25].
Important efforts are being made to characterize the phenomenon and some statistical models, and the development of physical theories to gain insight into clustering phenomena are currently being developed, along with predictive models that are yielding correct qualitative as well as quantitative predictions [382,381,73,131,157]; see [157] for a review. These models contain a number of assumptions and simplifications that limit their range of validity. Recently, Bragg and Collins [44] reviewed some of these formulations, recommending the one by Zaichik and Alipchenkov [381] as the most comprehensive and robust over a wide range of values (see below).
Nonetheless, there remain a number of open questions and inconsistencies in the literature that need to be addressed. For instance, there is a significant consensus that the strength of clustering peaks at [253,178], However Sumbekova et al. [328] concludes that the level of clustering, measured using a Voronoï tessellation of space, is most strongly related to the Taylor Reynolds number, less so to the average particle volume fraction and negligibly on the Stokes number. Moreover, Uhlmann and Chouippe [349] has found that clustering scaling seems to depend on the mode in which it is measured.
Understanding the physical mechanisms that govern the appearance of clustering is of extraordinary importance, not only due to its relevance to the assessment of the validity of the singleparticle theory (which includes the MRE), but also in understanding the statistical distribution of particles in space, statistics on the flow properties being sampled (of relevance to the simulation of chemical reactions, for example) and the prediction of the collision rates (of relevance to the study of initiation of rain and snow, for example). Of direct relevance to the present discussion is the study by Aliseda et al. [4], where the settling velocity of spheres was measured experimentally in a turbulent channel flow. A notable increase in the settling velocity was observed with respect to quiescent fluid conditions, especially for Stokes numbers on the order of unity. Furthermore, this velocity increased monotonically with the overall volume fraction, indicating the effect of threeway coupling. A phenomenological model based on the idea of particle clusters locally altering the average velocity seen by its constituent particles turned out to explain the observations very well. It is remarkable that the whole study was performed under volume fractions well below the limit given in Eq.~2.32, which clearly demonstrates the limitations of the singleparticle paradigm.
The estimates provided in Section 2.2.4 were derived under the assumption of an unbounded array of particles characterized by an average solid fraction under the hypothesis of having a homogeneous suspension. Therefore, the existence of significant inhomogeneities as those caused by preferential concentration in turbulent flows, invalidate the theory and the derived estimates. Aliseda et al. [4] proposed a phenomenological model to explain the effect of preferential concentration in a turbulent suspension of settling particles. They observed that the average settling velocity of the particles belonging to a cluster was reasonably well predicted by Eq.~2.33, with

(2.39) 
where is the volumetric shape factor [285], an orderone constant that depends on the shape of the cluster and is equal to one for spherical clusters and is the characteristic size of the cluster. This formula was derived under the hypothesis of quasisteady Stokesian flow around the cluster. That is, this work is concerned with the lowReynolds regime
Interestingly, this model was theoretically ratified by Jabin and Otto [181] under similar conditions, though no mention is made in their paper of the previous investigations by Aliseda et al. This remarkable work seems to have remained relatively unknown, perhaps because of this omission. This work refines the picture sketched by Aliseda et al. by

(2.40) 
where is an unknown constant and where is the average separation between the particles in the cluster. This expression is equivalent to Eq.~2.39 in the regime (), as it can be immediately shown using and .
Further support to this model can be found by using an argument due to Zaichik and Alipchenkov [381] (ZA henceforth). Note that this is the same work describing the model advocated by Bragg and Collins [44] for the description of clustering in turbulence, mentioned above. ZA used the theory of Batchelor [24], a generalization of that discussed in Section 2.2.4 allowing for nonuniform spatial distribution of particles, to provide an expression for settling coefficient in the much more general setting of isotropic turbulence. The general expression reads:

(2.41) 
where is the radial distance normalized by the particle radius ; , , and are the mobility functions of the particle pair [24] and is the radial distribution function (RDF), defined as the ratio of the probability density of particle pair relative to the same quantity in a uniform suspension, see, e.g., [380].
In simple words, the RDF gives, assuming a particle is found at a certain location, the average number of neighbours found at any given distance from the former, expressed as a proportion to the average number found at that distance in a homogeneous distribution. Note that this is an infinitessimal quantity. The latter can be approximated as

(2.42) 
where the number density is the number of particles per unit volume and is the ball of (nondimensional) radius centred at the particle.
The RDF is often employed in the fundamental study of turbulence dispersion [310,26]. This function would be exactly equal to one for a statistically homogeneous distribution. When a certain degree of inhomogeneity is present, one observes peaks in its value at distances of the order of the characteristic distance of the inhomogeneities (the diameter of the clusters, say). Typically, the RDF asymptotes to one at infinity, as homogeneity is often attained at large enough distances. This is roughly the case for many physical systems, including particle suspensions in turbulence [381] and the molecules in an ideal gas in thermal equilibrium [161] and even the distribution of matter in the universe [347].
The different terms in Eq.~2.41 can be identified as the effect of the particlepair interactions (), the effect of preferential concentration () and the effect of backflow (the resistance caused by the displaced fluid raising due to the settling of the particles cloud).
Therefore, Eq.~2.41 provides a means to calculate the effect of generic inhomogeneities on the settling velocity of the individual particles. Note that it can also be interpreted as the effect on the drag force (since in a timeaveraged sense both notions are onetoone related) of the presence of neighbours in remarkably general conditions, and in a timeaveraged sense. This is of interest to produce a first set of estimates for the effect of neighbouring particles in a wide range of situations involving a predominant drag force and arbitrary distribution of particles. We believe it can provide an orderofmagnitude guide for engineers that wonder about the magnitude of this effect. Nonetheless, the whole theory is subject to the assumptions underlying Batchelor's theory of sedimentation, and thus the restrictions mentioned at the end of Section 2.2.4, defining the Stokes regime for the whole cluster, apply here too.
The key is thus to provide an expression for to calculate the settling coefficient. ZA distinguish between three different regimes, valid for heavy particles ():
Vanishing inertia () with no turbulence In this case there is no clustering and thus , which leads to the same expression derived by Batchelor in [23], where a random, uniform distribution was assumed. In such case, the settling coefficient is given by Eq.~2.33.
Vanishing inertia () In general, in a turbulent flow field the particles interactions produce a drift velocity of the noninertial particles towards one another [382,73]. Taking this drift into account modifies the RDF (see [381]), leading to . That is, there is a small increase of the settling velocity that cannot however compensate the hindering effect completely, making the particles settle slower than isolated, by an amount that is about 4x10^{1} less than that corresponding to a laminar settling situation. In terms of orderofmagnitude arguments this is in any case a negligible correction.
Finite but small inertia () In this case, the RDF is derived from the equations presented in [381] for low inertia particles, based on a powerlaw model for the RDF:

(2.43) 
where is defined so as to guarantee continuity and and depend only on the Stokes number. can be interpreted as a characteristic size of a cluster, since beyond this length the RDF becomes uniform [381]. While this formulation is only strictly valid (within the ZA theory) at distances smaller than the Kolmogorov microlength scale , they use it anyway, warning of its exclusively qualitative value beyond that point. But since we are interested in distances of only up to the order of a few times this distance (the effect at larger distances would be more appropriately dealt with through a twoway coupled approach), and taking into account the excellent agreement of ZA's model with the results reported in [178], fitted from DNS results in the range , we will consider it reasonable for orderofmagnitude arguments.
Moreover, ZA point out that for heavy particles, with (although, judging from their [381], this argument could be extended to ), one has and so using Eq.~2.43 we can write

(2.44) 
which can be rewritten as

(2.45) 
where ( corresponds to in [381])

(2.46) 
The values of and are reasonably approximated in by [381]

(2.47) 
Large inertia () In this case the outweighs by a large margin, since the small disturbances caused by the surrounding neighbours have a relatively small effect on the sluggish, heavy particles. Therefore it is possible to take ^{3}

(2.48) 
where now the radial coordinate () is normalized by the Kolmogorov microscale. ZA provide a set of stochastic, probability conservation equations that can be numerically solved to determine .
The formulation by ZA is reviewed and compared to the alternative formulation by Bragg and Collins [44], who, while not addressing the issue of settling, compare it to the alternative formulation of Chun et al. [73] and give a synthetic review of the most important physical mechanisms that explain the development of inhomogeneities in particulate suspensions in turbulence. Furthermore, these authors propose to introduce a change in the original formulation by ZA to include the nonlocal diffusion effect, as had already been applied by [73]. This mechanism is due to the oversampling of strain vs. rotation regions by inertial particles. Based on an analogy between the two formulations, they were able to apply the same modification to the model by ZA, yielding more accurate results when compared to experimental results for all Stokes numbers, but especially for . For this range of values of the modification is reduced to modifying the value of in the model of Chun et al. by dividing it by the value . And since in this regime the values for for both models are extremely close, we can apply the same multiplicative factor to the value obtained from Eq.~2.47, obtaining a corrected approximation for small Stokes numbers. Doing so results in a much more perfect coincidence between the ZA theory and numerical results for at small values of (see Fig. 6 to appreciate the corrected level of matching).
Still, for larger values of the full ZA formulation must be solved numerically. This is cumbersome and does not allow us to obtain explicit analytical estimates as we would like. Furthermore, Ireland et al. [178] have shown the estimation of ZA to be accurate only up to , even when the nonlocal diffusion correction is applied due, they argue, to its poor prediction of the relative velocity statistics at .
We have thus taken a different path at . We use the fitted expression for given by the former authors in this regime. It is again based on the powerlaw model of Eq.~2.43, where the and are empirically derived. Evidence supporting the powerlaw scaling of the RDF has been presented for [288,178].
Once again, the coefficients and are functions of the only. This hypothesis can be seen to be reasonable up to and is independent (or only weakly dependent) of the Reynolds number, at least for the range of Reynolds numbers investigated so far [178]. A set of values of and , as a function of where obtained from DNS results and plotted in [178]. We provide a simple fit to this data points from the digitalized images of these plots, see Fig. 6. We found that a very convenient model for our fit was given by

(2.49) 
which is proportional to a lognormal distribution probability density function. The optimal parameters and were found to be

(2.50) 
(a)  (b) 
Figure 6: Fits for the and as a function of . The data were extracted from the digitalized Figure 22 in [178]. The continuous curves correspond to the fits in Eq.~2.49. The values of in parenthesis show the root mean square error of the approximation. 
Now, using Eq.~2.43 and stitching together the corrected version of ZA with the fits above for and we are able to produce an explicit analytic formula for the settling coefficient that can be used to establish the firsteffects of the influence of the neighbouring particles in a situation were nonuniform suspensions are to be expected, due to turbulence. This formulation is furthermore limited to heavy particles, but can be used as a conservative estimation for lighter particles while a more complete theory is still unavailable. The condition for neglecting the collective effects of neighbours then becomes, based on Eq.~2.33, or, assuming a maximum value for the acceptable influence at 1 change in the settling velocity, the maximum solid fraction can be expressed as

(2.51) 
where has the expression of Eq.~2.45 and where the expression for and can be taken from Eq.~2.47 (dividing by 4.2x10^{1}) for small values of and using Eqs.~2.49 and 2.50 for larger values. The behaviour of as a function of have been plotted in Fig. 7a for different values of the relative density.
(a) Settling coefficient as a function of the Stokes number using Eq.(2.45) with Eq.(2.47) for (with the nonlocal diffusion correction) and Eq.(2.50) for larger values.  (b) Estimation of the settling velocity using the same criterion with the data from [4]; the dots correspond to empirical data from [4], digitalized from Figure 17 in [381]. 
Figure 7: Effect of the Stokes number on the settling velocity based on the theory of Zaichik and Alipchenkov [381], with corrections from Bragg and Collins [44], Ireland et al. [178]. 
Note that the curves in Fig. 7a show a change in concavity at around toward a positive one. This is unexpected, as it has been empirically observed that tends to peak at one, and then monotonically decrease for larger Stokes numbers. This signals a clear flaw of the present model. Nonetheless, we believe the predictions are much more likely to be correct up to , as it is at around this point that the first signs of Reynolds number dependence in the and values are observed [44]. Therefore, Fig. 7a should be treated with caution and probably only for values of the Stokes number below one. Beyond this point, the value attained at one could be used as a conservative upper bound, as we know the maximum clustering is expected at around this point.
A final remark about this subject is the following. The powerlaw model for the RDF is supposed to be valid all the way down, below subKolmogorov scales, in accordance with DNS results for monodisperse suspensions. However, such model might suffer from the effects of a certain level of overidealization. Chun et al. [73] studied this issue theoretically for smallinertia particles and concluded that arbitrarily small differences between the diameters of the particle pairs lead to the appearance of a cutoff distance below which the RDF approaches a constant. This observation goes in the direction of moderating the amount of clustering expected at very small separations in realworld situations and should be kept in mind when interpreting our results, since our analysis has been limited to monodisperse dispersions.
Effect of collisions
In this work we are mainly interested in CFDDEM methods, where the MRE is enriched by the addition of contact forces to its RHS. In such a case the effect of collisions is naturally accounted for and the need to estimate the effects of its neglect diminishes. Nonetheless, in practice there are many situations in which it becomes advantageous to turn off interparticle interactions:
Therefore, it is interesting to study under what conditions interparticle interactions can be eliminated without incurring in serious error.
The evaluation of the first effects of collisions must be based on a statistical analysis because any useful characterization of a typical collision has to consider both the magnitude and angle of the approaching particles' velocities, which are stochastic variables by virtue of the chaotic nature of turbulent flows. Such characterization is challenging, as it must be derived from local and historydependent knowledge about the flow and the particles, including (but not limited to) the modelling of clustering that we discussed in the previous section.
The problem of estimating the rate of collision is central to the study of rain formation in clouds and, for a long time, it has been (and continuous to be) a topic of constant developments motivated by this problem and others [288, 123, 70].
For our purposes, though, we can use relatively crude estimates to bound the importance of these effects that can be useful for engineering considerations. In this regard, Loth [216] has established that interparticle collisions have a negligible effect to the overall trajectories of the particles if the average intercollision time is substantially larger than the particle relaxation time. This criterion naturally gives rise to the definition of a collisional Stokes number as

(2.52) 
where is the reciprocal of the average collision rate, . The above criterion thus requires an estimation of the collision rate, valid in as wide a range of regimes as possible.
In order to do so, one may use the following formula, that relates the collision rate of two generic neighbours, labelled , with their expected relative velocity:

(2.53) 
where

(2.54) 
where is the swept area formed by the projected silhouette of the two particles onto the plane orthogonal to , and are the respective radii and is the number density of the species of the particle . With this formula, the problem reduces to the determination of the expected number density (which can be obtained from the RDF) and the expected relative velocities, where the value of can be determined as a function of these values for spheres.
Loth [216] distinguishes between two sources of collisions, corresponding to different causes for the existence of the relative velocity: those due to turbulent fluctuations of the flow and those due to the different terminal velocities of particles of different terminal velocities under the influence of gravity.
For the case of gravitydriven settling, one can estimate the collision frequency experienced by a particle by calculating the swept volume per unit time of its neighbours as seen from a frame of reference moving with the average speed of the particles ensemble. Assuming monodispersity, the average number of collisions in a given time interval can be calculated as , corresponding to a cylindrical volume of base and height , where is an average relative velocity of a particle with respect to its surrounding ensemble. For small differences in the settling velocity, due to weak polydispersity this yields

(2.55) 
and the equation is valid to first order in , where is the range in diameters normalized by the average diameter.
The turbulencedriven collision frequency has been intensely studied [199,383,358]. Here we follow Loth [220] who bases the analysis on two limiting expressions are valid for two extreme regimes: the very low and very high Stokes numbers. For very low Staffman and Turner [301] derived the following expression, which has been validated in numerical simulations in the limit of very small Stokes numbers [331].

(2.56) 
where represents the usual Stokes number based on the Kolmogorov microscales.
At the opposite extreme, in the socalled uncorrelated regime, the following expression can be derived [220]

(2.57) 
where is the Stokes number based on the eddy turnover time and where is the lateral integral scale of the turbulent flows [280].
A simple approximation was proposed by Loth [220] to bridge both regimes in between (intermediate regime). This approximation is rough but sufficient to provide a first estimation of the importance of tracking collisions in dilute particleladen simulations. The result is

(2.58) 
For more accurate estimates see [331,383,71], where more advanced models for the collision frequency in the intermediate regime are discussed.
In order to unify most criteria, we can take, as a first approximation, the minimum of the two definitions Eqs.~2.55 and 2.58, so that the condition to be able to neglect collisions in turbulence becomes

(2.59) 
(^{1}) The asterisk has been added to avoid a symbolic clash with the settling coefficient, .
(^{2}) For particles lighter than the surrounding fluid, such as air bubbles in water, the inverse type of oversampling is observed [5], as expected.
(^{3}) Note that there appears in Equation 79 in [381] the lower limit of the integral is 0.0. Note however that this is inconsequential, since there is no contribution of the integral between 0.0 and because no particle centres can be found closer than .
In theory, the validity of the MRE should tend to be fully attained in the limit of a vanishingly small particle radius. Indeed, Eqs.~2.72.9 can all be fulfilled by making small enough. Small particles are thus the realm of the MRE. However, the implicit assumption of the applicability of the Navier–Stokes equations to the continuous phase in this limit is in fact flawed because the continuum hypothesis fails at sufficiently small sizes. In this section we look to determine the order of magnitude of the lower acceptable limits for the particle size for which the MRE still applies. Since the behaviour of a fluid near the referred limit is different for gases and liquids, we treat each one separately, keeping in mind mainly water and air as the archetypical examples.
Gas flow
At large scales, the evolution of gases is well described by the Navier–Stokes equations, which are based on the hypotheses of continuum and of pointwise thermodynamic equilibrium (or at least quasiequilibrium) that lead to the linear stressstrain rate relation [135]. The progressive deterioration of the validity of this assumption as the length scales shrink manifests itself as the socalled rarefaction effects. These effects arise due to the growing mean free path (the average distance that the molecules travel between successive impacts [313]) relative to the particle dimensions and the corresponding lower frequency rate of collisions ^{1}. The Knudsen number, , specifically quantifies this phenomenon and is defined as

(2.60) 
with the mean free path (MFP) and the representative physical dimension of the flow, which is usually taken as the particle radius in this context. In reducing the particle radius, the first rarefaction effects to be expected to appear are related to the breakdown of the noslip boundary condition [193], which is used to derive the MRE. Let us thus investigate the issue a bit further.
Upon impacting on a solid wall, the gas molecules exchange momentum in the tangential direction in other than idealized conditions. If all the reflections were merely specular (symmetric with respect to the surface average normal), there would be no average exchange of tangential momentum and one would have slip boundary conditions. For most materials, however, this is not the case [21], resulting in a net tangential momentum exchange. The main mechanism for such exchange is mediated by the molecularscale roughness on the wall, that randomizes the reflection angle of the particles, absorbing a part of the statistical bias in the tangential component of their momentum. The exchange is statistically systematic as long as the average relative velocity between gas and wall is nonzero.
As a result of the solid movement, near to the wall, the gas becomes statistically out of equilibrium, because there is a large number of molecules biased toward the wall velocity. This spatial shift in the mean velocity is rapidly homogenized over the neighbouring particles, moving away from the wall. This process is nonlinear (unlike the transmission of the macroscopic momentum, which, at close distances from the wall is linear) and spans a few mean free paths. Such outofequilibrium region is known as the Knudsen layer [193]. We can conclude that the assumption of noslip boundary condition at the solid wall must introduce an error of the order of , that is, proportional to .
While such an error is the manifestation of the onset of noncontinuum effects, the continuum model still has application for a range of Knudsen number values after these effects are first noticeable. This can be achieved through a generalization of the boundary conditions to partial slip (see [390] for a recent review) if is small enough. The set of low conditions for which this strategy is appropriate is called the slip regime. If one pretends to extend the applicability of the MRE to this regime, it is evident that one would need to modify the equation to adapt to the replacement of the noslip condition. Let us look into this.
At the very small particles limit that we are examining, the particle Reynolds number is expected to be extremely low, and so is the Stokes number (unless one is interested in high frequency oscillations, see [79]). This means that the prevailing force is the Stokes drag in this region (see Section 2.3). So to estimate the order of the error introduced by the use of the MRE in this range, we will first look at the drag force corrections for the lowKnudsen number limit (weak noncontinuum effects), neglecting for now the other forces.
As previously mentioned, the NavierStokes equations must be considered along slip boundary conditions. The simplest and most commonly used model is given by the Navier–Maxwell–Basset formula [390], that assumes the tangential slip velocity at any surface point on the sphere to be proportional to the local tangential stress. That is:

(2.61) 
where is the local tangential stress and is known as the Basset slip coefficient, in general ranging from zero (perfect slip) to infinity (noslip conditions).
The corresponding drag force on a sphere moving in uniform motion through an otherwise quiescent infinite fluid is (see, e.g. [283])

(2.62) 
since is a positive number.
In 1879, Maxwell derived such boundary condition from the kinetic theory of gases under isothermal conditions (we will assume isothermal conditions are a good approximation for the applications we are interested in). Those are given by the following expression, which relates the coefficient with the mean free path from the kinetic theory of monatomic gases in isothermal conditions.

(2.63) 
where is the tangential momentum accommodation coefficient ^{2}. This parameter is defined as the proportion of molecule collisions that result in a diffusive reflection (as opposed to specular) over the total. A molecule is said to undergo diffusive collision when its incidence angle (with respect to the wall's normal) is uncorrelated with the reflected angle, which is randomly distributed. On the contrary, specular reflection means that both incident and reflected angle coincide. The average of the tangential momentum exchange over a number of diffusive collisions equals the average incident momentum, since the average over the reflected momenta has null mean (noslip). In contrast, the averaged momentum exchange for a number of specular collision has zero mean, and the momentum transfer is null (full slip). The average in reality is somewhere in between these extremes. The value of depends on the gas and the surface material and finish, but experimental evidence shows to be between 2x10^{1} and 8x10^{1} (or ), with the former corresponding to specially treated, ultrasmooth surfaces and the latter to most practical surfaces [133]. Maxwell model is derived, by taking the first order (in ) effects into account, and thus its validity is bounded below by ; the socalled slip regime [193]. Higher order approximations have been derived, but this is unnecessary for our purposes.
Substituting these numbers into Eq.~2.62, we come up with the approximate requirement that , which imposing our accustomed one percent error bound reads

(2.64) 
for typical (most) surfaces. Taking the mean free path of atmospheric air as a reference, this means the particles radius should be larger than about 3 μm . Note that this requirement can become more strict for highly polished surfaces, dropping to for , putting the minimum radius at a ten times larger value.
While smallfrequency conditions are the most relevant for a majority of applications with low , it is still possible that higher frequencies arise, either due to the flow's own turbulent fluctuations or due to external forces, such as forced vibrations due to the workings of certain machinery etc. While an analysis analogous to that of Coimbra and Rangel is not available yet for small , a recent paper [376] addresses the case of imposed harmonic motion on the particle submerged in quiescent gas (rather than the perhaps more representative case of forced motion on the fluid and a free particle treated by the former authors) for small . We consider this case to be representative for unsteady motion in low conditions, at least for orderofmagnitude arguments. In the referred paper, expressions are provided for the asymptotic limits in which we are interested, that is, the limit of , where is a nondimensional frequency defined as (see [376])

(2.65) 
where is the molecular collision frequency, is the characteristic frequency of the relative motion and is a Stokes number ^{3}, defined as

(2.66) 
The first of these nondimensional numbers measures how the characteristic length scales of the system compare to the mean free path of the gas, which is related to the statistical convergence of the spatial homogenization, while the second number measures the ratio of the characteristic frequencies to the molecular collision frequency, which is related to the temporal homogenization.
Yap and Sader [376] provide the expression for the force in the slip regime, along with the force in the fully noslip limit (which can be found, for example, as a solved exercise in [202]). The ratio of the amplitudes of both forces (with the noslip case in the denominator), which we refer to as , provides an adequate measure of the distance from the MRE range of applicability, to first order in the . For conciseness, we do not reproduce its complicated expression here, but instead directly plot it in Fig. 13. In it, we highlight the one percent and ten percent error curves. Note how the force decreases both as and grow, as the momentum transfer at the particle's surface weakens.Figure 8: Ratio of force calculated as in [376], valid for small and , to the magnitude of the force calculated with the MRE, valid under the continuum, noslip conditions on the surface of the sphere. 
Note that the intersection of the graph with the plane recovers the steadystate asymptotic behaviour discussed above. For the case shown, , consistent with [376]. Such value comes as a result of refined gaskinetic simulations, which correct expression Eq.~2.63 and taking , see [21]. While the final value is about twice as big as the typical values expected following the above arguments, its order of magnitude is well within the mentioned limits and thus consistent for this matter.
In Fig. 8 we have limited the range of and so that is bounded above by 1x10^{1} since, as mentioned, the validity of this formulation is confined to small values of and . Nonetheless, within these restrictions, we have chosen the range of values for in Fig. 8 to cover all relative motions that could be driven by turbulence within the range of applicability of the MRE, but leaving out extremely large frequencies that would have to be excited by some external force. Indeed, the modified Stokes number can be related to Eq.~2.4 by using the estimate , with the characteristic time of the relative motion (which in turbulencedriven flows will be governed by the fluid's small scales motion). This gives

(2.67) 
The order of magnitude of the coefficient of can reach values of over 1x10^{3} for some materials (e.g.heavy metals suspended in air). Thus, roughly speaking, the graph covers , well beyond our range of interest (see Table 1).
In summary, at the smallsize limit that we are analysing it is unlikely to encounter large Stokes numbers due to the tiny inertia of the particles. Nonetheless, one should keep in mind that under unsteady circumstances, the low limitation is intensified, so particular attention must be paid to these situations if the MRE is to be applied.
Liquid flow
The molecules in a liquid are typically much closer than in a gas under the same conditions. Specifically, while the typical mean free path in atmospheric air is in the order of ten molecular lengths, the same distance is only in the order of a single molecule in water [134]. This explains the low compressibility of fluids with respect to gases. The result is that the concept of mean free path, and thus , is less useful for the theoretical characterization of liquids than for that of gases. Furthermore, the hypothesis of low density gas, where the molecular collisions can be regarded as binary and instantaneous, which is at the heart of the classical statistical treatment of gases does not apply. The result is that the molecular theory for liquids is much less developed than that of gases and, for example, a reliable analysis of the breakdown of the noslip boundary condition is not available yet [134]. Another consequence of the higher density of gases is the change in the order of appearance of smallsize effects (see [249]; compare Figures 1 and 4). Indeed, here noncontinuum effects come before thermodynamic nonequilibrium effects. That is, statistical fluctuations become significant before the noslip condition breaks down. These statistical fluctuations are due to a finite number of molecular collisions, leading to an erratic movement of the particles known as Brownian motion. Indeed, the assumption of a valid noslip boundary condition was employed by Einstein in his famous analysis of Brownian motion. Thus, the first smallsize effects in liquid suspensions are discussed in the following subsection, which deals with Brownian motion. The theory is however relevant to both gases and liquids.
Statistical fluctuations: Brownian motion
In order to study the onset of effect of Brownian fluctuations on the movement of the suspended particles, it is necessary to specify the length and time scales of interest. The situation is reminiscent of the study of turbulent dispersion of suspensions, where the averaged dispersion can be studied without having to resolve the smallest scales of the flow. This however means that their effect must be accounted for in an averaged sense. Similarly, Brownian motion presents a wide spectrum of time and length scales. Here too it is possible to avoid resolving all these scales and instead focus on their averaged effect on the particles motion, if the length scales of interest are significantly larger. Moreover, and in contrast with the turbulent dispersion problem, the statistical theory of Brownian motion has been well understood (at least for large times) for over a century, since the works of Einstein, Sutherland and Smoluchowski (see [39], where the relevant citations can be found, as well as a brief historical account of the subject). Let us thus begin by introducing the relevant time scales associated to Brownian motion.
Following Bian et al. [39] ^{4}

(2.68) 
where is the equivalent mass (particle mass plus added mass), is the local sound speed in the fluid, is the inverse of the mobility (that is, the constant quantity in the expression of drag coefficient ^{5}); and is the Stokes–Einstein–Smoluchowski diffusivity, given by the formula

(2.69) 
where is Boltzmann's constant and is the absolute temperature.
The meaning of the different time scales is as follows: is the time it takes a sound wave to travel through the fluid a distance equivalent to the particle's radius, ; is the time it takes for a particle to lose a fraction of its initial slip velocity, due to the action of the steady drag force and added mass forces alone; is the time it takes vorticity to diffuse over a length from the particle's surface, while the time scale is the longterm average time it takes for a Brownian particle to diffuse, again, to a distance (see the Langevin equation Eq.~2.79). Note that the definitions of (see beginning of this Section 2.2) and of (see Section 2.2.1) had already been introduced. We repeat them here for the reader's convenience.
Figure 9 shows the orders of magnitude of the different time scales for the same material combinations of Table 1 as a function of the particle radius. The individual combinations have not been identified for simplicity, but they can be inferred. Note that the values roughly reflect the succession in Eq.~2.68. That is, is much smaller than the rest, is typically smaller than for gases, while both are comparable for liquids and is much larger than the rest for all but the tiniest lengthscales, where the continuum hydrodynamic theory has started to break down anyway. Manifestly, such ordering is robust, at least for particles larger than 1x10^{2} nm.Figure 9: Time scales of Brownian motion normalized by the viscous diffusion time for the same combinations considered in Table 1. Dots indicate water as the suspending fluid, while crosses imply air. The size of the markers grows with the density of the particles. 
Sticking to the common theme of this section, we will now attempt to link this theory to the MRE. Our goal is to mark the onset of Brownian motion, providing a simple rule to assess whether Brownian effects can be safely neglected. We proceed as before, comparing these effects to that of other forces as the radius of the particle shrinks, since we know Brownian motion is mostly relevant at small sizes. The focus is placed on hydrodynamics, so we should compare Brownian motion to the motion caused by hydrodynamic forces. We would also like to take the buoyancy force into account.
Comparing different hydrodynamic forces is meaningful only when their characteristic times of application are the same, which is equivalent to comparing impulses. In turbulence, the natural reference is the microscopic Kolmogorov timescale, since it determines the time scale of the relative motion between fluid and particles (when turbulence is driving the relative motion). From hydrodynamic theory, we know that the Reynolds number represents the ratio of inertial versus diffusive momentum transport rates, that is . At the Kolmogorov scales the Reynolds number is by definition of order one, so

(2.70) 
since we have by assumption that . In fact, by the same argument, we have the even stronger

(2.71) 
This will have immediate consequence in our analysis, as we will see in the following lines.
On the other hand, Brownian motion presents a full spectrum of time and length scales for all time intervals below , due to its pseudofractal nature, down to the ballistic regime ^{6}, where the velocity of the particle becomes well defined as its trajectory starts to smooth out when observed at such small scales. The ballistic regime starts roughly below [39]. This kind of motion is more appropriately described in terms of statistical averages. Here we consider the mean squared displacement (MSD) or, rather, its square root as a measure of the displacement caused by Brownian motion:

(2.72) 
where

(2.73) 
That is, the ensemble average of the Euclidean norm of the displacement after a given time lapse . In this way, can be interpreted as the amplitude of the Brownian motion at the time scale . So that, instead of comparing forces, we can compare the characteristic displacements and to each other, where the latter is the characteristic displacement due to the hydrodynamic forces, still to be determined.
We expect the very small particles for which Brownian motion might be relevant to follow the background flow closely, due to their small inertia^{7}. So their velocity is expected to be overwhelmingly dominated by this contribution. However, it is the component of the motion relative to the fluid motion that determines many interesting phenomena in particulate flows, such as preferential concentration [112] or particle collisions [30], and so we will consider this component. After multiplying it by , we will obtain an adequate measure of displacement, due to the hydrodynamic forces.
Let us consider a turbulent flow. If the flow of interest were laminar, the same reasoning applies by replacing the Kolmogorov microscales with the typical scales of the problem at hand. If the main assumptions of the MRE are met and the particles are very small, so as to make Brownian motion potentially relevant, the scaling analysis of Balachandar [19] applies and the order of magnitude of the slip velocity can be estimated as:

(2.74) 
where

(2.75) 
and is the Kolmogorov time scale, while is the magnitude of the local fluid velocity. The estimate is derived under the hypothesis that the particle is most of the time close to steady state; that is, that the drag force dwarfs the other terms in the MRE most of the time. Such situation can be expected for very small Stokes numbers (for , according to Balachandar, as we will discuss in the following section^{8}).
The estimate above does not include the effect of buoyancy (and weight). Balachandar takes this into account separately by keeping track of the point at which the effect of gravity dominates the other forces. That is, when the settling velocity exceeds the slip velocity from equation Eq.~2.74, which happens beyond the point where the Kolmogorov acceleration (, with the turbulent energy dissipation rate) is smaller than . In order to derive our orderofmagnitude relation, it is sufficient to separately consider the steady settling velocity, , which is obtained by balancing the drag force with the submerged weight, giving

(2.76) 
Once the order of magnitude of the slip velocity has been identified, we can construct the following timedependent distance scale

(2.77) 
which is a measure of the inertial deviation from the perfect tracer particle trajectory of a given particle. Similarly, we define .
Now, the question is, in what circumstances is the effect of Brownian motion in deviating the particles' trajectories from the fluid stream lines of an order comparable to that caused by the particles' inertia (or buoyancy force)? We claim that an answer to this question can be given in terms of the quotient (or ); that is, the ratio of contributions to the deviation from a perfect tracer particle over one Kolmogorov eddy turnover. We propose to consider that the contribution of Brownian motion may be reasonably neglected whenever

(2.78) 
The Langevin theory [345] provides the right tool to calculate . This theory models the movement of the particle as a result of Newton's second law of motion, where the forces are given by a macroscopic contribution or resistance plus a random force originated from the molecular collisions which is modelled as Gaussian white noise ^{9}. In particular, the oftencalled modified Langevin equation has recently been experimentally verified [172,160] to be accurate all the way down to the ballistic regime. This version of the equation includes, apart from the steady drag, all the other terms in the MRE. From it, Clercx and Schram [75] have derived the exact expression of the MSD, which we have written as a function of the nondimensional time :

(2.79) 
where

(2.80) 
with

(2.81) 
The parameter is a small correction, relevant only for [277]. The function defined in Eq.~2.79 has two extreme regimes; those corresponding to (ballistic regime), and (diffusive regime). It is a direct calculation to check that their corresponding asymptotic relations are given by

(2.82) 
Eq.~2.79 can be used to compute . However, we would like to provide a simple rule that would serve for most situations. For that, it is enough to realize that the two asymptotic expressions Eq.~2.82 provide sufficient accuracy when used to estimate for and as the ballistic and diffusive regimes. By gluing together the two solutions at we generate the function , that approximates . Fig. 10 shows that the error made using such expression is at most of one order of magnitude, much less for the less extreme values of . Furthermore, such estimates are on the safe side, in the sense that they are upper bounds for the value of the MSD.
Fig. 10 (a). Evolution of 